[qfcc] Implement undual and use for regressive product

I'm not sure the regressive product is right (overall sign), but that's
actually partly a problem in the math itself (duals and the regressive
product still get poked at, so it may be just a matter of
interpretation).

tools/qfcc/include/algebra.h

@@ -115,6 +115,7 @@ pr_uint_t get_group_mask (const type_t *type, algebra_t *algebra) __attribute__(
 const expr_t *algebra_binary_expr (int op, const expr_t *e1, const expr_t *e2);
 const expr_t *algebra_negate (const expr_t *e);
 const expr_t *algebra_dual (const expr_t *e);
+const expr_t *algebra_undual (const expr_t *e);
 const expr_t *algebra_reverse (const expr_t *e);
 const expr_t *algebra_cast_expr (type_t *dstType, const expr_t *e);
 const expr_t *algebra_assign_expr (const expr_t *dst, const expr_t *src);

tools/qfcc/source/expr.c

@@ -2058,6 +2058,8 @@ bitnot_expr:
 			return algebra_reverse (e);
 		case DUAL:
 			return algebra_dual (e);
+		case UNDUAL:
+			return algebra_undual (e);
 	}
 	internal_error (e, 0);
 }

tools/qfcc/source/expr_algebra.c

@@ -2596,7 +2596,7 @@ regressive_product (const expr_t *e1, const expr_t *e2)
 	auto a = algebra_dual (e1);
 	auto b = algebra_dual (e2);
 	auto c = outer_product (a, b);
-	return algebra_dual (c);
+	return algebra_undual (c);
 }
 
 static const expr_t *
@@ -2762,12 +2762,13 @@ algebra_negate (const expr_t *e)
 	return mvec_gather (n, algebra);
 }
 
-const expr_t *
-algebra_dual (const expr_t *e)
+static const expr_t *
+hodge_dual (const expr_t *e, bool undual)
 {
 	auto algebra = algebra_context (get_type (e));
 	if (!algebra) {
-		return error (e, "cannot take the dual of a scalar without context");
+		return error (e, "cannot take the %s of a scalar without context",
+					  undual ? "undual" : "dual");
 	}
 	auto layout = &algebra->layout;
 
@@ -2776,6 +2777,7 @@ algebra_dual (const expr_t *e)
 	e = mvec_expr (e, algebra);
 	mvec_scatter (a, e, algebra);
 
+	int dim = algebra->dimension - algebra->zero;
 	pr_uint_t I_mask = (1u << algebra->dimension) - 1;
 	for (int i = 0; i < layout->count; i++) {
 		if (!a[i]) {
@@ -2784,21 +2786,39 @@ algebra_dual (const expr_t *e)
 		auto group = &layout->groups[i];
 		//FIXME assumes groups are mono-grade (either come up with something
 		//or reject mixed-grade groups)
-		pr_uint_t mask = I_mask ^ group->blades[0].mask;
-		int dual_ind = layout->group_map[layout->mask_map[mask]][0];
+		auto blade = group->blades[0];
+		pr_uint_t d_mask = I_mask ^ blade.mask;
+		int dual_ind = layout->group_map[layout->mask_map[d_mask]][0];
 		auto dual_group = &layout->groups[dual_ind];
 		auto dual_type = algebra_mvec_type (algebra, dual_group->group_mask);
 		auto dual = cast_expr (dual_type, a[i]);
-		if (algebra_count_flips (algebra, group->blades[0].mask, mask) & 1) {
-			b[dual_ind] = neg_expr (dual);
+		pr_uint_t mask = algebra_count_flips (algebra, blade.mask, d_mask);
+		if (mask & 1) {
+			dual = neg_expr (dual);
 		} else {
-			b[dual_ind] = dual;
+			dual = dual;
 		}
+		if (undual && ((dim * algebra_get_grade (get_type (a[i]))) & 1)) {
+			dual = neg_expr (dual);
+		}
+		b[dual_ind] = dual;
 	}
 
 	return mvec_gather (b, algebra);
 }
 
+const expr_t *
+algebra_dual (const expr_t *e)
+{
+	return hodge_dual (e, false);
+}
+
+const expr_t *
+algebra_undual (const expr_t *e)
+{
+	return hodge_dual (e, true);
+}
+
 static void
 set_sign (pr_type_t *val, int sign, const type_t *type)
 {

tools/qfcc/source/qc-lex.l

@@ -443,6 +443,8 @@ static keyword_t qf_keywords[] = {
 	{"@regressive",	REGRESSIVE,	},
 	{"@geometric",	GEOMETRIC,	},
 	{"@algebra",	ALGEBRA,	},
+	{"@dual",		DUAL,		},
+	{"@undual",		UNDUAL,		},
 };
 
 // These keywors are always available. Other than the @ keywords, they

tools/qfcc/source/qc-parse.y

@@ -145,7 +145,7 @@ int yylex (void);
 %left			'+' '-'
 %left			'*' '/' '%' MOD SCALE GEOMETRIC QMUL QVMUL VQMUL
 %left           HADAMARD CROSS DOT WEDGE REGRESSIVE
-%right	<op>	SIZEOF UNARY INCOP REVERSE STAR DUAL
+%right	<op>	SIZEOF UNARY INCOP REVERSE STAR DUAL UNDUAL
 %left			HYPERUNARY
 %left			'.' '(' '['
 
@@ -1730,6 +1730,7 @@ unary_expr
 	| unary_expr INCOP				{ $$ = incop_expr ($2, $1, 1); }
 	| unary_expr REVERSE			{ $$ = unary_expr (REVERSE, $1); }
 	| DUAL cast_expr %prec UNARY	{ $$ = unary_expr (DUAL, $2); }
+	| UNDUAL cast_expr %prec UNARY	{ $$ = unary_expr (UNDUAL, $2); }
 	| '+' cast_expr %prec UNARY	{ $$ = $2; }
 	| '-' cast_expr %prec UNARY	{ $$ = unary_expr ('-', $2); }
 	| '!' cast_expr %prec UNARY	{ $$ = unary_expr ('!', $2); }

tools/qfcc/test/pga2d.r

@@ -20,6 +20,96 @@ typedef union {
 	};
 } oddgrades_t;
 
+static int
+test_dual (void)
+{
+	@algebra (PGA) {
+		auto a = 1f;
+		auto a0 = e0;
+		auto a1 = e1;
+		auto a2 = e2;
+		auto a12 = e12;
+		auto a01 = e01;
+		auto a20 = e20;
+		auto a012 = e012;
+		#define TEST_DUAL(x, y) \
+			if (⋆x != y) { \
+				printf ("⋆" #x " != " #y "\n"); \
+				return 1; \
+			}
+		TEST_DUAL (a, e012);
+		TEST_DUAL (a0, e12);
+		TEST_DUAL (a1, e20);
+		TEST_DUAL (a2, e01);
+		TEST_DUAL (a20, e1);
+		TEST_DUAL (a01, e2);
+		TEST_DUAL (a12, e0);
+		TEST_DUAL (a012, 1);
+		#undef TEST_DUAL
+
+		#define TEST_DUAL(x) \
+			if (x * ⋆x != e012) { \
+				printf (#x " * ⋆" #x " != e012\n"); \
+				return 1; \
+			}
+		TEST_DUAL (a);
+		TEST_DUAL (a0);
+		TEST_DUAL (a1);
+		TEST_DUAL (a2);
+		TEST_DUAL (a20);
+		TEST_DUAL (a01);
+		TEST_DUAL (a12);
+		TEST_DUAL (a012);
+		#undef TEST_DUAL
+	}
+	return 0;
+}
+
+static int
+test_undual (void)
+{
+	@algebra (PGA) {
+		auto a = 1f;
+		auto a0 = e0;
+		auto a1 = e1;
+		auto a2 = e2;
+		auto a12 = e12;
+		auto a01 = e01;
+		auto a20 = e20;
+		auto a012 = e012;
+		#define TEST_UNDUAL(x, y) \
+			if (@undual x != y) { \
+				printf ("@undual " #x " != " #y "\n"); \
+				return 1; \
+			}
+		TEST_UNDUAL (a, e012);
+		TEST_UNDUAL (a0, e12);
+		TEST_UNDUAL (a1, e20);
+		TEST_UNDUAL (a2, e01);
+		TEST_UNDUAL (a20, e1);
+		TEST_UNDUAL (a01, e2);
+		TEST_UNDUAL (a12, e0);
+		TEST_UNDUAL (a012, 1);
+		#undef TEST_UNDUAL
+
+		#define TEST_UNDUAL(x) \
+			if (@undual x * x != e012) { \
+				printf ("@undual " #x " * " #x " != e012\n"); \
+				return 1; \
+			}
+		TEST_UNDUAL (a);
+		TEST_UNDUAL (a0);
+		TEST_UNDUAL (a1);
+		TEST_UNDUAL (a2);
+		TEST_UNDUAL (a20);
+		TEST_UNDUAL (a01);
+		TEST_UNDUAL (a12);
+		TEST_UNDUAL (a012);
+		#undef TEST_DUAL
+	}
+	return 0;
+}
+
 int
 main (void)
 {
@@ -269,7 +359,17 @@ main (void)
 	}
 	auto line = bvec ∨ bvecb;
 	if ((dvec3)line != '-11 8 34'd) {
-		printf ("bvec ∨ bvecb != '-11 8 34': %lv\n", line);
+		printf ("bvec(%lv) ∨ bvecb(%lv) != '-11 8 34': %lv\n", bvec, bvecb, line);
+		return 1;
+	}
+
+	if (test_dual ()) {
+		printf ("dual failed\n");
+		return 1;
+	}
+
+	if (test_undual ()) {
+		printf ("dual failed\n");
 		return 1;
 	}
 	return 0;		// to survive and prevail :)

tools/qfcc/test/pga3d.r

@@ -434,6 +434,76 @@ test_dual (void)
 		TEST_DUAL (a021);
 		TEST_DUAL (a123);
 		TEST_DUAL (a0123);
+		#undef TEST_DUAL
+	}
+	return 0;
+}
+
+static int
+test_undual (void)
+{
+	@algebra (PGA) {
+		auto a = 1f;
+		auto a0 = e0;
+		auto a1 = e1;
+		auto a2 = e2;
+		auto a3 = e3;
+		auto a23 = e23;
+		auto a31 = e31;
+		auto a12 = e12;
+		auto a01 = e01;
+		auto a02 = e02;
+		auto a03 = e03;
+		auto a032 = e032;
+		auto a013 = e013;
+		auto a021 = e021;
+		auto a123 = e123;
+		auto a0123 = e0123;
+		#define TEST_UNDUAL(x, y) \
+			if (@undual x != y) { \
+				printf ("@undual" #x " != " #y "\n"); \
+				return 1; \
+			}
+		TEST_UNDUAL (a, e0123);
+		TEST_UNDUAL (a0, -e123);
+		TEST_UNDUAL (a1, -e032);
+		TEST_UNDUAL (a2, -e013);
+		TEST_UNDUAL (a3, -e021);
+		TEST_UNDUAL (a23, e01);
+		TEST_UNDUAL (a31, e02);
+		TEST_UNDUAL (a12, e03);
+		TEST_UNDUAL (a01, e23);
+		TEST_UNDUAL (a02, e31);
+		TEST_UNDUAL (a03, e12);
+		TEST_UNDUAL (a032, e1);
+		TEST_UNDUAL (a013, e2);
+		TEST_UNDUAL (a021, e3);
+		TEST_UNDUAL (a123, e0);
+		TEST_UNDUAL (a0123, 1);
+		#undef TEST_UNDUAL
+
+		#define TEST_UNDUAL(x) \
+			if (@undual x * x != e0123) { \
+				printf ("@undual " #x " * " #x " != e0123\n"); \
+				return 1; \
+			}
+		TEST_UNDUAL (a);
+		TEST_UNDUAL (a0);
+		TEST_UNDUAL (a1);
+		TEST_UNDUAL (a2);
+		TEST_UNDUAL (a3);
+		TEST_UNDUAL (a23);
+		TEST_UNDUAL (a31);
+		TEST_UNDUAL (a12);
+		TEST_UNDUAL (a01);
+		TEST_UNDUAL (a02);
+		TEST_UNDUAL (a03);
+		TEST_UNDUAL (a032);
+		TEST_UNDUAL (a013);
+		TEST_UNDUAL (a021);
+		TEST_UNDUAL (a123);
+		TEST_UNDUAL (a0123);
+		#undef TEST_UNDUAL
 	}
 	return 0;
 }
@@ -549,6 +619,11 @@ main (void)
 		return 1;
 	}
 
+	if (test_undual ()) {
+		printf ("dual failed\n");
+		return 1;
+	}
+
 	if (test_basics ()) {
 		printf ("basics failed\n");
 		return 1;