Posted January 20, 2024 by Tomwol
#geometric algebra #PGA #XPBD
This is in no way an exhaustive list of all of the problems I ran into, but I wanted to list out some of the more memorable ones. Also when I was looking into using 4D PGA, I started off trying to implement it using the rust code generated at https://bivector.net/tools.html?p=4&q=0&r=1 and from what I knew about PGA from watching sudgylacmoe's videos and having implemented GA a little bit before. I wasn't really getting anywhere, until I found the resources at https://bivector.net/doc.html. Note: A lot of the problems I had probably could have been avoided by me properly learning this stuff before trying to implement it.
For XPBD I had followed https://johanhelsing.studio/posts/bevy-xpbd up to and including part 4, so I had, on my own, added an orientation component that used a Mat4, however there was still no rotational velocity. I had also generalized the shapes to use the same geometry builder as the manifold so I could write the contact solver more generally (although only hyperspheres can currently be dynamic).
Since position and rotation is combined, some of the stored quantities have different names:
fn solve_rate( mut query: Query<(&Pose, &mut Rate, &PrevRate, &Inertia, &Restitution), Without<isstatic>>, contacts: Res<contacts> ) { for (entity_a, entity_b, contact) in contacts.0.iter().cloned() { let ( (pose_a, mut rate_a, prev_rate_a, inertia_a, restitution_a), (pose_b, mut rate_b, prev_rate_b, inertia_b, restitution_b), ) = query.get_pair_mut(entity_a, entity_b).unwrap(); let n = contact.normal; let restitution = (restitution_a.0 + restitution_b.0) / 2.0; let normal_rate: Bivector = (n.x * e01 - n.y * e02 + n.z * e03 - n.w * e04).into(); let prev_relative_a: Bivector = (prev_rate_a.mv() - pose_a.inverse().apply(pose_b.0.apply(prev_rate_b.mv()))).into(); let prev_relative_b = (prev_rate_b.mv() - pose_b.inverse().apply(pose_a.0.apply(prev_rate_a.mv()))).into(); let relative_a = (rate_a.mv() - pose_a.inverse().apply(pose_b.0.apply(rate_b.mv()))).into(); let relative_b = (rate_b.mv() - pose_b.inverse().apply(pose_a.0.apply(rate_a.mv()))).into(); rate_a.0 += 0.5 * -normal_rate * (normal_rate.dot(relative_a) + (restitution * normal_rate.dot(prev_relative_a)).min(0.0)); rate_b.0 -= 0.5 * -normal_rate * (normal_rate.dot(relative_b) + (restitution * normal_rate.dot(prev_relative_b)).min(0.0)); } }
My plan, as of right now, is to eventually continue with 4D PGA after I work on some other features, hoping that when I go back to the code with a fresh mind I'll be able to fix it up. If that doesn't work out, I have a two alternatives in mind:
Currently the R401.rs file generated on bivector.net does not compile, the version that I edited and am now using is below. Although the current code generator in the ganja.js repo does work: codegen.
For the things that need to be fixed, it's easiest if you can do find/replace with regex:
(?<=[\+\-*])(\d)|(\d)(?=[\+\-*])
and substitute with $1.0
to fix this// Written by a generator written by enki. (edited by tomwol) #![allow(unused_imports)] #![allow(dead_code)] #![allow(non_upper_case_globals)] #![allow(non_snake_case)] #![allow(non_camel_case_types)] use std::fmt; use std::ops::{Index, IndexMut, Add, Sub, Mul, BitAnd, BitOr, BitXor, Not, Neg, Div}; type float_t = f32; // use std::f64::consts::PI; const PI: float_t = 3.14159265358979323846; const basis: &'static [&'static str] = &[ "1","e0","e1","e2","e3","e4","e01","e02","e03","e04","e12","e13","e14","e23","e24","e34","e012","e013","e014","e023","e024","e034","e123","e124","e134","e234","e0123","e0124","e0134","e0234","e1234","e01234" ]; const basis_count: usize = basis.len(); #[derive(Default,Debug,Clone,Copy,PartialEq)] pub struct R401 { mvec: [float_t; basis_count] } impl R401 { pub const fn zero() -> Self { Self { mvec: [0.0; basis_count] } } pub const fn new(f: float_t, idx: usize) -> Self { let mut ret = Self::zero(); ret.mvec[idx] = f; ret } pub const fn with(mut self, f: float_t, idx: usize) -> Self { self.mvec[idx] = f; self } } // basis vectors are available as global constants. pub const e0: R401 = R401::new(1.0, 1); pub const e1: R401 = R401::new(1.0, 2); pub const e2: R401 = R401::new(1.0, 3); pub const e3: R401 = R401::new(1.0, 4); pub const e4: R401 = R401::new(1.0, 5); pub const e01: R401 = R401::new(1.0, 6); pub const e02: R401 = R401::new(1.0, 7); pub const e03: R401 = R401::new(1.0, 8); pub const e04: R401 = R401::new(1.0, 9); pub const e12: R401 = R401::new(1.0, 10); pub const e13: R401 = R401::new(1.0, 11); pub const e14: R401 = R401::new(1.0, 12); pub const e23: R401 = R401::new(1.0, 13); pub const e24: R401 = R401::new(1.0, 14); pub const e34: R401 = R401::new(1.0, 15); pub const e012: R401 = R401::new(1.0, 16); pub const e013: R401 = R401::new(1.0, 17); pub const e014: R401 = R401::new(1.0, 18); pub const e023: R401 = R401::new(1.0, 19); pub const e024: R401 = R401::new(1.0, 20); pub const e034: R401 = R401::new(1.0, 21); pub const e123: R401 = R401::new(1.0, 22); pub const e124: R401 = R401::new(1.0, 23); pub const e134: R401 = R401::new(1.0, 24); pub const e234: R401 = R401::new(1.0, 25); pub const e0123: R401 = R401::new(1.0, 26); pub const e0124: R401 = R401::new(1.0, 27); pub const e0134: R401 = R401::new(1.0, 28); pub const e0234: R401 = R401::new(1.0, 29); pub const e1234: R401 = R401::new(1.0, 30); pub const e01234: R401 = R401::new(1.0, 31); impl Index<usize> for R401 { type Output = float_t; fn index<'a>(&'a self, index: usize) -> &'a Self::Output { &self.mvec[index] } } impl IndexMut<usize> for R401 { fn index_mut<'a>(&'a mut self, index: usize) -> &'a mut Self::Output { &mut self.mvec[index] } } impl fmt::Display for R401 { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { let mut n = 0; let ret = self.mvec.iter().enumerate().filter_map(|(i, &coeff)| { if coeff > 0.00001 || coeff < -0.00001 { n = 1; Some(format!("{}{}", format!("{:.*}", 7, coeff).trim_end_matches('0').trim_end_matches('.'), if i > 0 { basis[i] } else { "" } ) ) } else { None } }).collect::<Vec<String>>().join(" + "); if n==0 { write!(f,"0") } else { write!(f, "{}", ret) } } } impl From<float_t> for R401 { fn from(value: float_t) -> Self { R401::new(value, 0) } } impl From<R401> for float_t { fn from(value: R401) -> Self { value[0] } } // Reverse // Reverse the order of the basis blades. impl R401 { pub fn reverse(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[0]; res[1]=a[1]; res[2]=a[2]; res[3]=a[3]; res[4]=a[4]; res[5]=a[5]; res[6]=-a[6]; res[7]=-a[7]; res[8]=-a[8]; res[9]=-a[9]; res[10]=-a[10]; res[11]=-a[11]; res[12]=-a[12]; res[13]=-a[13]; res[14]=-a[14]; res[15]=-a[15]; res[16]=-a[16]; res[17]=-a[17]; res[18]=-a[18]; res[19]=-a[19]; res[20]=-a[20]; res[21]=-a[21]; res[22]=-a[22]; res[23]=-a[23]; res[24]=-a[24]; res[25]=-a[25]; res[26]=a[26]; res[27]=a[27]; res[28]=a[28]; res[29]=a[29]; res[30]=a[30]; res[31]=a[31]; res } } // Dual // Poincare duality operator. impl R401 { pub fn dual(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[31]; res[1]=a[30]; res[2]=a[29]; res[3]=a[28]; res[4]=a[27]; res[5]=a[26]; res[6]=a[25]; res[7]=a[24]; res[8]=a[23]; res[9]=a[22]; res[10]=a[21]; res[11]=a[20]; res[12]=a[19]; res[13]=a[18]; res[14]=a[17]; res[15]=a[16]; res[16]=a[15]; res[17]=a[14]; res[18]=a[13]; res[19]=a[12]; res[20]=a[11]; res[21]=a[10]; res[22]=a[9]; res[23]=a[8]; res[24]=a[7]; res[25]=a[6]; res[26]=a[5]; res[27]=a[4]; res[28]=a[3]; res[29]=a[2]; res[30]=a[1]; res[31]=a[0]; res } } // Hodge star impl R401 { pub fn dual_h(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[31]; res[1]=a[30]; res[2]=-a[29]; res[3]=a[28]; res[4]=-a[27]; res[5]=a[26]; res[6]=a[25]; res[7]=-a[24]; res[8]=a[23]; res[9]=-a[22]; res[10]=a[21]; res[11]=-a[20]; res[12]=a[19]; res[13]=a[18]; res[14]=-a[17]; res[15]=a[16]; res[16]=a[15]; res[17]=-a[14]; res[18]=a[13]; res[19]=a[12]; res[20]=-a[11]; res[21]=a[10]; res[22]=-a[9]; res[23]=a[8]; res[24]=-a[7]; res[25]=a[6]; res[26]=a[5]; res[27]=-a[4]; res[28]=a[3]; res[29]=-a[2]; res[30]=a[1]; res[31]=a[0]; res } pub fn undual_h(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[31]; res[1]=a[30]; res[2]=-a[29]; res[3]=a[28]; res[4]=-a[27]; res[5]=a[26]; res[6]=a[25]; res[7]=-a[24]; res[8]=a[23]; res[9]=-a[22]; res[10]=a[21]; res[11]=-a[20]; res[12]=a[19]; res[13]=a[18]; res[14]=-a[17]; res[15]=a[16]; res[16]=a[15]; res[17]=-a[14]; res[18]=a[13]; res[19]=a[12]; res[20]=-a[11]; res[21]=a[10]; res[22]=-a[9]; res[23]=a[8]; res[24]=-a[7]; res[25]=a[6]; res[26]=a[5]; res[27]=-a[4]; res[28]=a[3]; res[29]=-a[2]; res[30]=a[1]; res[31]=a[0]; res } } impl Neg for R401 { type Output = Self; fn neg(self) -> Self::Output { -1.0 * self } } impl Not for R401 { type Output = R401; fn not(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[31]; res[1]=a[30]; res[2]=a[29]; res[3]=a[28]; res[4]=a[27]; res[5]=a[26]; res[6]=a[25]; res[7]=a[24]; res[8]=a[23]; res[9]=a[22]; res[10]=a[21]; res[11]=a[20]; res[12]=a[19]; res[13]=a[18]; res[14]=a[17]; res[15]=a[16]; res[16]=a[15]; res[17]=a[14]; res[18]=a[13]; res[19]=a[12]; res[20]=a[11]; res[21]=a[10]; res[22]=a[9]; res[23]=a[8]; res[24]=a[7]; res[25]=a[6]; res[26]=a[5]; res[27]=a[4]; res[28]=a[3]; res[29]=a[2]; res[30]=a[1]; res[31]=a[0]; res } } // Conjugate // Clifford Conjugation impl R401 { pub fn conjugate(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[0]; res[1]=-a[1]; res[2]=-a[2]; res[3]=-a[3]; res[4]=-a[4]; res[5]=-a[5]; res[6]=-a[6]; res[7]=-a[7]; res[8]=-a[8]; res[9]=-a[9]; res[10]=-a[10]; res[11]=-a[11]; res[12]=-a[12]; res[13]=-a[13]; res[14]=-a[14]; res[15]=-a[15]; res[16]=a[16]; res[17]=a[17]; res[18]=a[18]; res[19]=a[19]; res[20]=a[20]; res[21]=a[21]; res[22]=a[22]; res[23]=a[23]; res[24]=a[24]; res[25]=a[25]; res[26]=a[26]; res[27]=a[27]; res[28]=a[28]; res[29]=a[29]; res[30]=a[30]; res[31]=-a[31]; res } } // Involute // Main involution impl R401 { pub fn involute(self: Self) -> R401 { let mut res = R401::zero(); let a = self; res[0]=a[0]; res[1]=-a[1]; res[2]=-a[2]; res[3]=-a[3]; res[4]=-a[4]; res[5]=-a[5]; res[6]=a[6]; res[7]=a[7]; res[8]=a[8]; res[9]=a[9]; res[10]=a[10]; res[11]=a[11]; res[12]=a[12]; res[13]=a[13]; res[14]=a[14]; res[15]=a[15]; res[16]=-a[16]; res[17]=-a[17]; res[18]=-a[18]; res[19]=-a[19]; res[20]=-a[20]; res[21]=-a[21]; res[22]=-a[22]; res[23]=-a[23]; res[24]=-a[24]; res[25]=-a[25]; res[26]=a[26]; res[27]=a[27]; res[28]=a[28]; res[29]=a[29]; res[30]=a[30]; res[31]=-a[31]; res } } // Mul // The geometric product. impl Mul for R401 { type Output = R401; fn mul(self: R401, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0]=b[0]*a[0]+b[2]*a[2]+b[3]*a[3]+b[4]*a[4]+b[5]*a[5]-b[10]*a[10]-b[11]*a[11]-b[12]*a[12]-b[13]*a[13]-b[14]*a[14]-b[15]*a[15]-b[22]*a[22]-b[23]*a[23]-b[24]*a[24]-b[25]*a[25]+b[30]*a[30]; res[1]=b[1]*a[0]+b[0]*a[1]-b[6]*a[2]-b[7]*a[3]-b[8]*a[4]-b[9]*a[5]+b[2]*a[6]+b[3]*a[7]+b[4]*a[8]+b[5]*a[9]-b[16]*a[10]-b[17]*a[11]-b[18]*a[12]-b[19]*a[13]-b[20]*a[14]-b[21]*a[15]-b[10]*a[16]-b[11]*a[17]-b[12]*a[18]-b[13]*a[19]-b[14]*a[20]-b[15]*a[21]+b[26]*a[22]+b[27]*a[23]+b[28]*a[24]+b[29]*a[25]-b[22]*a[26]-b[23]*a[27]-b[24]*a[28]-b[25]*a[29]+b[31]*a[30]+b[30]*a[31]; res[2]=b[2]*a[0]+b[0]*a[2]-b[10]*a[3]-b[11]*a[4]-b[12]*a[5]+b[3]*a[10]+b[4]*a[11]+b[5]*a[12]-b[22]*a[13]-b[23]*a[14]-b[24]*a[15]-b[13]*a[22]-b[14]*a[23]-b[15]*a[24]+b[30]*a[25]-b[25]*a[30]; res[3]=b[3]*a[0]+b[10]*a[2]+b[0]*a[3]-b[13]*a[4]-b[14]*a[5]-b[2]*a[10]+b[22]*a[11]+b[23]*a[12]+b[4]*a[13]+b[5]*a[14]-b[25]*a[15]+b[11]*a[22]+b[12]*a[23]-b[30]*a[24]-b[15]*a[25]+b[24]*a[30]; res[4]=b[4]*a[0]+b[11]*a[2]+b[13]*a[3]+b[0]*a[4]-b[15]*a[5]-b[22]*a[10]-b[2]*a[11]+b[24]*a[12]-b[3]*a[13]+b[25]*a[14]+b[5]*a[15]-b[10]*a[22]+b[30]*a[23]+b[12]*a[24]+b[14]*a[25]-b[23]*a[30]; res[5]=b[5]*a[0]+b[12]*a[2]+b[14]*a[3]+b[15]*a[4]+b[0]*a[5]-b[23]*a[10]-b[24]*a[11]-b[2]*a[12]-b[25]*a[13]-b[3]*a[14]-b[4]*a[15]-b[30]*a[22]-b[10]*a[23]-b[11]*a[24]-b[13]*a[25]+b[22]*a[30]; res[6]=b[6]*a[0]+b[2]*a[1]-b[1]*a[2]+b[16]*a[3]+b[17]*a[4]+b[18]*a[5]+b[0]*a[6]-b[10]*a[7]-b[11]*a[8]-b[12]*a[9]+b[7]*a[10]+b[8]*a[11]+b[9]*a[12]-b[26]*a[13]-b[27]*a[14]-b[28]*a[15]+b[3]*a[16]+b[4]*a[17]+b[5]*a[18]-b[22]*a[19]-b[23]*a[20]-b[24]*a[21]+b[19]*a[22]+b[20]*a[23]+b[21]*a[24]-b[31]*a[25]-b[13]*a[26]-b[14]*a[27]-b[15]*a[28]+b[30]*a[29]-b[29]*a[30]-b[25]*a[31]; res[7]=b[7]*a[0]+b[3]*a[1]-b[16]*a[2]-b[1]*a[3]+b[19]*a[4]+b[20]*a[5]+b[10]*a[6]+b[0]*a[7]-b[13]*a[8]-b[14]*a[9]-b[6]*a[10]+b[26]*a[11]+b[27]*a[12]+b[8]*a[13]+b[9]*a[14]-b[29]*a[15]-b[2]*a[16]+b[22]*a[17]+b[23]*a[18]+b[4]*a[19]+b[5]*a[20]-b[25]*a[21]-b[17]*a[22]-b[18]*a[23]+b[31]*a[24]+b[21]*a[25]+b[11]*a[26]+b[12]*a[27]-b[30]*a[28]-b[15]*a[29]+b[28]*a[30]+b[24]*a[31]; res[8]=b[8]*a[0]+b[4]*a[1]-b[17]*a[2]-b[19]*a[3]-b[1]*a[4]+b[21]*a[5]+b[11]*a[6]+b[13]*a[7]+b[0]*a[8]-b[15]*a[9]-b[26]*a[10]-b[6]*a[11]+b[28]*a[12]-b[7]*a[13]+b[29]*a[14]+b[9]*a[15]-b[22]*a[16]-b[2]*a[17]+b[24]*a[18]-b[3]*a[19]+b[25]*a[20]+b[5]*a[21]+b[16]*a[22]-b[31]*a[23]-b[18]*a[24]-b[20]*a[25]-b[10]*a[26]+b[30]*a[27]+b[12]*a[28]+b[14]*a[29]-b[27]*a[30]-b[23]*a[31]; res[9]=b[9]*a[0]+b[5]*a[1]-b[18]*a[2]-b[20]*a[3]-b[21]*a[4]-b[1]*a[5]+b[12]*a[6]+b[14]*a[7]+b[15]*a[8]+b[0]*a[9]-b[27]*a[10]-b[28]*a[11]-b[6]*a[12]-b[29]*a[13]-b[7]*a[14]-b[8]*a[15]-b[23]*a[16]-b[24]*a[17]-b[2]*a[18]-b[25]*a[19]-b[3]*a[20]-b[4]*a[21]+b[31]*a[22]+b[16]*a[23]+b[17]*a[24]+b[19]*a[25]-b[30]*a[26]-b[10]*a[27]-b[11]*a[28]-b[13]*a[29]+b[26]*a[30]+b[22]*a[31]; res[10]=b[10]*a[0]+b[3]*a[2]-b[2]*a[3]+b[22]*a[4]+b[23]*a[5]+b[0]*a[10]-b[13]*a[11]-b[14]*a[12]+b[11]*a[13]+b[12]*a[14]-b[30]*a[15]+b[4]*a[22]+b[5]*a[23]-b[25]*a[24]+b[24]*a[25]-b[15]*a[30]; res[11]=b[11]*a[0]+b[4]*a[2]-b[22]*a[3]-b[2]*a[4]+b[24]*a[5]+b[13]*a[10]+b[0]*a[11]-b[15]*a[12]-b[10]*a[13]+b[30]*a[14]+b[12]*a[15]-b[3]*a[22]+b[25]*a[23]+b[5]*a[24]-b[23]*a[25]+b[14]*a[30]; res[12]=b[12]*a[0]+b[5]*a[2]-b[23]*a[3]-b[24]*a[4]-b[2]*a[5]+b[14]*a[10]+b[15]*a[11]+b[0]*a[12]-b[30]*a[13]-b[10]*a[14]-b[11]*a[15]-b[25]*a[22]-b[3]*a[23]-b[4]*a[24]+b[22]*a[25]-b[13]*a[30]; res[13]=b[13]*a[0]+b[22]*a[2]+b[4]*a[3]-b[3]*a[4]+b[25]*a[5]-b[11]*a[10]+b[10]*a[11]-b[30]*a[12]+b[0]*a[13]-b[15]*a[14]+b[14]*a[15]+b[2]*a[22]-b[24]*a[23]+b[23]*a[24]+b[5]*a[25]-b[12]*a[30]; res[14]=b[14]*a[0]+b[23]*a[2]+b[5]*a[3]-b[25]*a[4]-b[3]*a[5]-b[12]*a[10]+b[30]*a[11]+b[10]*a[12]+b[15]*a[13]+b[0]*a[14]-b[13]*a[15]+b[24]*a[22]+b[2]*a[23]-b[22]*a[24]-b[4]*a[25]+b[11]*a[30]; res[15]=b[15]*a[0]+b[24]*a[2]+b[25]*a[3]+b[5]*a[4]-b[4]*a[5]-b[30]*a[10]-b[12]*a[11]+b[11]*a[12]-b[14]*a[13]+b[13]*a[14]+b[0]*a[15]-b[23]*a[22]+b[22]*a[23]+b[2]*a[24]+b[3]*a[25]-b[10]*a[30]; res[16]=b[16]*a[0]+b[10]*a[1]-b[7]*a[2]+b[6]*a[3]-b[26]*a[4]-b[27]*a[5]+b[3]*a[6]-b[2]*a[7]+b[22]*a[8]+b[23]*a[9]+b[1]*a[10]-b[19]*a[11]-b[20]*a[12]+b[17]*a[13]+b[18]*a[14]-b[31]*a[15]+b[0]*a[16]-b[13]*a[17]-b[14]*a[18]+b[11]*a[19]+b[12]*a[20]-b[30]*a[21]-b[8]*a[22]-b[9]*a[23]+b[29]*a[24]-b[28]*a[25]+b[4]*a[26]+b[5]*a[27]-b[25]*a[28]+b[24]*a[29]-b[21]*a[30]-b[15]*a[31]; res[17]=b[17]*a[0]+b[11]*a[1]-b[8]*a[2]+b[26]*a[3]+b[6]*a[4]-b[28]*a[5]+b[4]*a[6]-b[22]*a[7]-b[2]*a[8]+b[24]*a[9]+b[19]*a[10]+b[1]*a[11]-b[21]*a[12]-b[16]*a[13]+b[31]*a[14]+b[18]*a[15]+b[13]*a[16]+b[0]*a[17]-b[15]*a[18]-b[10]*a[19]+b[30]*a[20]+b[12]*a[21]+b[7]*a[22]-b[29]*a[23]-b[9]*a[24]+b[27]*a[25]-b[3]*a[26]+b[25]*a[27]+b[5]*a[28]-b[23]*a[29]+b[20]*a[30]+b[14]*a[31]; res[18]=b[18]*a[0]+b[12]*a[1]-b[9]*a[2]+b[27]*a[3]+b[28]*a[4]+b[6]*a[5]+b[5]*a[6]-b[23]*a[7]-b[24]*a[8]-b[2]*a[9]+b[20]*a[10]+b[21]*a[11]+b[1]*a[12]-b[31]*a[13]-b[16]*a[14]-b[17]*a[15]+b[14]*a[16]+b[15]*a[17]+b[0]*a[18]-b[30]*a[19]-b[10]*a[20]-b[11]*a[21]+b[29]*a[22]+b[7]*a[23]+b[8]*a[24]-b[26]*a[25]-b[25]*a[26]-b[3]*a[27]-b[4]*a[28]+b[22]*a[29]-b[19]*a[30]-b[13]*a[31]; res[19]=b[19]*a[0]+b[13]*a[1]-b[26]*a[2]-b[8]*a[3]+b[7]*a[4]-b[29]*a[5]+b[22]*a[6]+b[4]*a[7]-b[3]*a[8]+b[25]*a[9]-b[17]*a[10]+b[16]*a[11]-b[31]*a[12]+b[1]*a[13]-b[21]*a[14]+b[20]*a[15]-b[11]*a[16]+b[10]*a[17]-b[30]*a[18]+b[0]*a[19]-b[15]*a[20]+b[14]*a[21]-b[6]*a[22]+b[28]*a[23]-b[27]*a[24]-b[9]*a[25]+b[2]*a[26]-b[24]*a[27]+b[23]*a[28]+b[5]*a[29]-b[18]*a[30]-b[12]*a[31]; res[20]=b[20]*a[0]+b[14]*a[1]-b[27]*a[2]-b[9]*a[3]+b[29]*a[4]+b[7]*a[5]+b[23]*a[6]+b[5]*a[7]-b[25]*a[8]-b[3]*a[9]-b[18]*a[10]+b[31]*a[11]+b[16]*a[12]+b[21]*a[13]+b[1]*a[14]-b[19]*a[15]-b[12]*a[16]+b[30]*a[17]+b[10]*a[18]+b[15]*a[19]+b[0]*a[20]-b[13]*a[21]-b[28]*a[22]-b[6]*a[23]+b[26]*a[24]+b[8]*a[25]+b[24]*a[26]+b[2]*a[27]-b[22]*a[28]-b[4]*a[29]+b[17]*a[30]+b[11]*a[31]; res[21]=b[21]*a[0]+b[15]*a[1]-b[28]*a[2]-b[29]*a[3]-b[9]*a[4]+b[8]*a[5]+b[24]*a[6]+b[25]*a[7]+b[5]*a[8]-b[4]*a[9]-b[31]*a[10]-b[18]*a[11]+b[17]*a[12]-b[20]*a[13]+b[19]*a[14]+b[1]*a[15]-b[30]*a[16]-b[12]*a[17]+b[11]*a[18]-b[14]*a[19]+b[13]*a[20]+b[0]*a[21]+b[27]*a[22]-b[26]*a[23]-b[6]*a[24]-b[7]*a[25]-b[23]*a[26]+b[22]*a[27]+b[2]*a[28]+b[3]*a[29]-b[16]*a[30]-b[10]*a[31]; res[22]=b[22]*a[0]+b[13]*a[2]-b[11]*a[3]+b[10]*a[4]-b[30]*a[5]+b[4]*a[10]-b[3]*a[11]+b[25]*a[12]+b[2]*a[13]-b[24]*a[14]+b[23]*a[15]+b[0]*a[22]-b[15]*a[23]+b[14]*a[24]-b[12]*a[25]+b[5]*a[30]; res[23]=b[23]*a[0]+b[14]*a[2]-b[12]*a[3]+b[30]*a[4]+b[10]*a[5]+b[5]*a[10]-b[25]*a[11]-b[3]*a[12]+b[24]*a[13]+b[2]*a[14]-b[22]*a[15]+b[15]*a[22]+b[0]*a[23]-b[13]*a[24]+b[11]*a[25]-b[4]*a[30]; res[24]=b[24]*a[0]+b[15]*a[2]-b[30]*a[3]-b[12]*a[4]+b[11]*a[5]+b[25]*a[10]+b[5]*a[11]-b[4]*a[12]-b[23]*a[13]+b[22]*a[14]+b[2]*a[15]-b[14]*a[22]+b[13]*a[23]+b[0]*a[24]-b[10]*a[25]+b[3]*a[30]; res[25]=b[25]*a[0]+b[30]*a[2]+b[15]*a[3]-b[14]*a[4]+b[13]*a[5]-b[24]*a[10]+b[23]*a[11]-b[22]*a[12]+b[5]*a[13]-b[4]*a[14]+b[3]*a[15]+b[12]*a[22]-b[11]*a[23]+b[10]*a[24]+b[0]*a[25]-b[2]*a[30]; res[26]=b[26]*a[0]+b[22]*a[1]-b[19]*a[2]+b[17]*a[3]-b[16]*a[4]+b[31]*a[5]+b[13]*a[6]-b[11]*a[7]+b[10]*a[8]-b[30]*a[9]+b[8]*a[10]-b[7]*a[11]+b[29]*a[12]+b[6]*a[13]-b[28]*a[14]+b[27]*a[15]+b[4]*a[16]-b[3]*a[17]+b[25]*a[18]+b[2]*a[19]-b[24]*a[20]+b[23]*a[21]-b[1]*a[22]+b[21]*a[23]-b[20]*a[24]+b[18]*a[25]+b[0]*a[26]-b[15]*a[27]+b[14]*a[28]-b[12]*a[29]+b[9]*a[30]+b[5]*a[31]; res[27]=b[27]*a[0]+b[23]*a[1]-b[20]*a[2]+b[18]*a[3]-b[31]*a[4]-b[16]*a[5]+b[14]*a[6]-b[12]*a[7]+b[30]*a[8]+b[10]*a[9]+b[9]*a[10]-b[29]*a[11]-b[7]*a[12]+b[28]*a[13]+b[6]*a[14]-b[26]*a[15]+b[5]*a[16]-b[25]*a[17]-b[3]*a[18]+b[24]*a[19]+b[2]*a[20]-b[22]*a[21]-b[21]*a[22]-b[1]*a[23]+b[19]*a[24]-b[17]*a[25]+b[15]*a[26]+b[0]*a[27]-b[13]*a[28]+b[11]*a[29]-b[8]*a[30]-b[4]*a[31]; res[28]=b[28]*a[0]+b[24]*a[1]-b[21]*a[2]+b[31]*a[3]+b[18]*a[4]-b[17]*a[5]+b[15]*a[6]-b[30]*a[7]-b[12]*a[8]+b[11]*a[9]+b[29]*a[10]+b[9]*a[11]-b[8]*a[12]-b[27]*a[13]+b[26]*a[14]+b[6]*a[15]+b[25]*a[16]+b[5]*a[17]-b[4]*a[18]-b[23]*a[19]+b[22]*a[20]+b[2]*a[21]+b[20]*a[22]-b[19]*a[23]-b[1]*a[24]+b[16]*a[25]-b[14]*a[26]+b[13]*a[27]+b[0]*a[28]-b[10]*a[29]+b[7]*a[30]+b[3]*a[31]; res[29]=b[29]*a[0]+b[25]*a[1]-b[31]*a[2]-b[21]*a[3]+b[20]*a[4]-b[19]*a[5]+b[30]*a[6]+b[15]*a[7]-b[14]*a[8]+b[13]*a[9]-b[28]*a[10]+b[27]*a[11]-b[26]*a[12]+b[9]*a[13]-b[8]*a[14]+b[7]*a[15]-b[24]*a[16]+b[23]*a[17]-b[22]*a[18]+b[5]*a[19]-b[4]*a[20]+b[3]*a[21]-b[18]*a[22]+b[17]*a[23]-b[16]*a[24]-b[1]*a[25]+b[12]*a[26]-b[11]*a[27]+b[10]*a[28]+b[0]*a[29]-b[6]*a[30]-b[2]*a[31]; res[30]=b[30]*a[0]+b[25]*a[2]-b[24]*a[3]+b[23]*a[4]-b[22]*a[5]+b[15]*a[10]-b[14]*a[11]+b[13]*a[12]+b[12]*a[13]-b[11]*a[14]+b[10]*a[15]+b[5]*a[22]-b[4]*a[23]+b[3]*a[24]-b[2]*a[25]+b[0]*a[30]; res[31]=b[31]*a[0]+b[30]*a[1]-b[29]*a[2]+b[28]*a[3]-b[27]*a[4]+b[26]*a[5]+b[25]*a[6]-b[24]*a[7]+b[23]*a[8]-b[22]*a[9]+b[21]*a[10]-b[20]*a[11]+b[19]*a[12]+b[18]*a[13]-b[17]*a[14]+b[16]*a[15]+b[15]*a[16]-b[14]*a[17]+b[13]*a[18]+b[12]*a[19]-b[11]*a[20]+b[10]*a[21]-b[9]*a[22]+b[8]*a[23]-b[7]*a[24]+b[6]*a[25]+b[5]*a[26]-b[4]*a[27]+b[3]*a[28]-b[2]*a[29]+b[1]*a[30]+b[0]*a[31]; res } } // Wedge // The outer product. (MEET) impl BitXor for R401 { type Output = R401; fn bitxor(self: R401, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0]=b[0]*a[0]; res[1]=b[1]*a[0]+b[0]*a[1]; res[2]=b[2]*a[0]+b[0]*a[2]; res[3]=b[3]*a[0]+b[0]*a[3]; res[4]=b[4]*a[0]+b[0]*a[4]; res[5]=b[5]*a[0]+b[0]*a[5]; res[6]=b[6]*a[0]+b[2]*a[1]-b[1]*a[2]+b[0]*a[6]; res[7]=b[7]*a[0]+b[3]*a[1]-b[1]*a[3]+b[0]*a[7]; res[8]=b[8]*a[0]+b[4]*a[1]-b[1]*a[4]+b[0]*a[8]; res[9]=b[9]*a[0]+b[5]*a[1]-b[1]*a[5]+b[0]*a[9]; res[10]=b[10]*a[0]+b[3]*a[2]-b[2]*a[3]+b[0]*a[10]; res[11]=b[11]*a[0]+b[4]*a[2]-b[2]*a[4]+b[0]*a[11]; res[12]=b[12]*a[0]+b[5]*a[2]-b[2]*a[5]+b[0]*a[12]; res[13]=b[13]*a[0]+b[4]*a[3]-b[3]*a[4]+b[0]*a[13]; res[14]=b[14]*a[0]+b[5]*a[3]-b[3]*a[5]+b[0]*a[14]; res[15]=b[15]*a[0]+b[5]*a[4]-b[4]*a[5]+b[0]*a[15]; res[16]=b[16]*a[0]+b[10]*a[1]-b[7]*a[2]+b[6]*a[3]+b[3]*a[6]-b[2]*a[7]+b[1]*a[10]+b[0]*a[16]; res[17]=b[17]*a[0]+b[11]*a[1]-b[8]*a[2]+b[6]*a[4]+b[4]*a[6]-b[2]*a[8]+b[1]*a[11]+b[0]*a[17]; res[18]=b[18]*a[0]+b[12]*a[1]-b[9]*a[2]+b[6]*a[5]+b[5]*a[6]-b[2]*a[9]+b[1]*a[12]+b[0]*a[18]; res[19]=b[19]*a[0]+b[13]*a[1]-b[8]*a[3]+b[7]*a[4]+b[4]*a[7]-b[3]*a[8]+b[1]*a[13]+b[0]*a[19]; res[20]=b[20]*a[0]+b[14]*a[1]-b[9]*a[3]+b[7]*a[5]+b[5]*a[7]-b[3]*a[9]+b[1]*a[14]+b[0]*a[20]; res[21]=b[21]*a[0]+b[15]*a[1]-b[9]*a[4]+b[8]*a[5]+b[5]*a[8]-b[4]*a[9]+b[1]*a[15]+b[0]*a[21]; res[22]=b[22]*a[0]+b[13]*a[2]-b[11]*a[3]+b[10]*a[4]+b[4]*a[10]-b[3]*a[11]+b[2]*a[13]+b[0]*a[22]; res[23]=b[23]*a[0]+b[14]*a[2]-b[12]*a[3]+b[10]*a[5]+b[5]*a[10]-b[3]*a[12]+b[2]*a[14]+b[0]*a[23]; res[24]=b[24]*a[0]+b[15]*a[2]-b[12]*a[4]+b[11]*a[5]+b[5]*a[11]-b[4]*a[12]+b[2]*a[15]+b[0]*a[24]; res[25]=b[25]*a[0]+b[15]*a[3]-b[14]*a[4]+b[13]*a[5]+b[5]*a[13]-b[4]*a[14]+b[3]*a[15]+b[0]*a[25]; res[26]=b[26]*a[0]+b[22]*a[1]-b[19]*a[2]+b[17]*a[3]-b[16]*a[4]+b[13]*a[6]-b[11]*a[7]+b[10]*a[8]+b[8]*a[10]-b[7]*a[11]+b[6]*a[13]+b[4]*a[16]-b[3]*a[17]+b[2]*a[19]-b[1]*a[22]+b[0]*a[26]; res[27]=b[27]*a[0]+b[23]*a[1]-b[20]*a[2]+b[18]*a[3]-b[16]*a[5]+b[14]*a[6]-b[12]*a[7]+b[10]*a[9]+b[9]*a[10]-b[7]*a[12]+b[6]*a[14]+b[5]*a[16]-b[3]*a[18]+b[2]*a[20]-b[1]*a[23]+b[0]*a[27]; res[28]=b[28]*a[0]+b[24]*a[1]-b[21]*a[2]+b[18]*a[4]-b[17]*a[5]+b[15]*a[6]-b[12]*a[8]+b[11]*a[9]+b[9]*a[11]-b[8]*a[12]+b[6]*a[15]+b[5]*a[17]-b[4]*a[18]+b[2]*a[21]-b[1]*a[24]+b[0]*a[28]; res[29]=b[29]*a[0]+b[25]*a[1]-b[21]*a[3]+b[20]*a[4]-b[19]*a[5]+b[15]*a[7]-b[14]*a[8]+b[13]*a[9]+b[9]*a[13]-b[8]*a[14]+b[7]*a[15]+b[5]*a[19]-b[4]*a[20]+b[3]*a[21]-b[1]*a[25]+b[0]*a[29]; res[30]=b[30]*a[0]+b[25]*a[2]-b[24]*a[3]+b[23]*a[4]-b[22]*a[5]+b[15]*a[10]-b[14]*a[11]+b[13]*a[12]+b[12]*a[13]-b[11]*a[14]+b[10]*a[15]+b[5]*a[22]-b[4]*a[23]+b[3]*a[24]-b[2]*a[25]+b[0]*a[30]; res[31]=b[31]*a[0]+b[30]*a[1]-b[29]*a[2]+b[28]*a[3]-b[27]*a[4]+b[26]*a[5]+b[25]*a[6]-b[24]*a[7]+b[23]*a[8]-b[22]*a[9]+b[21]*a[10]-b[20]*a[11]+b[19]*a[12]+b[18]*a[13]-b[17]*a[14]+b[16]*a[15]+b[15]*a[16]-b[14]*a[17]+b[13]*a[18]+b[12]*a[19]-b[11]*a[20]+b[10]*a[21]-b[9]*a[22]+b[8]*a[23]-b[7]*a[24]+b[6]*a[25]+b[5]*a[26]-b[4]*a[27]+b[3]*a[28]-b[2]*a[29]+b[1]*a[30]+b[0]*a[31]; res } } // Vee // The regressive product. (JOIN) impl BitAnd for R401 { type Output = R401; fn bitand(self: R401, b: R401) -> R401 { // return (self.dual_h() ^ b.dual_h()).undual_h(); let mut res = R401::zero(); let a = self; res[31]=1.0*(a[31]*b[31]); res[30]=1.0*(a[30]*b[31]+a[31]*b[30]); res[29]=-1.0*(a[29]*-1.0*b[31]+a[31]*b[29]*-1.0); res[28]=1.0*(a[28]*b[31]+a[31]*b[28]); res[27]=-1.0*(a[27]*-1.0*b[31]+a[31]*b[27]*-1.0); res[26]=1.0*(a[26]*b[31]+a[31]*b[26]); res[25]=1.0*(a[25]*b[31]+a[29]*-1.0*b[30]-a[30]*b[29]*-1.0+a[31]*b[25]); res[24]=-1.0*(a[24]*-1.0*b[31]+a[28]*b[30]-a[30]*b[28]+a[31]*b[24]*-1.0); res[23]=1.0*(a[23]*b[31]+a[27]*-1.0*b[30]-a[30]*b[27]*-1.0+a[31]*b[23]); res[22]=-1.0*(a[22]*-1.0*b[31]+a[26]*b[30]-a[30]*b[26]+a[31]*b[22]*-1.0); res[21]=1.0*(a[21]*b[31]+a[28]*b[29]*-1.0-a[29]*-1.0*b[28]+a[31]*b[21]); res[20]=-1.0*(a[20]*-1.0*b[31]+a[27]*-1.0*b[29]*-1.0-a[29]*-1.0*b[27]*-1.0+a[31]*b[20]*-1.0); res[19]=1.0*(a[19]*b[31]+a[26]*b[29]*-1.0-a[29]*-1.0*b[26]+a[31]*b[19]); res[18]=1.0*(a[18]*b[31]+a[27]*-1.0*b[28]-a[28]*b[27]*-1.0+a[31]*b[18]); res[17]=-1.0*(a[17]*-1.0*b[31]+a[26]*b[28]-a[28]*b[26]+a[31]*b[17]*-1.0); res[16]=1.0*(a[16]*b[31]+a[26]*b[27]*-1.0-a[27]*-1.0*b[26]+a[31]*b[16]); res[15]=1.0*(a[15]*b[31]+a[21]*b[30]-a[24]*-1.0*b[29]*-1.0+a[25]*b[28]+a[28]*b[25]-a[29]*-1.0*b[24]*-1.0+a[30]*b[21]+a[31]*b[15]); res[14]=-1.0*(a[14]*-1.0*b[31]+a[20]*-1.0*b[30]-a[23]*b[29]*-1.0+a[25]*b[27]*-1.0+a[27]*-1.0*b[25]-a[29]*-1.0*b[23]+a[30]*b[20]*-1.0+a[31]*b[14]*-1.0); res[13]=1.0*(a[13]*b[31]+a[19]*b[30]-a[22]*-1.0*b[29]*-1.0+a[25]*b[26]+a[26]*b[25]-a[29]*-1.0*b[22]*-1.0+a[30]*b[19]+a[31]*b[13]); res[12]=1.0*(a[12]*b[31]+a[18]*b[30]-a[23]*b[28]+a[24]*-1.0*b[27]*-1.0+a[27]*-1.0*b[24]*-1.0-a[28]*b[23]+a[30]*b[18]+a[31]*b[12]); res[11]=-1.0*(a[11]*-1.0*b[31]+a[17]*-1.0*b[30]-a[22]*-1.0*b[28]+a[24]*-1.0*b[26]+a[26]*b[24]*-1.0-a[28]*b[22]*-1.0+a[30]*b[17]*-1.0+a[31]*b[11]*-1.0); res[10]=1.0*(a[10]*b[31]+a[16]*b[30]-a[22]*-1.0*b[27]*-1.0+a[23]*b[26]+a[26]*b[23]-a[27]*-1.0*b[22]*-1.0+a[30]*b[16]+a[31]*b[10]); res[9]=-1.0*(a[9]*-1.0*b[31]+a[18]*b[29]*-1.0-a[20]*-1.0*b[28]+a[21]*b[27]*-1.0+a[27]*-1.0*b[21]-a[28]*b[20]*-1.0+a[29]*-1.0*b[18]+a[31]*b[9]*-1.0); res[8]=1.0*(a[8]*b[31]+a[17]*-1.0*b[29]*-1.0-a[19]*b[28]+a[21]*b[26]+a[26]*b[21]-a[28]*b[19]+a[29]*-1.0*b[17]*-1.0+a[31]*b[8]); res[7]=-1.0*(a[7]*-1.0*b[31]+a[16]*b[29]*-1.0-a[19]*b[27]*-1.0+a[20]*-1.0*b[26]+a[26]*b[20]*-1.0-a[27]*-1.0*b[19]+a[29]*-1.0*b[16]+a[31]*b[7]*-1.0); res[6]=1.0*(a[6]*b[31]+a[16]*b[28]-a[17]*-1.0*b[27]*-1.0+a[18]*b[26]+a[26]*b[18]-a[27]*-1.0*b[17]*-1.0+a[28]*b[16]+a[31]*b[6]); res[5]=1.0*(a[5]*b[31]+a[9]*-1.0*b[30]-a[12]*b[29]*-1.0+a[14]*-1.0*b[28]-a[15]*b[27]*-1.0+a[18]*b[25]-a[20]*-1.0*b[24]*-1.0+a[21]*b[23]+a[23]*b[21]-a[24]*-1.0*b[20]*-1.0+a[25]*b[18]+a[27]*-1.0*b[15]-a[28]*b[14]*-1.0+a[29]*-1.0*b[12]-a[30]*b[9]*-1.0+a[31]*b[5]); res[4]=-1.0*(a[4]*-1.0*b[31]+a[8]*b[30]-a[11]*-1.0*b[29]*-1.0+a[13]*b[28]-a[15]*b[26]+a[17]*-1.0*b[25]-a[19]*b[24]*-1.0+a[21]*b[22]*-1.0+a[22]*-1.0*b[21]-a[24]*-1.0*b[19]+a[25]*b[17]*-1.0+a[26]*b[15]-a[28]*b[13]+a[29]*-1.0*b[11]*-1.0-a[30]*b[8]+a[31]*b[4]*-1.0); res[3]=1.0*(a[3]*b[31]+a[7]*-1.0*b[30]-a[10]*b[29]*-1.0+a[13]*b[27]*-1.0-a[14]*-1.0*b[26]+a[16]*b[25]-a[19]*b[23]+a[20]*-1.0*b[22]*-1.0+a[22]*-1.0*b[20]*-1.0-a[23]*b[19]+a[25]*b[16]+a[26]*b[14]*-1.0-a[27]*-1.0*b[13]+a[29]*-1.0*b[10]-a[30]*b[7]*-1.0+a[31]*b[3]); res[2]=-1.0*(a[2]*-1.0*b[31]+a[6]*b[30]-a[10]*b[28]+a[11]*-1.0*b[27]*-1.0-a[12]*b[26]+a[16]*b[24]*-1.0-a[17]*-1.0*b[23]+a[18]*b[22]*-1.0+a[22]*-1.0*b[18]-a[23]*b[17]*-1.0+a[24]*-1.0*b[16]+a[26]*b[12]-a[27]*-1.0*b[11]*-1.0+a[28]*b[10]-a[30]*b[6]+a[31]*b[2]*-1.0); res[1]=1.0*(a[1]*b[31]+a[6]*b[29]*-1.0-a[7]*-1.0*b[28]+a[8]*b[27]*-1.0-a[9]*-1.0*b[26]+a[16]*b[21]-a[17]*-1.0*b[20]*-1.0+a[18]*b[19]+a[19]*b[18]-a[20]*-1.0*b[17]*-1.0+a[21]*b[16]+a[26]*b[9]*-1.0-a[27]*-1.0*b[8]+a[28]*b[7]*-1.0-a[29]*-1.0*b[6]+a[31]*b[1]); res[0]=1.0*(a[0]*b[31]+a[1]*b[30]-a[2]*-1.0*b[29]*-1.0+a[3]*b[28]-a[4]*-1.0*b[27]*-1.0+a[5]*b[26]+a[6]*b[25]-a[7]*-1.0*b[24]*-1.0+a[8]*b[23]-a[9]*-1.0*b[22]*-1.0+a[10]*b[21]-a[11]*-1.0*b[20]*-1.0+a[12]*b[19]+a[13]*b[18]-a[14]*-1.0*b[17]*-1.0+a[15]*b[16]+a[16]*b[15]-a[17]*-1.0*b[14]*-1.0+a[18]*b[13]+a[19]*b[12]-a[20]*-1.0*b[11]*-1.0+a[21]*b[10]-a[22]*-1.0*b[9]*-1.0+a[23]*b[8]-a[24]*-1.0*b[7]*-1.0+a[25]*b[6]+a[26]*b[5]-a[27]*-1.0*b[4]*-1.0+a[28]*b[3]-a[29]*-1.0*b[2]*-1.0+a[30]*b[1]+a[31]*b[0]); res } } // Dot // The inner product. impl BitOr for R401 { type Output = R401; fn bitor(self: R401, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0]=b[0]*a[0]+b[2]*a[2]+b[3]*a[3]+b[4]*a[4]+b[5]*a[5]-b[10]*a[10]-b[11]*a[11]-b[12]*a[12]-b[13]*a[13]-b[14]*a[14]-b[15]*a[15]-b[22]*a[22]-b[23]*a[23]-b[24]*a[24]-b[25]*a[25]+b[30]*a[30]; res[1]=b[1]*a[0]+b[0]*a[1]-b[6]*a[2]-b[7]*a[3]-b[8]*a[4]-b[9]*a[5]+b[2]*a[6]+b[3]*a[7]+b[4]*a[8]+b[5]*a[9]-b[16]*a[10]-b[17]*a[11]-b[18]*a[12]-b[19]*a[13]-b[20]*a[14]-b[21]*a[15]-b[10]*a[16]-b[11]*a[17]-b[12]*a[18]-b[13]*a[19]-b[14]*a[20]-b[15]*a[21]+b[26]*a[22]+b[27]*a[23]+b[28]*a[24]+b[29]*a[25]-b[22]*a[26]-b[23]*a[27]-b[24]*a[28]-b[25]*a[29]+b[31]*a[30]+b[30]*a[31]; res[2]=b[2]*a[0]+b[0]*a[2]-b[10]*a[3]-b[11]*a[4]-b[12]*a[5]+b[3]*a[10]+b[4]*a[11]+b[5]*a[12]-b[22]*a[13]-b[23]*a[14]-b[24]*a[15]-b[13]*a[22]-b[14]*a[23]-b[15]*a[24]+b[30]*a[25]-b[25]*a[30]; res[3]=b[3]*a[0]+b[10]*a[2]+b[0]*a[3]-b[13]*a[4]-b[14]*a[5]-b[2]*a[10]+b[22]*a[11]+b[23]*a[12]+b[4]*a[13]+b[5]*a[14]-b[25]*a[15]+b[11]*a[22]+b[12]*a[23]-b[30]*a[24]-b[15]*a[25]+b[24]*a[30]; res[4]=b[4]*a[0]+b[11]*a[2]+b[13]*a[3]+b[0]*a[4]-b[15]*a[5]-b[22]*a[10]-b[2]*a[11]+b[24]*a[12]-b[3]*a[13]+b[25]*a[14]+b[5]*a[15]-b[10]*a[22]+b[30]*a[23]+b[12]*a[24]+b[14]*a[25]-b[23]*a[30]; res[5]=b[5]*a[0]+b[12]*a[2]+b[14]*a[3]+b[15]*a[4]+b[0]*a[5]-b[23]*a[10]-b[24]*a[11]-b[2]*a[12]-b[25]*a[13]-b[3]*a[14]-b[4]*a[15]-b[30]*a[22]-b[10]*a[23]-b[11]*a[24]-b[13]*a[25]+b[22]*a[30]; res[6]=b[6]*a[0]+b[16]*a[3]+b[17]*a[4]+b[18]*a[5]+b[0]*a[6]-b[26]*a[13]-b[27]*a[14]-b[28]*a[15]+b[3]*a[16]+b[4]*a[17]+b[5]*a[18]-b[31]*a[25]-b[13]*a[26]-b[14]*a[27]-b[15]*a[28]-b[25]*a[31]; res[7]=b[7]*a[0]-b[16]*a[2]+b[19]*a[4]+b[20]*a[5]+b[0]*a[7]+b[26]*a[11]+b[27]*a[12]-b[29]*a[15]-b[2]*a[16]+b[4]*a[19]+b[5]*a[20]+b[31]*a[24]+b[11]*a[26]+b[12]*a[27]-b[15]*a[29]+b[24]*a[31]; res[8]=b[8]*a[0]-b[17]*a[2]-b[19]*a[3]+b[21]*a[5]+b[0]*a[8]-b[26]*a[10]+b[28]*a[12]+b[29]*a[14]-b[2]*a[17]-b[3]*a[19]+b[5]*a[21]-b[31]*a[23]-b[10]*a[26]+b[12]*a[28]+b[14]*a[29]-b[23]*a[31]; res[9]=b[9]*a[0]-b[18]*a[2]-b[20]*a[3]-b[21]*a[4]+b[0]*a[9]-b[27]*a[10]-b[28]*a[11]-b[29]*a[13]-b[2]*a[18]-b[3]*a[20]-b[4]*a[21]+b[31]*a[22]-b[10]*a[27]-b[11]*a[28]-b[13]*a[29]+b[22]*a[31]; res[10]=b[10]*a[0]+b[22]*a[4]+b[23]*a[5]+b[0]*a[10]-b[30]*a[15]+b[4]*a[22]+b[5]*a[23]-b[15]*a[30]; res[11]=b[11]*a[0]-b[22]*a[3]+b[24]*a[5]+b[0]*a[11]+b[30]*a[14]-b[3]*a[22]+b[5]*a[24]+b[14]*a[30]; res[12]=b[12]*a[0]-b[23]*a[3]-b[24]*a[4]+b[0]*a[12]-b[30]*a[13]-b[3]*a[23]-b[4]*a[24]-b[13]*a[30]; res[13]=b[13]*a[0]+b[22]*a[2]+b[25]*a[5]-b[30]*a[12]+b[0]*a[13]+b[2]*a[22]+b[5]*a[25]-b[12]*a[30]; res[14]=b[14]*a[0]+b[23]*a[2]-b[25]*a[4]+b[30]*a[11]+b[0]*a[14]+b[2]*a[23]-b[4]*a[25]+b[11]*a[30]; res[15]=b[15]*a[0]+b[24]*a[2]+b[25]*a[3]-b[30]*a[10]+b[0]*a[15]+b[2]*a[24]+b[3]*a[25]-b[10]*a[30]; res[16]=b[16]*a[0]-b[26]*a[4]-b[27]*a[5]-b[31]*a[15]+b[0]*a[16]+b[4]*a[26]+b[5]*a[27]-b[15]*a[31]; res[17]=b[17]*a[0]+b[26]*a[3]-b[28]*a[5]+b[31]*a[14]+b[0]*a[17]-b[3]*a[26]+b[5]*a[28]+b[14]*a[31]; res[18]=b[18]*a[0]+b[27]*a[3]+b[28]*a[4]-b[31]*a[13]+b[0]*a[18]-b[3]*a[27]-b[4]*a[28]-b[13]*a[31]; res[19]=b[19]*a[0]-b[26]*a[2]-b[29]*a[5]-b[31]*a[12]+b[0]*a[19]+b[2]*a[26]+b[5]*a[29]-b[12]*a[31]; res[20]=b[20]*a[0]-b[27]*a[2]+b[29]*a[4]+b[31]*a[11]+b[0]*a[20]+b[2]*a[27]-b[4]*a[29]+b[11]*a[31]; res[21]=b[21]*a[0]-b[28]*a[2]-b[29]*a[3]-b[31]*a[10]+b[0]*a[21]+b[2]*a[28]+b[3]*a[29]-b[10]*a[31]; res[22]=b[22]*a[0]-b[30]*a[5]+b[0]*a[22]+b[5]*a[30]; res[23]=b[23]*a[0]+b[30]*a[4]+b[0]*a[23]-b[4]*a[30]; res[24]=b[24]*a[0]-b[30]*a[3]+b[0]*a[24]+b[3]*a[30]; res[25]=b[25]*a[0]+b[30]*a[2]+b[0]*a[25]-b[2]*a[30]; res[26]=b[26]*a[0]+b[31]*a[5]+b[0]*a[26]+b[5]*a[31]; res[27]=b[27]*a[0]-b[31]*a[4]+b[0]*a[27]-b[4]*a[31]; res[28]=b[28]*a[0]+b[31]*a[3]+b[0]*a[28]+b[3]*a[31]; res[29]=b[29]*a[0]-b[31]*a[2]+b[0]*a[29]-b[2]*a[31]; res[30]=b[30]*a[0]+b[0]*a[30]; res[31]=b[31]*a[0]+b[0]*a[31]; res } } // Add // Multivector addition impl Add for R401 { type Output = R401; fn add(self: R401, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]+b[0]; res[1] = a[1]+b[1]; res[2] = a[2]+b[2]; res[3] = a[3]+b[3]; res[4] = a[4]+b[4]; res[5] = a[5]+b[5]; res[6] = a[6]+b[6]; res[7] = a[7]+b[7]; res[8] = a[8]+b[8]; res[9] = a[9]+b[9]; res[10] = a[10]+b[10]; res[11] = a[11]+b[11]; res[12] = a[12]+b[12]; res[13] = a[13]+b[13]; res[14] = a[14]+b[14]; res[15] = a[15]+b[15]; res[16] = a[16]+b[16]; res[17] = a[17]+b[17]; res[18] = a[18]+b[18]; res[19] = a[19]+b[19]; res[20] = a[20]+b[20]; res[21] = a[21]+b[21]; res[22] = a[22]+b[22]; res[23] = a[23]+b[23]; res[24] = a[24]+b[24]; res[25] = a[25]+b[25]; res[26] = a[26]+b[26]; res[27] = a[27]+b[27]; res[28] = a[28]+b[28]; res[29] = a[29]+b[29]; res[30] = a[30]+b[30]; res[31] = a[31]+b[31]; res } } // Sub // Multivector subtraction impl Sub for R401 { type Output = R401; fn sub(self: R401, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]-b[0]; res[1] = a[1]-b[1]; res[2] = a[2]-b[2]; res[3] = a[3]-b[3]; res[4] = a[4]-b[4]; res[5] = a[5]-b[5]; res[6] = a[6]-b[6]; res[7] = a[7]-b[7]; res[8] = a[8]-b[8]; res[9] = a[9]-b[9]; res[10] = a[10]-b[10]; res[11] = a[11]-b[11]; res[12] = a[12]-b[12]; res[13] = a[13]-b[13]; res[14] = a[14]-b[14]; res[15] = a[15]-b[15]; res[16] = a[16]-b[16]; res[17] = a[17]-b[17]; res[18] = a[18]-b[18]; res[19] = a[19]-b[19]; res[20] = a[20]-b[20]; res[21] = a[21]-b[21]; res[22] = a[22]-b[22]; res[23] = a[23]-b[23]; res[24] = a[24]-b[24]; res[25] = a[25]-b[25]; res[26] = a[26]-b[26]; res[27] = a[27]-b[27]; res[28] = a[28]-b[28]; res[29] = a[29]-b[29]; res[30] = a[30]-b[30]; res[31] = a[31]-b[31]; res } } // smul // scalar/multivector multiplication impl Mul<R401> for float_t { type Output = R401; fn mul(self: float_t, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a*b[0]; res[1] = a*b[1]; res[2] = a*b[2]; res[3] = a*b[3]; res[4] = a*b[4]; res[5] = a*b[5]; res[6] = a*b[6]; res[7] = a*b[7]; res[8] = a*b[8]; res[9] = a*b[9]; res[10] = a*b[10]; res[11] = a*b[11]; res[12] = a*b[12]; res[13] = a*b[13]; res[14] = a*b[14]; res[15] = a*b[15]; res[16] = a*b[16]; res[17] = a*b[17]; res[18] = a*b[18]; res[19] = a*b[19]; res[20] = a*b[20]; res[21] = a*b[21]; res[22] = a*b[22]; res[23] = a*b[23]; res[24] = a*b[24]; res[25] = a*b[25]; res[26] = a*b[26]; res[27] = a*b[27]; res[28] = a*b[28]; res[29] = a*b[29]; res[30] = a*b[30]; res[31] = a*b[31]; res } } // muls // multivector/scalar multiplication impl Mul<float_t> for R401 { type Output = R401; fn mul(self: R401, b: float_t) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]*b; res[1] = a[1]*b; res[2] = a[2]*b; res[3] = a[3]*b; res[4] = a[4]*b; res[5] = a[5]*b; res[6] = a[6]*b; res[7] = a[7]*b; res[8] = a[8]*b; res[9] = a[9]*b; res[10] = a[10]*b; res[11] = a[11]*b; res[12] = a[12]*b; res[13] = a[13]*b; res[14] = a[14]*b; res[15] = a[15]*b; res[16] = a[16]*b; res[17] = a[17]*b; res[18] = a[18]*b; res[19] = a[19]*b; res[20] = a[20]*b; res[21] = a[21]*b; res[22] = a[22]*b; res[23] = a[23]*b; res[24] = a[24]*b; res[25] = a[25]*b; res[26] = a[26]*b; res[27] = a[27]*b; res[28] = a[28]*b; res[29] = a[29]*b; res[30] = a[30]*b; res[31] = a[31]*b; res } } impl Div<float_t> for R401 { type Output = R401; fn div(self: R401, b: float_t) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]/b; res[1] = a[1]/b; res[2] = a[2]/b; res[3] = a[3]/b; res[4] = a[4]/b; res[5] = a[5]/b; res[6] = a[6]/b; res[7] = a[7]/b; res[8] = a[8]/b; res[9] = a[9]/b; res[10] = a[10]/b; res[11] = a[11]/b; res[12] = a[12]/b; res[13] = a[13]/b; res[14] = a[14]/b; res[15] = a[15]/b; res[16] = a[16]/b; res[17] = a[17]/b; res[18] = a[18]/b; res[19] = a[19]/b; res[20] = a[20]/b; res[21] = a[21]/b; res[22] = a[22]/b; res[23] = a[23]/b; res[24] = a[24]/b; res[25] = a[25]/b; res[26] = a[26]/b; res[27] = a[27]/b; res[28] = a[28]/b; res[29] = a[29]/b; res[30] = a[30]/b; res[31] = a[31]/b; res } } // sadd // scalar/multivector addition impl Add<R401> for float_t { type Output = R401; fn add(self: float_t, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a+b[0]; res[1] = b[1]; res[2] = b[2]; res[3] = b[3]; res[4] = b[4]; res[5] = b[5]; res[6] = b[6]; res[7] = b[7]; res[8] = b[8]; res[9] = b[9]; res[10] = b[10]; res[11] = b[11]; res[12] = b[12]; res[13] = b[13]; res[14] = b[14]; res[15] = b[15]; res[16] = b[16]; res[17] = b[17]; res[18] = b[18]; res[19] = b[19]; res[20] = b[20]; res[21] = b[21]; res[22] = b[22]; res[23] = b[23]; res[24] = b[24]; res[25] = b[25]; res[26] = b[26]; res[27] = b[27]; res[28] = b[28]; res[29] = b[29]; res[30] = b[30]; res[31] = b[31]; res } } // adds // multivector/scalar addition impl Add<float_t> for R401 { type Output = R401; fn add(self: R401, b: float_t) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]+b; res[1] = a[1]; res[2] = a[2]; res[3] = a[3]; res[4] = a[4]; res[5] = a[5]; res[6] = a[6]; res[7] = a[7]; res[8] = a[8]; res[9] = a[9]; res[10] = a[10]; res[11] = a[11]; res[12] = a[12]; res[13] = a[13]; res[14] = a[14]; res[15] = a[15]; res[16] = a[16]; res[17] = a[17]; res[18] = a[18]; res[19] = a[19]; res[20] = a[20]; res[21] = a[21]; res[22] = a[22]; res[23] = a[23]; res[24] = a[24]; res[25] = a[25]; res[26] = a[26]; res[27] = a[27]; res[28] = a[28]; res[29] = a[29]; res[30] = a[30]; res[31] = a[31]; res } } // ssub // scalar/multivector subtraction impl Sub<R401> for float_t { type Output = R401; fn sub(self: float_t, b: R401) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a-b[0]; res[1] = -b[1]; res[2] = -b[2]; res[3] = -b[3]; res[4] = -b[4]; res[5] = -b[5]; res[6] = -b[6]; res[7] = -b[7]; res[8] = -b[8]; res[9] = -b[9]; res[10] = -b[10]; res[11] = -b[11]; res[12] = -b[12]; res[13] = -b[13]; res[14] = -b[14]; res[15] = -b[15]; res[16] = -b[16]; res[17] = -b[17]; res[18] = -b[18]; res[19] = -b[19]; res[20] = -b[20]; res[21] = -b[21]; res[22] = -b[22]; res[23] = -b[23]; res[24] = -b[24]; res[25] = -b[25]; res[26] = -b[26]; res[27] = -b[27]; res[28] = -b[28]; res[29] = -b[29]; res[30] = -b[30]; res[31] = -b[31]; res } } // subs // multivector/scalar subtraction impl Sub<float_t> for R401 { type Output = R401; fn sub(self, b: float_t) -> R401 { let mut res = R401::zero(); let a = self; res[0] = a[0]-b; res[1] = a[1]; res[2] = a[2]; res[3] = a[3]; res[4] = a[4]; res[5] = a[5]; res[6] = a[6]; res[7] = a[7]; res[8] = a[8]; res[9] = a[9]; res[10] = a[10]; res[11] = a[11]; res[12] = a[12]; res[13] = a[13]; res[14] = a[14]; res[15] = a[15]; res[16] = a[16]; res[17] = a[17]; res[18] = a[18]; res[19] = a[19]; res[20] = a[20]; res[21] = a[21]; res[22] = a[22]; res[23] = a[23]; res[24] = a[24]; res[25] = a[25]; res[26] = a[26]; res[27] = a[27]; res[28] = a[28]; res[29] = a[29]; res[30] = a[30]; res[31] = a[31]; res } } impl R401 { pub fn norm(self: Self) -> float_t { let scalar_part = (self * self.conjugate())[0]; scalar_part.abs().sqrt() } pub fn inorm(self: Self) -> float_t { self.dual().norm() } pub fn normalized(self: Self) -> Self { self * (1.0 / self.norm()) } pub fn commutator(self, mv: R401) -> Self { 0.5 * (self * mv - mv * self) } } fn main() { println!("e0*e0 : {}", e0 * e0); println!("pss : {}", e01234); println!("pss*pss : {}", e01234*e01234); }