From 1cad35536664c1996fe69ac7a630ea4c749d0739 Mon Sep 17 00:00:00 2001 From: Bill Currie Date: Wed, 30 Aug 2023 15:00:27 +0900 Subject: [PATCH] [qfcc] Write tests for 3d PGA geometric products And fix the flood of errors. --- tools/qfcc/source/expr_algebra.c | 88 ++++++++++++------ tools/qfcc/test/pga3d.r | 155 ++++++++++++++++++++++++++++++- 2 files changed, 210 insertions(+), 33 deletions(-) diff --git a/tools/qfcc/source/expr_algebra.c b/tools/qfcc/source/expr_algebra.c index 2ee35774a..031d94918 100644 --- a/tools/qfcc/source/expr_algebra.c +++ b/tools/qfcc/source/expr_algebra.c @@ -977,12 +977,12 @@ pga3_x_y_z_w_geom_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto va = offset_cast (vtype, a, 0); auto sa = offset_cast (stype, a, 3); - auto cv = scale_expr (geom_type, b, sa); + auto cv = scale_expr (vtype, b, sa); auto cs = dot_expr (stype, va, b); - c[0] = cross_expr (algebra_mvec_type (alg, 0x01), b, va); - c[5] = sum_expr (geom_type, '+', new_extend_expr (cv, geom_type, 0, false), - new_extend_expr (cs, geom_type, 0, true)); - c[5]->e.expr.type = geom_type; + c[0] = new_extend_expr (cross_expr (vtype, b, va), + algebra_mvec_type (alg, 0x01), 0, false); + c[5] = sum_expr (geom_type, '-', new_extend_expr (cs, geom_type, 0, true), + new_extend_expr (cv, geom_type, 0, false)); } static void @@ -993,8 +993,8 @@ pga3_x_y_z_w_geom_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto va = offset_cast (vtype, a, 0); auto cs = unary_expr ('-', dot_expr (stype, va, b)); - c[0] = cross_expr (algebra_mvec_type (alg, 0x01), va, b); - c[5] = new_extend_expr (cs, geom_type, 0, true); + c[5] = new_extend_expr (cross_expr (vtype, va, b), geom_type, 0, false); + c[0] = new_extend_expr (cs, algebra_mvec_type (alg, 0x01), 0, true); } static void @@ -1005,8 +1005,8 @@ pga3_x_y_z_w_geom_wxyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto va = offset_cast (vtype, a, 0); auto sb = offset_cast (stype, b, 0); - auto cv = scale_expr (geom_type, va, sb); - c[5] = new_extend_expr (cv, geom_type, 0, true); + auto cv = scale_expr (vtype, va, sb); + c[5] = new_extend_expr (cv, geom_type, 0, false); } static void @@ -1020,7 +1020,7 @@ pga3_x_y_z_w_geom_wzy_wxz_wyx_xyz (expr_t **c, expr_t *a, expr_t *b, auto vb = offset_cast (vtype, b, 0); auto sb = offset_cast (stype, b, 3); c[1] = scale_expr (algebra_mvec_type (alg, 0x02), va, sb); - c[3] = cross_expr (algebra_mvec_type (alg, 0x80), vb, va); + c[3] = cross_expr (algebra_mvec_type (alg, 0x08), vb, va); c[4] = dot_expr (geom_type, a, b); } @@ -1032,9 +1032,10 @@ pga3_yz_zx_xy_geom_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto vb = offset_cast (vtype, b, 0); auto sb = offset_cast (stype, b, 3); - auto cv = scale_expr (geom_type, a, sb); + auto cv = scale_expr (vtype, a, sb); auto cs = dot_expr (stype, vb, a); - c[0] = cross_expr (algebra_mvec_type (alg, 0x01), a, vb); + c[0] = new_extend_expr (cross_expr (vtype, vb, a), + algebra_mvec_type (alg, 0x01), 0, false); c[5] = sum_expr (geom_type, '-', new_extend_expr (cs, geom_type, 0, true), new_extend_expr (cv, geom_type, 0, false)); } @@ -1053,7 +1054,6 @@ pga3_yz_zx_xy_geom_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) c[4] = dot_expr (algebra_mvec_type (alg, 0x10), b, a); } -#define pga3_wxyz_geom_yz_zx_xy pga3_yz_zx_xy_geom_wxyz static void pga3_yz_zx_xy_geom_wxyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) { @@ -1062,7 +1062,6 @@ pga3_yz_zx_xy_geom_wxyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) c[3] = scale_expr (algebra_mvec_type (alg, 0x08), a, unary_expr ('-', sb)); } -#define pga3_wzy_wxz_wyx_xyz_geom_yz_zx_xy pga3_yz_zx_xy_geom_wzy_wxz_wyx_xyz static void pga3_yz_zx_xy_geom_wzy_wxz_wyx_xyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) @@ -1072,11 +1071,12 @@ pga3_yz_zx_xy_geom_wzy_wxz_wyx_xyz (expr_t **c, expr_t *a, expr_t *b, auto geom_type = algebra_mvec_type (alg, 0x01); auto vb = offset_cast (vtype, b, 0); auto sb = offset_cast (stype, b, 3); - auto cv = scale_expr (geom_type, a, sb); + auto cv = scale_expr (vtype, a, sb); auto cs = dot_expr (stype, vb, a); c[0] = sum_expr (geom_type, '-', new_extend_expr (cs, geom_type, 0, true), new_extend_expr (cv, geom_type, 0, false)); - c[5] = cross_expr (algebra_mvec_type (alg, 0x20), vb, a); + c[5] = new_extend_expr (cross_expr (vtype, vb, a), + algebra_mvec_type (alg, 0x20), 0, false); } static void @@ -1087,14 +1087,15 @@ pga3_wx_wy_wz_geom_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto vb = offset_cast (vtype, b, 0); auto cs = dot_expr (stype, vb, a); - c[0] = cross_expr (algebra_mvec_type (alg, 0x01), vb, a); - c[5] = new_extend_expr (cs, geom_type, 0, true); + c[5] = new_extend_expr (cross_expr (vtype, vb, a), + geom_type, 0, false); + c[0] = new_extend_expr (cs, algebra_mvec_type (alg, 0x01), 0, true); } static void pga3_wx_wy_wz_geom_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) { - c[3] = cross_expr (algebra_mvec_type (alg, 0x08), a, b); + c[3] = cross_expr (algebra_mvec_type (alg, 0x08), b, a); c[4] = dot_expr (algebra_mvec_type (alg, 0x10), a, b); } @@ -1103,10 +1104,11 @@ pga3_wx_wy_wz_geom_wzy_wxz_wyx_xyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) { auto stype = alg->type; + auto vtype = vector_type (stype, 3); auto geom_type = algebra_mvec_type (alg, 0x20); auto vs = offset_cast (stype, b, 3); - auto cv = scale_expr (geom_type, a, unary_expr ('-', vs)); - c[5] = new_extend_expr (cv, geom_type, 0, true); + auto cv = scale_expr (vtype, a, unary_expr ('-', vs)); + c[5] = new_extend_expr (cv, geom_type, 0, false); } static void @@ -1117,8 +1119,16 @@ pga3_wxyz_geom_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) auto geom_type = algebra_mvec_type (alg, 0x20); auto sa = offset_cast (stype, a, 0); auto vb = offset_cast (vtype, b, 0); - auto cv = scale_expr (geom_type, vb, unary_expr ('-', sa)); - c[5] = new_extend_expr (cv, geom_type, 0, true); + auto cv = scale_expr (vtype, vb, unary_expr ('-', sa)); + c[5] = new_extend_expr (cv, geom_type, 0, false); +} + +static void +pga3_wxyz_geom_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) +{ + auto stype = alg->type; + auto sa = offset_cast (stype, a, 0); + c[3] = scale_expr (algebra_mvec_type (alg, 0x08), b, unary_expr ('-', sa)); } static void @@ -1144,8 +1154,25 @@ pga3_wzy_wxz_wyx_xyz_geom_x_y_z_w (expr_t **c, expr_t *a, expr_t *b, auto sa = offset_cast (stype, a, 3); auto vb = offset_cast (vtype, b, 0); c[1] = scale_expr (algebra_mvec_type (alg, 0x02), vb, sa); - c[3] = cross_expr (algebra_mvec_type (alg, 0x08), vb, va); - c[0] = unary_expr ('-', dot_expr (geom_type, b, a)); + c[3] = cross_expr (algebra_mvec_type (alg, 0x08), va, vb); + c[4] = unary_expr ('-', dot_expr (geom_type, b, a)); +} + +static void +pga3_wzy_wxz_wyx_xyz_geom_yz_zx_xy (expr_t **c, expr_t *a, expr_t *b, + algebra_t *alg) +{ + auto stype = alg->type; + auto vtype = vector_type (stype, 3); + auto geom_type = algebra_mvec_type (alg, 0x01); + auto va = offset_cast (vtype, a, 0); + auto sa = offset_cast (stype, a, 3); + auto cv = scale_expr (vtype, b, sa); + auto cs = dot_expr (stype, va, b); + c[0] = sum_expr (geom_type, '-', new_extend_expr (cs, geom_type, 0, true), + new_extend_expr (cv, geom_type, 0, false)); + c[5] = new_extend_expr (cross_expr (vtype, b, va), + algebra_mvec_type (alg, 0x20), 0, false); } static void @@ -1153,10 +1180,11 @@ pga3_wzy_wxz_wyx_xyz_geom_wx_wy_wz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) { auto stype = alg->type; + auto vtype = vector_type (stype, 3); auto geom_type = algebra_mvec_type (alg, 0x20); - auto vs = offset_cast (stype, b, 3); - auto cv = scale_expr (geom_type, a, vs); - c[5] = new_extend_expr (cv, geom_type, 0, true); + auto sa = offset_cast (stype, a, 3); + auto cv = scale_expr (vtype, b, sa); + c[5] = new_extend_expr (cv, geom_type, 0, false); } static void @@ -1164,8 +1192,8 @@ pga3_wzy_wxz_wyx_xyz_geom_wxyz (expr_t **c, expr_t *a, expr_t *b, algebra_t *alg) { auto stype = alg->type; - auto sa = offset_cast (stype, a, 0); - auto sb = offset_cast (stype, b, 3); + auto sa = offset_cast (stype, a, 3); + auto sb = offset_cast (stype, b, 0); auto geom_type = algebra_mvec_type (alg, 0x01); auto cs = scale_expr (stype, sa, sb); c[0] = new_extend_expr (cs, geom_type, 0, true); diff --git a/tools/qfcc/test/pga3d.r b/tools/qfcc/test/pga3d.r index 52a7d2d5b..cbe5e8347 100644 --- a/tools/qfcc/test/pga3d.r +++ b/tools/qfcc/test/pga3d.r @@ -14,9 +14,12 @@ typedef union { } bivector_t; typedef union { PGA.group_mask(0x1e) mvec; + bivector_t bvec; struct { - bivector_t bvec; + PGA.group_mask(0x02) dir; scalar_t scalar; + PGA.group_mask(0x08) mom; + quadvector_t qvec; }; } evengrades_t; typedef union { @@ -204,8 +207,7 @@ test_dot (void) return 1; } - b.dir = nil;//FIXME allow bvec = mom - b.mom = qvec • bvec.bvec; + b.bvec = qvec • bvec.bvec; if ((vec3)b.dir != '0 0 0' || (vec3)b.mom != '-88 32 -16') { printf ("qvec • bvec != '0 0 0' '-88 32 -16': %v %v\n", b.dir, b.mom); return 1; @@ -238,6 +240,136 @@ test_dot (void) return 0; } +static int +test_geom (void) +{ + scalar_t scalar; + 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 + 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 (qvec * qvecb != 0) { + printf ("(qvec * qvecb) != 0': %g\n", qvec * qvecb); + return 1; + } + + auto a = tvec * qvec; + if ((vec4)a != '0 0 0 -16') { + printf ("vec * vecb != '0 0 0 -16': %v\n", a); + return 1; + } + + bivector_t b = { .mom = bvec.bvec * qvec }; + if ((vec3)b.dir != '0 0 0' || (vec3)b.mom != '-88 32 -16') { + printf ("bvec * qvec != '0 0 0' '-88 32 -16': %v %v\n", b.dir, b.mom); + return 1; + } + + auto c = vec * qvec; + if ((vec4)c != '24 -16 32 0') { + printf ("vec * qvec != '24 -16 32 0': %v\n", c); + return 1; + } + + oddgrades_t d = { .mvec = bvec.bvec * tvec }; + if ((vec4)d.vec != '22 -8 4 9' || (vec4)d.tvec != '30 85 -34 0') { + printf ("bvec * tvec != '22 -8 4 9' '30 85 -34 0': %q %q\n", + d.vec, d.tvec); + return 1; + } + + evengrades_t e = { .mvec = vec * tvec }; + if ((vec3)e.dir != '-6 4 -8' || (vec3)e.mom != '30 17 -14' + || e.scalar || (scalar_t)e.qvec != 21) { + printf ("vec * tvec != 0 '-6 4 -8' '30 17 -14' 21: %g %v %v %g\n", + e.scalar, e.bvec.dir, e.bvec.mom, (scalar_t) e.qvec); + return 1; + } + + e.mvec = bvec.bvec * bvec.bvec; + if (e.scalar != -141 || (vec3)e.bvec.dir || (vec3)e.bvec.mom + || (scalar_t)e.qvec != -78) { + printf ("bvec * bvec != -141 '0 0 0' '0 0 0' -78: %g %v %v %g\n", + e.scalar, e.bvec.dir, e.bvec.mom, (scalar_t)e.qvec); + return 1; + } + + d = { .mvec = vec * bvec.bvec }; + if ((vec4)d.vec != '-12 -38 -10 -9' || (vec4)d.tvec != '-45 -29 7 49') { + printf ("vec * bvec != '-12 -38 -10 -9' '-45 -29 7 49': %v %v\n", + d.vec, d.tvec); + return 1; + } + + e.mvec = vec * vecb; + if (e.scalar != -9 || (scalar_t)e.qvec + || (vec3)e.bvec.dir != '-48 20 46' || (vec3)e.bvec.mom != '44 -14 52') { + printf ("(vec * vecb) != -9 '-48 20 46' '44 -14 52' 0': %g %v %v %g\n", + e.scalar, e.bvec.dir, e.bvec.mom, (scalar_t) e.qvec); + return 1; + } + + a = qvec * tvec; + if ((vec4)a != '0 0 0 16') { + printf ("vec * vecb != '0 0 0 16': %v\n", a); + return 1; + } + + e.mvec = qvec * bvec.bvec; + if (e.scalar || (scalar_t)e.qvec + || (vec3)e.bvec.dir || (vec3)e.bvec.mom != '-88 32 -16') { + printf ("(vec * vecb) != 0 '0 0 0' '-88 32 -16' 0': %g %v %v %g\n", + e.scalar, e.bvec.dir, e.bvec.mom, (scalar_t) e.qvec); + return 1; + } + + c = qvec * vec; + if ((vec4)c != '-24 16 -32 0') { + printf ("qvec * vec != '-24 16 -32 0': %q\n", c); + return 1; + } + + d = { .mvec = tvec * bvec.bvec }; + if ((vec4)d.vec != '22 -8 4 9' || (vec4)d.tvec != '-30 -85 34 0') { + printf ("vec * bvec != '22 -8 4 9' '-30 -85 34 0': %q %q\n", + d.vec, d.tvec); + return 1; + } + + e.mvec = tvec * vec; + if (e.scalar || (scalar_t)e.qvec != -21 + || (vec3)e.bvec.dir != '-6 4 -8' || (vec3)e.bvec.mom != '30 17 -14') { + printf ("tvec * vec != 0 '-6 4 -8' '30 17 -14' -21: %g %v %v %g\n", + e.scalar, e.bvec.dir, e.bvec.mom, (scalar_t) e.qvec); + return 1; + } + + d.mvec = bvec.bvec * vec; + if ((vec4)d.vec != '12 38 10 9' || (vec4)d.tvec != '-45 -29 7 49') { + printf ("vec * bvec != '12 38 10 9' '-45 -29 7 49': %q %q\n", + d.vec, d.tvec); + return 1; + } + + return 0; +} + int main (void) { @@ -267,6 +399,18 @@ main (void) printf ("quadvector has wrong size: %d\n", sizeof (quadvector_t)); return 1; } + if (sizeof (evengrades_t) != sizeof (bivector_t)) { + evengrades_t e; + #define offsetof(s,f) ((int)(&((s*)0).f)) + printf ("evengrades has wrong size: %d\n", sizeof (evengrades_t)); + printf ("mvec: %d, bvec: %d, scalar: %d, qvec: %d\n", + sizeof (e.mvec), sizeof (e.bvec), + sizeof (e.scalar), sizeof (e.qvec)); + printf ("mvec: %d, bvec: %d, scalar: %d, qvec: %d\n", + offsetof (evengrades_t, mvec), offsetof (evengrades_t, bvec), + offsetof (evengrades_t, scalar), offsetof (evengrades_t, qvec)); + return 1; + } if (test_wedge ()) { printf ("wedge products failed\n"); @@ -278,5 +422,10 @@ main (void) return 1; } + if (test_geom ()) { + printf ("geometric products failed\n"); + return 1; + } + return 0; // to survive and prevail :) }