fix: mitigate floating point issues (#21)

* feat: modulo angles internally

* fix: properly adapt to size

* feat: use f64s to store values

* feat: apply correction after every full rotation

* fix: correct bad generation of complex rotation matrices

feat: add more vec operations

chore: remove ncube correction as it won't work on complex rotations

closes #10

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ncube"
-version = "2.3.0"
+version = "2.4.0"
 edition = "2021"
 authors = ["Nuno David <email@ndavd.com>"]
 license = "MIT"

src/camera.rs

@@ -1,3 +1,4 @@
+use crate::resources::SIZE;
 use crate::vec::{SphericalCoordinate, SphericalCoordinateSystem};
 use bevy::prelude::*;
 use bevy::render::camera::ScalingMode;
@@ -16,12 +17,12 @@ const WHEEL_FACTOR: f32 = if cfg!(target_family = "wasm") {
     0.001
 } else {
     0.3
-};
+} * SIZE;
 const MOTION_FACTOR: f32 = 0.3;
 
 pub fn get_default_camera_transform() -> Transform {
     Transform::from_translation(Vec3::from_spherical(SphericalCoordinate::new(
-        4.0,
+        4.0 * SIZE,
         std::f32::consts::PI / 2.0,
         0.0,
     )))
@@ -55,20 +56,24 @@ fn spawn_camera(mut commands: Commands) {
 }
 
 fn spawn_lighting(mut commands: Commands) {
+    let intensity = 1500.0 * SIZE.powi(2);
+    let range = 20.0 * SIZE;
     commands.spawn(PointLightBundle {
         point_light: PointLight {
-            intensity: 1500.0,
+            intensity,
+            range,
             ..default()
         },
-        transform: Transform::from_xyz(3.0, 8.0, 4.0),
+        transform: Transform::from_translation(Vec3::new(3.0, 8.0, 4.0) * SIZE),
         ..default()
     });
     commands.spawn(PointLightBundle {
         point_light: PointLight {
-            intensity: 1500.0,
+            intensity,
+            range,
             ..default()
         },
-        transform: Transform::from_xyz(-4.0, 8.0, -4.0),
+        transform: Transform::from_translation(Vec3::new(-4.0, 8.0, -4.0) * SIZE),
         ..default()
     });
 }
@@ -122,8 +127,8 @@ fn update_camera(
             curr_pos_spherical.r,
             (curr_pos_spherical.theta + event.delta.y.to_radians() * -MOTION_FACTOR)
                 // Prevent flipping effect
-                .min(std::f32::consts::PI * 0.99999)
-                .max(0.00001),
+                .min(std::f32::consts::PI * (1.0 - std::f32::EPSILON))
+                .max(std::f32::EPSILON),
             curr_pos_spherical.phi + event.delta.x.to_radians() * -MOTION_FACTOR,
         ));
         *camera_transform = camera_transform.looking_at(Vec3::ZERO, Vec3::Y);

src/main.rs

@@ -68,6 +68,7 @@ fn spawn_hypercube(
     mut ncube_vertices_3d: ResMut<NCubeVertices3D>,
     ncube_unlit: Res<NCubeUnlit>,
     ncube_edge_color: Res<NCubeEdgeColor>,
+    ncube_edge_thickness: Res<NCubeEdgeThickness>,
     ncube_face_color: Res<NCubeFaceColor>,
     q_ncube_entities: Query<Entity, With<NCubeMesh>>,
 ) {
@@ -81,14 +82,14 @@ fn spawn_hypercube(
             commands.entity(entity).despawn();
         });
 
-        **ncube = ncube::NCube::new(**ncube_dimension, 1.0);
+        **ncube = ncube::NCube::new(**ncube_dimension, ncube.size);
         let planes_of_rotation = usize::pair_permutations(0, **ncube_dimension - 1);
-        let mut rotations: HashMap<(usize, usize), (f32, f32)> = HashMap::new();
+        let mut rotations: HashMap<(usize, usize), (f64, f64)> = HashMap::new();
         let mut angles = Vec::new();
         for plane in &planes_of_rotation {
             let v = match ncube_rotations.get(plane) {
                 Some(v) => *v,
-                None => (0.0_f32, 0.0_f32),
+                None => (0.0, 0.0),
             };
             rotations.insert(*plane, v);
             angles.push(v.0);
@@ -113,7 +114,7 @@ fn spawn_hypercube(
                     ..default()
                 }),
                 transform: edge::Edge::transform(
-                    0.01,
+                    **ncube_edge_thickness,
                     ncube_vertices_3d[*i],
                     ncube_vertices_3d[*j],
                 ),
@@ -234,23 +235,26 @@ fn update_ncube_meshes(
 
 fn rotate_ncube(
     time: Res<Time>,
+    ncube_planes_of_rotation: Res<NCubePlanesOfRotation>,
+    ncube_is_paused: Res<NCubeIsPaused>,
     mut ncube: ResMut<NCube>,
     mut ncube_rotations: ResMut<NCubeRotations>,
-    ncube_planes_of_rotation: Res<NCubePlanesOfRotation>,
     mut ncube_vertices_3d: ResMut<NCubeVertices3D>,
-    ncube_is_paused: Res<NCubeIsPaused>,
 ) {
     if **ncube_is_paused {
         return;
     }
-    let dt = time.delta_seconds();
+    let dt: f64 = time.delta_seconds().into();
     let mut das = Vec::new();
     for i in 0..ncube_planes_of_rotation.len() {
         let plane = ncube_planes_of_rotation[i];
         let (angle, vel) = *ncube_rotations.get(&plane).unwrap();
         let da = dt * vel;
         das.push(da);
-        ncube_rotations.insert((plane.0, plane.1), (angle + da, vel));
+        ncube_rotations.insert(
+            (plane.0, plane.1),
+            ((angle + da) % std::f64::consts::TAU, vel),
+        );
     }
     ncube.rotate(&ncube_planes_of_rotation, &das);
     **ncube_vertices_3d = ncube.perspective_project_vertices();

src/mat.rs

@@ -9,12 +9,12 @@ macro_rules! emat {
 pub struct Mat {
     pub rows: usize,
     pub cols: usize,
-    matrix: Vec<f32>,
+    matrix: Vec<f64>,
 }
 
 #[allow(dead_code)]
 impl Mat {
-    pub fn new(mat: &[&[f32]]) -> Self {
+    pub fn new(mat: &[&[f64]]) -> Self {
         let row_len = mat[0].len();
         assert!(mat.iter().all(|r| row_len == r.len()));
         Self {
@@ -40,29 +40,38 @@ impl Mat {
         }
     }
 
-    pub fn fill(n: f32, rows: usize, cols: usize) -> Self {
+    pub fn fill(n: f64, rows: usize, cols: usize) -> Self {
         Self {
             rows,
             cols,
             matrix: vec![n; rows * cols],
         }
     }
 
-    pub fn rotation(
+    pub fn rotation(rows: usize, cols: usize, plane: (usize, usize), theta: f64) -> Self {
+        let mut m = Self::identity(rows, cols);
+        let assign_element = |element: &mut f64, v: f64| {
+            *element = if *element == 0.0 { v } else { *element * v };
+        };
+        assign_element(&mut emat!(m[plane.0][plane.0]), theta.cos());
+        assign_element(&mut emat!(m[plane.0][plane.1]), -theta.sin());
+        assign_element(&mut emat!(m[plane.1][plane.0]), theta.sin());
+        assign_element(&mut emat!(m[plane.1][plane.1]), theta.cos());
+        m
+    }
+
+    pub fn from_rotations(
         rows: usize,
         cols: usize,
         planes: &Vec<(usize, usize)>,
-        thetas: &Vec<f32>,
+        thetas: &Vec<f64>,
     ) -> Self {
         let mut m = Self::identity(rows, cols);
-        let assign_element = |element: &mut f32, v: f32| {
-            *element = if *element == 0.0 { v } else { *element * v };
-        };
         for i in 0..planes.len() {
-            assign_element(&mut emat!(m[planes[i].0][planes[i].0]), thetas[i].cos());
-            assign_element(&mut emat!(m[planes[i].0][planes[i].1]), -thetas[i].sin());
-            assign_element(&mut emat!(m[planes[i].1][planes[i].0]), thetas[i].sin());
-            assign_element(&mut emat!(m[planes[i].1][planes[i].1]), thetas[i].cos());
+            if thetas[i] == 0.0 {
+                continue;
+            }
+            m = m * Self::rotation(rows, cols, planes[i], thetas[i])
         }
         m
     }
@@ -91,7 +100,7 @@ impl std::ops::Mul for Mat {
     type Output = Self;
     fn mul(self, rhs: Self) -> Self::Output {
         assert_eq!(self.cols, rhs.rows);
-        let mut m: Vec<f32> = Vec::with_capacity(self.rows * rhs.cols);
+        let mut m: Vec<f64> = Vec::with_capacity(self.rows * rhs.cols);
         for i in 0..self.rows {
             for j in 0..rhs.cols {
                 m.push(
@@ -107,9 +116,9 @@ impl std::ops::Mul for Mat {
     }
 }
 
-impl std::ops::Mul<Vec<f32>> for Mat {
-    type Output = Vec<f32>;
-    fn mul(self, mut rhs: Vec<f32>) -> Self::Output {
+impl std::ops::Mul<Vec<f64>> for Mat {
+    type Output = Vec<f64>;
+    fn mul(self, mut rhs: Vec<f64>) -> Self::Output {
         assert_eq!(rhs.len(), self.cols);
         let v = rhs.clone();
         rhs.truncate(self.rows);
@@ -124,9 +133,9 @@ impl std::ops::Mul<Vec<f32>> for Mat {
     }
 }
 
-impl std::ops::Mul<f32> for Mat {
+impl std::ops::Mul<f64> for Mat {
     type Output = Self;
-    fn mul(self, rhs: f32) -> Self::Output {
+    fn mul(self, rhs: f64) -> Self::Output {
         let mut m = self.clone();
         for i in 0..self.rows {
             for j in 0..self.cols {
@@ -136,9 +145,9 @@ impl std::ops::Mul<f32> for Mat {
         m
     }
 }
-impl std::ops::Add<f32> for Mat {
+impl std::ops::Add<f64> for Mat {
     type Output = Self;
-    fn add(self, rhs: f32) -> Self::Output {
+    fn add(self, rhs: f64) -> Self::Output {
         let mut m = self.clone();
         for i in 0..self.rows {
             for j in 0..self.cols {
@@ -148,9 +157,9 @@ impl std::ops::Add<f32> for Mat {
         m
     }
 }
-impl std::ops::Sub<f32> for Mat {
+impl std::ops::Sub<f64> for Mat {
     type Output = Self;
-    fn sub(self, rhs: f32) -> Self::Output {
+    fn sub(self, rhs: f64) -> Self::Output {
         let mut m = self.clone();
         for i in 0..self.rows {
             for j in 0..self.cols {
@@ -204,8 +213,8 @@ mod tests {
             &[2.0, 1.0, 2.0, 0.0],
             &[3.0, 0.0, 1.0, 2.0],
         ]);
-        let b: Vec<f32> = Vec::from([1.0, 2.0, 3.0, 4.0]);
-        let c: Vec<f32> = Vec::from([22.0, 10.0, 14.0]);
+        let b: Vec<f64> = Vec::from([1.0, 2.0, 3.0, 4.0]);
+        let c: Vec<f64> = Vec::from([22.0, 10.0, 14.0]);
         assert_eq!(a * b, c);
     }
 }

src/ncube.rs

@@ -55,7 +55,7 @@ where
 #[derive(Debug, Clone)]
 pub struct NCube {
     pub dimensions: usize,
-    pub size: f32,
+    pub size: f64,
     /// Cartesian coordinates of the vertices of the hypercube.
     pub vertices: NVertices,
     /// Vertex indices of the edges of the hypercube.
@@ -65,7 +65,7 @@ pub struct NCube {
 }
 
 #[derive(Debug, Clone)]
-pub struct NVertices(pub Vec<Vec<f32>>);
+pub struct NVertices(pub Vec<Vec<f64>>);
 
 #[derive(Debug, Clone)]
 /// Each edge is composed of 2 vertices (index)
@@ -102,7 +102,7 @@ impl std::fmt::Display for NVertices {
 #[allow(dead_code)]
 impl NCube {
     /// Creates an `n` dimensional hypercube of size `s`.
-    pub fn new(n: usize, s: f32) -> Self {
+    pub fn new(n: usize, s: f64) -> Self {
         let vertices = Self::_vertices(n, s);
         Self {
             dimensions: n,
@@ -118,11 +118,16 @@ impl NCube {
         Self::_face_count(self.dimensions, m)
     }
 
+    /// Computes the diagonal of the hypercube
+    pub fn diagonal_length(&self) -> f64 {
+        self.size * (self.dimensions as f64).sqrt()
+    }
+
     fn _face_count(n: usize, m: usize) -> usize {
         2_usize.pow((n - m).try_into().unwrap()) * n.permute(m)
     }
 
-    fn _vertices(n: usize, s: f32) -> NVertices {
+    fn _vertices(n: usize, s: f64) -> NVertices {
         let s = s / 2.0;
         let v_count = Self::_face_count(n, 0);
         let vertices = (0..v_count)
@@ -131,7 +136,7 @@ impl NCube {
                     .map(|j| {
                         let direction =
                             -1 + 2 * ((i / 2_usize.pow(j.try_into().unwrap())) % 2 == 0) as i8;
-                        s * direction as f32
+                        s * direction as f64
                     })
                     .collect::<Vec<_>>()
             })
@@ -163,7 +168,7 @@ impl NCube {
 
     // NOTE: This was not trivial
     fn _faces(vertices: &NVertices, n: usize) -> NFaces {
-        let extract_faces = |vertices: Vec<(usize, &Vec<f32>)>| {
+        let extract_faces = |vertices: Vec<(usize, &Vec<f64>)>| {
             vertices
                 .permute_four()
                 .iter()
@@ -211,18 +216,19 @@ impl NCube {
         NFaces(faces)
     }
 
-    pub fn rotate(&mut self, planes: &Vec<(usize, usize)>, theta_rads: &Vec<f32>) -> &mut Self {
+    pub fn rotate(&mut self, planes: &Vec<(usize, usize)>, theta_rads: &Vec<f64>) -> &mut Self {
         for vertex in &mut self.vertices.0 {
-            *vertex = Mat::rotation(self.dimensions, self.dimensions, planes, theta_rads)
+            *vertex = Mat::from_rotations(self.dimensions, self.dimensions, planes, theta_rads)
                 * vertex.clone();
         }
         self
     }
 
     pub fn perspective_project_vertices(&self) -> Vec<Vec3> {
         let projection_count = self.dimensions - 3;
-        let proj_m = |from_d: usize, to_d: usize, q: f32| {
-            Mat::identity(to_d, from_d) * (1.0 / (self.size * 1.5 - q))
+        let proj_m = |from_d: usize, to_d: usize, q: f64| {
+            let f = self.size / (self.size * 1.5 - q);
+            Mat::identity(to_d, from_d) * f
         };
         let mut v = self.vertices.0.clone();
         for i in 0..projection_count {
@@ -233,7 +239,9 @@ impl NCube {
                 v[v_index] = m * v[v_index].clone();
             }
         }
-        v.iter().map(|x| Vec3::new(x[0], x[1], x[2])).collect()
+        v.iter()
+            .map(|x| Vec3::new(x[0] as f32, x[1] as f32, x[2] as f32))
+            .collect()
     }
 }