Fun with benchmarking. Reverted to using built-in trigonometric functions!

src/grid.rs

@@ -9,7 +9,7 @@ pub struct Grid {
     height: usize,
     data: Vec<f32>,
 
-    // The scratch space for the blur operation.
+    // Scratch space for the blur operation.
     buf: Vec<f32>,
     blur: Blur,
 }

src/lib.rs

@@ -0,0 +1,4 @@
+mod blur;
+mod grid;
+pub mod model;
+pub mod trig;

src/main.rs

@@ -1,7 +1,4 @@
-mod blur;
-mod grid;
-mod model;
-mod trig;
+use physarum::model;
 
 fn main() {
     let model = model::Model::new(4, 4, 20, 1);

src/model.rs

@@ -1,4 +1,4 @@
-use crate::{grid::Grid, trig};
+use crate::grid::Grid;
 
 use rand::{thread_rng, Rng};
 
@@ -97,17 +97,18 @@ impl Model {
         }
     }
 
+    /// Perform a single simulation step.
     pub fn step(&mut self) {
         let sensor_distance = self.config.sensor_distance;
         let sensor_angle = self.config.sensor_angle;
 
         for agent in self.agents.iter_mut() {
-            let xc = agent.x + trig::cos(agent.angle) * sensor_distance;
-            let yc = agent.y + trig::sin(agent.angle) * sensor_distance;
-            let xl = agent.x + trig::cos(agent.angle - sensor_angle) * sensor_distance;
-            let yl = agent.y + trig::sin(agent.angle - sensor_angle) * sensor_distance;
-            let xr = agent.x + trig::cos(agent.angle + sensor_angle) * sensor_distance;
-            let yr = agent.y + trig::sin(agent.angle + sensor_angle) * sensor_distance;
+            let xc = agent.x + agent.angle.cos() * sensor_distance;
+            let yc = agent.y + agent.angle.sin() * sensor_distance;
+            let xl = agent.x + (agent.angle - sensor_angle).cos() * sensor_distance;
+            let yl = agent.y + (agent.angle - sensor_angle).sin() * sensor_distance;
+            let xr = agent.x + (agent.angle + sensor_angle).cos() * sensor_distance;
+            let yr = agent.y + (agent.angle + sensor_angle).sin() * sensor_distance;
 
             // Sense
             let trail_c = self.grid.get(xc, yc);

src/trig.rs

@@ -1,10 +1,12 @@
 /// From https://bits.stephan-brumme.com/absFloat.html
-pub(crate) fn abs(x: f32) -> f32 {
+#[inline(always)]
+fn abs(x: f32) -> f32 {
     f32::from_bits(x.to_bits() & 0x7FFF_FFFF)
 }
 
 /// Branchless floor implementation
-pub(crate) fn floor(x: f32) -> f32 {
+#[inline(always)]
+fn floor(x: f32) -> f32 {
     let mut x_trunc = (x as i32) as f32;
     x_trunc -= (x < x_trunc) as i32 as f32;
     x_trunc
@@ -12,7 +14,7 @@ pub(crate) fn floor(x: f32) -> f32 {
 
 /// Approximates `cos(x)` in radians with the maximum error of `0.002`
 /// https://stackoverflow.com/posts/28050328/revisions
-pub(crate) fn cos(mut x: f32) -> f32 {
+pub fn cos(mut x: f32) -> f32 {
     const ALPHA: f32 = 0.5 * std::f32::consts::FRAC_1_PI;
     x *= ALPHA;
     x -= 0.25_f32 + floor(x + 0.25_f32);
@@ -22,6 +24,32 @@ pub(crate) fn cos(mut x: f32) -> f32 {
 }
 
 /// Approximates `sin(x)` in radians with the maximum error of `0.002`
-pub(crate) fn sin(x: f32) -> f32 {
+pub fn sin(x: f32) -> f32 {
     cos(x - std::f32::consts::FRAC_PI_2)
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use itertools::repeat_n;
+
+    #[test]
+    fn test_cos() {
+        let n_points = 1000;
+        let x: Vec<f32> = repeat_n(std::f32::consts::TAU / (n_points - 1) as f32, n_points)
+            .enumerate()
+            .map(|(i, delta)| i as f32 * delta)
+            .collect();
+
+        let exact: Vec<f32> = x.iter().map(|v| v.cos()).collect();
+        let appr: Vec<f32> = x.iter().map(|v| cos(*v)).collect();
+
+        let mut max_error = 0.0_f32;
+        for (y_exact, y_appr) in exact.iter().zip(&appr) {
+            max_error = (y_exact - y_appr).abs().max(max_error);
+        }
+
+        // The error bound is even better!
+        assert!(max_error <= 0.0011);
+    }
+}