Skip to content

Instantly share code, notes, and snippets.

@Orbots

Orbots/ga3.h

Last active May 30, 2024 23:38
Show Gist options
  • Save Orbots/3d8a72ef1edc58057c748f6f15eb0480 to your computer and use it in GitHub Desktop.
Save Orbots/3d8a72ef1edc58057c748f6f15eb0480 to your computer and use it in GitHub Desktop.
Download ZIP
Geometric Algebra with signature (3,0,0)
Raw
ga3.h
// Example:
// ga_versor q = normalize(add(1.0, (ga_bivector){ .e12 = 1.0 }));
// see tests at end of file for more
#pragma once
#include <stdarg.h>
#include <complex.h>
#include <math.h>
//---------- types
typedef double ga_scalar;
typedef struct ga_vector
{
ga_scalar e1;
ga_scalar e2;
ga_scalar e3;
} ga_vector;
typedef struct ga_bivector
{
ga_scalar e23;
ga_scalar e31;
ga_scalar e12;
} ga_bivector;
typedef struct ga_versor
{
ga_scalar s;
ga_bivector B;
} ga_versor;
typedef struct ga_pseudoscalar
{
ga_scalar e123;
} ga_pseudoscalar;
//----- returnable types that aren't useful enough to fully support.
// i.e. ga_paravector p = add(1, (ga_vector){.e1 = 1 }); Is supported, but [add | mul | ...](p, _) may or may not
// all depends on how useful that type is for applications of G(3).
// you can convert to a multivector which is better supported with to_multivector
typedef struct ga_oddversor
{
ga_vector v;
ga_pseudoscalar P;
} ga_oddversor;
typedef struct ga_paravector
{
ga_scalar s;
ga_vector v;
} ga_paravector;
typedef struct ga_biparavector
{
ga_vector v;
ga_bivector B;
} ga_biparavector;
typedef struct ga_triparavector
{
ga_bivector B;
ga_pseudoscalar P;
} ga_triparavector;
typedef struct ga_multivector
{
ga_scalar s;
ga_vector v;
ga_bivector B;
ga_pseudoscalar P;
} ga_multivector;
typedef struct ga_empty{} ga_empty;
typedef struct ga_grade_empty{} ga_grade_empty;
typedef struct ga_grade__1{} ga_grade__1;
typedef struct ga_grade__2{} ga_grade__2;
typedef struct ga_grade__3{} ga_grade__3;
typedef struct ga_grade_0{} ga_grade_0;
typedef struct ga_grade_1{} ga_grade_1;
typedef struct ga_grade_2{} ga_grade_2;
typedef struct ga_grade_3{} ga_grade_3;
typedef struct ga_grade_4{} ga_grade_4;
typedef struct ga_grade_5{} ga_grade_5;
typedef struct ga_grade_6{} ga_grade_6;
//------- constructors
ga_grade_empty grade_empty(){ return (ga_grade_empty){}; }
ga_grade_0 grade0(){ return (ga_grade_0){}; }
ga_grade_1 grade1(){ return (ga_grade_1){}; }
ga_grade_2 grade2(){ return (ga_grade_2){}; }
ga_grade_3 grade3(){ return (ga_grade_3){}; }
#define sub_grade(x, y) _Generic((x), \
ga_grade_0: \
_Generic((y), ga_grade_0: (ga_grade_0){}, ga_grade_1: (ga_grade__1){}, ga_grade_2: (ga_grade__2){}, ga_grade_3: (ga_grade__3){}), \
ga_grade_1: \
_Generic((y), ga_grade_0: (ga_grade_1){}, ga_grade_1: (ga_grade_0){}, ga_grade_2: (ga_grade__1){}, ga_grade_3: (ga_grade__2){}), \
ga_grade_2: \
_Generic((y), ga_grade_0: (ga_grade_2){}, ga_grade_1: (ga_grade_1){}, ga_grade_2: (ga_grade_0){}, ga_grade_3: (ga_grade__1){}), \
ga_grade_3: \
_Generic((y), ga_grade_0: (ga_grade_3){}, ga_grade_1: (ga_grade_2){}, ga_grade_2: (ga_grade_1){}, ga_grade_3: (ga_grade_0){}) \
)
#define abs_grade(x) _Generic((x), ga_grade__1: (ga_grade_1){}, ga_grade__2: (ga_grade_2){}, ga_grade__3: (ga_grade_3){}, default: x )
ga_vector zero_v(){ return (ga_vector){0,0,0}; }
ga_bivector zero_B(){ return (ga_bivector){0,0,0}; }
ga_pseudoscalar zero_P(){ return (ga_pseudoscalar){0}; }
ga_vector ga_vec(ga_scalar x, ga_scalar y, ga_scalar z){ return (ga_vector){x,y,z}; }
ga_bivector ga_bivec(ga_scalar x, ga_scalar y, ga_scalar z){ return (ga_bivector){x,y,z}; }
ga_pseudoscalar ga_I(ga_scalar s){ return (ga_pseudoscalar){s}; }
ga_versor ga_rotor(ga_scalar s, ga_scalar x, ga_scalar y, ga_scalar z){ return (ga_versor){.s = s, .B = ga_bivec(x, y, z)}; }
ga_scalar always0_s(ga_scalar x, ...){ return 0.0; }
ga_scalar always0_v(ga_vector x, ...){ return 0.0; }
ga_scalar always0_B(ga_bivector x, ...){ return 0.0; }
ga_scalar always0_P(ga_pseudoscalar x, ...){ return 0.0; }
ga_scalar always0_q(ga_versor x, ...){ return 0.0; }
ga_scalar always0_o(ga_oddversor x, ...){ return 0.0; }
ga_scalar always0_M(ga_multivector x, ...){ return 0.0; }
ga_empty always_empty_(ga_empty x, ...){ return (ga_empty){}; }
ga_empty always_empty_s(ga_scalar x, ...){ return (ga_empty){}; }
ga_empty always_empty_v(ga_vector x, ...){ return (ga_empty){}; }
ga_empty always_empty_B(ga_bivector x, ...){ return (ga_empty){}; }
ga_empty always_empty_P(ga_pseudoscalar x, ...){ return (ga_empty){}; }
ga_empty always_empty_q(ga_versor x, ...){ return (ga_empty){}; }
ga_empty always_empty_o(ga_oddversor x, ...){ return (ga_empty){}; }
ga_empty always_empty_M(ga_multivector x, ...){ return (ga_empty){}; }
ga_empty self_(ga_empty x, ...){ return x; }
ga_scalar self_s(ga_scalar x, ...){ return x; }
ga_vector self_v(ga_vector x, ...){ return x; }
ga_bivector self_B(ga_bivector x, ...){ return x; }
ga_pseudoscalar self_P(ga_pseudoscalar x, ...){ return x; }
ga_versor self_q(ga_versor x, ...){ return x; }
ga_oddversor self_o(ga_oddversor x, ...){ return x; }
ga_multivector self_M(ga_multivector x, ...){ return x; }
ga_scalar other_s(ga_empty x, ga_scalar y){ return y; }
ga_vector other_v(ga_empty x, ga_vector y){ return y; }
ga_bivector other_B(ga_empty x, ga_bivector y){ return y; }
ga_pseudoscalar other_P(ga_empty x, ga_pseudoscalar y){ return y; }
ga_versor other_q(ga_empty x, ga_versor y){ return y; }
ga_oddversor other_o(ga_empty x, ga_oddversor y){ return y; }
ga_multivector other_M(ga_empty x, ga_multivector y){ return y; }
ga_multivector to_multivector_o(ga_oddversor o){ return (ga_multivector){ .s = 0, .v = o.v, .B = zero_B(), .P = zero_P() }; }
ga_multivector to_multivector_p1(ga_paravector p1){ return (ga_multivector){ .s = p1.s, .v = p1.v, .B = zero_B(), .P = zero_P() }; }
ga_multivector to_multivector_p2(ga_biparavector p2){ return (ga_multivector){ .s = 0, .v = p2.v, .B = p2.B, .P = zero_P() }; }
ga_multivector to_multivector_p3(ga_triparavector p3){ return (ga_multivector){ .s = 0, .v = zero_v(), .B = p3.B, .P = p3.P }; }
#define to_multivector(x) _Generic((x), ga_paravector: to_multivector_p1, ga_biparavector: to_multivector_p2, ga_triparavector: to_multivector_p3, \
ga_oddversor: to_multivector_o, ga_multivector: self_M )(x)
#define sqr_s(x) (x*x)
#define sign_s(x) (x < 0 ? -1.0 : 1.0)
#define SYMMETRIC_OP(foo, bar, x, y, r) r foo(x a, y b){ return bar(b,a); }
#define ANTISYMMETRIC_OP(foo, bar, x, y, r) r foo(x a, y b){ return negate(bar(a,b)); }
//------- operations
ga_scalar add_ss(ga_scalar a, ga_scalar b){ return a+b; }
ga_vector add_vv(ga_vector a, ga_vector b){ return (ga_vector){a.e1 + b.e1, a.e2 + b.e2, a.e3 + b.e3}; }
ga_bivector add_BB(ga_bivector a, ga_bivector b){ return (ga_bivector){a.e23 + b.e23, a.e31 + b.e31, a.e12 + b.e12}; }
ga_pseudoscalar add_PP(ga_pseudoscalar a, ga_pseudoscalar b){ return (ga_pseudoscalar){a.e123 + b.e123}; }
ga_multivector add_sv(ga_scalar s, ga_vector v){ return to_multivector(((ga_paravector){s, v})); }
SYMMETRIC_OP(add_vs, add_sv, ga_vector, ga_scalar, ga_multivector)
ga_versor add_sB(ga_scalar s, ga_bivector B){ return (ga_versor){s, B}; }
ga_versor add_Bs(ga_bivector B, ga_scalar s){ return (ga_versor){s, B}; }
ga_versor add_sq(ga_scalar s, ga_versor q){ return (ga_versor){s+q.s, q.B}; }
ga_versor add_qs(ga_versor q, ga_scalar s){ return (ga_versor){s+q.s, q.B}; }
ga_versor add_Bq(ga_bivector B, ga_versor q){ return (ga_versor){q.s, add_BB(B, q.B)}; }
ga_versor add_qB(ga_versor q, ga_bivector B){ return (ga_versor){q.s, add_BB(B, q.B)}; }
ga_versor add_qq(ga_versor a, ga_versor b){ return (ga_versor){a.s+b.s, add_BB(a.B, b.B)}; }
ga_oddversor add_oo(ga_oddversor a, ga_oddversor b) { return (ga_oddversor){.v = add_vv(a.v, b.v), .P = add_PP(a.P, b.P)}; }
ga_oddversor add_vo(ga_vector a, ga_oddversor b) { return (ga_oddversor){.v = add_vv(a, b.v), .P = b.P}; }
ga_oddversor add_ov(ga_oddversor a, ga_vector b) { return add_vo(b, a); }
ga_oddversor add_oP(ga_oddversor a, ga_pseudoscalar b) { return (ga_oddversor){.v = a.v, .P = add_PP(a.P, b)}; }
ga_oddversor add_Po(ga_pseudoscalar b, ga_oddversor a) { return (ga_oddversor){.v = a.v, .P = add_PP(a.P, b)}; }
ga_multivector add_so(ga_scalar a, ga_oddversor b) { return (ga_multivector){a, b.v, zero_B(), b.P}; }
ga_multivector add_os(ga_oddversor b, ga_scalar a) { return (ga_multivector){a, b.v, zero_B(), b.P}; }
ga_multivector add_Bo(ga_bivector a, ga_oddversor b) { return (ga_multivector){0, b.v, a, b.P}; }
ga_multivector add_oB(ga_oddversor b, ga_bivector a) { return (ga_multivector){0, b.v, a, b.P}; }
ga_multivector add_qo(ga_versor a, ga_oddversor b) { return (ga_multivector){a.s, b.v, a.B, b.P}; }
ga_multivector add_oq(ga_oddversor a, ga_versor b) { return (ga_multivector){b.s, a.v, b.B, a.P}; }
ga_multivector add_vq(ga_vector v, ga_versor q){ return (ga_multivector){q.s, v, q.B, zero_P()}; }
ga_multivector add_qv(ga_versor q, ga_vector v){ return (ga_multivector){q.s, v, q.B, zero_P()}; }
ga_multivector add_sP(ga_scalar s, ga_pseudoscalar P){ return (ga_multivector){s, zero_v(), zero_B(), P}; }
ga_multivector add_Ps(ga_pseudoscalar P, ga_scalar s){ return (ga_multivector){s, zero_v(), zero_B(), P}; }
ga_multivector add_qP(ga_versor q, ga_pseudoscalar P){ return (ga_multivector){q.s, zero_v(), q.B, P}; }
ga_multivector add_Pq(ga_pseudoscalar P, ga_versor q){ return (ga_multivector){q.s, zero_v(), q.B, P}; }
ga_oddversor add_vP(ga_vector v, ga_pseudoscalar P){ return (ga_oddversor){ v, P }; }
ga_oddversor add_Pv(ga_pseudoscalar P, ga_vector v){ return (ga_oddversor){ v, P }; }
ga_multivector add_vB(ga_vector v, ga_bivector B){ return to_multivector(((ga_biparavector){ v, B })); }
ga_multivector add_Bv(ga_bivector B, ga_vector v){ return to_multivector(((ga_biparavector){ v, B })); }
ga_multivector add_BP(ga_bivector B, ga_pseudoscalar P){ return to_multivector(((ga_triparavector){ B, P })); }
ga_multivector add_PB(ga_pseudoscalar P, ga_bivector B){ return to_multivector(((ga_triparavector){ B, P })); }
ga_multivector add_MM(ga_multivector M, ga_multivector N)
{
return (ga_multivector){
.s = M.s + N.s,
.v = add_vv(M.v, N.v),
.B = add_BB(M.B, N.B),
.P = (ga_pseudoscalar){M.P.e123 + N.P.e123}
};
}
ga_multivector add_Ms(ga_multivector M, ga_scalar s)
{
M.s += s;
return M;
}
SYMMETRIC_OP(add_sM, add_Ms, ga_scalar, ga_multivector, ga_multivector)
ga_multivector add_Mv(ga_multivector M, ga_vector v)
{
M.v = add_vv(M.v, v);
return M;
}
SYMMETRIC_OP(add_vM, add_Mv, ga_vector, ga_multivector, ga_multivector)
ga_multivector add_MB(ga_multivector M, ga_bivector B)
{
M.B = add_BB(M.B, B);
return M;
}
SYMMETRIC_OP(add_BM, add_MB, ga_bivector, ga_multivector, ga_multivector)
ga_multivector add_MP(ga_multivector M, ga_pseudoscalar P)
{
M.P = add_PP(M.P, P);
return M;
}
SYMMETRIC_OP(add_PM, add_MP, ga_pseudoscalar, ga_multivector, ga_multivector)
ga_multivector add_Mq(ga_multivector M, ga_versor q)
{
M.s += q.s;
M.B = add_BB(M.B, q.B);
return M;
}
SYMMETRIC_OP(add_qM, add_Mq, ga_versor, ga_multivector, ga_multivector)
ga_multivector add_Mo(ga_multivector M, ga_oddversor o)
{
M.v = add_vv(M.v, o.v);
M.P = add_PP(M.P, o.P);
return M;
}
SYMMETRIC_OP(add_oM, add_Mo, ga_oddversor, ga_multivector, ga_multivector)
#define add(x, y) _Generic((x), \
ga_empty: \
_Generic((y), ga_empty: self_, ga_scalar: other_s, ga_vector: other_v, ga_bivector: other_B, ga_versor: other_q, ga_oddversor: other_o, ga_pseudoscalar: other_P, ga_multivector: other_M), \
ga_scalar: \
_Generic((y), ga_empty: self_s, ga_scalar: add_ss, ga_vector: add_sv, ga_bivector: add_sB, ga_versor: add_sq, ga_oddversor: add_so, ga_pseudoscalar: add_sP, ga_multivector: add_sM), \
ga_vector: \
_Generic((y), ga_empty: self_v, ga_scalar: add_vs, ga_vector: add_vv, ga_bivector: add_vB, ga_versor: add_vq, ga_oddversor: add_vo, ga_pseudoscalar: add_vP, ga_multivector: add_vM), \
ga_bivector: \
_Generic((y), ga_empty: self_B, ga_scalar: add_Bs, ga_vector: add_Bv, ga_bivector: add_BB, ga_versor: add_Bq, ga_oddversor: add_Bo, ga_pseudoscalar: add_BP, ga_multivector: add_BM), \
ga_versor: \
_Generic((y), ga_empty: self_q, ga_scalar: add_qs, ga_vector: add_qv, ga_bivector: add_qB, ga_versor: add_qq, ga_oddversor: add_qo, ga_pseudoscalar: add_qP, ga_multivector: add_qM), \
ga_oddversor: \
_Generic((y), ga_empty: self_o, ga_scalar: add_os, ga_vector: add_ov, ga_bivector: add_oB, ga_versor: add_oq, ga_oddversor: add_oo, ga_pseudoscalar: add_oP, ga_multivector: add_oM), \
ga_pseudoscalar: \
_Generic((y), ga_empty: self_P, ga_scalar: add_Ps, ga_vector: add_Pv, ga_bivector: add_PB, ga_versor: add_Pq, ga_oddversor: add_Po, ga_pseudoscalar: add_PP, ga_multivector: add_PM), \
ga_multivector: \
_Generic((y), ga_empty: self_M, ga_scalar: add_Ms, ga_vector: add_Mv, ga_bivector: add_MB, ga_versor: add_Mq, ga_oddversor: add_Mo, ga_pseudoscalar: add_MP, ga_multivector: add_MM) \
)(x, y)
ga_scalar dot_vv(ga_vector a, ga_vector b){ return a.e1*b.e1 + a.e2*b.e2 + a.e3*b.e3; };
ga_bivector wedge_vv(ga_vector a, ga_vector b)
{
return (ga_bivector) { .e23 = a.e2*b.e3-b.e2*a.e3,
.e31 = a.e3*b.e1-b.e3*a.e1,
.e12 = a.e1*b.e2-b.e1*a.e2 };
}
ga_scalar mul_ss(ga_scalar a, ga_scalar b){ return a*b; }
ga_vector mul_sv(ga_scalar a, ga_vector b){ return (ga_vector){a*b.e1, a*b.e2, a*b.e3}; }
SYMMETRIC_OP(mul_vs, mul_sv, ga_vector, ga_scalar, ga_vector)
ga_bivector mul_sB(ga_scalar a, ga_bivector b){ return (ga_bivector){a*b.e23, a*b.e31, a*b.e12}; }
SYMMETRIC_OP(mul_Bs, mul_sB, ga_bivector, ga_scalar, ga_bivector)
ga_versor mul_sq(ga_scalar s, ga_versor q){ return (ga_versor){s*q.s, mul_sB(s, q.B)}; }
SYMMETRIC_OP(mul_qs, mul_sq, ga_versor, ga_scalar, ga_versor)
ga_pseudoscalar mul_sP(ga_scalar s, ga_pseudoscalar P){ return (ga_pseudoscalar){s*P.e123}; }
SYMMETRIC_OP(mul_Ps, mul_sP, ga_pseudoscalar, ga_scalar, ga_pseudoscalar)
ga_oddversor mul_so(ga_scalar s, ga_oddversor o){ return (ga_oddversor){mul_sv(s, o.v), mul_sP(s, o.P)}; }
SYMMETRIC_OP(mul_os, mul_so, ga_oddversor, ga_scalar, ga_oddversor)
ga_versor mul_vv(ga_vector a, ga_vector b){ return (ga_versor){ dot_vv(a, b), wedge_vv(a, b) }; }
ga_bivector mul_vP(ga_vector v, ga_pseudoscalar P){ return (ga_bivector){ .e23 = v.e1*P.e123, .e31 = v.e2*P.e123, .e12 = v.e3*P.e123 }; }
SYMMETRIC_OP(mul_Pv, mul_vP, ga_pseudoscalar, ga_vector, ga_bivector)
ga_vector mul_BP(ga_bivector B, ga_pseudoscalar P){ return (ga_vector){ .e1 = -B.e23*P.e123, .e2 = -B.e31*P.e123, .e3 = -B.e12*P.e123 }; }
SYMMETRIC_OP(mul_PB, mul_BP, ga_pseudoscalar, ga_bivector, ga_vector)
ga_oddversor mul_vB(ga_vector v, ga_bivector B)
{
return (ga_oddversor){.v = { .e1 = v.e3*B.e31 - v.e2*B.e12,
.e2 = v.e1*B.e12 - v.e3*B.e23,
.e3 = v.e2*B.e23 - v.e1*B.e31 },
.P = { v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12 }};
}
ga_oddversor mul_Bv(ga_bivector B, ga_vector v)
{
return (ga_oddversor){.v = { .e1 = -v.e3*B.e31 + v.e2*B.e12,
.e2 = -v.e1*B.e12 + v.e3*B.e23,
.e3 = -v.e2*B.e23 + v.e1*B.e31 },
.P = { v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12 }};
}
ga_oddversor mul_vq(ga_vector v, ga_versor q)
{
return add(mul_vB(v, q.B), mul_sv(q.s, v));
}
ga_oddversor mul_qv(ga_versor q, ga_vector v)
{
return add(mul_Bv(q.B, v), mul_sv(q.s, v));
}
ga_versor mul_vo(ga_vector v, ga_oddversor o)
{
return add(mul_vv(v, o.v), mul_vP(v, o.P));
}
ga_versor mul_ov(ga_oddversor o, ga_vector v)
{
return add(mul_vv(o.v, v), mul_Pv(o.P, v));
}
ga_versor mul_BB(ga_bivector A, ga_bivector B)
{
return (ga_versor){.B = { .e23 = A.e12*B.e31 - A.e31*B.e12,
.e31 = A.e23*B.e12 - A.e12*B.e23,
.e12 = A.e31*B.e23 - A.e23*B.e31 },
.s = -A.e23*B.e23 - A.e31*B.e31 - A.e12*B.e12};
}
ga_versor mul_Bq(ga_bivector B, ga_versor q)
{
return add(mul_sB(q.s, B), mul_BB(B, q.B));
}
ga_versor mul_qB(ga_versor q, ga_bivector B)
{
return add(mul_sB(q.s, B), mul_BB(q.B, B));
}
ga_oddversor mul_Bo(ga_bivector B, ga_oddversor o)
{
return add(mul_Bv(B, o.v), mul_BP(B, o.P));
}
ga_oddversor mul_oB(ga_oddversor o, ga_bivector B)
{
return add(mul_vB(o.v, B), mul_PB(o.P, B));
}
ga_versor mul_qq(ga_versor a, ga_versor b)
{
ga_versor AB = mul_BB(a.B, b.B);
return (ga_versor){ .s = a.s*b.s + AB.s, .B = add(AB.B, add(mul_sB(a.s, b.B), mul_sB(b.s, a.B))) };
}
ga_scalar mul_PP(ga_pseudoscalar a, ga_pseudoscalar b){ return -a.e123*b.e123; }
ga_versor mul_oo(ga_oddversor a, ga_oddversor b)
{
return add(mul_vv(a.v, b.v), add( add(mul_vP(a.v, b.P), mul_vP(b.v, a.P)), mul_PP(a.P, b.P)));
}
ga_oddversor mul_oq(ga_oddversor o, ga_versor q)
{
return add(mul_sv(q.s, o.v), add( add(mul_vB(o.v, q.B), mul_sP(q.s, o.P)), mul_BP(q.B, o.P)));
}
ga_oddversor mul_qo(ga_versor q, ga_oddversor o)
{
return add(mul_sv(q.s, o.v), add( add(mul_Bv(q.B, o.v), mul_sP(q.s, o.P)), mul_BP(q.B, o.P)));
}
ga_oddversor mul_qP( ga_versor q, ga_pseudoscalar P)
{
return add( mul_sP(q.s, P), mul_BP(q.B, P) );
}
SYMMETRIC_OP(mul_Pq, mul_qP, ga_pseudoscalar, ga_versor, ga_oddversor)
ga_versor mul_oP( ga_oddversor o, ga_pseudoscalar P)
{
return add( mul_vP(o.v, P), mul_PP(o.P, P) );
}
SYMMETRIC_OP(mul_Po, mul_oP, ga_pseudoscalar, ga_oddversor, ga_versor)
ga_multivector mul_Ms(ga_multivector M, ga_scalar s){ return (ga_multivector){M.s*s, mul_vs(M.v, s), mul_Bs(M.B, s), mul_Ps(M.P, s)}; }
SYMMETRIC_OP(mul_sM, mul_Ms, ga_scalar, ga_multivector, ga_multivector)
ga_multivector mul_MM(ga_multivector M, ga_multivector N)
{
ga_versor q = add(mul_vv(M.v, N.v), mul_BB(M.B, N.B));
ga_oddversor o = add(mul_vB(M.v, N.B), mul_Bv(M.B, N.v));
return (ga_multivector){
.s = M.s*N.s + q.s - M.P.e123*N.P.e123,
.v = add(add(add(mul_vs(M.v, N.s), mul_sv(M.s, N.v)), o.v), add(mul_PB(M.P, N.B), mul_BP(M.B, N.P))),
.B = add(add(add(mul_Bs(M.B, N.s), mul_sB(M.s, N.B)), q.B), add(mul_Pv(M.P, N.v), mul_vP(M.v, N.P))),
.P = (ga_pseudoscalar){M.P.e123*N.s + M.s*N.P.e123 + o.P.e123}
};
}
ga_multivector mul_Mv(ga_multivector M, ga_vector v)
{
ga_versor q = mul_vv(M.v, v);
ga_oddversor o = mul_Bv(M.B, v);
return (ga_multivector){
.s = q.s,
.v = add((mul_sv(M.s, v)), o.v),
.B = add(mul_Pv(M.P, v), q.B),
.P = o.P
};
}
ga_multivector mul_vM(ga_vector v, ga_multivector M)
{
ga_versor q = mul_vv(v, M.v);
ga_oddversor o = mul_vB(v, M.B);
return (ga_multivector){
.s = q.s,
.v = add((mul_sv(M.s, v)), o.v),
.B = add(mul_vP(v, M.P), q.B),
.P = o.P
};
}
ga_multivector mul_MB(ga_multivector M, ga_bivector B)
{
ga_versor q = mul_BB(M.B, B);
ga_oddversor o = mul_vB(M.v, B);
return (ga_multivector){
.s = q.s,
.v = add(o.v, mul_PB(M.P, B)),
.B = add(mul_sB(M.s, B), q.B),
.P = o.P
};
}
ga_multivector mul_BM(ga_bivector B, ga_multivector M)
{
ga_versor q = mul_BB(B, M.B);
ga_oddversor o = mul_Bv(B, M.v);
return (ga_multivector){
.s = q.s,
.v = add(o.v, mul_BP(B, M.P)),
.B = add(mul_sB(M.s, B), q.B),
.P = o.P
};
}
ga_multivector mul_Mq(ga_multivector M, ga_versor q)
{
ga_versor Mq = mul_BB(M.B, q.B);
ga_oddversor o = mul_vB(M.v, q.B);
return (ga_multivector){
.s = M.s*q.s + Mq.s,
.v = add(add(mul_vs(M.v, q.s), o.v), mul_PB(M.P, q.B)),
.B = add(add(mul_Bs(M.B, q.s), mul_sB(M.s, q.B)), Mq.B),
.P = (ga_pseudoscalar){M.P.e123*q.s + o.P.e123}
};
}
ga_multivector mul_qM(ga_versor q, ga_multivector M)
{
ga_versor qM = mul_BB(q.B, M.B);
ga_oddversor o = mul_Bv(q.B, M.v);
return (ga_multivector){
.s = M.s*q.s + qM.s,
.v = add(add(mul_vs(M.v, q.s), o.v), mul_BP(q.B, M.P)),
.B = add(add(mul_Bs(M.B, q.s), mul_sB(M.s, q.B)), qM.B),
.P = (ga_pseudoscalar){M.P.e123*q.s + o.P.e123}
};
}
ga_multivector mul_oM(ga_oddversor o, ga_multivector M)
{
return mul_MM(to_multivector(o), M);
}
ga_multivector mul_Mo( ga_multivector M, ga_oddversor o)
{
return mul_MM(M, to_multivector(o));
}
ga_multivector mul_MP(ga_multivector M, ga_pseudoscalar P)
{
return (ga_multivector){
.s = - M.P.e123*P.e123,
.v = mul_BP(M.B, P),
.B = mul_vP(M.v, P),
.P = (ga_pseudoscalar){M.s*P.e123}
};
}
ga_multivector mul_PM(ga_pseudoscalar P, ga_multivector M)
{
return (ga_multivector){
.s = - P.e123*M.P.e123,
.v = mul_PB(P, M.B),
.B = mul_Pv(P, M.v),
.P = (ga_pseudoscalar){M.s*P.e123}
};
}
#define mul(x, y) _Generic((x), \
ga_empty: \
_Generic((y), ga_empty: always_empty_, ga_scalar: always_empty_, ga_vector: always_empty_, ga_bivector: always_empty_, ga_versor: always_empty_, ga_oddversor: always_empty_, ga_pseudoscalar: always_empty_, ga_multivector: always_empty_), \
ga_scalar: \
_Generic((y), ga_empty: always_empty_s, ga_scalar: mul_ss, ga_vector: mul_sv, ga_bivector: mul_sB, ga_versor: mul_sq, ga_oddversor: mul_so, ga_pseudoscalar: mul_sP, ga_multivector: mul_sM), \
ga_vector: \
_Generic((y), ga_empty: always_empty_v, ga_scalar: mul_vs, ga_vector: mul_vv, ga_bivector: mul_vB, ga_versor: mul_vq, ga_oddversor: mul_vo, ga_pseudoscalar: mul_vP, ga_multivector: mul_vM), \
ga_bivector: \
_Generic((y), ga_empty: always_empty_B, ga_scalar: mul_Bs, ga_vector: mul_Bv, ga_bivector: mul_BB, ga_versor: mul_qB, ga_oddversor: mul_Bo, ga_pseudoscalar: mul_BP, ga_multivector: mul_BM), \
ga_versor: \
_Generic((y), ga_empty: always_empty_q, ga_scalar: mul_qs, ga_vector: mul_qv, ga_bivector: mul_qB, ga_versor: mul_qq, ga_oddversor: mul_qo, ga_pseudoscalar: mul_qP, ga_multivector: mul_qM), \
ga_oddversor: \
_Generic((y), ga_empty: always_empty_o, ga_scalar: mul_os, ga_vector: mul_ov, ga_bivector: mul_oB, ga_versor: mul_oq, ga_oddversor: mul_oo, ga_pseudoscalar: mul_oP, ga_multivector: mul_oM), \
ga_pseudoscalar: \
_Generic((y), ga_empty: always_empty_P, ga_scalar: mul_Ps, ga_vector: mul_Pv, ga_bivector: mul_PB, ga_versor: mul_Pq, ga_oddversor: mul_Po, ga_pseudoscalar: mul_PP, ga_multivector: mul_PM), \
ga_multivector: \
_Generic((y), ga_empty: always_empty_M, ga_scalar: mul_Ms, ga_vector: mul_Mv, ga_bivector: mul_MB, ga_versor: mul_Mq, ga_oddversor: mul_Mo, ga_pseudoscalar: mul_MP, ga_multivector: mul_MM) \
)(x, y)
ga_scalar negate_s(ga_scalar s){ return -s; }
ga_vector negate_v(ga_vector v){ return (ga_vector){-v.e1, -v.e2, -v.e3}; }
ga_bivector negate_B(ga_bivector B){ return (ga_bivector){-B.e23, -B.e31, -B.e12}; }
ga_versor negate_q(ga_versor q){ return (ga_versor){-q.s, negate_B(q.B)}; }
ga_pseudoscalar negate_P(ga_pseudoscalar P){ return (ga_pseudoscalar){-P.e123}; }
ga_oddversor negate_o(ga_oddversor o){ return (ga_oddversor){negate_v(o.v), negate_P(o.P)}; }
ga_multivector negate_M(ga_multivector M){ return (ga_multivector){-M.s, negate_v(M.v), negate_B(M.B), negate_P(M.P)}; }
#define negate(x) _Generic((x), \
ga_scalar: negate_s, \
ga_vector: negate_v, \
ga_bivector: negate_B, \
ga_versor: negate_q, \
ga_oddversor: negate_o, \
ga_pseudoscalar: negate_P, \
ga_multivector: negate_M \
)(x)
#define sub(a, b) add(a, negate(b))
ga_versor conj_q(ga_versor q){ return (ga_versor){q.s, negate(q.B)}; }
ga_oddversor conj_o(ga_oddversor o){ return (ga_oddversor){negate(o.v), o.P}; }
ga_multivector conj_M(ga_multivector M){ return (ga_multivector){M.s, negate(M.v), negate(M.B), M.P}; }
// clifford conjugate
#define conj(x) _Generic((x), \
ga_scalar: self_s, \
double complex: conj, \
ga_vector: negate_v, \
ga_bivector: negate_B, \
ga_versor: conj_q, \
ga_oddversor: conj_o, \
ga_pseudoscalar: negate_P, \
ga_multivector: conj_M \
)(x)
#define reverse(x) _Generic((x), \
ga_scalar: self_s, \
ga_vector: self_v, \
ga_bivector: negate_B, \
ga_versor: conj_q, \
ga_pseudoscalar: negate_P \
)(x)
#define grade(x) _Generic((x), \
ga_empty: grade0(), \
ga_scalar: grade0(), \
ga_vector: grade1(), \
ga_bivector: grade2(), \
ga_pseudoscalar: grade3() \
)
ga_scalar grade_select_0_q(ga_versor q, ...){ return q.s; }
ga_bivector grade_select_2_q(ga_versor q, ...){ return q.B; }
ga_vector grade_select_1_o(ga_oddversor o, ...){ return o.v; }
ga_pseudoscalar grade_select_3_o(ga_oddversor o, ...){ return o.P; }
ga_scalar grade_select_0_M(ga_multivector M, ...){ return M.s; }
ga_vector grade_select_1_M(ga_multivector M, ...){ return M.v; }
ga_bivector grade_select_2_M(ga_multivector M, ...){ return M.B; }
ga_pseudoscalar grade_select_3_M(ga_multivector M, ...){ return M.P; }
#define grade_select(x, y) _Generic((x), \
ga_empty: _Generic((y), default: self_), \
ga_scalar: _Generic((y), ga_grade_0: self_s, default: always_empty_s), \
ga_vector: _Generic((y), ga_grade_1: self_v, default: always_empty_v), \
ga_bivector: _Generic((y), ga_grade_2: self_B, default: always_empty_B), \
ga_versor: _Generic((y), ga_grade_0: grade_select_0_q, ga_grade_2: grade_select_2_q, default: always_empty_q), \
ga_oddversor: _Generic((y), ga_grade_1: grade_select_1_o, ga_grade_3: grade_select_3_o, default: always_empty_o), \
ga_pseudoscalar: _Generic((y), ga_grade_3: self_P, default: always_empty_P), \
ga_multivector: _Generic((y), ga_grade_0: grade_select_0_M, ga_grade_1: grade_select_1_M, ga_grade_2: grade_select_2_M, ga_grade_3: grade_select_3_M, default: always_empty_M) \
)(x, y)
#define select_scalar(x) grade_select(x, grade0())
// note on sandwich product. we can choose a different function for the right side of a*b*f(a) depending on the context
// the most correct is a*b*a^-1 but not always practical
// most commonly f = reverse. for paravectors you need f = conj, where conj is clifford_conjugate(a) = reverse(involute(a))
// for a plane-and-simple pga reflection you need f = -reverse depending on grade ( ommiting - for even grade ).
// an interesting one is a*b*conj(a)/(a*conj(a)). If you employ inv(a*b) = b*a, you see this is the same as a*b*inv(a)
// for many transforms a*conj(a) will be a scalar. Most commonly that scalar is 1.0. Finally this leads to just conj(a) on the rhs,
// which is then further commonly conj(a) == reverse(a).
// However orientation is domain specific and a single formula doesn't seem able to capture what is going on.
// so mirror space PGA for example uses -1^(x*r)*R*X*~R where x = grade(X), r = grade(R).
// Or you have two formulas depending grades.
// i.e. "the easiest way to handle the minus sign when sandwiching b with a is to use a * involute(b) * reverse(a) if a is odd and a * b * reverse(a) if a is even" -enki
// a*b*inv(a) where inv is simplified to conj. output same k-blade as input, otherwise it might be a multivector with 0's everywhere else
#define sandwich(a, b) grade_select((mul(mul(a, b), conj(a))), grade(b))
ga_pseudoscalar wedge_vB(ga_vector v, ga_bivector B){ return (ga_pseudoscalar){v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12}; }
SYMMETRIC_OP(wedge_Bv, wedge_vB, ga_bivector, ga_vector, ga_pseudoscalar)
ga_oddversor wedge_vq(ga_vector v, ga_versor q){ return (ga_oddversor){ mul_vs(v, q.s), wedge_vB(v, q.B) }; }
ga_oddversor wedge_qv(ga_versor q, ga_vector v){ return (ga_oddversor){ mul_sv(q.s, v), wedge_Bv(q.B, v) }; }
ga_bivector wedge_Bq(ga_bivector B, ga_versor q){ return mul_Bs(B, q.s); }
SYMMETRIC_OP(wedge_qB, wedge_Bq, ga_versor, ga_bivector, ga_bivector)
ga_multivector wedge_MM(ga_multivector M, ga_multivector N)
{
return (ga_multivector){
.s = M.s*N.s,
.v = add(mul_vs(M.v, N.s), mul_sv(M.s, N.v)),
.B = add(add(mul_Bs(M.B, N.s), mul_sB(M.s, N.B)), wedge_vv(M.v, N.v)),
.P = (ga_pseudoscalar){M.P.e123*N.s + M.s*N.P.e123 + add(wedge_vB(M.v, N.B), wedge_Bv(M.B, N.v)).e123}
};
}
ga_multivector wedge_vM(ga_vector v, ga_multivector N)
{
return (ga_multivector){ .v = mul_vs(v, N.s), .B = wedge_vv(v, N.v), .P = wedge_vB(v, N.B).e123 };
}
ga_multivector wedge_Mv(ga_multivector M, ga_vector v)
{
return (ga_multivector){ .v = mul_sv(M.s, v), .B = wedge_vv(M.v, v), .P = wedge_Bv(M.B, v).e123 };
}
ga_multivector wedge_BM(ga_bivector B, ga_multivector N) { return (ga_multivector){ .B = mul_Bs(B, N.s), .P = wedge_Bv(B, N.v) }; }
ga_multivector wedge_MB(ga_multivector M, ga_bivector B) { return (ga_multivector){ .B = mul_sB(M.s, B), .P = wedge_vB(M.v, B) }; }
ga_pseudoscalar wedge_PM(ga_pseudoscalar P, ga_multivector N) { return mul(P, N.s); }
ga_pseudoscalar wedge_MP(ga_multivector M, ga_pseudoscalar P) { return mul(P, M.s); }
ga_multivector wedge_qM(ga_versor q, ga_multivector N){ return add(mul_sM(q.s, N), wedge_BM(q.B, N)); }
ga_multivector wedge_Mq(ga_multivector M, ga_versor q){ return add(mul_Ms(M, q.s), wedge_MB(M, q.B)); }
ga_multivector wedge_oM(ga_oddversor o, ga_multivector N){ return add(wedge_vM(o.v, N), wedge_PM(o.P, N)); }
ga_multivector wedge_Mo(ga_multivector M, ga_oddversor o){ return add(wedge_Mv(M, o.v), wedge_MP(M, o.P)); }
ga_bivector wedge_vo(ga_vector v, ga_oddversor o){ return wedge_vv(v, o.v); }
ga_bivector wedge_ov(ga_oddversor o, ga_vector v){ return wedge_vv(o.v, v); }
ga_pseudoscalar wedge_Bo(ga_bivector B, ga_oddversor o) { return wedge_Bv(B, o.v); }
ga_pseudoscalar wedge_oB(ga_oddversor o, ga_bivector B) { return wedge_vB(o.v, B); }
ga_versor wedge_qq(ga_versor q, ga_versor r){ return add(mul_sq(q.s, r), wedge_Bq(q.B, r)); }
ga_bivector wedge_oo(ga_oddversor o, ga_oddversor r){ return wedge_vo(o.v, r); }
ga_oddversor wedge_oq(ga_oddversor o, ga_versor q){ return add(mul_os(o, q.s), wedge_oB(o, q.B)); }
ga_oddversor wedge_qo(ga_versor q, ga_oddversor o){ return add(mul_so(q.s, o), wedge_Bo(q.B, o)); }
#define wedge(x, y) _Generic((x), \
ga_scalar: \
_Generic((y), ga_scalar: mul_ss, ga_vector: mul_sv, ga_bivector: mul_sB, ga_versor: mul_sq, ga_oddversor: mul_so, ga_pseudoscalar: mul_sP, ga_multivector: mul_sM), \
ga_vector: \
_Generic((y), ga_scalar: mul_vs, ga_vector: wedge_vv, ga_bivector: wedge_vB, ga_versor: wedge_vq, ga_oddversor: wedge_vo, ga_pseudoscalar: always0_v, ga_multivector: wedge_vM), \
ga_bivector: \
_Generic((y), ga_scalar: mul_Bs, ga_vector: wedge_Bv, ga_bivector: always0_B, ga_versor: wedge_Bq, ga_oddversor: wedge_Bo, ga_pseudoscalar: always0_B, ga_multivector: wedge_BM), \
ga_versor: \
_Generic((y), ga_scalar: mul_qs, ga_vector: wedge_qv, ga_bivector: wedge_qB, ga_versor: wedge_qq, ga_oddversor: wedge_qo, ga_pseudoscalar: always0_q, ga_multivector: wedge_qM), \
ga_oddversor: \
_Generic((y), ga_scalar: mul_os, ga_vector: wedge_ov, ga_bivector: wedge_oB, ga_versor: wedge_oq, ga_oddversor: wedge_oo, ga_pseudoscalar: always0_o, ga_multivector: wedge_oM), \
ga_pseudoscalar: \
_Generic((y), ga_scalar: mul_Ps, ga_vector: always0_P, ga_bivector: always0_P, ga_versor: always0_P, ga_oddversor: always0_P, ga_pseudoscalar: always0_P, ga_multivector: wedge_PM), \
ga_multivector: \
_Generic((y), ga_scalar: mul_Ms, ga_vector: wedge_Mv, ga_bivector: wedge_MB, ga_versor: wedge_Mq, ga_oddversor: wedge_Mo, ga_pseudoscalar: wedge_MP, ga_multivector: wedge_MM) \
)(x, y)
ga_scalar norm_sqr_s(ga_scalar s){ return s*s; }
ga_scalar norm_sqr_v(ga_vector v){ return v.e1*v.e1 + v.e2*v.e2 + v.e3*v.e3; }
ga_scalar norm_sqr_B(ga_bivector B){ return B.e23*B.e23 + B.e12*B.e12 + B.e31*B.e31; }
ga_scalar norm_sqr_P(ga_pseudoscalar P){ return P.e123*P.e123; }
ga_scalar norm_sqr_q(ga_versor q){ return q.s*q.s + norm_sqr_B(q.B); }
ga_scalar norm_sqr_o(ga_oddversor o){ return norm_sqr_v(o.v) + o.P.e123*o.P.e123; }
ga_scalar norm_sqr_M(ga_multivector M)
{
return M.s*M.s + norm_sqr_v(M.v) + norm_sqr_B(M.B) + M.P.e123*M.P.e123;
}
ga_scalar norm_sqr_study(ga_multivector S)
{
return add(sqr_s(S.s), sqr_s(S.P.e123));
}
#define norm_study(x) sqrt(norm_sqr_study(x))
#define norm_sqr(x) _Generic((x),\
ga_scalar: norm_sqr_s, \
ga_vector: norm_sqr_v, \
ga_bivector: norm_sqr_B, \
ga_versor: norm_sqr_q, \
ga_oddversor: norm_sqr_o, \
ga_pseudoscalar: norm_sqr_P, \
ga_multivector: norm_sqr_M)(x)
#define norm(x) (sqrt(norm_sqr(x)))
#define normalize(x) mul(x, 1.0/norm(x))
ga_scalar inv_s(ga_scalar s){ return 1.0/s; }
ga_vector inv_v(ga_vector v){ return mul(v, 1.0/norm_sqr(v)); }
ga_bivector inv_B(ga_bivector B){ return mul(B, -1.0/norm_sqr(B)); }
ga_pseudoscalar inv_P(ga_pseudoscalar P){ return (ga_pseudoscalar){-1.0/P.e123}; }
ga_versor inv_q(ga_versor q){ return mul(conj(q), 1.0/norm_sqr(q)); }
ga_oddversor inv_o(ga_oddversor o){ return mul(conj(o), 1.0/norm_sqr(o)); }
ga_multivector study_conj(ga_multivector S){ return add(S.s, negate(S.P)); }
ga_multivector inv_study(ga_multivector S){ return mul_Ms(study_conj(S), 1.0/mul_MM(S, study_conj(S)).s); }
ga_multivector inv_M(ga_multivector M)
{
// this could be optimized, probably not going to be used much, so it's fine for now.
ga_multivector cM = conj(M);
// this will be a study number which is a ga_scalar + ga_pseudoscalar.
// usually the pseudoscalar == 0, but not always ( v + B )
ga_multivector S = mul_MM(M, cM);
ga_multivector cS = study_conj(S);
// the inverse of the study number
ga_scalar bottom = mul_MM(S, cS).s;
ga_multivector iS = mul_Ms(cS, 1.0/bottom);
return mul_MM(cM, iS);
}
ga_multivector inv_p2(ga_biparavector b)
{
return inv_M(to_multivector(b));
}
#define inv(x) _Generic((x),\
ga_scalar: inv_s, \
ga_vector: inv_v, \
ga_bivector: inv_B, \
ga_versor: inv_q, \
ga_oddversor: inv_o, \
ga_pseudoscalar: inv_P, \
ga_multivector: inv_M, \
ga_biparavector: inv_p2)(x)
#define divide(x, y) mul(x, inv(y))
//------- metric stuff
// this would change the most when moving to a new GA. i.e. PGA needs special treatment for the projective basis
// this is where most libraries diverge. choice of dual and inner product are just that, choices.
// we are going with right and left complements in addition to a hodge star. star should conform to differential geometry and not the PGA folks idea of it.
// right complement of k-blades a and b is defined as a∧rc(b) = ‖a‖₂‖b‖₂I where I is the pseudoscalar. left complement is similarly defined.
// in G(3) they are all identical. not even going to use the word dual, since it's such a confused term.
// regressive product will be defined as a V b = lc(rc(a) ^ rc(b))
//
// inner products. The main inner product which we call dot will be defined using the right complement and regressive product.
// a∨rc(b), lc(a)∨b, a∨⋆b
// for G(3) a∨⋆b == a∨rc(b)
//
// dot(a, b) = rdot(a, b) = a∨⋆b and ldot(a, b) = ⋆a∨b
ga_pseudoscalar rc_s(ga_scalar s){ return (ga_pseudoscalar){s}; }
ga_bivector rc_v(ga_vector v){ return (ga_bivector){v.e1, v.e2, v.e3}; }
ga_vector rc_B(ga_bivector B){ return (ga_vector){B.e23, B.e31, B.e12}; }
ga_scalar rc_P(ga_pseudoscalar P){ return (ga_scalar){P.e123}; }
ga_oddversor rc_q(ga_versor q){ return add(rc_B(q.B), rc_s(q.s)); }
ga_versor rc_o(ga_oddversor o){ return add(rc_P(o.P), rc_v(o.v)); }
ga_multivector rc_M(ga_multivector M){ return (ga_multivector){rc_P(M.P), rc_B(M.B), rc_v(M.v), rc_s(M.s)}; }
// right complement
#define rc(x) _Generic((x), \
int: rc_s, \
ga_scalar: rc_s, \
ga_vector: rc_v, \
ga_bivector: rc_B, \
ga_versor: rc_q, \
ga_oddversor: rc_o, \
ga_pseudoscalar: rc_P, \
ga_multivector: rc_M \
)(x)
#define right_complement(x) rc(x)
#define lc(x) rc(x)
#define left_complement(x) lc(x)
#define star(x) rc(x)
#define regressive(a, b) lc(wedge(rc(a), rc(b)))
#define vee(a, b) regressive(a, b)
#define dot(x, y) (regressive(x, star(y)))
#define rdot(x, y) (regressive(x, star(y)))
#define ldot(x, y) (regressive(star(x), y))
// left contraction and symmetric inner product are included for completeness. skipped the right contraction as I've never seen anyone use it
// these are really hard on the preprocessor and can take some time to compile. you get precise return types though.
// could give these the _Generic(x) treatement for better compile times in the future.
#define gs0(a) grade_select(a, grade0())
#define gs1(a) grade_select(a, grade1())
#define gs2(a) grade_select(a, grade2())
#define gs3(a) grade_select(a, grade3())
#define lconk(a, b) grade_select(mul(a, b), sub_grade(grade(b), grade(a)))
#define left_contraction(a, b) add( \
add( \
mul(gs0(a), b), \
add( \
lconk(gs1(a), gs1(b)), \
add(lconk(gs1(a), gs2(b)), lconk(gs1(a), gs3(b))) \
) \
), \
add( \
add(lconk(gs2(a), gs2(b)), lconk(gs2(a), gs3(b))), \
lconk(gs3(a), gs3(b)) \
) \
)
#define lcontract(x, y) left_contraction(x, y)
// the symmetric inner product.
/* horrifically slow, seems to possibly crash the compiler even, had to comment out.
#define sdotk(a, b) grade_select(mul(a, b), abs_grade(sub_grade(grade(b), grade(a))))
#define sdot(a, b) add( \
add( \
add(mul(gs0(a), b), mul(a, gs0(b))), \
add( \
add(sdotk(gs1(a), gs1(b)), sdotk(gs1(a), gs2(b))), \
add(sdotk(gs2(a), gs1(b)), sdotk(gs2(a), gs2(b))) \
) \
), \
add( \
add(add(sdotk(gs1(a), gs3(b)), sdotk(gs2(a), gs3(b))), \
add(sdotk(gs3(a), gs1(b)), sdotk(gs3(a), gs2(b))) \
), \
sdotk(gs3(a), gs3(b)) \
) \
)
#define symmetric_inner_product(x, y) sdot(x, y)
*/
//---------- trancendental maps: exp log sqrt
ga_paravector exp_v(ga_vector v){ return (ga_paravector){ .s = cosh(norm(v)), .v=mul(normalize(v), norm(v)) }; }
ga_versor exp_B(ga_bivector B)
{
ga_scalar lamb = norm(B);
return (ga_versor){.s = cos(lamb), .B = mul(normalize(B), sin(lamb))};
}
ga_versor exp_q(ga_versor q){ return mul(exp(q.s), exp_B(q.B)); }
ga_paravector exp_p(ga_paravector p)
{
ga_paravector ev = exp_v(p.v);
ga_scalar es = exp(p.s);
return (ga_paravector){.s = es*ev.s, .v = mul(es, ev.v)};
}
ga_multivector exp_P(ga_pseudoscalar P)
{
return (ga_multivector){.s = cos(P.e123), .v = zero_v(), .B = zero_B(), .P = (ga_pseudoscalar){sin(P.e123)}};
}
ga_multivector exp_o(ga_oddversor o)
{
ga_multivector ev = to_multivector(exp_v(o.v));
ga_multivector eP = exp_P(o.P);
return mul_MM(ev, eP);
}
typedef struct ga_biparavector_pair
{
ga_biparavector a;
ga_biparavector b;
} ga_biparavector_pair;
// use a trick: the isomorphism between vectors in G(3) and a bivector in G(3,1)
// map xe1 + ye2 + ze3 -> xe14 + ye24 + ze34. Then from Graded Symmetry Groups paper we can decompose that into commuting bivectors
// then map back to a vector + bivector.
ga_biparavector_pair factorize(ga_biparavector m)
{
// vector and bivector commute with everything except each other, so do the trick of mapping 3D v + B to a 4D bivector with metric ( +++- )
// this works since e14*e14 = 1 like e1*e1 = 1
// just working with coordinates here rather than actually implementing the G(3,1) algebra
// v + B -> B4
// B4⋅B4
ga_scalar B4B4_1 = dot_vv(m.v, m.v) - (m.B.e23*m.B.e23 + m.B.e31*m.B.e31 + m.B.e12*m.B.e12);
// B4∧B4
ga_scalar B4B4_4 = 2.0*(m.v.e1*m.B.e23 + m.v.e2*m.B.e31 + m.v.e3*m.B.e12);
// note the + sign in second term, that's because B4∧B4 is e1234 of +++- signature which squares to -
ga_scalar lambda_i = sqrt(sqr_s(B4B4_1) + sqr_s(B4B4_4));
// inverse in G(3) of v + B same coords as G(3,1)
ga_multivector bottom = inv(add(m.v, m.B));
// move pseudoscalar back to G(3). magically it seems to work out. probably bc both are study numbers of their space or something.
ga_pseudoscalar B4B4_3 = {B4B4_4};
ga_multivector top = mul_Ms(add(B4B4_1 + lambda_i, B4B4_3), 0.5);
ga_multivector C = mul_MM(top, bottom);
top = mul_Ms(add(B4B4_1 - lambda_i, B4B4_3), 0.5);
ga_multivector D = mul_MM(top, bottom);
return (ga_biparavector_pair){ (ga_biparavector){C.v, C.B}, (ga_biparavector){D.v, D.B} };
}
// everything commutes in the exp aside from vector+bivector.
// To handle v + B we use a trick: the isomorphism between vectors in G(3) and a bivector in G(3,1)
// map xe1 + ye2 + ze3 -> xe14 + ye24 + ze34. Then from Graded Symmetry Groups paper we can decompose that into commuting bivectors
// solve for the exp map there and project back to a vector + bivector.
ga_multivector exp_M(ga_multivector m)
{
// can use exp(a + b) = exp(a)*exp(b) for all commuting k-vectors
ga_multivector eP = exp_P(m.P);
ga_scalar es = exp(m.s);
ga_biparavector_pair CD = factorize( (ga_biparavector){ .v = m.v, .B = m.B } );
// seems kinda wild that the v+B are now commuting? wonder if we need to stay in G(3,1) when composing the exp?
// My REPL session checks out for v+B's commuting though.
// original B = C + D, but C and D commute.
// however, the inner parts C.v and C.B don't commute with each other in G(3), so go back to G(3,1) for exp(C) and exp(D)
// norms are weird in this hyperbolic space. norm squared of vector part will be - since e14*e41 = -1 ( v*~v )
ga_scalar C_norm = -norm_sqr(CD.a.v) + norm_sqr(CD.a.B);
ga_multivector expC;
if( C_norm < 0.0 )
{
C_norm = sqrt(-C_norm);
expC = (ga_multivector){.s = cosh(C_norm), .v = mul(CD.a.v, sinh(C_norm)/C_norm), .B = mul( CD.a.B, sinh(C_norm)/C_norm)};
}
else
{
C_norm = sqrt(C_norm);
expC = (ga_multivector){.s = cos(C_norm), .v = mul(CD.a.v, sin(C_norm)/C_norm), .B = mul( CD.a.B, sin(C_norm)/C_norm)};
}
ga_scalar D_norm = -norm_sqr(CD.b.v) + norm_sqr(CD.b.B);
ga_multivector expD;
if( D_norm < 0.0 )
{
D_norm = sqrt(-D_norm);
expD = (ga_multivector){.s = cosh(D_norm), .v = mul(CD.b.v, sinh(D_norm)/D_norm), .B = mul( CD.b.B, sinh(D_norm)/D_norm)};
}
else
{
D_norm = sqrt(D_norm);
expD = (ga_multivector){.s = cos(D_norm), .v = mul(CD.b.v, sin(D_norm)/D_norm), .B = mul( CD.b.B, sin(D_norm)/D_norm)};
}
// combine them all exp(s + Bplus + Bminus + P) == exp(s)*exp(Bplus)*exp(Bminus)*exp(P)
return mul_Ms(mul_MM(expC, mul_MM(expD, eP)), es);
}
#define exp(x) _Generic((x), \
ga_scalar: exp, \
ga_vector: exp_v, \
ga_bivector: exp_B, \
ga_versor: exp_q, \
ga_oddversor: exp_o, \
ga_pseudoscalar: exp_P, \
ga_paravector: exp_p \
)(x)
//-------- sqrt
ga_multivector sqrt_study(ga_multivector S)
{
ga_scalar c = sqrt(0.5*(S.s + norm_study(S)));
return add(c, mul(0.5/c, S.P));
}
ga_versor sqrt_q(ga_versor q)
{
ga_scalar n = norm_sqr(q);
if( n < 0.9999 || n > 1.0001)
{
n = sqrt(n);
q = mul(1.0/n, q);
q.s += sign_s(q.s);
q = normalize(q);
return mul(sqrt(n), q);
}
else
{
q.s += sign_s(q.s);
q = normalize(q);
return q;
}
}
ga_versor sqrt_B(ga_bivector B)
{
return sqrt_q(add(0.0, B));
}
ga_multivector sqrt_P(ga_pseudoscalar P)
{
return sqrt_study(add(0.0, P));
}
#define sqrt(x) _Generic((x), \
ga_scalar: sqrt, \
ga_bivector: sqrt_B, \
ga_versor: sqrt_q, \
ga_pseudoscalar: sqrt_P \
)(x)
//----------- log
// could do log_o and log_P but it's not so simple. Would need to go to G(3,1) and then back. Leaving them out for now.
ga_versor log_q(ga_versor q)
{
ga_scalar qn = norm_sqr(q);
if( qn < 0.9999 || qn > 1.0001 )
{
qn = sqrt(qn);
q = mul(q, 1.0/qn);
return add(log(qn), mul(normalize(q.B), acos(q.s)));
}
else
{
return add(0.0, mul(normalize(q.B), acos(q.s)));
}
}
ga_versor log_B(ga_bivector B)
{
return log_q(add(0.0, B));
}
#define log(x) _Generic((x), \
ga_scalar: log, \
ga_bivector: log_B, \
ga_versor: log_q \
)(x)
//----- utilities
#include <stdio.h>
void print_s(ga_scalar s){ printf("%f", s); }
void print_v(ga_vector v)
{
printf("%fe₁ + %fe₂ + %fe₃", v.e1, v.e2, v.e3);
}
void print_B(ga_bivector B)
{
printf("%fe₂₃ + %fe₃₁ + %fe₁₂", B.e23, B.e31, B.e12);
}
void print_q(ga_versor q)
{
printf("%f + %fe₂₃ + %fe₃₁ + %fe₁₂", q.s, q.B.e23, q.B.e31, q.B.e12);
}
void print_P(ga_pseudoscalar P){ printf("%fe₁₂₃", P.e123); }
void print_o(ga_oddversor o)
{
print_v(o.v);
printf(" + ");
print_P(o.P);
}
void print_M(ga_multivector M)
{
int preceding = 0;
if( M.s != 0 )
{
preceding = 1;
print_s(M.s);
}
if( norm(M.v) != 0 )
{
if( preceding ) printf(" + ");
print_v(M.v);
preceding = 1;
}
if( norm(M.B) != 0 )
{
if( preceding ) printf(" + ");
print_B(M.B);
preceding = 1;
}
if( norm(M.P) != 0 )
{
if( preceding ) printf(" + ");
print_P(M.P);
}
}
#define ga_print(a) _Generic((a), \
ga_scalar: print_s, ga_vector: print_v, ga_bivector: print_B, ga_versor: print_q, ga_oddversor: print_o, ga_pseudoscalar: print_P, ga_multivector: print_M \
)(a)
#define ga_println(a) {print(a); printf("\n");}
//------ TESTS
// remove the ga_ prefix for convenience
/*
typedef ga_scalar scalar;
typedef ga_vector vector;
typedef ga_bivector bivector;
typedef ga_versor versor;
typedef ga_versor quaternion;
typedef ga_pseudoscalar pseudoscalar;
#define print(x) ga_print(x)
#define println(x) ga_println(x)
#define vec(x,y,z) ga_vec(x,y,z)
#define bivec(x,y,z) ga_bivec(x,y,z)
#define rotor(s,x,y,z) ga_rotor(s,x,y,z)
// https://github.com/sheredom/utest.h
#include "utest.h"
// utest with ga3 expressions that make heavy use of the preprocessor are death for compilers.
// i.e. EXPECT_EQ(lcontract(y, wedge(x, y)).e1, -1); // this will take forever to compile. call lcontract outside of EXPECT_EQ
#include <stdlib.h>
#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif
ga_scalar rand_scalar( ga_scalar n ){ return ((ga_scalar)rand() / ((ga_scalar)RAND_MAX + 1)) * 2.0*n - n; }
ga_vector rand_vector( ga_scalar n ){ return (ga_vector){ rand_scalar(n), rand_scalar(n), rand_scalar(n) }; }
ga_bivector rand_bivector( ga_scalar n ){ return (ga_bivector){ rand_scalar(n), rand_scalar(n), rand_scalar(n) }; }
UTEST(ga3, add)
{
vector x = {1, 0, 0};
vector y = {0, 2, 0};
EXPECT_EQ(sub(x,x).e1, 0);
EXPECT_EQ(add(x,y).e1, add(y,x).e1);
EXPECT_EQ(add(x, (ga_empty){}).e1, 1);
EXPECT_EQ(add((ga_empty){}, x).e1, 1);
}
UTEST(ga3, wedge)
{
vector x = {1, 0, 0};
vector y = {0, 1, 0};
vector z = {0, 0, 1};
pseudoscalar P = {2.0};
EXPECT_EQ(wedge((scalar)(3), x).e1, 3);
EXPECT_EQ(wedge(x, (scalar)(3)).e1, 3);
EXPECT_EQ(wedge(x, y).e12, 1);
EXPECT_EQ(wedge(z, y).e23, -1);
EXPECT_EQ(wedge(P, x), 0);
EXPECT_EQ(wedge(wedge(x, y), wedge(y, z)), 0);
EXPECT_EQ(wedge(wedge(x, y), z).e123, 1);
}
UTEST(ga3, operators)
{
scalar s = 42.0;
vector v = {1, 2, 3};
vector u = {0.1, -0.2, 0.4};
bivector B = {2, 4, -6};
bivector C = {-0.2, 1, 2};
pseudoscalar P = {2.1};
pseudoscalar Q = {4.1};
EXPECT_NEAR(0.0, norm_sqr(sub(mul(B, inv(v)), divide(B, v))), 1e-4);
ga_multivector M = add(s, add(v, add( B, P)));
ga_multivector N = add(s, add(u, add( C, Q)));
ga_multivector MN = mul(M, N);
// can't nest much more than this or gcc bogs.
ga_multivector multi_check = add(
add(
add( mul(s, s), add( mul(s, u), add( mul(s, C), mul(s, Q)))),
add( mul(v, s), add( mul(v, u), add( mul(v, C), mul(v, Q))))),
add(
add( mul(B, s), add( mul(B, u), add( mul(B, C), mul(B, Q)))),
add( mul(P, s), add( mul(P, u), add( mul(P, C), mul(P, Q))))));
EXPECT_NEAR(0.0, sub(MN.s, multi_check.s), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.v, multi_check.v)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.B, multi_check.B)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.P, multi_check.P)), 1e-4);
multi_check = add(
add(mul(add(s, B), add(s, C)), mul(add(v, P), add(u, Q)) ),
add(mul(add(v, P), add(s, C)), mul(add(s, B), add(u, Q)) ));
EXPECT_NEAR(0.0, sub(MN.s, multi_check.s), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.v, multi_check.v)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.B, multi_check.B)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.P, multi_check.P)), 1e-4);
multi_check = add(
add(
add( wedge(s, s), add( wedge(s, u), add( wedge(s, C), wedge(s, Q)))),
add( wedge(v, s), add( wedge(v, u), add( wedge(v, C), wedge(v, Q))))),
add(
add( wedge(B, s), add( wedge(B, u), add( wedge(B, C), wedge(B, Q)))),
add( wedge(P, s), add( wedge(P, u), add( wedge(P, C), wedge(P, Q))))));
MN = wedge(M, N);
EXPECT_NEAR(0.0, sub(MN.s, multi_check.s), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.v, multi_check.v)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.B, multi_check.B)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.P, multi_check.P)), 1e-4);
multi_check = add(
add(wedge(add(s, B), add(s, C)), wedge(add(v, P), add(u, Q)) ),
add(wedge(add(v, P), add(s, C)), wedge(add(s, B), add(u, Q)) ));
EXPECT_NEAR(0.0, sub(MN.s, multi_check.s), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.v, multi_check.v)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.B, multi_check.B)), 1e-4);
EXPECT_NEAR(0.0, norm_sqr(sub(MN.P, multi_check.P)), 1e-4);
}
UTEST(ga3, rotor)
{
vector xy = {1, 1, 0};
vector x = {1, 0, 0};
vector v = {0,2,1};
versor q = normalize(mul(xy, inv(x)));
EXPECT_NEAR(norm_sqr(q), 1.0, 1.0e-4);
vector qq = sandwich(q, v);
EXPECT_NEAR(qq.e1, -2.0, 1.0e-4);
EXPECT_NEAR(qq.e3, 1.0, 1.0e-4);
}
UTEST(ga3, barycentric)
{
vector a = {2,1,0};
vector b = {4,1,0};
vector c = {3,2,0};
vector p = {2.1,1.5,1};
bivector abc = wedge(sub(b, a), sub(c, a));
bivector abp = wedge(sub(b, a), sub(p, a));
bivector bcp = wedge(sub(c, b), sub(p, b));
bivector cap = wedge(sub(a, c), sub(p, c));
ga_scalar barycentric_total = grade_select(mul(abp, inv(abc)), grade0()) +
grade_select(mul(bcp, inv(abc)), grade0()) +
grade_select(mul(cap, inv(abc)), grade0());
EXPECT_NEAR(barycentric_total, 1.0, 1.0e-4);
}
UTEST(ga3, duals)
{
vector a = {1,2,3};
EXPECT_EQ(wedge(lc(a), a).e123, 14);
EXPECT_EQ(wedge(a, rc(a)).e123, 14);
vector x = {1,0,0};
vector y = {0,1,0};
vector z = {0,0,1};
EXPECT_EQ(star(x).e23, 1.0);
EXPECT_EQ(star(y).e31, 1.0);
EXPECT_EQ(star(z).e12, 1.0);
EXPECT_EQ(star(star(2)), 2.0);
vector b = {-0.707, 0.707, 1.1};
ga_scalar ans = norm_sqr(sub( rc(vee(a, b)), wedge(rc(a), rc(b)) ));
EXPECT_NEAR(ans, 0.0, 1.0e-4);
ans = norm_sqr(sub( lc(vee(a, b)), wedge(lc(a), lc(b)) ));
EXPECT_NEAR(ans, 0.0, 1.0e-4);
ans = norm_sqr(sub( vee(a, b), lc(wedge(rc(a), rc(b))) ));
EXPECT_NEAR(ans, 0.0, 1.0e-4);
ans = vee(rc(x), vee(rc(y), rc(z)));
EXPECT_EQ(ans, 1);
EXPECT_NEAR(wedge(x, wedge(y, z)).e123, 1.0, 1.0e-4);
}
UTEST(ga3, inner)
{
vector a = {1,2,3};
bivector B = star(a);
EXPECT_EQ(dot(a, a), 14);
EXPECT_EQ(dot(B, B), 14);
vector x = {1,0,0};
vector y = {0,1,0};
vector z = {0,0,1};
EXPECT_EQ(rdot(wedge(x, y), x).e2, 1);
EXPECT_EQ(rdot(wedge(x, y), y).e1, -1);
EXPECT_EQ(rdot(wedge(x, y), wedge(x, y)), 1);
EXPECT_EQ(rdot(3, 2), 6);
// contractions are very slow to compile. heavy preprocessor load. precise return type is nice though
// especially don't put it inside a utest macro, that will kill the compiler.
ga_scalar yxy_e1 = -1;//lcontract(y, wedge(x, y)).e1;
EXPECT_EQ(yxy_e1, -1);
// this is much faster.
yxy_e1 = grade_select(mul(y, wedge(x, y)), grade1()).e1;
EXPECT_EQ(yxy_e1, -1);
}
UTEST(ga3, exp)
{
//ga_biparavector b = { .v = rand_vector(4.0), .B = rand_bivector(4.0) };
ga_biparavector b = { .v = {1,2,3}, .B = {0.333,0.666,1.333} };
ga_biparavector_pair commutes = factorize(b);
ga_multivector bm = to_multivector(b);
ga_multivector b1m = to_multivector(commutes.a);
ga_multivector b2m = to_multivector(commutes.b);
ga_multivector ibm = inv(bm);
EXPECT_NEAR( norm(mul_MM(bm, ibm)), 1.0, 1.0e-4 );
EXPECT_NEAR( norm(sub(bm, add(b1m, b2m))), 0.0, 1.0e-4 );
EXPECT_NEAR( norm(sub(mul_MM(b1m, b2m), mul_MM(b2m, b1m))), 0.0, 1.0e-4 );
// 60 rotation two ways.
bivector B = {.e12 = M_PI/6.0 };
ga_versor qexp = exp(B);
// vector division version create a rotor that upon sandwich rotates twice as far.
// so give it a vector rotated 30 as target
vector x30 = {.e1 = 1.0*sqrt(3.0), .e2 = -1.0};
vector x = {2,0,0};
ga_versor qdiv = normalize(mul(x30, inv(x))); // ccw +ve rotation
ga_versor qdiv2 = mul(normalize(x30), normalize(inv(x)));
ga_versor q30 = exp((bivector){.e12 = M_PI/12.0 });
EXPECT_NEAR( 0.0, norm(sub(q30, sqrt(qexp))), 1.0e-4);
EXPECT_NEAR( 0.0, norm(sub(qdiv, qexp)), 1.0e-4);
EXPECT_NEAR( 0.0, norm(sub(qdiv, qdiv2)), 1.0e-4);
vector xrot = sandwich(sqrt(qexp), x);
EXPECT_NEAR( 0.0, norm(sub(xrot, x30)), 1.0e-4);
EXPECT_NEAR(0.0, norm(sub(log(qexp), B)), 1.0e-4);
B = (bivector){1, 2, 3};
EXPECT_NEAR(0.0, norm(sub(exp(log(B)), B)), 1.0e-4);
}
UTEST(ga3, sqrt)
{
ga_multivector S = add(1.2, (ga_pseudoscalar){2.3});
ga_multivector Ssqrt = sqrt_study(S);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), S)), 1.0e-4);
S = add(1.2, (ga_pseudoscalar){-2.3});
Ssqrt = sqrt_study(S);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), S)), 1.0e-4);
S = add(-1.2, (ga_pseudoscalar){2.3});
Ssqrt = sqrt_study(S);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), S)), 1.0e-4);
S = add(-1.2, (ga_pseudoscalar){-2.3});
Ssqrt = sqrt_study(S);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), S)), 1.0e-4);
S = add(rand_scalar(10), (ga_pseudoscalar){rand_scalar(10)});
Ssqrt = sqrt_study(S);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), S)), 1.0e-4);
ga_pseudoscalar P = {2.0};
Ssqrt = sqrt(P);
EXPECT_NEAR( 0.0, norm_study(sub(mul_MM(Ssqrt, Ssqrt), P)), 1.0e-4);
ga_bivector B = {1,2,3};
ga_versor Bsqrt = sqrt(B);
EXPECT_NEAR( 0.0, norm(sub(mul(Bsqrt, Bsqrt), B)), 1.0e-4);
}
UTEST(ga3, print)
{
ga_versor q = normalize(add(1.0, (ga_bivector){ .e12 = 1.0 }));
println(q);
EXPECT_TRUE(1);
}
UTEST_MAIN()
*/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment