feat: algebra

Cargo.toml

@@ -10,9 +10,4 @@ readme = "README.md"
 rust-version = "1.80.0"
 
 [dependencies]
-extensor-macros = { path = "macros/", version = "0.1.0" }
-
-
-[dev-dependencies]
-log = "0.4.21"
-env_logger = "0.11.3"
+const-default = "1.0"

src/algebra/mod.rs

@@ -0,0 +1,39 @@
+use core::marker::PhantomData;
+
+use module::{Module, Vector};
+use tensor::Tensor;
+
+use super::*;
+
+pub mod quadratic_form;
+
+// #[derive(Clone, Copy)]
+// pub struct TensorAlgebra<M: Module>(PhantomData<M>);
+
+// impl<M:Module> Module for TensorAlgebra<M> {
+//     type Ring = M::Ring;
+// }
+
+// impl<M:Module> Add for TensorAlgebra<M> {
+//     type Output = TensorAlgebra<M>;
+
+//     fn add(self, rhs: Self) -> Self::Output {
+//         todo!()
+//     }
+// }
+
+// impl<M:Module> Neg for TensorAlgebra<M> {
+//     type Output = TensorAlgebra<M>;
+
+//     fn neg(self) -> Self::Output {
+//         todo!()
+//     }
+// }
+
+// impl<M:Module> Mul<M::Ring> for TensorAlgebra<M> {
+//     type Output = TensorAlgebra<M>;
+
+//     fn mul(self, rhs: M::Ring) -> Self::Output {
+//         todo!()
+//     }
+// }

src/algebra/quadratic_form.rs

@@ -0,0 +1,63 @@
+use core::ops::AddAssign;
+
+use module::{Module, Vector};
+
+use super::*;
+
+pub struct SymmetricIndex<const D: usize>
+where
+    [(); D * (D + 1) / 2]:,
+{
+    // Store only upper triangular indices
+    // Length is (D * (D + 1)) / 2
+    indices: [(usize, usize); (D * (D + 1)) / 2],
+}
+
+pub struct QuadraticForm<const D: usize, F> {
+    eigenbasis: [Vector<D, F>; D],
+    eigenvalues: [F; D],
+}
+
+impl<const D: usize, F: Copy + ConstDefault + Mul<Output = F> + AddAssign> QuadraticForm<D, F> {
+    pub const fn new_diagonal(eigenvalues: [F; D]) -> Self {
+        let mut eigenbasis = [Vector::DEFAULT; D];
+        let mut i = 0;
+        while i < D {
+            let mut j = 0;
+            while j < D {
+                if i == j {
+                    eigenbasis[i].0[i] = eigenvalues[i];
+                }
+                j += 1;
+            }
+            i += 1;
+        }
+
+        Self {
+            eigenbasis,
+            eigenvalues,
+        }
+    }
+
+    // TODO: Make `const` if we ever get a const `Mul`
+    pub fn eval(&self, lhs: Vector<D, F>, rhs: Vector<D, F>) -> F {
+        let mut sum = F::DEFAULT;
+        let mut i = 0;
+        while i < D {
+            // Project vectors onto eigenbasis
+            let mut lhs_comp = F::DEFAULT;
+            let mut rhs_comp = F::DEFAULT;
+            let mut j = 0;
+            while j < D {
+                lhs_comp += lhs.0[j] * self.eigenbasis[i].0[j];
+                rhs_comp += rhs.0[j] * self.eigenbasis[i].0[j];
+                j += 1;
+            }
+
+            // Multiply by eigenvalue and add to sum
+            sum += lhs_comp * rhs_comp * self.eigenvalues[i];
+            i += 1;
+        }
+        sum
+    }
+}