namespace reps { TD struct mvector { array<typename D::Number, D::Dim> values; typename D::Number& operator [] (int i) { return values[i]; } const typename D::Number& operator [] (int i) const { return values[i]; } mvector operator + (const mvector& M) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] + M[i]; return result; } mvector operator - (const mvector& M) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] - M[i]; return result; } mvector operator * (const typename D::Number& x) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] * x; return result; } mvector operator / (const typename D::Number& x) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] / x; return result; } mvector operator * (int x) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] * x; return result; } mvector operator / (int x) const { mvector result; for(int i=0; i<D::Dim; i++) result[i] = self[i] / x; return result; } }; TD struct matrix { array<array<typename D::Number, D::Dim>, D::Dim> values; array<typename D::Number, D::Dim>& operator [] (int i) { return values[i]; } const array<typename D::Number, D::Dim>& operator [] (int i) const { return values[i]; } matrix operator * (const matrix& M) const { matrix result; for(int i=0; i<D::Dim; i++) for(int k=0; k<D::Dim; k++) { result[i][k] = typename D::Number(0); for(int j=0; j<D::Dim; j++) result[i][k] += self[i][j] * M[j][k]; } return result; } mvector<D> operator * (const mvector<D>& V) const { mvector<D> result; for(int i=0; i<D::Dim; i++) { result[i] = typename D::Number(0); for(int j=0; j<D::Dim; j++) result[i] += self[i][j] * V[j]; } return result; } matrix operator * (const typename D::Number& x) const { matrix result; for(int i=0; i<D::Dim; i++) for(int j=0; j<D::Dim; j++) result[i][j] = self[i][j] * x; return result; } matrix operator / (const typename D::Number& x) const { matrix result; for(int i=0; i<D::Dim; i++) for(int j=0; j<D::Dim; j++) result[i][j] = self[i][j] / x; return result; } matrix operator * (int x) const { matrix result; for(int i=0; i<D::Dim; i++) for(int j=0; j<D::Dim; j++) result[i][j] = self[i][j] * x; return result; } matrix operator / (int x) const { matrix result; for(int i=0; i<D::Dim; i++) for(int j=0; j<D::Dim; j++) result[i][j] = self[i][j] / x; return result; } }; TD constexpr mvector<D> zero_vector() { mvector<D> result; for(int i=0; i<D::Dim; i++) result[i] = typename D::Number(0); return result; } TD constexpr mvector<D> unit_vector(int id) { mvector<D> result; for(int i=0; i<D::Dim; i++) result[i] = typename D::Number(0); result[id] = typename D::Number(1); return result; } TD struct multivector_data { using Number = typename D::Number; static constexpr int Dim = 1<<D::Dim; static constexpr int Flipped = -1; }; TD using multivector = mvector<multivector_data<D>>; TD std::string nz(const multivector<D>& a) { constexpr int mdim = 1<<D::Dim; using Number = typename D::Number; hr::shstream str; for(int i=0; i<mdim; i++) if(abs(a[i]) > Number(1e-10)) { if(str.s != "") print(str, " "); if(a[i] > Number(0)) print(str, "+"); print(str, a[i]); for(int u=0; u<D::Dim; u++) if(i & (1<<u)) print(str, hr::s0 + char('A'+u)); } if(str.s == "") return "0"; return str.s; } TD constexpr multivector<D> unit(const typename D::Number& a) { auto res = zero_vector<multivector_data<D>>(); res[0] = a; return res; } TD constexpr multivector<D> embed(const mvector<D>& a) { auto res = zero_vector<multivector_data<D>>(); for(int i=0; i<D::Dim; i++) res[1<<i] = a[i]; return res; } TD constexpr mvector<D> unembed(const multivector<D>& a) { mvector<D> res; for(int i=0; i<D::Dim; i++) res[i] = a[1<<i]; return res; } /* for clarity */ using mvindex = int; using signtype = int; /* mvindex decimal 10 (binary 1010) corresponds to unit_vector(1) * unit_vector(3) */ TD constexpr signtype conj_sign(mvindex mvid) { int b = __builtin_popcount(mvid); b = b * (b+1) / 2; return (b&1) ? -1 : 1; } TD constexpr signtype tra_sign(mvindex mvid) { int b = __builtin_popcount(mvid); b = b * (b-1) / 2; return (b&1) ? -1 : 1; } TD constexpr signtype mul_sign(mvindex a, mvindex b) { int flips = 0; for(int i=0; i<D::Dim; i++) if(b & (1<<i)) { // we will need to swap it with that many 1-bits of a flips += __builtin_popcount(a & ((1<<i)-1)); if((i == D::Flipped) && (a & (1<<i))) flips++; } return (flips&1) ? -1 : 1; } TD struct all { static constexpr bool check(mvindex a) { return true; } static constexpr bool isflat(mvindex a) { return false; } }; TD struct even { static constexpr bool check(mvindex a) { return __builtin_popcount(a) % 2 == 0; } static constexpr bool isflat(mvindex a) { return false; } }; TD struct flat_even { static constexpr bool check(mvindex a) { return __builtin_popcount(a) % 2 == 0; } static constexpr bool isflat(mvindex a) { return nm == nmFlatten && a == 0; } }; TD struct odd { static constexpr bool check(mvindex a) { return __builtin_popcount(a) % 2 == 1; } static constexpr bool isflat(mvindex a) { return false; } }; TD struct units { static constexpr bool check(mvindex a) { return a == 0; } static constexpr bool isflat(mvindex a) { return false; } }; TD struct rotational { static constexpr bool check(mvindex a) { return __builtin_popcount(a) % 2 == 0 && a < (1<<(D::Dim-1)); } static constexpr bool isflat(mvindex a) { return false; } }; TD struct underling { static constexpr bool check(mvindex a) { return __builtin_popcount(a) == 1; } static constexpr bool isflat(mvindex a) { return false; } }; TD struct flat_underling { static constexpr bool check(mvindex a) { return __builtin_popcount(a) == 1; } static constexpr bool isflat(mvindex a) { return nm == nmFlatten && a == 1<<(D::Dim-1); } }; TD struct poincare { static constexpr bool check(mvindex a) { return __builtin_popcount(a ^ (1<<(D::Dim-1))) == 1; } static constexpr bool isflat(mvindex a) { return false; } }; TD multivector<D> multimul(const multivector<D>& a, const multivector<D>& b) { constexpr int mdim = 1<<D::Dim; auto res = zero_vector<multivector_data<D>>(); for(mvindex i=0; i<mdim; i++) for(mvindex j=0; j<mdim; j++) { res[i^j] += a[i] * b[j] * mul_sign<D>(i, j); } return res; } template<class A, class B, class C, class D> multivector<D> chkmul(const multivector<D>& a, const multivector<D>& b) { constexpr int mdim = 1<<D::Dim; auto res = zero_vector<multivector_data<D>>(); /* we initialize with 0s and then add stuff, so one add per component is not necessary */ for(mvindex i=0; i<mdim; i++) if(C::check(i)) cbc[cbcAdd]--; for(mvindex i=0; i<mdim; i++) if(A::check(i)) for(mvindex j=0; j<mdim; j++) if(B::check(j) && C::check(i^j)) { if(A::isflat(i)) res[i^j] += b[j] * mul_sign<D>(i, j); else if(B::isflat(j)) res[i^j] += a[i] * mul_sign<D>(i, j); else res[i^j] += a[i] * b[j] * mul_sign<D>(i, j); } return res; } TD multivector<D> conjugate(const multivector<D>& a) { constexpr int mdim = 1<<D::Dim; auto res = a; for(int i=0; i<mdim; i++) res[i] *= conj_sign<D>(i); return res; } TD multivector<D> transpose(const multivector<D>& a) { constexpr int mdim = 1<<D::Dim; auto res = a; for(int i=0; i<mdim; i++) res[i] *= tra_sign<D>(i); return res; } template<class C, class D> multivector<D> apply_nm(multivector<D> a); TD using poincare_rotation = std::pair<multivector<D>, multivector<D>>; /** decompose o into the poincare part and the rotational component */ TD poincare_rotation<D> despin2(const multivector<D>& a) { auto p = a; for(int i=(1<<(D::Dim-1)); i<(1<<(D::Dim)); i++) p[i] = typename D::Number(0); p = p * pow(chkmul<rotational<D>,rotational<D>,units<D>>(p, conjugate(p))[0], -0.5); auto p1 = chkmul<even<D>,rotational<D>,poincare<D>>(a, conjugate(p)); return {apply_nm<poincare<D>, D>(p1), p}; } /** remove the rotational component of a, leaving only the poincare part */ TD multivector<D> despin(const multivector<D>& a) { auto p = a; for(int i=(1<<(D::Dim-1)); i<(1<<(D::Dim)); i++) p[i] = typename D::Number(0); auto p1 = chkmul<even<D>,rotational<D>,poincare<D>>(a, conjugate(p)); if(nm == nmInvariant) return p1 * pow(chkmul<rotational<D>,rotational<D>,units<D>>(p, conjugate(p))[0], -0.5); return apply_nm<poincare<D>, D>(p1); } TD std::string nzv(const mvector<D>& a) { return "vector(" + nz(embed(a)) + ")"; } TD std::string nzv(const matrix<D>& a) { return "<matrix>"; } template<class C, class D> typename D::Number sqnorm(multivector<D> a) { using N = typename D::Number; auto res = chkmul<C, C, units<D>>(a, conjugate(a))[0]; if(res <= N(0) || isinf(res) || isnan(res)) res = N(1); return res; } TD typename D::Number sqnorm(mvector<D> a) { using N = typename D::Number; N res(0); for(int i=0; i<D::Dim; i++) res += a[i] * a[i] * (i == D::Flipped ? -1:1); if(D::Flipped != -1) res = -res; if(nm ==nmWeak && (res <= N(0) || isinf(res) || isnan(res))) res = N(1); return res; } /** if nm is set to nmFlatten or nmForced or nmBinary, apply the requested operation */ template<class C, class D> multivector<D> flatten(multivector<D> a) { using N = typename D::Number; auto divby = a[0]; a[0] = N(1); for(int i=1; i<(1<<D::Dim); i++) if(C::check(i)) a[i] /= divby; return a; } template<class C, class D> multivector<D> apply_nm(multivector<D> a) { if(nm == nmFlatten) return flatten<C>(a); if(nm == nmForced || nm == nmWeak) return a * pow(sqnorm<C,D>(a), -0.5); if(nm == nmBinary) { while(a[0] >= 2) { a = a / 2; } while(a[0] > 0 && a[0] < 0.5) { a = a * 2; } } return a; } TD mvector<D> apply_nm(mvector<D> a) { if(nm == nmFlatten) { cbc[cbcDiv]--; return a / a[D::Dim-1]; } if(nm == nmForced || nm == nmWeak) return a * pow(sqnorm<D>(a), -0.5); if(nm == nmBinary) { while(a[D::Dim-1] >= 2) { a = a / 2; } while(a[D::Dim-1] > 0 && a[D::Dim-1] < 0.5) { a = a * 2; } } return a; } /** get b which is a coordinate of a, but in normalized form. That is, if a is normalized simply return b, otherwise, multiply b appropriately */ template<class C, class D, class E> E get_normalized(multivector<D> a, E b) { if(nm != nmInvariant && nm != nmForced) return b * pow(sqnorm<C,D>(a), -0.5); return b; } template<class D, class E> E get_normalized(mvector<D> a, E b) { if(nm != nmInvariant && nm != nmForced) return b * pow(sqnorm<D>(a), -0.5); return b; } }