[qfcc] Write tests for 3d PGA wedge products

And fix the discovered bugs.
This commit is contained in:
Bill Currie 2023-08-29 19:42:00 +09:00
parent 33295a8ad9
commit c42aaea060
2 changed files with 92 additions and 18 deletions

View file

@ -715,7 +715,6 @@ pga3_x_y_z_w_wedge_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
scale_expr (wedge_type, va, sb));
}
#define pga3_yz_zx_xy_wedge_x_y_z_w pga3_x_y_z_w_wedge_yz_zx_xy
static void
pga3_x_y_z_w_wedge_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
{
@ -724,7 +723,7 @@ pga3_x_y_z_w_wedge_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
auto wedge_type = algebra_mvec_type (alg, 0x20);
auto va = offset_cast (vtype, a, 0);
auto sa = offset_cast (stype, a, 3);
auto cv = scale_expr (get_type (b), b, sa);
auto cv = scale_expr (vtype, b, sa);
auto cs = dot_expr (stype, va, b);
c[5] = sum_expr (wedge_type, '-',
@ -732,8 +731,6 @@ pga3_x_y_z_w_wedge_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
new_extend_expr (cv, wedge_type, 0, false));
}
// vector-bivector wedge is commutative
#define pga3_wx_wy_wz_wedge_x_y_z_w pga3_x_y_z_w_wedge_wx_wy_wz
static void
pga3_x_y_z_w_wedge_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
{
@ -742,10 +739,7 @@ pga3_x_y_z_w_wedge_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
auto wedge_type = algebra_mvec_type (alg, 0x20);
auto va = offset_cast (vtype, a, 0);
auto cv = cross_expr (vtype, va, b);
auto cs = new_zero_expr (stype);
c[5] = sum_expr (wedge_type, '+',
new_extend_expr (cs, wedge_type, 0, true),
new_extend_expr (cv, wedge_type, 0, false));
c[5] = new_extend_expr (cv, wedge_type, 0, false);
}
static void
@ -755,6 +749,22 @@ pga3_x_y_z_w_wedge_wzy_wxz_wyx_xyz (expr_t **c, expr_t *a, expr_t *b,
c[4] = dot_expr (algebra_mvec_type (alg, 0x10), a, b);
}
static void
pga3_yz_zx_xy_wedge_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
{
auto stype = alg->type;
auto vtype = vector_type (stype, 3);
auto wedge_type = algebra_mvec_type (alg, 0x20);
auto vb = offset_cast (vtype, b, 0);
auto sb = offset_cast (stype, b, 3);
auto cv = scale_expr (vtype, a, sb);
auto cs = dot_expr (stype, vb, a);
c[5] = sum_expr (wedge_type, '-',
new_extend_expr (cs, wedge_type, 0, true),
new_extend_expr (cv, wedge_type, 0, false));
}
// bivector-bivector wedge is commutative
#define pga3_wx_wy_wz_wedge_yz_zx_xy pga3_yz_zx_xy_wedge_wx_wy_wz
static void
@ -763,6 +773,17 @@ pga3_yz_zx_xy_wedge_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
c[4] = dot_expr (algebra_mvec_type (alg, 0x10), a, b);
}
static void
pga3_wx_wy_wz_wedge_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg)
{
auto stype = alg->type;
auto vtype = vector_type (stype, 3);
auto wedge_type = algebra_mvec_type (alg, 0x20);
auto vb = offset_cast (vtype, b, 0);
auto cv = cross_expr (vtype, vb, a);
c[5] = new_extend_expr (cv, wedge_type, 0, false);
}
static void
pga3_wzy_wxz_wyx_xyz_wedge_x_y_z_w (expr_t **c, expr_t *a, expr_t *b,
algebra_t *alg)

View file

@ -58,42 +58,95 @@ main (void)
}
scalar_t scalar;
vector_t vec, vecb;
vector_t vec, vecb, vecc, vecd;
bivector_t bvec, bvecb;
trivector_t tvec, tvecb;
quadvector_t qvec, qvecb;
@algebra (PGA) {
scalar = 42;
vec = 3*e1 - 2*e2 + e0;
bvec.bvec = 4*e20 - 3*e01 + 2*e12;
tvec = 7*e012;
vec = 3*e1 - 2*e2 + 4*e3 + e0;
bvec.bvec = -3*e01 + 5*e02 + 7*e03 + 11*e23 - 4*e31 + 2*e12;
tvec = 1*e032 + 4*e013 + 7*e021 - 2*e123;
qvec = 8*e0123;
vecb = 5*e1 + 12*e2 - 13*e0;
bvecb.bvec = 6*e20 + 4*e01 + 1*e12;
tvecb = 3*e032;
qvecb = 1*e0123;
vecc = 5*e1 + 4*e2 - 3*e3 + 4*e0;
vecd = e1 + e2 + e3 + e0;
}
if (scalar != 42) {
printf ("scalar != 42: %g\n", scalar);
return 1;
}
if ((vec4)vec != '3 -2 0 1') {
printf ("vec != '3 -2 0 1': %q\n", vec);
if ((vec4)vec != '3 -2 4 1') {
printf ("vec != '3 -2 4 1': %q\n", vec);
return 1;
}
if ((vec3)bvec.dir != '0 0 2' || (vec3)bvec.mom != '-3 -4 0') {
printf ("bvec != '0 0 2' '-3 -4 0': %v %v\n", bvec.dir, bvec.mom);
if ((vec3)bvec.dir != '11 -4 2' || (vec3)bvec.mom != '-3 5 7') {
printf ("bvec != '11 -4 2' '-3 5 7': %v %v\n", bvec.dir, bvec.mom);
return 1;
}
if ((vec4)tvec != '0 0 -7 0') {
printf ("tvec != '0 0 -7': %g\n", tvec);
if ((vec4)tvec != '1 4 7 -2') {
printf ("tvec != '1 4 7 -2': %g\n", tvec);
return 1;
}
if ((scalar_t)qvec != 8) {
printf ("tvec != 8: %g\n", (scalar_t) qvec);
return 1;
}
bivector_t a = { .bvec = vecvecb };
if ((vec3)a.dir != '-48 20 46' || (vec3)a.mom != '44 -14 52') {
printf ("vec ∧ vecb != '-48 20 46' '44 -14 52': %v %v\n", a.dir, a.mom);
return 1;
}
auto b = a.bvecvecc;
if ((vec4)b != '358 -472 -430 -298') {
printf ("a ∧ vecc != '358 -472 -430 -298': %v\n", b);
return 1;
}
auto c = bvecd;
if ((scalar_t)c != 842) {
printf ("b ∧ vecd != 742': %g\n", (scalar_t) c);
return 1;
}
a.bvec = vecbvec;
if ((vec3)a.dir != '48 -20 -46' || (vec3)a.mom != '-44 14 -52') {
printf ("vecb ∧ vec != '48 -20 -46' '-44 14 -52': %v %v\n",
a.dir, a.mom);
return 1;
}
b = vecca.bvec;
if ((vec4)b != '-358 472 430 298') {
printf ("vecc ∧ a != '-358 472 430 298': %v\n", b);
return 1;
}
c = vecdb;
if ((scalar_t)c != 842) {
printf ("vecd ^ b != 842': %g\n", (scalar_t) c);
return 1;
}
c = a.bvec(veccvecd);
if ((scalar_t)c != -842) {
printf ("a ∧ (vecc ∧ vecd) != -742': %g\n", (scalar_t) c);
return 1;
}
c = (vecdvecc)a.bvec;
if ((scalar_t)c != 842) {
printf ("(vecd ∧ vecc) ∧ a != 742': %g\n", (scalar_t) c);
return 1;
}
return 0; // to survive and prevail :)
}