[gatest] Add basic geometric algebra test

This is for developing methods of implementing geometric algebra and
eventually playing with it visually.

ruamoko/Makemodule.am

@@ -6,5 +6,6 @@ include ruamoko/lib/Makemodule.am
 include ruamoko/game/Makemodule.am
 include ruamoko/gui/Makemodule.am
 include ruamoko/cl_menu/Makemodule.am
+include ruamoko/gatest/Makemodule.am
 include ruamoko/scheme/Makemodule.am
 include ruamoko/qwaq/Makemodule.am

ruamoko/gatest/Makemodule.am

@@ -0,0 +1,35 @@
+noinst_PROGRAMS += ruamoko/gatest/gatest.dat$(EXEEXT)
+
+gatest_dat_src= \
+	ruamoko/gatest/algebra.r \
+	ruamoko/gatest/basisblade.r \
+	ruamoko/gatest/basisgroup.r \
+	ruamoko/gatest/basislayout.r \
+	ruamoko/gatest/main.r \
+	ruamoko/gatest/metric.r \
+	ruamoko/gatest/multivector.r \
+	ruamoko/gatest/util.r
+
+ruamoko_gatest_gatest_dat_SOURCES=$(gatest_dat_src)
+ruamoko_gatest_gatest_obj=$(ruamoko_gatest_gatest_dat_SOURCES:.r=.o)
+ruamoko_gatest_gatest_dep=$(call qcautodep,$(ruamoko_gatest_gatest_dat_SOURCES:.o=.Qo))
+ruamoko/gatest/gatest.dat$(EXEEXT): $(ruamoko_gatest_gatest_obj) $(QFCC_DEP) $(libui) ruamoko/lib/libr.a
+	$(V_QFCCLD)$(QLINK) -o $@ $(ruamoko_gatest_gatest_obj) $(libui) -lr
+include $(ruamoko_gatest_gatest_dep) # am--include-marker
+r_depfiles_remade += $(ruamoko_gatest_gatest_dep)
+
+EXTRA_PROGRAMS += \
+	ruamoko/gatest/gatest
+
+EXTRA_DIST += \
+	$(gatest_dat_src) \
+	ruamoko/gatest/algebra.h \
+	ruamoko/gatest/basisblade.h \
+	ruamoko/gatest/basisgroup.h \
+	ruamoko/gatest/basislayout.h \
+	ruamoko/gatest/metric.h \
+	ruamoko/gatest/multivector.h \
+	ruamoko/gatest/util.h
+CLEANFILES += \
+	ruamoko/gatest/*.dat \
+	ruamoko/gatest/*.sym

ruamoko/gatest/algebra.h

@@ -0,0 +1,28 @@
+#ifndef __algebra_h
+#define __algebra_h
+
+#include "Object.h"
+
+@class Metric;
+@class BasisGroup;
+@class BasisLayout;
+
+@interface Algebra : Object
+{
+	Metric     *metric;
+	BasisGroup **grades;
+	BasisLayout *layout;
+	int         num_components;
+	int         dimension;
+}
++(Algebra *) R:(int)p, int m, int z;
++(Algebra *) PGA:(int)n;
+-(void) print;
+-(BasisGroup *)grade:(int)k;
+-(BasisLayout *)layout;
+-(Metric *) metric;
+-(int)count;
+-(int)dimension;
+@end
+
+#endif//__algebra_h

ruamoko/gatest/algebra.r

@@ -0,0 +1,135 @@
+#include <stdlib.h>
+
+#include "algebra.h"
+#include "metric.h"
+#include "basisblade.h"
+#include "basisgroup.h"
+#include "basislayout.h"
+#include "util.h"
+
+@implementation Algebra
++(Algebra *) R:(int)p, int m, int z
+{
+	Algebra    *a = [[[Algebra alloc] init] autorelease];
+	a.metric = [[Metric R:p, m, z] retain];
+
+	int         d = p + m + z;
+	a.dimension = d;
+	a.num_components = 1 << d;
+	BasisBlade **blades = obj_malloc (a.num_components * sizeof (BasisBlade *));
+	int        *counts = binomial (d);
+	int        *indices = obj_malloc ((d + 1) * sizeof (int));
+	
+	indices[0] = 0;
+	for (int i = 0; i < d; i++) {
+		indices[i + 1] = counts[i];
+	}
+	prefixsum (indices, d + 1);
+	for (int i = 0; i < a.num_components; i++) {
+		int         grade = count_bits (i);
+		int         ind = indices[grade]++;
+		blades[ind] = [BasisBlade basis:i];
+	}
+
+	a.grades = obj_malloc ((d + 1) * sizeof (BasisGroup *));
+	for (int i = 0; i < d + 1; i++) {
+		int         c = counts[i];
+		int         ind = indices[i];
+		a.grades[i] = [BasisGroup new:c basis:blades + ind - c];
+	}
+
+	if (p == 3 && m == 0 && z == 1) {
+		// 3d PGA (w squares to 0, x y z square to +1):
+		// : x   y   z   w
+		// : yz  zx  xy  1
+		// : xw  yw  zw  xyzw
+		// : zyw xzw yxw xyz
+		BasisBlade *pga_blades[16] = {
+			blades[1],  blades[2],  blades[3],  blades[4],
+			blades[7],  blades[6],  blades[5],  blades[0],
+			blades[8],  blades[9],  blades[10], blades[15],
+			blades[14], blades[13], blades[12], blades[11],
+		};
+		BasisGroup *pga_groups[4] = {
+			[BasisGroup new:4 basis:pga_blades +  0],
+			[BasisGroup new:4 basis:pga_blades +  4],
+			[BasisGroup new:4 basis:pga_blades +  8],
+			[BasisGroup new:4 basis:pga_blades + 12],
+		};
+		a.layout = [[BasisLayout new:4 groups: pga_groups] retain];
+	} else if (p == 2 && m == 0 && z == 1) {
+		// 2d PGA (w squares to 0, x y square to +1):
+		// : x   y   w  1
+		// : yw  wx  xy xyw
+		BasisBlade *pga_blades[8] = {
+			blades[1], blades[2], blades[3], blades[0],
+			blades[6], blades[5], blades[4], blades[7],
+		};
+		BasisGroup *pga_groups[2] = {
+			[BasisGroup new:4 basis:pga_blades + 0],
+			[BasisGroup new:4 basis:pga_blades + 4],
+		};
+		a.layout = [[BasisLayout new:2 groups: pga_groups] retain];
+	} else {
+		// just use the grades as the default layout
+		a.layout = [[BasisLayout new:d + 1 groups: a.grades] retain];
+	}
+
+	obj_free (indices);
+	obj_free (counts);
+	obj_free (blades);
+	return a;
+}
+
++(Algebra *) PGA:(int)n
+{
+	return [Algebra R:n, 0, 1];
+}
+
+-(void)dealloc
+{
+	obj_free (grades);
+	[metric release];
+	[layout release];
+	[super dealloc];
+}
+
+-(BasisGroup *)grade:(int)k
+{
+	return grades[k];
+}
+
+-(BasisLayout *)layout
+{
+	return layout;
+}
+
+-(Metric *) metric
+{
+	return metric;
+}
+
+-(int)count
+{
+	return num_components;
+}
+
+-(int)dimension
+{
+	return dimension;
+}
+
+-(void) print
+{
+	int count = [layout count];
+	for (int i = 0; i < count; i++) {
+		BasisGroup *g = [layout group:i];
+		int c = [g count];
+		printf ("%d %d %@\n", i, c, [g set]);
+		for (int j = 0; j < c; j++) {
+			printf ("    %@\n", [g bladeAt:j]);
+		}
+	}
+}
+
+@end

ruamoko/gatest/basisblade.h

@@ -0,0 +1,26 @@
+#ifndef __basisblade_h
+#define __basisblade_h
+#include <Object.h>
+
+@class Metric;
+
+@interface BasisBlade : Object
+{
+	unsigned    mask;
+	double      scale;
+}
++(BasisBlade *) scalar:(double) scale;
++(BasisBlade *) zero;
++(BasisBlade *) basis:(unsigned) mask;
++(BasisBlade *) basis:(unsigned) mask scale:(double) scale;
+-(BasisBlade *) product:(BasisBlade *) b isOuter:(int)outer metric:(Metric *) m;
+-(BasisBlade *) outerProduct:(BasisBlade *) b;
+-(BasisBlade *) geometricProduct:(BasisBlade *) b metric:(Metric *) m;
+-(BasisBlade *) geometricProduct:(BasisBlade *) b;
+-(int) grade;
+-(unsigned) mask;
+-(double) scale;
+-(string) name;
+@end
+
+#endif//__basisblade_h

ruamoko/gatest/basisblade.r

@@ -0,0 +1,106 @@
+#include <string.h>
+#include "basisblade.h"
+#include "metric.h"
+#include "util.h"
+
+@implementation BasisBlade : Object
+-(BasisBlade *) initWithMask:(unsigned) mask scale:(double) scale
+{
+	if (!(self = [super init])) {
+		return self;
+	}
+	self.mask = mask;
+	self.scale = scale;
+	return self;
+}
+
++(BasisBlade *) scalar:(double) scale
+{
+	return [[[BasisBlade alloc] initWithMask:0 scale:scale] autorelease];
+}
+
++(BasisBlade *) zero
+{
+	return [[[BasisBlade alloc] initWithMask:0 scale:0] autorelease];
+}
+
++(BasisBlade *) basis:(unsigned) mask
+{
+	return [[[BasisBlade alloc] initWithMask:mask scale:1] autorelease];
+}
+
++(BasisBlade *) basis:(unsigned) mask scale:(double) scale
+{
+	return [[[BasisBlade alloc] initWithMask:mask scale:scale] autorelease];
+}
+
+-(BasisBlade *) product:(BasisBlade *) b isOuter:(int)outer metric:(Metric *)m
+{
+	if (outer && (mask & b.mask)) {
+		// the two blades share at least one basis vector
+		return [BasisBlade zero];
+	}
+	if (!scale || !b.scale) {
+		return [BasisBlade zero];
+	}
+	int         sign = 1 - (-(count_flips (mask, b.mask) & 1) & 2);
+	double      s = scale * b.scale * sign;
+	if (m) {
+		s *= [m apply: mask, b.mask];
+		if (!s) {
+			return [BasisBlade zero];
+		}
+	}
+	return [BasisBlade basis:(mask ^ b.mask) scale:s];
+}
+
+-(BasisBlade *) outerProduct:(BasisBlade *) b
+{
+	return [self product:b isOuter:1 metric:nil];
+}
+
+-(BasisBlade *) geometricProduct:(BasisBlade *) b metric:(Metric *) m
+{
+	return [self product:b isOuter:0 metric:m];
+}
+
+-(BasisBlade *) geometricProduct:(BasisBlade *) b
+{
+	return [self product:b isOuter:0 metric:nil];
+}
+
+-(int) grade
+{
+	return count_bits (mask);
+}
+
+-(unsigned) mask
+{
+	return mask;
+}
+
+-(double) scale
+{
+	return scale;
+}
+
+-(string) name
+{
+	string      basis = "";
+
+	for (int i = 0; i < 32; i++) {
+		if (mask & (1 << i)) {
+			basis += sprintf("%x", i + 1);
+		}
+	}
+	if (basis) {
+		basis = "e" + basis;
+	}
+	return basis;
+}
+
+-(string) describe
+{
+	return sprintf ("%g%s", scale, [self name]);
+}
+@end

ruamoko/gatest/basisgroup.h

@@ -0,0 +1,25 @@
+#ifndef __basisgroup_h
+#define __basisgroup_h
+#include <Object.h>
+
+@class Metric;
+@class BasisBlade;
+@class Set;
+
+@interface BasisGroup : Object
+{
+	int         count;
+	uivec2      range;
+	BasisBlade **blades;
+	int        *map;
+	Set        *set;
+}
++(BasisGroup *) new:(int) count basis:(BasisBlade **)blades;
+-(int)count;
+-(uivec2)blade_range;
+-(BasisBlade *) bladeAt:(int) ind;
+-(BasisBlade *) blade:(unsigned) mask;
+-(Set *) set;
+@end
+
+#endif//__basisgroup_h

ruamoko/gatest/basisgroup.r

@@ -0,0 +1,69 @@
+#include <Set.h>
+#include "basisblade.h"
+#include "basisgroup.h"
+
+@implementation BasisGroup
++(BasisGroup *) new:(int) count basis:(BasisBlade **)blades
+{
+	BasisGroup *group = [[[BasisGroup alloc] init] autorelease];
+	group.blades = obj_malloc (count * sizeof (BasisBlade *));
+	group.set = [[Set set] retain];
+	group.range = { ~0u, 0 };
+	for (int i = 0; i < count; i++) {
+		group.blades[i] = [blades[i] retain];
+		unsigned    m = [blades[i] mask];
+		if (m < group.range[0]) {
+			group.range[0] = m;
+		}
+		if (m > group.range[1]) {
+			group.range[1] = m;
+		}
+		[group.set add:m];
+	}
+	int         num = group.range[1] - group.range[0] + 1;
+	group.map = obj_malloc (num * sizeof (int));
+	for (int i = 0; i < count; i++) {
+		group.map[[blades[i] mask] - group.range[0]] = i;
+	}
+	group.count = count;
+	return group;
+}
+
+-(void)dealloc
+{
+	[set release];
+	for (int i = 0; i < count; i++) {
+		[blades[i] release];
+	}
+	obj_free (blades);
+	[super dealloc];
+}
+
+-(int)count
+{
+	return count;
+}
+
+-(uivec2)blade_range
+{
+	return range;
+}
+
+-(BasisBlade *) bladeAt:(int) ind
+{
+	return blades[ind];
+}
+
+-(BasisBlade *) blade:(unsigned) mask
+{
+	if (![set is_member:mask]) {
+		return nil;
+	}
+	return blades[map[mask - range[0]]];
+}
+
+-(Set *) set
+{
+	return set;
+}
+@end

ruamoko/gatest/basislayout.h

@@ -0,0 +1,28 @@
+#ifndef __basislayout_h
+#define __basislayout_h
+#include <Object.h>
+
+@class BasisGroup;
+@class Set;
+
+@interface BasisLayout : Object
+{
+	int         count;
+	uivec2      range;
+	BasisGroup **groups;
+	ivec3      *group_map;
+	int        *mask_map;
+	int         blade_count;
+	Set        *set;
+}
++(BasisLayout *) new:(int) count groups:(BasisGroup **)groups;
+-(int)count;
+-(int)num_components;
+-(int)blade_count;
+-(BasisGroup *) group:(int) ind;
+-(BasisBlade *) bladeAt:(int) ind;
+-(int) bladeIndex:(unsigned) mask;
+-(Set *) set;
+@end
+
+#endif//__basislayout_h

ruamoko/gatest/basislayout.r

@@ -0,0 +1,111 @@
+#include <stdlib.h>
+#include <Set.h>
+#include "basisblade.h"
+#include "basisgroup.h"
+#include "basislayout.h"
+#include "util.h"
+
+@implementation BasisLayout
++(BasisLayout *) new:(int) count groups:(BasisGroup **)groups
+{
+	BasisLayout *layout = [[[BasisLayout alloc] init] autorelease];
+	layout.count = count;
+	layout.groups = obj_malloc (count * sizeof (BasisGroup *));
+	layout.set = [[Set set] retain];
+	layout.range = { ~0u, 0 };
+	int        *group_base = obj_malloc ((count + 1) * sizeof (int));
+	group_base[0] = 0;
+	int         num_blades = 0;
+	for (int i = 0; i < count; i++) {
+		layout.groups[i] = [groups[i] retain];
+		[layout.set union:[groups[i] set]];
+		group_base[i + 1] = [groups[i] count];
+		num_blades += group_base[i + 1];
+
+		uivec2      r = [groups[i] blade_range];
+		if (r[0] < layout.range[0]) {
+			layout.range[0] = r[0];
+		}
+		if (r[1] > layout.range[1]) {
+			layout.range[1] = r[1];
+		}
+	}
+	prefixsum (group_base, count);
+	layout.blade_count = num_blades;
+	layout.group_map = obj_malloc (num_blades * sizeof (ivec3));
+
+	int         num = layout.range[1] - layout.range[0] + 1;
+	layout.mask_map = obj_malloc (num * sizeof (int));
+	int *group_inds = obj_malloc ((count + 1) * sizeof (int));
+	for (int i = 0; i < count; i++) {
+		BasisGroup *g = groups[i];
+		group_inds[i] = 0;
+		for (int j = 0; j < [g count]; j++) {
+			BasisBlade *b = [g bladeAt:j];
+			layout.mask_map[[b mask] - layout.range[0]] = group_inds[count];
+			layout.group_map[group_inds[count]][0] = i;
+			layout.group_map[group_inds[count]][1] = group_inds[i]++;
+			layout.group_map[group_inds[count]][2] = group_base[i];
+			group_inds[count]++;
+		}
+	}
+	obj_free (group_inds);
+	return layout;
+}
+
+-(void)dealloc
+{
+	[set release];
+	for (int i = 0; i < count; i++) {
+		[groups[i] release];
+	}
+	obj_free (groups);
+	obj_free (group_map);
+	[super dealloc];
+}
+
+-(int)count
+{
+	return count;
+}
+
+-(int)num_components
+{
+	int         num_components = 0;
+	// assumes there is no overlap
+	for (int i = 0; i < count; i++) {
+		num_components += [groups[i] count];
+	}
+	return num_components;
+}
+
+-(int)blade_count
+{
+	return blade_count;
+}
+
+-(BasisGroup *) group:(int) ind
+{
+	return groups[ind];
+}
+
+-(BasisBlade *) bladeAt:(int) ind
+{
+	ivec3       gm = group_map[ind];
+	BasisGroup *g = groups[gm[0]];
+	return [g bladeAt:ind - gm[2]];
+}
+
+-(int) bladeIndex:(unsigned) mask
+{
+	if (![set is_member:mask]) {
+		return 0;
+	}
+	return mask_map[mask - range[0]];
+}
+
+-(Set *) set
+{
+	return set;
+}
+@end

ruamoko/gatest/main.r

@@ -0,0 +1,98 @@
+#include <AutoreleasePool.h>
+#include "algebra.h"
+#include "basisblade.h"
+#include "basisgroup.h"
+#include "metric.h"
+#include "multivector.h"
+#include "util.h"
+
+@static AutoreleasePool *autorelease_pool;
+@static void
+arp_start (void)
+{
+	autorelease_pool = [[AutoreleasePool alloc] init];
+}
+
+@static void
+arp_end (void)
+{
+	[autorelease_pool release];
+	autorelease_pool = nil;
+}
+
+int
+main ()
+{
+	arp_start ();
+
+	BasisBlade *a = [[BasisBlade basis:1] retain];
+	BasisBlade *b = [[BasisBlade basis:2] retain];
+	BasisBlade *c = [[BasisBlade basis:4] retain];
+	BasisBlade *d = [[BasisBlade basis:8] retain];
+	BasisBlade *blades[] = {a, b, c, d};
+	static string names[] = {"a", "b", "c", "d"};
+
+//	printf ("a: %@\n", a);
+//	printf ("b: %@\n", b);
+//	printf ("c: %@\n", c);
+//	printf ("d: %@\n", d);
+
+	arp_end ();
+#if 0
+	arp_start ();
+	for (int i = 0; i < 4; i++) {
+		arp_end ();
+		arp_start ();
+		BasisBlade *vec = blades[i];
+		printf ("%s: %@\n", names[i], vec);
+		for (int j = 0; j < 4; j++) {
+			BasisBlade *bvec = [vec outerProduct:blades[j]];
+			if (![bvec scale]) {
+				continue;
+			}
+			printf ("%s^%s: %@\n", names[i], names[j], bvec);
+			for (int k = 0; k < 4; k++) {
+				BasisBlade *tvec = [bvec outerProduct:blades[k]];
+				if (![tvec scale]) {
+					continue;
+				}
+				printf ("%s^%s^%s: %@\n", names[i], names[j], names[k],
+						tvec);
+				for (int l = 0; l < 4; l++) {
+					BasisBlade *qvec = [tvec outerProduct:blades[l]];
+					if (![qvec scale]) {
+						continue;
+					}
+					printf ("%s^%s^%s^%s: %@\n",
+							names[i], names[j], names[k], names[l],
+							qvec);
+				}
+			}
+		}
+	}
+	arp_end ();
+#endif
+	arp_start ();
+
+	Metric *m = [Metric R:3,0,1];
+	BasisBlade *ad = [a geometricProduct:d metric:m];
+	BasisBlade *prod = [ad geometricProduct:ad metric:m];
+	printf ("%s%s %s%s: %@\n",
+			names[0], names[3], names[0], names[3], prod);
+
+	Algebra *alg = [Algebra R:3, 0, 1];
+	double vals1[32] = {1, 0, 0, 8};
+	static double vals2[32] = {0, 1, 0, 8};//FIXME qfcc bug (static)
+	static double origin_vals[32] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1};
+	MultiVector *plane1 = [MultiVector new:alg values:vals1];
+	MultiVector *plane2 = [MultiVector new:alg values:vals2];
+	MultiVector *origin = [MultiVector new:alg values:origin_vals];
+
+	MultiVector *line = [plane1 wedge:plane2];
+	MultiVector *point = [[line dot:origin] product:line];
+	printf ("plane1:%@\nplane2:%@\nline:%@\norigin:%@\n", plane1, plane2, line, origin);
+	printf ("point:%@\n", point);
+
+	arp_end ();
+	return 0;
+}

ruamoko/gatest/metric.h

@@ -0,0 +1,16 @@
+#ifndef __metric_h
+#define __metric_h
+
+#include <Object.h>
+
+@interface Metric : Object
+{
+	unsigned    plus;	// mask of elements that square to +1
+	unsigned    minus;	// mask of elements that square to -1
+	unsigned    zero;	// mask of elements that square to 0
+}
++(Metric *)R:(int)p, int m, int z;
+-(double)apply:(unsigned) a, unsigned b;
+@end
+
+#endif//__metric_h

ruamoko/gatest/metric.r

@@ -0,0 +1,37 @@
+#include "metric.h"
+
+@implementation Metric
++(Metric *)R:(int)p, int m, int z
+{
+	Metric     *metric = [[[Metric alloc] init] autorelease];
+	metric.plus = (1 << p) - 1;
+	metric.minus = ((1 << m) - 1) << p;
+	metric.zero = ((1 << z) - 1) << (p + m);
+	return metric;
+}
+
+static double
+count_minus (unsigned minus)
+{
+	double      s = 1;
+	while (minus) {
+		if (minus & 1) {
+			s = -s;
+		}
+		minus >>= 1;
+	}
+	return s;
+}
+
+-(double)apply:(unsigned) a, unsigned b
+{
+	// find all the squared elements
+	unsigned    c = a & b;
+	// any elements that square to 0 result in 0
+	if (c & zero) {
+		return 0;
+	}
+	return count_minus (c & minus);
+}
+
+@end

ruamoko/gatest/multivector.h

@@ -0,0 +1,28 @@
+#ifndef __multivector_h
+#define __multivector_h
+#include <Object.h>
+
+@class Algebra;
+@class BasisLayout;
+
+@interface MultiVector : Object
+{
+	double     *components;
+	Algebra    *algebra;
+	BasisLayout *layout;
+	int         num_components;
+}
++(MultiVector *) new:(Algebra *) algebra;
+// NOTE: values must have the same layout as algebra
++(MultiVector *) new:(Algebra *) algebra values:(double *) values;
++(MultiVector *) new:(Algebra *) algebra group:(BasisGroup *) group;
+// NOTE: values must have the same layout as group
++(MultiVector *) new:(Algebra *) algebra group:(BasisGroup *) group values:(double *) values;
++(MultiVector *) copy:(MultiVector *) src;
+-(MultiVector *) product:(MultiVector *) rhs;
+-(MultiVector *) wedge:(MultiVector *) rhs;
+-(MultiVector *) dot:(MultiVector *) rhs;
+-(MultiVector *) dual;
+@end
+
+#endif//__multivector_h

ruamoko/gatest/multivector.r

@@ -0,0 +1,201 @@
+#include <string.h>
+#include "algebra.h"
+#include "basisblade.h"
+#include "basislayout.h"
+#include "multivector.h"
+#include "util.h"
+
+@implementation MultiVector
+static MultiVector *new_mv (Algebra *algebra, BasisLayout *layout)
+{
+	MultiVector *mv = [[[MultiVector alloc] init] autorelease];
+	mv.algebra = [algebra retain];
+	layout = layout ? layout : [algebra layout];
+	mv.layout = [layout retain];
+	mv.num_components = [layout num_components];
+	mv.components = obj_malloc (mv.num_components * sizeof (double));
+	return mv;
+}
+
++(MultiVector *) new:(Algebra *) algebra
+{
+	MultiVector *mv = new_mv (algebra, nil);
+	for (int i = 0; i < mv.num_components; i++) {
+		mv.components[i] = 0;
+	}
+	return mv;
+}
+
++(MultiVector *) new:(Algebra *) algebra values:(double *) values
+{
+	MultiVector *mv = new_mv (algebra, nil);
+	for (int i = 0; i < mv.num_components; i++) {
+		mv.components[i] = values[i];
+	}
+	return mv;
+}
+
++(MultiVector *) new:(Algebra *) algebra group:(BasisGroup *) group
+{
+	BasisLayout *layout = [BasisLayout new:1 groups:&group];
+	MultiVector *mv = new_mv (algebra, layout);
+	for (int i = 0; i < mv.num_components; i++) {
+		mv.components[i] = 0;
+	}
+	return mv;
+}
+
++(MultiVector *) new:(Algebra *) algebra group:(BasisGroup *) group values:(double *) values
+{
+	BasisLayout *layout = [BasisLayout new:1 groups:&group];
+	MultiVector *mv = new_mv (algebra, layout);
+	for (int i = 0; i < mv.num_components; i++) {
+		mv.components[i] = values[i];
+	}
+	return mv;
+}
+
++(MultiVector *) copy:(MultiVector *) src
+{
+	MultiVector *mv = new_mv (src.algebra, src.layout);
+	for (int i = 0; i < mv.num_components; i++) {
+		mv.components[i] = src.components[i];
+	}
+	return mv;
+}
+
+-(void)dealloc
+{
+	[algebra release];
+	[layout release];
+	obj_free (components);
+	[super dealloc];
+}
+
+-(string)describe
+{
+	string str = nil;
+
+	for (int i = 0; i < num_components; i++) {
+		if (!components[i]) {
+			continue;
+		}
+		BasisBlade *b = [layout bladeAt:i];
+		str = sprintf ("%s%s%g%s", str, str ? " + " : "", components[i], [b name]);
+	}
+	if (!str) {
+		str = "0";
+	}
+	return str;
+}
+
+-(MultiVector *) product:(MultiVector *) rhs
+{
+	MultiVector *prod = new_mv (algebra, nil);
+	for (int i = 0; i < num_components; i++) {
+		if (!components[i]) {
+			continue;
+		}
+		double lc = components[i];
+		BasisBlade *lb = [layout bladeAt:i];
+		for (int j = 0; j < rhs.num_components; j++) {
+			if (!rhs.components[j]) {
+				continue;
+			}
+			double rc = rhs.components[j];
+			BasisBlade *rb = [rhs.layout bladeAt:j];
+			BasisBlade *b = [lb geometricProduct:rb metric:[algebra metric]];
+			double s = [b scale];
+			if (!s) {
+				continue;
+			}
+			int ind = [rhs.layout bladeIndex:[b mask]];
+			prod.components[ind] += s * lc * rc;
+		}
+	}
+	return prod;
+}
+
+-(MultiVector *) wedge:(MultiVector *) rhs
+{
+	MultiVector *prod = new_mv (algebra, nil);
+	for (int i = 0; i < num_components; i++) {
+		if (!components[i]) {
+			continue;
+		}
+		double lc = components[i];
+		BasisBlade *lb = [layout bladeAt:i];
+		for (int j = 0; j < rhs.num_components; j++) {
+			if (!rhs.components[j]) {
+				continue;
+			}
+			double rc = rhs.components[j];
+			BasisBlade *rb = [rhs.layout bladeAt:j];
+			BasisBlade *b = [lb outerProduct:rb];
+			double s = [b scale];
+			if (!s) {
+				continue;
+			}
+			int ind = [rhs.layout bladeIndex:[b mask]];
+			prod.components[ind] += s * lc * rc;
+		}
+	}
+	return prod;
+}
+
+-(MultiVector *) dot:(MultiVector *) rhs
+{
+	MultiVector *prod = new_mv (algebra, nil);
+	for (int i = 0; i < num_components; i++) {
+		if (!components[i]) {
+			continue;
+		}
+		double lc = components[i];
+		BasisBlade *lb = [layout bladeAt:i];
+		int    lg = [lb grade];
+		for (int j = 0; j < rhs.num_components; j++) {
+			if (!rhs.components[j]) {
+				continue;
+			}
+			double rc = rhs.components[j];
+			BasisBlade *rb = [rhs.layout bladeAt:j];
+			int         rg = [rb grade];
+			BasisBlade *b = [lb geometricProduct:rb];
+			int         g = [b grade];
+			if ((lg <= rg) && (g != rg - lg)) {
+				continue;
+			}
+			if ((lg > rg) && (g != lg - rg)) {
+				continue;
+			}
+			double s = [b scale];
+			if (!s) {
+				continue;
+			}
+			int ind = [rhs.layout bladeIndex:[b mask]];
+			prod.components[ind] += s * lc * rc;
+		}
+	}
+	return prod;
+}
+
+-(MultiVector *) dual
+{
+	MultiVector *dual = new_mv (algebra, nil);
+	unsigned dual_mask = (1 << [algebra dimension]) - 1;
+	printf ("dual: %x %d\n", dual_mask, [algebra dimension]);
+	for (int i = 0; i < num_components; i++) {
+		if (!components[i]) {
+			continue;
+		}
+		double lc = components[i];
+		BasisBlade *lb = [layout bladeAt:i];
+		unsigned mask = [lb mask] ^ dual_mask;
+		double s = 1;
+		int ind = [layout bladeIndex:mask];
+		printf ("    : %x %d\n", mask, ind);
+		dual.components[ind] += s * lc;
+	}
+	return dual;
+}
+@end

ruamoko/gatest/util.h

@@ -0,0 +1,11 @@
+#ifndef __util_h
+#define __util_h
+
+void printf(string fmt, ...);
+void traceon(void);
+void traceoff(void);
+int *binomial (int n);
+int count_bits (unsigned v);
+int count_flips (unsigned a, unsigned b);
+
+#endif//__util_h

ruamoko/gatest/util.r

@@ -0,0 +1,40 @@
+#include <runtime.h>
+#include "util.h"
+
+void printf(string fmt, ...) = #0;
+void traceon(void) = #0;
+void traceoff(void) = #0;
+
+int *
+binomial (int n)
+{
+	int        *coef = obj_malloc ((n + 1) * sizeof (int));
+	int         c = 1;
+	for (int i = 0; i < n + 1; i++) {
+		coef[i] = c;
+		c = (c * (n - i)) / (i + 1);
+	}
+	return coef;
+}
+
+int
+count_bits (unsigned v)
+{
+	int         c = 0;
+	for (; v; c++) {
+		v &= v - 1u; //XXX bug in qfcc: just 1 results in extra temp
+	}
+	return c;
+}
+
+int
+count_flips (unsigned a, unsigned b)
+{
+	int         c = 0;
+	a >>= 1;
+	while (a) {
+		c += count_bits (a & b);
+		a >>= 1;
+	}
+	return c;
+}