Barry Schwartz avatar Barry Schwartz committed 2290745

A new arrangement -- include the framework for 5-dimensional conformal space, but priority is on implementing plane euclidean space.

Comments (0)

Files changed (10)

 CPPFLAGS += -I"${top_srcdir}/c" -I"${top_builddir}/c"
 
 lib_LTLIBRARIES = libebeno.la
-libebeno_la_SOURCES = c/nan.c ${srcdir}/c/e2ga.c					\
+libebeno_la_SOURCES = c/nan.c ${srcdir}/c/ga.c						\
 	fortran/pascals_triangle.f90 fortran/sorting.f90				\
 	fortran/basic_algebra.f90 fortran/linear_algebra.f90			\
 	fortran/geometric_algebra.f90 fortran/basic_geometry.f90		\
 pkgconfigdir="${libdir}/pkgconfig"
 nodist_pkgconfig_DATA = ebeno.pc
 
-${srcdir}/c/e2ga.h: ${srcdir}/c/e2ga.xml
+${srcdir}/c/ga.h: ${srcdir}/c/ga.xml
 	$(G25) $<
-#	-$(MKDIR_P) ${builddir}/c
-	mv e2ga.c e2ga.h ${srcdir}/c
-${srcdir}/c/e2ga.c: ${srcdir}/c/e2ga.h ${srcdir}/c/e2ga.xml
+	mv ga.c ga.h ${srcdir}/c
+${srcdir}/c/ga.c: ${srcdir}/c/ga.h ${srcdir}/c/ga.xml
 	@$(force_rebuild)
 
-EXTRA_DIST = ${srcdir}/c/e2ga.xml ${srcdir}/c/e2ga.c ${srcdir}/c/e2ga.h
+EXTRA_DIST = ${srcdir}/c/ga.xml ${srcdir}/c/ga.c ${srcdir}/c/ga.h
 
 #--------------------------------------------------------------------------
 
 EXTRA_DIST += COPYING INSTALL
 MOSTLYCLEANFILES = $(foreach f, $(MODFILES), $(call modfile,${f}))	\
 	Doxyfile
-MAINTAINERCLEANFILES = ${srcdir}/c/e2ga.c ${srcdir}/c/e2ga.h			\
+MAINTAINERCLEANFILES = ${srcdir}/c/ga.c ${srcdir}/c/ga.h				\
 	pure/ebeno/basic_algebra_h.pure pure/ebeno/basic_geometry_h.pure	\
 	pure/ebeno/convex_hull_h.pure pure/ebeno/curve_h.pure				\
 	pure/ebeno/geometric_algebra_h.pure

c/e2ga.c

-/*
-Copyright (c) 2012 Barry Schwartz
-*/
-/*
-
-    Copyright (c) 2012 Barry Schwartz
-
-    Permission is hereby granted, free of charge, to any person obtaining a copy
-    of this software and associated documentation files (the "Software"), to deal
-    in the Software without restriction, including without limitation the rights
-    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-    copies of the Software, and to permit persons to whom the Software is
-    furnished to do so, subject to the following conditions:
-
-    The above copyright notice and this permission notice shall be included in
-    all copies or substantial portions of the Software.
-
-    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-    THE SOFTWARE.
-  
-*/
-#include <stdio.h>
-#include "e2ga.h"
-
-const int e2ga_spaceDim = 2;
-const int e2ga_nbGroups = 3;
-const int e2ga_metricEuclidean = 1;
-const char *e2ga_basisVectorNames[2] = {
-	"e1", "e2"
-};
-const int e2ga_grades[] = {GRADE_0, GRADE_1, GRADE_2, 0, 0, 0};
-const int e2ga_groups[] = {GROUP_0, GROUP_1, GROUP_2};
-const int e2ga_groupSize[3] = {
-	1, 2, 1
-};
-const int e2ga_mvSize[8] = {
-	0, 1, 2, 3, 1, 2, 3, 4};
-const int e2ga_basisElements[4][3] = {
-	{-1},
-	{0, -1},
-	{1, -1},
-	{0, 1, -1}
-};
-const double e2ga_basisElementSignByIndex[4] =
-	{1, 1, 1, 1};
-const double e2ga_basisElementSignByBitmap[4] =
-	{1, 1, 1, 1};
-const int e2ga_basisElementIndexByBitmap[4] =
-	{0, 1, 2, 3};
-const int e2ga_basisElementBitmapByIndex[4] =
-	{0, 1, 2, 3};
-const int e2ga_basisElementGradeByBitmap[4] =
-	{0, 1, 1, 2};
-const int e2ga_basisElementGroupByBitmap[4] =
-	{0, 1, 1, 2};
-
-const char *g_e2gaTypenames[] = 
-{
-	"There are no specialized types defined"
-};
-
-
-void e2ga_double_zero_1(double *dst) {
-	dst[0]=0.0;
-}
-void e2ga_double_copy_1(double *dst, const double *src) {
-	dst[0] = src[0];
-}
-void e2ga_double_zero_2(double *dst) {
-	dst[0]=dst[1]=0.0;
-}
-void e2ga_double_copy_2(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-}
-void e2ga_double_zero_3(double *dst) {
-	dst[0]=dst[1]=dst[2]=0.0;
-}
-void e2ga_double_copy_3(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-}
-void e2ga_double_zero_4(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=0.0;
-}
-void e2ga_double_copy_4(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-}
-void e2ga_double_zero_5(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=0.0;
-}
-void e2ga_double_copy_5(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-}
-void e2ga_double_zero_6(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=0.0;
-}
-void e2ga_double_copy_6(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-}
-void e2ga_double_zero_7(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=0.0;
-}
-void e2ga_double_copy_7(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-}
-void e2ga_double_zero_8(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=0.0;
-}
-void e2ga_double_copy_8(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-}
-void e2ga_double_zero_9(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=0.0;
-}
-void e2ga_double_copy_9(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-}
-void e2ga_double_zero_10(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=0.0;
-}
-void e2ga_double_copy_10(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-}
-void e2ga_double_zero_11(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=0.0;
-}
-void e2ga_double_copy_11(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-}
-void e2ga_double_zero_12(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=0.0;
-}
-void e2ga_double_copy_12(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-	dst[11] = src[11];
-}
-void e2ga_double_zero_13(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=0.0;
-}
-void e2ga_double_copy_13(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-	dst[11] = src[11];
-	dst[12] = src[12];
-}
-void e2ga_double_zero_14(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=0.0;
-}
-void e2ga_double_copy_14(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-	dst[11] = src[11];
-	dst[12] = src[12];
-	dst[13] = src[13];
-}
-void e2ga_double_zero_15(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=dst[14]=0.0;
-}
-void e2ga_double_copy_15(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-	dst[11] = src[11];
-	dst[12] = src[12];
-	dst[13] = src[13];
-	dst[14] = src[14];
-}
-void e2ga_double_zero_16(double *dst) {
-	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=dst[14]=dst[15]=0.0;
-}
-void e2ga_double_copy_16(double *dst, const double *src) {
-	dst[0] = src[0];
-	dst[1] = src[1];
-	dst[2] = src[2];
-	dst[3] = src[3];
-	dst[4] = src[4];
-	dst[5] = src[5];
-	dst[6] = src[6];
-	dst[7] = src[7];
-	dst[8] = src[8];
-	dst[9] = src[9];
-	dst[10] = src[10];
-	dst[11] = src[11];
-	dst[12] = src[12];
-	dst[13] = src[13];
-	dst[14] = src[14];
-	dst[15] = src[15];
-}
-/** Sets N doubles to zero */
-void e2ga_double_zero_N(double *dst, int N) {
-	int i = 0;
-	while ((N-i) > 16) {
-		e2ga_double_zero_16(dst + i);
-		i += 16;
-	}
-	for (; i < N; i++)
-		dst[i] = 0.0;
-}
-/** Copies N doubles from 'src' to 'dst' */
-void e2ga_double_copy_N(double *dst, const double *src, int N) {
-	int i = 0;
-	while ((N-i) > 16) {
-		e2ga_double_copy_16(dst + i, src + i);
-		i += 16;
-	}
-	for (; i < N; i++)
-		dst[i] = src[i];
-}
-
-/** This function is not for external use. It compresses arrays of floats for storage in multivectors. */
-void compress(const double *c, double *cc, int *cgu, double epsilon, int gu) {
-	int i, j, ia = 0, ib = 0, f, s;
-	*cgu = 0;
-
-	/* for all grade parts... */
-	for (i = 0; i <= 3; i++) {
-		/* check if grade part has memory use: */
-		if (!(gu & (1 << i))) continue;
-
-		/* check if abs coordinates of grade part are all < epsilon */
-		s = e2ga_groupSize[i];
-		j = ia + s;
-		f = 0;
-		for (; ia < j; ia++)
-			if (fabs(c[ia]) > epsilon) {f = 1; break;}
-		ia = j;
-		if (f) {
-						e2ga_double_copy_N(cc + ib, c + ia - s, s);
-			ib += s;
-			*cgu |= (1 << i);
-		}
-	}
-}
-
-/** This function is not for external use. It decompresses the coordinates stored in this */
-void expand(const double *ptrs[3], const mv *src) {
-	const double *c = src->c;
-	
-	if (src->gu & 1) {
-		ptrs[0] =  c;
-		c += 1;
-	}
-	else ptrs[0] = NULL;	
-	if (src->gu & 2) {
-		ptrs[1] =  c;
-		c += 2;
-	}
-	else ptrs[1] = NULL;	
-	if (src->gu & 4) {
-		ptrs[2] =  c;
-	}
-	else ptrs[2] = NULL;	
-}
-
-
-void swapPointer(void **P1, void **P2)
-{
-	void *tmp = *P1;
-	*P1 = *P2;
-	*P2 = tmp;
-}
-
-/* 
-These strings determine how the output of string() is formatted.
-You can alter them at runtime using e2ga_setStringFormat().
-*/
- 
-const char *e2ga_string_fp = "%2.2f"; 
-const char *e2ga_string_start = ""; 
-const char *e2ga_string_end = ""; 
-const char *e2ga_string_mul = "*"; 
-const char *e2ga_string_wedge = "^"; 
-const char *e2ga_string_plus = " + "; 
-const char *e2ga_string_minus = " - "; 
-
-void e2ga_setStringFormat(const char *what, const char *format) {
- 
-	if (!strcmp(what, "fp")) 
-		e2ga_string_fp = (format) ? format : "%2.2f"; 
-	else if (!strcmp(what, "start")) 
-		e2ga_string_start = (format) ? format : ""; 
-	else if (!strcmp(what, "end")) 
-		e2ga_string_end = (format) ? format : ""; 
-	else if (!strcmp(what, "mul")) 
-		e2ga_string_mul = (format) ? format : "*"; 
-	else if (!strcmp(what, "wedge")) 
-		e2ga_string_wedge = (format) ? format : "^"; 
-	else if (!strcmp(what, "plus")) 
-		e2ga_string_plus = (format) ? format : " + "; 
-	else if (!strcmp(what, "minus")) 
-		e2ga_string_minus = (format) ? format : " - ";
-}
-
-
-#ifdef WIN32
-#define snprintf _snprintf
-#pragma warning(disable:4996) /* quit your whining already */
-#endif /* WIN32 */
-const char *toString_mv(const mv *V, char *str, int maxLength, const char *fp) 
-{
-	int dummyArg = 0; /* prevents compiler warning on some platforms */
-	int stdIdx = 0, l;
-	char tmpBuf[256], tmpFloatBuf[256]; 
-	int i, j, k = 0, bei, ia = 0, s = e2ga_mvSize[V->gu], p = 0, cnt = 0;
-
-	/* set up the floating point precision */
-	if (fp == NULL) fp = e2ga_string_fp;
-
-	/* start the string */
-	l = snprintf(tmpBuf, 256, "%s", e2ga_string_start);
-	if (stdIdx + l <= maxLength) {
-		strcpy(str + stdIdx, tmpBuf);
-		stdIdx += l;
-	}
-	else {
-		snprintf(str, maxLength, "toString_mv: buffer too small");
-		return str;
-	}
-
-	/* print all coordinates */
-	for (i = 0; i <= 3; i++) {
-		if (V->gu & (1 << i)) {
-			for (j = 0; j < e2ga_groupSize[i]; j++) {
-				double coord = (double)e2ga_basisElementSignByIndex[ia] * V->c[k];
-				/* goal: print [+|-]V->c[k][* basisVector1 ^ ... ^ basisVectorN] */			
-				snprintf(tmpFloatBuf, 256, fp, fabs(coord));
-				if (atof(tmpFloatBuf) != 0.0) {
-					l = 0;
-
-					/* print [+|-] */
-					l += snprintf(tmpBuf + l, 256-l, "%s", (coord >= 0.0) 
-						? (cnt ? e2ga_string_plus : "")
-						: e2ga_string_minus);
-						
-					/* print obj.m_c[k] */
-					l += snprintf(tmpBuf + l, 256-l, tmpFloatBuf, dummyArg);
-
-					if (i) { /* if not grade 0, print [* basisVector1 ^ ... ^ basisVectorN] */
-						l += snprintf(tmpBuf + l, 256-l, "%s", e2ga_string_mul);
-
-						/* print all basis vectors */
-						bei = 0;
-						while (e2ga_basisElements[ia][bei] >= 0) {
-							l += snprintf(tmpBuf + l, 256-l, "%s%s", (bei) ? e2ga_string_wedge : "", 
-							 e2ga_basisVectorNames[e2ga_basisElements[ia][bei]]);
-							 bei++;
-						}
-					}
-
-					/* copy all to 'str' */
-					if (stdIdx + l <= maxLength) {
-						strcpy(str + stdIdx, tmpBuf);
-						stdIdx += l;
-					}
-					else {
-						snprintf(str, maxLength, "toString_mv: buffer too small");
-						return str;
-					}
-					cnt++;
-				}
-				k++; ia++;
-			}
-		}
-		else ia += e2ga_groupSize[i];
-	}
-
-    /* if no coordinates printed: 0 */
-	l = 0;
-	if (cnt == 0) {
-		l += snprintf(tmpBuf + l, 256-l, "0");
-	}
-
-    /* end the string */
-	l += snprintf(tmpBuf + l, 256-l, "%s", e2ga_string_end);
-	if (stdIdx + l <= maxLength) {
-		strcpy(str + stdIdx, tmpBuf);
-		stdIdx += l;
-		return str;
-	}
-	else {
-		snprintf(str, maxLength, "toString_mv: buffer too small\n");
-		return str;
-	}
-}
-
-/**
- * Computes the partial geometric product of two multivectors (group 0  x  group 0 -> group 0)
- */
-void gp_default_0_0_0(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 0  x  group 1 -> group 1)
- */
-void gp_default_0_1_1(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 0  x  group 2 -> group 2)
- */
-void gp_default_0_2_2(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 1  x  group 0 -> group 1)
- */
-void gp_default_1_0_1(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 1  x  group 1 -> group 0)
- */
-void gp_default_1_1_0(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 1  x  group 1 -> group 2)
- */
-void gp_default_1_1_2(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 1  x  group 2 -> group 1)
- */
-void gp_default_1_2_1(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 2  x  group 0 -> group 2)
- */
-void gp_default_2_0_2(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 2  x  group 1 -> group 1)
- */
-void gp_default_2_1_1(const double *A, const double *B, double *C);
-/**
- * Computes the partial geometric product of two multivectors (group 2  x  group 2 -> group 0)
- */
-void gp_default_2_2_0(const double *A, const double *B, double *C);
-/**
- * copies coordinates of group 0
- */
-void copyGroup_0(const double *A, double *C);
-/**
- * copies and multiplies (by s) coordinates of group 0
- */
-void copyMul_0(const double *A, double *C, double s);
-/**
- * copies and divides (by s) coordinates of group 0
- */
-void copyDiv_0(const double *A, double *C, double s);
-/**
- * adds coordinates of group 0 from variable A to C
- */
-void add_0(const double *A, double *C);
-/**
- * subtracts coordinates of group 0 in variable A from C
- */
-void sub_0(const double *A, double *C);
-/**
- * negate coordinates of group 0 of variable A
- */
-void neg_0(const double *A, double *C);
-/**
- * adds coordinates of group 0 of variables A and B
- */
-void add2_0_0(const double *A, const double *B, double *C);
-/**
- * subtracts coordinates of group 0 of variables A from B
- */
-void sub2_0_0(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise multiplication of coordinates of group 0 of variables A and B
- */
-void hp_0_0(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise division of coordinates of group 0 of variables A and B
- * (no checks for divide by zero are made)
- */
-void ihp_0_0(const double *A, const double *B, double *C);
-/**
- * check for equality up to eps of coordinates of group 0 of variables A and B
- */
-int equals_0_0(const double *A, const double *B, double eps);
-/**
- * checks if coordinates of group 0 of variable A are zero up to eps
- */
-int zeroGroup_0(const double *A, double eps);
-/**
- * copies coordinates of group 1
- */
-void copyGroup_1(const double *A, double *C);
-/**
- * copies and multiplies (by s) coordinates of group 1
- */
-void copyMul_1(const double *A, double *C, double s);
-/**
- * copies and divides (by s) coordinates of group 1
- */
-void copyDiv_1(const double *A, double *C, double s);
-/**
- * adds coordinates of group 1 from variable A to C
- */
-void add_1(const double *A, double *C);
-/**
- * subtracts coordinates of group 1 in variable A from C
- */
-void sub_1(const double *A, double *C);
-/**
- * negate coordinates of group 1 of variable A
- */
-void neg_1(const double *A, double *C);
-/**
- * adds coordinates of group 1 of variables A and B
- */
-void add2_1_1(const double *A, const double *B, double *C);
-/**
- * subtracts coordinates of group 1 of variables A from B
- */
-void sub2_1_1(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise multiplication of coordinates of group 1 of variables A and B
- */
-void hp_1_1(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise division of coordinates of group 1 of variables A and B
- * (no checks for divide by zero are made)
- */
-void ihp_1_1(const double *A, const double *B, double *C);
-/**
- * check for equality up to eps of coordinates of group 1 of variables A and B
- */
-int equals_1_1(const double *A, const double *B, double eps);
-/**
- * checks if coordinates of group 1 of variable A are zero up to eps
- */
-int zeroGroup_1(const double *A, double eps);
-/**
- * copies coordinates of group 2
- */
-void copyGroup_2(const double *A, double *C);
-/**
- * copies and multiplies (by s) coordinates of group 2
- */
-void copyMul_2(const double *A, double *C, double s);
-/**
- * copies and divides (by s) coordinates of group 2
- */
-void copyDiv_2(const double *A, double *C, double s);
-/**
- * adds coordinates of group 2 from variable A to C
- */
-void add_2(const double *A, double *C);
-/**
- * subtracts coordinates of group 2 in variable A from C
- */
-void sub_2(const double *A, double *C);
-/**
- * negate coordinates of group 2 of variable A
- */
-void neg_2(const double *A, double *C);
-/**
- * adds coordinates of group 2 of variables A and B
- */
-void add2_2_2(const double *A, const double *B, double *C);
-/**
- * subtracts coordinates of group 2 of variables A from B
- */
-void sub2_2_2(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise multiplication of coordinates of group 2 of variables A and B
- */
-void hp_2_2(const double *A, const double *B, double *C);
-/**
- * performs coordinate-wise division of coordinates of group 2 of variables A and B
- * (no checks for divide by zero are made)
- */
-void ihp_2_2(const double *A, const double *B, double *C);
-/**
- * check for equality up to eps of coordinates of group 2 of variables A and B
- */
-int equals_2_2(const double *A, const double *B, double eps);
-/**
- * checks if coordinates of group 2 of variable A are zero up to eps
- */
-int zeroGroup_2(const double *A, double eps);
-/**
- * Computes the partial dual (w.r.t. full space) of a multivector.
- */
-void dual_default_0_2(const double *A, double *C);
-/**
- * Computes the partial undual (w.r.t. full space) of a multivector.
- */
-void undual_default_0_2(const double *A, double *C);
-/**
- * Computes the partial dual (w.r.t. full space) of a multivector.
- */
-void dual_default_1_1(const double *A, double *C);
-/**
- * Computes the partial undual (w.r.t. full space) of a multivector.
- */
-void undual_default_1_1(const double *A, double *C);
-/**
- * Computes the partial dual (w.r.t. full space) of a multivector.
- */
-void dual_default_2_0(const double *A, double *C);
-/**
- * Computes the partial undual (w.r.t. full space) of a multivector.
- */
-void undual_default_2_0(const double *A, double *C);
-void gp_default_0_0_0(const double *A, const double *B, double *C) {
-	C[0] += A[0]*B[0];
-}
-void gp_default_0_1_1(const double *A, const double *B, double *C) {
-	C[0] += A[0]*B[0];
-	C[1] += A[0]*B[1];
-}
-void gp_default_0_2_2(const double *A, const double *B, double *C) {
-	gp_default_0_0_0(A, B, C);
-}
-void gp_default_1_0_1(const double *A, const double *B, double *C) {
-	C[0] += A[0]*B[0];
-	C[1] += A[1]*B[0];
-}
-void gp_default_1_1_0(const double *A, const double *B, double *C) {
-	C[0] += (A[0]*B[0]+A[1]*B[1]);
-}
-void gp_default_1_1_2(const double *A, const double *B, double *C) {
-	C[0] += (A[0]*B[1]-A[1]*B[0]);
-}
-void gp_default_1_2_1(const double *A, const double *B, double *C) {
-	C[0] += -A[1]*B[0];
-	C[1] += A[0]*B[0];
-}
-void gp_default_2_0_2(const double *A, const double *B, double *C) {
-	gp_default_0_0_0(A, B, C);
-}
-void gp_default_2_1_1(const double *A, const double *B, double *C) {
-	C[0] += A[0]*B[1];
-	C[1] += -A[0]*B[0];
-}
-void gp_default_2_2_0(const double *A, const double *B, double *C) {
-	C[0] += -A[0]*B[0];
-}
-void copyGroup_0(const double *A, double *C) {
-	C[0] = A[0];
-}
-void copyMul_0(const double *A, double *C, double s) {
-	C[0] = A[0]*s;
-}
-void copyDiv_0(const double *A, double *C, double s) {
-	C[0] = A[0]/s;
-}
-void add_0(const double *A, double *C) {
-	C[0] += A[0];
-}
-void sub_0(const double *A, double *C) {
-	C[0] -= A[0];
-}
-void neg_0(const double *A, double *C) {
-	C[0] = -A[0];
-}
-void add2_0_0(const double *A, const double *B, double *C) {
-	C[0] = (A[0]+B[0]);
-}
-void sub2_0_0(const double *A, const double *B, double *C) {
-	C[0] = (A[0]-B[0]);
-}
-void hp_0_0(const double *A, const double *B, double *C) {
-	C[0] = A[0]*B[0];
-}
-void ihp_0_0(const double *A, const double *B, double *C) {
-	C[0] = A[0]/((B[0]));
-}
-int equals_0_0(const double *A, const double *B, double eps) {
-		if (((A[0] - B[0]) < -eps) || ((A[0] - B[0]) > eps)) return 0;
-	return 1;
-}
-int zeroGroup_0(const double *A, double eps) {
-		if ((A[0] < -eps) || (A[0] > eps)) return 0;
-		return 1;
-}
-void copyGroup_1(const double *A, double *C) {
-	C[0] = A[0];
-	C[1] = A[1];
-}
-void copyMul_1(const double *A, double *C, double s) {
-	C[0] = A[0]*s;
-	C[1] = A[1]*s;
-}
-void copyDiv_1(const double *A, double *C, double s) {
-	C[0] = A[0]/s;
-	C[1] = A[1]/s;
-}
-void add_1(const double *A, double *C) {
-	C[0] += A[0];
-	C[1] += A[1];
-}
-void sub_1(const double *A, double *C) {
-	C[0] -= A[0];
-	C[1] -= A[1];
-}
-void neg_1(const double *A, double *C) {
-	C[0] = -A[0];
-	C[1] = -A[1];
-}
-void add2_1_1(const double *A, const double *B, double *C) {
-	C[0] = (A[0]+B[0]);
-	C[1] = (A[1]+B[1]);
-}
-void sub2_1_1(const double *A, const double *B, double *C) {
-	C[0] = (A[0]-B[0]);
-	C[1] = (A[1]-B[1]);
-}
-void hp_1_1(const double *A, const double *B, double *C) {
-	C[0] = A[0]*B[0];
-	C[1] = A[1]*B[1];
-}
-void ihp_1_1(const double *A, const double *B, double *C) {
-	C[0] = A[0]/((B[0]));
-	C[1] = A[1]/((B[1]));
-}
-int equals_1_1(const double *A, const double *B, double eps) {
-		if (((A[0] - B[0]) < -eps) || ((A[0] - B[0]) > eps)) return 0;
-		if (((A[1] - B[1]) < -eps) || ((A[1] - B[1]) > eps)) return 0;
-	return 1;
-}
-int zeroGroup_1(const double *A, double eps) {
-		if ((A[0] < -eps) || (A[0] > eps)) return 0;
-		if ((A[1] < -eps) || (A[1] > eps)) return 0;
-		return 1;
-}
-void copyGroup_2(const double *A, double *C) {
-	copyGroup_0(A, C);
-}
-void copyMul_2(const double *A, double *C, double s) {
-	copyMul_0(A, C, s);
-}
-void copyDiv_2(const double *A, double *C, double s) {
-	copyDiv_0(A, C, s);
-}
-void add_2(const double *A, double *C) {
-	add_0(A, C);
-}
-void sub_2(const double *A, double *C) {
-	sub_0(A, C);
-}
-void neg_2(const double *A, double *C) {
-	neg_0(A, C);
-}
-void add2_2_2(const double *A, const double *B, double *C) {
-	add2_0_0(A, B, C);
-}
-void sub2_2_2(const double *A, const double *B, double *C) {
-	sub2_0_0(A, B, C);
-}
-void hp_2_2(const double *A, const double *B, double *C) {
-	hp_0_0(A, B, C);
-}
-void ihp_2_2(const double *A, const double *B, double *C) {
-	ihp_0_0(A, B, C);
-}
-int equals_2_2(const double *A, const double *B, double eps) {
-	return equals_0_0(A, B, eps);
-}
-int zeroGroup_2(const double *A, double eps) {
-	return zeroGroup_0(A, eps);
-}
-void dual_default_0_2(const double *A, double *C) {
-	C[0] = -A[0];
-}
-void undual_default_0_2(const double *A, double *C) {
-	C[0] = A[0];
-}
-void dual_default_1_1(const double *A, double *C) {
-	C[0] = A[1];
-	C[1] = -A[0];
-}
-void undual_default_1_1(const double *A, double *C) {
-	C[0] = -A[1];
-	C[1] = A[0];
-}
-void dual_default_2_0(const double *A, double *C) {
-	undual_default_0_2(A, C);
-}
-void undual_default_2_0(const double *A, double *C) {
-	dual_default_0_2(A, C);
-}
-void mv_setZero(mv *M) {
-	M->gu = 0;
-}
-void mv_setScalar(mv *M, double val) {
-	M->gu = 1;
-	M->c[0] = val;
-}
-void mv_setArray(mv *M, int gu, const double *arr) {
-	M->gu = gu;
-	e2ga_double_copy_N(M->c, arr, e2ga_mvSize[gu]);
-
-}
-void mv_copy(mv *dst, const mv *src) {
-	int i;
-	dst->gu = src->gu;
-	for (i = 0; i < e2ga_mvSize[src->gu]; i++)
-		dst->c[i] = (double)src->c[i];
-}
-
-void mv_reserveGroup_0(mv *A) {
-	int groupUsageBelow, groupUsageAbove, newGroupUsage, newGroupUsageBelowNextGroup;
-	int i;
-	double *dst, *src;
-	if ((A->gu & 1) == 0) {
-
-		groupUsageBelow = A->gu & 0;
-		groupUsageAbove = A->gu ^ groupUsageBelow;
-		newGroupUsage = A->gu | 1;
-		newGroupUsageBelowNextGroup = newGroupUsage & 1;
-
-		dst = A->c + e2ga_mvSize[newGroupUsageBelowNextGroup];
-		src = A->c + e2ga_mvSize[groupUsageBelow];
-		for (i = e2ga_mvSize[groupUsageAbove]-1; i >= 0; i--) /* work from end to start of array to avoid overwriting (dst is always beyond src) */
-			dst[i] = src[i];
-		e2ga_double_zero_1(A->c);
-
-		A->gu = newGroupUsage;
-	}
-}
-
-void mv_reserveGroup_1(mv *A) {
-	int groupUsageBelow, groupUsageAbove, newGroupUsage, newGroupUsageBelowNextGroup;
-	int i;
-	double *dst, *src;
-	if ((A->gu & 2) == 0) {
-
-		groupUsageBelow = A->gu & 1;
-		groupUsageAbove = A->gu ^ groupUsageBelow;
-		newGroupUsage = A->gu | 2;
-		newGroupUsageBelowNextGroup = newGroupUsage & 3;
-
-		dst = A->c + e2ga_mvSize[newGroupUsageBelowNextGroup];
-		src = A->c + e2ga_mvSize[groupUsageBelow];
-		for (i = e2ga_mvSize[groupUsageAbove]-1; i >= 0; i--) /* work from end to start of array to avoid overwriting (dst is always beyond src) */
-			dst[i] = src[i];
-		e2ga_double_zero_2(A->c);
-
-		A->gu = newGroupUsage;
-	}
-}
-
-void mv_reserveGroup_2(mv *A) {
-	int groupUsageBelow, groupUsageAbove, newGroupUsage, newGroupUsageBelowNextGroup;
-	if ((A->gu & 4) == 0) {
-
-		groupUsageBelow = A->gu & 3;
-		groupUsageAbove = A->gu ^ groupUsageBelow;
-		newGroupUsage = A->gu | 4;
-		newGroupUsageBelowNextGroup = newGroupUsage & 7;
-
-		e2ga_double_zero_1(A->c);
-
-		A->gu = newGroupUsage;
-	}
-}
-
-double mv_scalar(const mv *A) {
-	return (A->gu & 1) ? A->c[e2ga_mvSize[A->gu & 0] + 0] : 0.0;
-}
-double mv_double(const mv *A) {
-	return mv_scalar(A);
-}
-
-void mv_set_scalar(mv *A, double scalar_coord) {
-	mv_reserveGroup_0(A);
-	A->c[e2ga_mvSize[A->gu & 0] + 0] = scalar_coord;
-}
-
-double mv_e1(const mv *A) {
-	return (A->gu & 2) ? A->c[e2ga_mvSize[A->gu & 1] + 0] : 0.0;
-}
-
-void mv_set_e1(mv *A, double e1_coord) {
-	mv_reserveGroup_1(A);
-	A->c[e2ga_mvSize[A->gu & 1] + 0] = e1_coord;
-}
-
-double mv_e2(const mv *A) {
-	return (A->gu & 2) ? A->c[e2ga_mvSize[A->gu & 1] + 1] : 0.0;
-}
-
-void mv_set_e2(mv *A, double e2_coord) {
-	mv_reserveGroup_1(A);
-	A->c[e2ga_mvSize[A->gu & 1] + 1] = e2_coord;
-}
-
-double mv_e1_e2(const mv *A) {
-	return (A->gu & 4) ? A->c[e2ga_mvSize[A->gu & 3] + 0] : 0.0;
-}
-
-void mv_set_e1_e2(mv *A, double e1_e2_coord) {
-	mv_reserveGroup_2(A);
-	A->c[e2ga_mvSize[A->gu & 3] + 0] = e1_e2_coord;
-}
-
-double mv_largestCoordinate(const mv *x) {
-	double maxValue = 0.0;
-	int nbC = e2ga_mvSize[x->gu], i;
-	for (i = 0; i < nbC; i++)
-		if (fabs(x->c[i]) > maxValue) maxValue = fabs(x->c[i]);
-	return maxValue;
-}
-
-double mv_largestBasisBlade(const mv *x, unsigned int *bm) {
-	int nc = e2ga_mvSize[x->gu];
-	double maxC = -1.0, C;
-
-	int idx = 0;
-	int group = 0;
-	int i = 0, j;
-	*bm = 0;
-
-	while (i < nc) {
-		if (x->gu & (1 << group)) {
-			for (j = 0; j < e2ga_groupSize[group]; j++) {
-				C = fabs(x->c[i]);
-				if (C > maxC) {
-					maxC = C;
-					*bm = e2ga_basisElementBitmapByIndex[idx];
-				}
-				idx++;
-				i++;
-			}
-		}
-		else idx += e2ga_groupSize[group];
-		group++;
-	}
-
-	return maxC;
-} /* end of mv::largestBasisBlade() */
-
-
-
-
-
-
-
-
-
-
-void add_mv_mv(mv *_dst, const mv *a, const mv *b)
-{
-	int aidx = 0, bidx = 0, cidx = 0;
-	_dst->gu = a->gu | b->gu;
-	
-	if (a->gu & 1) {
-		if (b->gu & 1) {
-			add2_0_0(a->c + aidx, b->c + bidx, _dst->c + cidx);
-			bidx += 1;
-		}
-		else copyGroup_0(a->c + aidx, _dst->c + cidx);
-		aidx += 1;
-		cidx += 1;
-	}
-	else if (b->gu & 1) {
-		copyGroup_0(b->c + bidx, _dst->c + cidx);
-		bidx += 1;
-		cidx += 1;
-	}
-	
-	if (a->gu & 2) {
-		if (b->gu & 2) {
-			add2_1_1(a->c + aidx, b->c + bidx, _dst->c + cidx);
-			bidx += 2;
-		}
-		else copyGroup_1(a->c + aidx, _dst->c + cidx);
-		aidx += 2;
-		cidx += 2;
-	}
-	else if (b->gu & 2) {
-		copyGroup_1(b->c + bidx, _dst->c + cidx);
-		bidx += 2;
-		cidx += 2;
-	}
-	
-	if (a->gu & 4) {
-		if (b->gu & 4) {
-			add2_2_2(a->c + aidx, b->c + bidx, _dst->c + cidx);
-		}
-		else copyGroup_2(a->c + aidx, _dst->c + cidx);
-		cidx += 1;
-	}
-	else if (b->gu & 4) {
-		copyGroup_2(b->c + bidx, _dst->c + cidx);
-		cidx += 1;
-	}
-}
-void applyUnitVersor_mv_mv(mv *_dst, const mv *a, const mv *b)
-{
-	mv rev; /* temp space for reverse */
-	mv tmp, tmp2; /* temp variables */
-	reverse_mv(&rev, a); /* compute inverse or reverse */
-	gp_mv_mv(&tmp, a, b); /* compute geometric product a b */
-	gp_mv_mv(&tmp2, &tmp, &rev); /* compute geometric product (a b) &rev */
-	
-	extractGrade_mv(_dst, &tmp2, b->gu); /* ditch grade parts which were not in b */
-}
-void applyVersor_mv_mv(mv *_dst, const mv *a, const mv *b)
-{
-	mv inv; /* temp space for reverse */
-	mv tmp, tmp2; /* temp variables */
-	versorInverse_mv(&inv, a); /* compute inverse or reverse */
-	gp_mv_mv(&tmp, a, b); /* compute geometric product a b */
-	gp_mv_mv(&tmp2, &tmp, &inv); /* compute geometric product (a b) &inv */
-	
-	extractGrade_mv(_dst, &tmp2, b->gu); /* ditch grade parts which were not in b */
-}
-void applyVersorWI_mv_mv_mv(mv *_dst, const mv *a, const mv *b, const mv *c)
-{
-	mv tmp, tmp2; /* temp variables */
-	gp_mv_mv(&tmp, a, b); /* compute geometric product a b */
-	gp_mv_mv(&tmp2, &tmp, c); /* compute geometric product (a b) c */
-	
-	extractGrade_mv(_dst, &tmp2, b->gu); /* ditch grade parts which were not in b */
-}
-void div_mv_double(mv *_dst, const mv *a, const double b)
-{
-	int idx = 0;
-	_dst->gu = a->gu;
-	
-	if (a->gu & 1) {
-		copyDiv_0(a->c + idx, _dst->c + idx, b);
-		idx += 1;
-	}
-	
-	if (a->gu & 2) {
-		copyDiv_1(a->c + idx, _dst->c + idx, b);
-		idx += 2;
-	}
-	
-	if (a->gu & 4) {
-		copyDiv_2(a->c + idx, _dst->c + idx, b);
-	}
-}
-void dual_mv(mv *_dst, const mv *a)
-{
-	int idx = 0;
-	double c[4];
-	e2ga_double_zero_4(c);
-	if (a->gu & 1) {
-		dual_default_0_2(a->c + idx, c + 3);
-		idx += 1;
-	}
-	
-	if (a->gu & 2) {
-		dual_default_1_1(a->c + idx, c + 1);
-		idx += 2;
-	}
-	
-	if (a->gu & 4) {
-		dual_default_2_0(a->c + idx, c + 0);
-	}
-	
-	compress(c, _dst->c, &(_dst->gu), 0.0, 7);
-}
-int equals_mv_mv_double(const mv *a, const mv *b, const double c)
-{
-	int aidx = 0, bidx = 0;
-	
-	if (a->gu & 1) {
-		if (b->gu & 1) {
-			if (!equals_0_0(a->c + aidx, b->c + bidx, c)) return 0;
-			bidx += 1;
-		}
-		else if (!zeroGroup_0(a->c + bidx, c)) return 0;
-		aidx += 1;
-	}
-	else if (b->gu & 1) {
-		if (!zeroGroup_0(b->c + bidx, c)) return 0;
-		bidx += 1;
-	}
-	
-	if (a->gu & 2) {
-		if (b->gu & 2) {
-			if (!equals_1_1(a->c + aidx, b->c + bidx, c)) return 0;
-			bidx += 2;
-		}
-		else if (!zeroGroup_1(a->c + bidx, c)) return 0;
-		aidx += 2;
-	}
-	else if (b->gu & 2) {
-		if (!zeroGroup_1(b->c + bidx, c)) return 0;
-		bidx += 2;
-	}
-	
-	if (a->gu & 4) {
-		if (b->gu & 4) {
-			if (!equals_2_2(a->c + aidx, b->c + bidx, c)) return 0;
-		}
-		else if (!zeroGroup_2(a->c + bidx, c)) return 0;
-	}
-	else if (b->gu & 4) {
-		if (!zeroGroup_2(b->c + bidx, c)) return 0;
-	}
-	return 1;
-}
-void extractGrade_mv(mv *_dst, const mv *a, int groupBitmap)
-{
-	int aidx = 0, cidx = 0;
-	_dst->gu = a->gu & groupBitmap;
-	
-	if (a->gu & 1) {
-		if (groupBitmap & 1) {
-			copyGroup_0(a->c + aidx, _dst->c + cidx);
-			cidx += 1;
-		}
-		aidx += 1;
-	}
-	
-	if (a->gu & 2) {
-		if (groupBitmap & 2) {
-			copyGroup_1(a->c + aidx, _dst->c + cidx);
-			cidx += 2;
-		}
-		aidx += 2;
-	}
-	
-	if (a->gu & 4) {
-		if (groupBitmap & 4) {
-			copyGroup_2(a->c + aidx, _dst->c + cidx);
-		}
-	}
-}
-void extractGrade0_mv(mv *_dst, const mv *a)
-{
-	extractGrade_mv(_dst, a, 1);
-}
-void extractGrade1_mv(mv *_dst, const mv *a)
-{
-	extractGrade_mv(_dst, a, 2);
-}
-void extractGrade2_mv(mv *_dst, const mv *a)
-{
-	extractGrade_mv(_dst, a, 4);
-}
-void gp_mv_mv(mv *_dst, const mv *a, const mv *b)
-{
-	double c[4];
-	const double* _a[3];
-	const double* _b[3];
-	expand(_a, a);
-	expand(_b, b);
-	e2ga_double_zero_4(c);
-	if (a->gu & 1) {
-		if (b->gu & 1) {
-			gp_default_0_0_0(_a[0], _b[0], c + 0);
-		}
-		if (b->gu & 2) {
-			gp_default_0_1_1(_a[0], _b[1], c + 1);
-		}
-		if (b->gu & 4) {
-			gp_default_0_2_2(_a[0], _b[2], c + 3);
-		}
-	}
-	if (a->gu & 2) {
-		if (b->gu & 1) {
-			gp_default_1_0_1(_a[1], _b[0], c + 1);
-		}
-		if (b->gu & 2) {
-			gp_default_1_1_0(_a[1], _b[1], c + 0);
-			gp_default_1_1_2(_a[1], _b[1], c + 3);
-		}
-		if (b->gu & 4) {
-			gp_default_1_2_1(_a[1], _b[2], c + 1);
-		}
-	}
-	if (a->gu & 4) {
-		if (b->gu & 1) {
-			gp_default_2_0_2(_a[2], _b[0], c + 3);
-		}
-		if (b->gu & 2) {
-			gp_default_2_1_1(_a[2], _b[1], c + 1);
-		}
-		if (b->gu & 4) {
-			gp_default_2_2_0(_a[2], _b[2], c + 0);
-		}
-	}
-	compress(c, _dst->c, &(_dst->gu), 0.0, 7);
-}
-double norm2_mv(const mv *a)
-{
-	double c[1];
-	double n2 = 0.0;
-	int idx = 0;
-	
-	if (a->gu & 1) { /* group 0 (grade 0) */
-		c[0] = 0.0;
-		gp_default_0_0_0(a->c + idx, a->c + idx, c);
-		n2 += c[0];
-		idx += 1;
-	}
-	
-	if (a->gu & 2) { /* group 1 (grade 1) */
-		c[0] = 0.0;
-		gp_default_1_1_0(a->c + idx, a->c + idx, c);
-		n2 += c[0];
-		idx += 2;
-	}
-	
-	if (a->gu & 4) { /* group 2 (grade 2) */
-		c[0] = 0.0;
-		gp_default_2_2_0(a->c + idx, a->c + idx, c);
-		n2 -= c[0];
-	}
-	return n2;
-}
-double norm2_mv_returns_scalar(const mv *a) {
-	return norm2_mv(a);
-}
-void reverse_mv(mv *_dst, const mv *a)
-{
-	int idx = 0;
-	_dst->gu = a->gu;
-	
-	if (a->gu & 1) {
-		copyGroup_0(a->c + idx, _dst->c + idx);
-		idx += 1;
-	}
-	
-	if (a->gu & 2) {
-		copyGroup_1(a->c + idx, _dst->c + idx);
-		idx += 2;
-	}
-	
-	if (a->gu & 4) {
-		neg_2(a->c + idx, _dst->c + idx);
-	}
-}
-void subtract_mv_mv(mv *_dst, const mv *a, const mv *b)
-{
-	int aidx = 0, bidx = 0, cidx = 0;
-	_dst->gu = a->gu | b->gu;
-	
-	if (a->gu & 1) {
-		if (b->gu & 1) {
-			sub2_0_0(a->c + aidx, b->c + bidx, _dst->c + cidx);
-			bidx += 1;
-		}
-		else copyGroup_0(a->c + aidx, _dst->c + cidx);
-		aidx += 1;
-		cidx += 1;
-	}
-	else if (b->gu & 1) {
-		neg_0(b->c + bidx, _dst->c + cidx);
-		bidx += 1;
-		cidx += 1;
-	}
-	
-	if (a->gu & 2) {
-		if (b->gu & 2) {
-			sub2_1_1(a->c + aidx, b->c + bidx, _dst->c + cidx);
-			bidx += 2;
-		}
-		else copyGroup_1(a->c + aidx, _dst->c + cidx);
-		aidx += 2;
-		cidx += 2;
-	}
-	else if (b->gu & 2) {
-		neg_1(b->c + bidx, _dst->c + cidx);
-		bidx += 2;
-		cidx += 2;
-	}
-	
-	if (a->gu & 4) {
-		if (b->gu & 4) {
-			sub2_2_2(a->c + aidx, b->c + bidx, _dst->c + cidx);
-		}
-		else copyGroup_2(a->c + aidx, _dst->c + cidx);
-		cidx += 1;
-	}
-	else if (b->gu & 4) {
-		neg_2(b->c + bidx, _dst->c + cidx);
-		cidx += 1;
-	}
-}
-void undual_mv(mv *_dst, const mv *a)
-{
-	int idx = 0;
-	double c[4];
-	e2ga_double_zero_4(c);
-	if (a->gu & 1) {
-		undual_default_0_2(a->c + idx, c + 3);
-		idx += 1;
-	}
-	
-	if (a->gu & 2) {
-		undual_default_1_1(a->c + idx, c + 1);
-		idx += 2;
-	}
-	
-	if (a->gu & 4) {
-		undual_default_2_0(a->c + idx, c + 0);
-	}
-	
-	compress(c, _dst->c, &(_dst->gu), 0.0, 7);
-}
-void versorInverse_mv(mv *_dst, const mv *a)
-{
-	int idx = 0;
-	double n2 = norm2_mv_returns_scalar(a);
-	_dst->gu = a->gu;
-	
-	if (a->gu & 1) {
-		copyDiv_0(a->c + idx, _dst->c + idx, n2);
-		idx += 1;
-	}
-	
-	if (a->gu & 2) {
-		copyDiv_1(a->c + idx, _dst->c + idx, n2);
-		idx += 2;
-	}
-	
-	if (a->gu & 4) {
-		copyDiv_2(a->c + idx, _dst->c + idx, -n2);
-	}
-}

c/e2ga.h

-/*
-Copyright (c) 2012 Barry Schwartz
-*/
-/*
-
-    Copyright (c) 2012 Barry Schwartz
-
-    Permission is hereby granted, free of charge, to any person obtaining a copy
-    of this software and associated documentation files (the "Software"), to deal
-    in the Software without restriction, including without limitation the rights
-    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-    copies of the Software, and to permit persons to whom the Software is
-    furnished to do so, subject to the following conditions:
-
-    The above copyright notice and this permission notice shall be included in
-    all copies or substantial portions of the Software.
-
-    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-    THE SOFTWARE.
-  
-*/
-
-/*! \mainpage e2ga documentation
- *
- * e2ga implementation generated by Gaigen 2.5. 
- * 
- * 
- * License: 
-
-    Copyright (c) 2012 Barry Schwartz
-
-    Permission is hereby granted, free of charge, to any person obtaining a copy
-    of this software and associated documentation files (the "Software"), to deal
-    in the Software without restriction, including without limitation the rights
-    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-    copies of the Software, and to permit persons to whom the Software is
-    furnished to do so, subject to the following conditions:
-
-    The above copyright notice and this permission notice shall be included in
-    all copies or substantial portions of the Software.
-
-    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-    THE SOFTWARE.
-    
- * 
- * \section intro_sec Introduction
- *
- * Todo
- * 
- */
-#ifndef _E2GA_H_
-#define _E2GA_H_
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-
-/* group: 1*/
-#define GROUP_0 1
-/* group: e1, e2*/
-#define GROUP_1 2
-/* group: e1^e2*/
-#define GROUP_2 4
-#define GRADE_0 1
-#define GRADE_1 2
-#define GRADE_2 4
-
-
-/** The dimension of the space. */
-extern const int e2ga_spaceDim;
-
-/** Number of groups/grades of coordinates in a multivector. */
-extern const int e2ga_nbGroups;
-
-/** The constants for the grades in an array. */
-extern const int e2ga_grades[];
-
-/** The constants for the groups in an array. */
-extern const int e2ga_groups[];
-
-/** Is the metric of the space Euclidean? (0 or 1). */
-extern const int e2ga_metricEuclidean;
-
-/** This array can be used to lookup the number of coordinates for a group part of a general multivector. */
-extern const int e2ga_groupSize[3];
-
-/** This array can be used to lookup the number of coordinates based on a group usage bitmap. */
-extern const int e2ga_mvSize[8];
-
-/** This array of ASCIIZ strings contains the names of the basis vectors. */
-extern const char *e2ga_basisVectorNames[2];
-
-/** This array of integers contains the order of basis elements in the general multivector
-  * Use it to answer: 'what basis vectors are in the basis element at position [x]'? */
-extern const int e2ga_basisElements[4][3];
-
-/** This array of integers contains the 'sign' (even/odd permutation of canonical order) of basis elements in the general multivector
-  * Use it to answer 'what is the permutation of the coordinate at index [x]'? */
-extern const double e2ga_basisElementSignByIndex[4];
-
-/** This array of integers contains the 'sign' (even/odd permutation of canonical order) of basis elements in the general multivector
-  * Use it to answer 'what is the permutation of the coordinate of bitmap [x]'? */
-extern const double e2ga_basisElementSignByBitmap[4];
-
-/** This array of integers contains the order of basis elements in the general multivector
-   * Use it to answer: 'at what index do I find basis element [x] (x = basis vector bitmap)?' */
-extern const int e2ga_basisElementIndexByBitmap[4];
-
-/** This array of integers contains the indices of basis elements in the general multivector
-  * Use it to answer: 'what basis element do I find at index [x]'? */
-extern const int e2ga_basisElementBitmapByIndex[4];
-
-/** This array of grade of each basis elements in the general multivector
-   * Use it to answer: 'what is the grade of basis element bitmap [x]'? */
-extern const int e2ga_basisElementGradeByBitmap[4];
-
-/** This array of group of each basis elements in the general multivector
-  * Use it to answer: 'what is the group of basis element bitmap [x]'? */
-extern const int e2ga_basisElementGroupByBitmap[4];
-
-/**
- * This struct can hold a general multivector.
- * 
- * The coordinates are stored in type double.
- * 
- * There are 3 coordinate groups:
- * group 0:1  (grade 0).
- * group 1:e1, e2  (grade 1).
- * group 2:e1^e2  (grade 2).
- * 
- * 4 doubles are allocated inside the struct.
- */
-typedef struct 
-{
-	/** group/grade usage (a bitmap which specifies which groups/grades are stored in 'c', below). */
-	int gu;
-	/** The coordinates (full). */
-	double c[4];
-} mv;
-
-/**
-This function alters the formatting of 'string()'.
-'format' = NULL will give you back the default.
-*/
-void e2ga_setStringFormat(const char *what, const char *format);
-
-extern const char *e2ga_string_fp; /* = \"%2.2f\" */
-extern const char *e2ga_string_start; /* = \"\" */
-extern const char *e2ga_string_end; /* = \"\" */
-extern const char *e2ga_string_mul; /* = \"*\" */
-extern const char *e2ga_string_wedge; /* = \"^\" */
-extern const char *e2ga_string_plus; /* = \" + \" */
-extern const char *e2ga_string_minus; /* = \" - \" */
-
-
-/** Writes value of 'V' to 'str' using float point precision 'fp' (e.g. %f). 'maxLength' is the length of 'str'. 'str' is returned. */
-const char *toString_mv(const mv *V, char *str, int maxLength, const char *fp);
-
-
-
-
-
-
-/** Sets 1 double to zero */
-void e2ga_double_zero_1(double *dst);
-/** Copies 1 double from 'src' to 'dst' */
-void e2ga_double_copy_1(double *dst, const double *src);
-/** Sets 2 doubles to zero */
-void e2ga_double_zero_2(double *dst);
-/** Copies 2 doubles from 'src' to 'dst' */
-void e2ga_double_copy_2(double *dst, const double *src);
-/** Sets 3 doubles to zero */
-void e2ga_double_zero_3(double *dst);
-/** Copies 3 doubles from 'src' to 'dst' */
-void e2ga_double_copy_3(double *dst, const double *src);
-/** Sets 4 doubles to zero */
-void e2ga_double_zero_4(double *dst);
-/** Copies 4 doubles from 'src' to 'dst' */
-void e2ga_double_copy_4(double *dst, const double *src);
-/** Sets 5 doubles to zero */
-void e2ga_double_zero_5(double *dst);
-/** Copies 5 doubles from 'src' to 'dst' */
-void e2ga_double_copy_5(double *dst, const double *src);
-/** Sets 6 doubles to zero */
-void e2ga_double_zero_6(double *dst);
-/** Copies 6 doubles from 'src' to 'dst' */
-void e2ga_double_copy_6(double *dst, const double *src);
-/** Sets 7 doubles to zero */
-void e2ga_double_zero_7(double *dst);
-/** Copies 7 doubles from 'src' to 'dst' */
-void e2ga_double_copy_7(double *dst, const double *src);
-/** Sets 8 doubles to zero */
-void e2ga_double_zero_8(double *dst);
-/** Copies 8 doubles from 'src' to 'dst' */
-void e2ga_double_copy_8(double *dst, const double *src);
-/** Sets 9 doubles to zero */
-void e2ga_double_zero_9(double *dst);
-/** Copies 9 doubles from 'src' to 'dst' */
-void e2ga_double_copy_9(double *dst, const double *src);
-/** Sets 10 doubles to zero */
-void e2ga_double_zero_10(double *dst);
-/** Copies 10 doubles from 'src' to 'dst' */
-void e2ga_double_copy_10(double *dst, const double *src);
-/** Sets 11 doubles to zero */
-void e2ga_double_zero_11(double *dst);
-/** Copies 11 doubles from 'src' to 'dst' */
-void e2ga_double_copy_11(double *dst, const double *src);
-/** Sets 12 doubles to zero */
-void e2ga_double_zero_12(double *dst);
-/** Copies 12 doubles from 'src' to 'dst' */
-void e2ga_double_copy_12(double *dst, const double *src);
-/** Sets 13 doubles to zero */
-void e2ga_double_zero_13(double *dst);
-/** Copies 13 doubles from 'src' to 'dst' */
-void e2ga_double_copy_13(double *dst, const double *src);
-/** Sets 14 doubles to zero */
-void e2ga_double_zero_14(double *dst);
-/** Copies 14 doubles from 'src' to 'dst' */
-void e2ga_double_copy_14(double *dst, const double *src);
-/** Sets 15 doubles to zero */
-void e2ga_double_zero_15(double *dst);
-/** Copies 15 doubles from 'src' to 'dst' */
-void e2ga_double_copy_15(double *dst, const double *src);
-/** Sets 16 doubles to zero */
-void e2ga_double_zero_16(double *dst);
-/** Copies 16 doubles from 'src' to 'dst' */
-void e2ga_double_copy_16(double *dst, const double *src);
-/** Sets N doubles to zero */
-void e2ga_double_zero_N(double *dst, int N);
-/** Copies N doubles from 'src' to 'dst' */
-void e2ga_double_copy_N(double *dst, const double *src, int N);
-/** Sets a mv to zero */
-void mv_setZero(mv *M);
-/** Sets a mv to a scalar value */
-void mv_setScalar(mv *M, double val);
-/** Sets a mv to the value in the array. 'gu' is a group usage bitmap. */
-void mv_setArray(mv *M, int gu, const double *arr);
-/** Copies a mv */
-void mv_copy(mv *dst, const mv *src);
-/** Allocates memory to store coordinate group 0 */
-void mv_reserveGroup_0(mv *A);
-/** Allocates memory to store coordinate group 1 */
-void mv_reserveGroup_1(mv *A);
-/** Allocates memory to store coordinate group 2 */
-void mv_reserveGroup_2(mv *A);
-/** Returns the 1 coordinate of 'A' */
-double mv_scalar(const mv *A);
-/** Returns the 1 coordinate of 'A' */
-double mv_double(const mv *A);
-/** Sets the 1 coordinate of 'A' */
-void mv_set_scalar(mv *A, double scalar_coord);
-/** Returns the e1 coordinate of 'A' */
-double mv_e1(const mv *A);
-/** Sets the e1 coordinate of 'A' */
-void mv_set_e1(mv *A, double e1_coord);
-/** Returns the e2 coordinate of 'A' */
-double mv_e2(const mv *A);
-/** Sets the e2 coordinate of 'A' */
-void mv_set_e2(mv *A, double e2_coord);
-/** Returns the e1^e2 coordinate of 'A' */
-double mv_e1_e2(const mv *A);
-/** Sets the e1^e2 coordinate of 'A' */
-void mv_set_e1_e2(mv *A, double e1_e2_coord);
-
-/** Returns absolute largest coordinate in mv. */
-double mv_largestCoordinate(const mv *x);
-/** Returns absolute largest coordinate in mv, and the bitmap of the corresponding basis blade (in 'bm'). */
-double mv_largestBasisBlade(const mv *x, unsigned int *bm);
-
-
-
-
-
-
-
-
-/**
- * Returns mv + mv.
- */
-void add_mv_mv(mv *_dst, const mv *a, const mv *b);
-/**
- * Returns a * b * reverse(a) using default metric. Only gives the correct result when the versor has a positive squared norm.
- * 
- */
-void applyUnitVersor_mv_mv(mv *_dst, const mv *a, const mv *b);
-/**
- * Returns a * b * inverse(a) using default metric.
- */
-void applyVersor_mv_mv(mv *_dst, const mv *a, const mv *b);
-/**
- * Returns a * b * reverse(a) using default metric. Only gives the correct result when the versor has a positive squared norm.
- * 
- */
-void applyVersorWI_mv_mv_mv(mv *_dst, const mv *a, const mv *b, const mv *c);
-/**
- * Returns a / b
- */
-void div_mv_double(mv *_dst, const mv *a, const double b);
-/**
- * Returns dual of mv using default metric.
- */
-void dual_mv(mv *_dst, const mv *a);
-/**
- * Returns whether input multivectors are equal up to an epsilon c.
- */
-int equals_mv_mv_double(const mv *a, const mv *b, const double c);
-/**
- * Returns grade groupBitmap of  mv.
- */
-void extractGrade_mv(mv *_dst, const mv *a, int groupBitmap);
-/**
- * Returns grade 0 of  mv.
- */
-void extractGrade0_mv(mv *_dst, const mv *a);
-/**
- * Returns grade 1 of  mv.
- */
-void extractGrade1_mv(mv *_dst, const mv *a);
-/**
- * Returns grade 2 of  mv.
- */
-void extractGrade2_mv(mv *_dst, const mv *a);
-/**
- * Returns geometric product of mv and mv.
- */
-void gp_mv_mv(mv *_dst, const mv *a, const mv *b);
-/**
- * Returns norm2 of mv using default metric.
- */
-double norm2_mv(const mv *a);
-/**
- * internal conversion function
- */
-double norm2_mv_returns_scalar(const mv *a);
-/**
- * Returns reverse of mv.
- */
-void reverse_mv(mv *_dst, const mv *a);
-/**
- * Returns mv - mv.
- */
-void subtract_mv_mv(mv *_dst, const mv *a, const mv *b);
-/**
- * Returns undual of mv using default metric.
- */
-void undual_mv(mv *_dst, const mv *a);
-/**
- * Returns versor inverse of a using default metric.
- */
-void versorInverse_mv(mv *_dst, const mv *a);
-#endif /* _E2GA_H_ */

c/e2ga.xml

-<?xml version="1.0" encoding="utf-8" ?>
-
-<g25spec license="custom" language="c" namespace="e2ga"
-         coordStorage="array" defaultOperatorBindings="false"
-         dimension="2" testSuite="false" reportUsage="false" 
-         gmvCode="expand" parser="none"
-         copyright="Copyright (c) 2012 Barry Schwartz">
-
-  <customLicense>
-    Copyright (c) 2012 Barry Schwartz
-
-    Permission is hereby granted, free of charge, to any person obtaining a copy
-    of this software and associated documentation files (the "Software"), to deal
-    in the Software without restriction, including without limitation the rights
-    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-    copies of the Software, and to permit persons to whom the Software is
-    furnished to do so, subject to the following conditions:
-
-    The above copyright notice and this permission notice shall be included in
-    all copies or substantial portions of the Software.
-
-    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-    THE SOFTWARE.
-  </customLicense>
-
-  <inline constructors="false" set="false" assign="false"
-          operators="false" functions="false"/>
-  <floatType type="double"/>
-
-  <basisVectorNames name1="e1" name2="e2"/>
-  <metric name="default"> e1.e1 = e2.e2 = 1 </metric>
-
-  <mv name="mv" compress="byGroup" coordinateOrder="custom" memAlloc="full">
-    <group> scalar </group>
-    <group> e1 e2  </group>
-    <group> e1^e2  </group>
-  </mv>
-
-  <function name="add" arg1="mv" arg2="mv"/>
-  <function name="applyUnitVersor" arg1="mv" arg2="mv"/>
-  <function name="applyVersor" arg1="mv" arg2="mv"/>
-  <function name="applyVersorWI" arg1="mv" arg2="mv" arg3="mv"/>
-  <function name="div" arg1="mv" arg2="double"/>
-  <function name="dual" arg1="mv"/>
-  <function name="equals" arg1="mv" arg2="mv" arg3="double"/>
-  <function name="extractGrade" arg1="mv"/>
-  <function name="extractGrade0" arg1="mv"/>
-  <function name="extractGrade1" arg1="mv"/>
-  <function name="extractGrade2" arg1="mv"/>
-  <function name="gp" arg1="mv" arg2="mv"/>
-  <function name="norm2" arg1="mv"/>
-  <function name="reverse" arg1="mv"/>
-  <function name="subtract" arg1="mv" arg2="mv"/>
-  <function name="undual" arg1="mv"/>
-  <function name="versorInverse" arg1="mv"/>
-
-</g25spec>
+/*
+Copyright (c) 2012 Barry Schwartz
+*/
+/*
+
+    Copyright (c) 2012 Barry Schwartz
+
+    Permission is hereby granted, free of charge, to any person obtaining a copy
+    of this software and associated documentation files (the "Software"), to deal
+    in the Software without restriction, including without limitation the rights
+    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+    copies of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included in
+    all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+    THE SOFTWARE.
+  
+*/
+#include <stdio.h>
+#include "ga.h"
+
+const int ga_spaceDim = 5;
+const int ga_nbGroups = 6;
+const int ga_metricEuclidean = 0;
+const char *ga_basisVectorNames[5] = {
+	"no", "e1", "e2", "e3", "ni"
+};
+const int ga_grades[] = {GRADE_0, GRADE_1, GRADE_2, GRADE_3, GRADE_4, GRADE_5, 0, 0, 0, 0, 0, 0};
+const int ga_groups[] = {GROUP_0, GROUP_1, GROUP_2, GROUP_3, GROUP_4, GROUP_5};
+const int ga_groupSize[6] = {
+	1, 5, 10, 10, 5, 1
+};
+const int ga_mvSize[64] = {
+	0, 1, 5, 6, 10, 11, 15, 16, 10, 11, 15, 16, 20, 21, 25, 26, 5, 6, 10, 11, 
+	15, 16, 20, 21, 15, 16, 20, 21, 25, 26, 30, 31, 1, 2, 6, 7, 11, 12, 16, 17, 
+	11, 12, 16, 17, 21, 22, 26, 27, 6, 7, 11, 12, 16, 17, 21, 22, 16, 17, 21, 22, 
+	26, 27, 31, 32};
+const int ga_basisElements[32][6] = {
+	{-1},
+	{0, -1},
+	{1, -1},
+	{2, -1},
+	{3, -1},
+	{4, -1},
+	{0, 1, -1},
+	{0, 2, -1},
+	{1, 2, -1},
+	{0, 3, -1},
+	{1, 3, -1},
+	{2, 3, -1},
+	{0, 4, -1},
+	{1, 4, -1},
+	{2, 4, -1},
+	{3, 4, -1},
+	{0, 1, 2, -1},
+	{0, 1, 3, -1},
+	{0, 2, 3, -1},
+	{1, 2, 3, -1},
+	{0, 1, 4, -1},
+	{0, 2, 4, -1},
+	{1, 2, 4, -1},
+	{0, 3, 4, -1},
+	{1, 3, 4, -1},
+	{2, 3, 4, -1},
+	{0, 1, 2, 3, -1},
+	{0, 1, 2, 4, -1},
+	{0, 1, 3, 4, -1},
+	{0, 2, 3, 4, -1},
+	{1, 2, 3, 4, -1},
+	{0, 1, 2, 3, 4, -1}
+};
+const double ga_basisElementSignByIndex[32] =
+	{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+const double ga_basisElementSignByBitmap[32] =
+	{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+const int ga_basisElementIndexByBitmap[32] =
+	{0, 1, 2, 6, 3, 7, 8, 16, 4, 9, 10, 17, 11, 18, 19, 26, 5, 12, 13, 20, 14, 21, 22, 27, 15, 23, 24, 28, 25, 29, 30, 31};
+const int ga_basisElementBitmapByIndex[32] =
+	{0, 1, 2, 4, 8, 16, 3, 5, 6, 9, 10, 12, 17, 18, 20, 24, 7, 11, 13, 14, 19, 21, 22, 25, 26, 28, 15, 23, 27, 29, 30, 31};
+const int ga_basisElementGradeByBitmap[32] =
+	{0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5};
+const int ga_basisElementGroupByBitmap[32] =
+	{0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5};
+
+const char *g_gaTypenames[] = 
+{
+	"mvec2",
+	"vec2",
+	"bvec2"
+};
+
+
+void ga_double_zero_1(double *dst) {
+	dst[0]=0.0;
+}
+void ga_double_copy_1(double *dst, const double *src) {
+	dst[0] = src[0];
+}
+void ga_double_zero_2(double *dst) {
+	dst[0]=dst[1]=0.0;
+}
+void ga_double_copy_2(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+}
+void ga_double_zero_3(double *dst) {
+	dst[0]=dst[1]=dst[2]=0.0;
+}
+void ga_double_copy_3(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+}
+void ga_double_zero_4(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=0.0;
+}
+void ga_double_copy_4(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+}
+void ga_double_zero_5(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=0.0;
+}
+void ga_double_copy_5(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+}
+void ga_double_zero_6(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=0.0;
+}
+void ga_double_copy_6(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+}
+void ga_double_zero_7(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=0.0;
+}
+void ga_double_copy_7(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+}
+void ga_double_zero_8(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=0.0;
+}
+void ga_double_copy_8(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+}
+void ga_double_zero_9(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=0.0;
+}
+void ga_double_copy_9(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+}
+void ga_double_zero_10(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=0.0;
+}
+void ga_double_copy_10(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+}
+void ga_double_zero_11(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=0.0;
+}
+void ga_double_copy_11(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+}
+void ga_double_zero_12(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=0.0;
+}
+void ga_double_copy_12(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+	dst[11] = src[11];
+}
+void ga_double_zero_13(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=0.0;
+}
+void ga_double_copy_13(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+	dst[11] = src[11];
+	dst[12] = src[12];
+}
+void ga_double_zero_14(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=0.0;
+}
+void ga_double_copy_14(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+	dst[11] = src[11];
+	dst[12] = src[12];
+	dst[13] = src[13];
+}
+void ga_double_zero_15(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=dst[14]=0.0;
+}
+void ga_double_copy_15(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+	dst[11] = src[11];
+	dst[12] = src[12];
+	dst[13] = src[13];
+	dst[14] = src[14];
+}
+void ga_double_zero_16(double *dst) {
+	dst[0]=dst[1]=dst[2]=dst[3]=dst[4]=dst[5]=dst[6]=dst[7]=dst[8]=dst[9]=dst[10]=dst[11]=dst[12]=dst[13]=dst[14]=dst[15]=0.0;
+}
+void ga_double_copy_16(double *dst, const double *src) {
+	dst[0] = src[0];
+	dst[1] = src[1];
+	dst[2] = src[2];
+	dst[3] = src[3];
+	dst[4] = src[4];
+	dst[5] = src[5];
+	dst[6] = src[6];
+	dst[7] = src[7];
+	dst[8] = src[8];
+	dst[9] = src[9];
+	dst[10] = src[10];
+	dst[11] = src[11];
+	dst[12] = src[12];
+	dst[13] = src[13];
+	dst[14] = src[14];
+	dst[15] = src[15];
+}
+/** Sets N doubles to zero */
+void ga_double_zero_N(double *dst, int N) {
+	int i = 0;
+	while ((N-i) > 16) {
+		ga_double_zero_16(dst + i);
+		i += 16;
+	}
+	for (; i < N; i++)
+		dst[i] = 0.0;
+}
+/** Copies N doubles from 'src' to 'dst' */
+void ga_double_copy_N(double *dst, const double *src, int N) {
+	int i = 0;
+	while ((N-i) > 16) {
+		ga_double_copy_16(dst + i, src + i);
+		i += 16;
+	}
+	for (; i < N; i++)
+		dst[i] = src[i];
+}
+
+/** This function is not for external use. It compresses arrays of floats for storage in multivectors. */
+void compress(const double *c, double *cc, int *cgu, double epsilon, int gu) {
+	int i, j, ia = 0, ib = 0, f, s;
+	*cgu = 0;
+
+	/* for all grade parts... */
+	for (i = 0; i <= 6; i++) {
+		/* check if grade part has memory use: */
+		if (!(gu & (1 << i))) continue;
+
+		/* check if abs coordinates of grade part are all < epsilon */
+		s = ga_groupSize[i];
+		j = ia + s;
+		f = 0;
+		for (; ia < j; ia++)
+			if (fabs(c[ia]) > epsilon) {f = 1; break;}
+		ia = j;
+		if (f) {
+						ga_double_copy_N(cc + ib, c + ia - s, s);
+			ib += s;
+			*cgu |= (1 << i);
+		}
+	}
+}
+
+/** This function is not for external use. It decompresses the coordinates stored in this */
+void expand(const double *ptrs[6], const mvec3c *src) {
+	const double *c = src->c;
+	
+	if (src->gu & 1) {
+		ptrs[0] =  c;
+		c += 1;
+	}
+	else ptrs[0] = NULL;	
+	if (src->gu & 2) {
+		ptrs[1] =  c;
+		c += 5;
+	}
+	else ptrs[1] = NULL;	
+	if (src->gu & 4) {
+		ptrs[2] =  c;
+		c += 10;
+	}
+	else ptrs[2] = NULL;	
+	if (src->gu & 8) {
+		ptrs[3] =  c;
+		c += 10;
+	}
+	else ptrs[3] = NULL;	
+	if (src->gu & 16) {
+		ptrs[4] =  c;
+		c += 5;
+	}
+	else ptrs[4] = NULL;	
+	if (src->gu & 32) {
+		ptrs[5] =  c;
+	}
+	else ptrs[5] = NULL;	
+}
+
+
+void swapPointer(void **P1, void **P2)
+{
+	void *tmp = *P1;
+	*P1 = *P2;
+	*P2 = tmp;
+}
+
+/* 
+These strings determine how the output of string() is formatted.
+You can alter them at runtime using ga_setStringFormat().
+*/
+ 
+const char *ga_string_fp = "%2.2f"; 
+const char *ga_string_start = ""; 
+const char *ga_string_end = ""; 
+const char *ga_string_mul = "*"; 
+const char *ga_string_wedge = "^"; 
+const char *ga_string_plus = " + "; 
+const char *ga_string_minus = " - "; 
+
+void ga_setStringFormat(const char *what, const char *format) {
+ 
+	if (!strcmp(what, "fp")) 
+		ga_string_fp = (format) ? format : "%2.2f"; 
+	else if (!strcmp(what, "start")) 
+		ga_string_start = (format) ? format : ""; 
+	else if (!strcmp(what, "end")) 
+		ga_string_end = (format) ? format : ""; 
+	else if (!strcmp(what, "mul")) 
+		ga_string_mul = (format) ? format : "*"; 
+	else if (!strcmp(what, "wedge")) 
+		ga_string_wedge = (format) ? format : "^"; 
+	else if (!strcmp(what, "plus")) 
+		ga_string_plus = (format) ? format : " + "; 
+	else if (!strcmp(what, "minus")) 
+		ga_string_minus = (format) ? format : " - ";
+}
+
+const char *toString_mvec2(const mvec2 *V, char *str, int maxLength, const char *fp)
+{
+	mvec3c tmp;
+	mvec2_to_mvec3c(&tmp,V);
+	return toString_mvec3c(&tmp, str, maxLength, fp);
+}
+const char *toString_vec2(const vec2 *V, char *str, int maxLength, const char *fp)
+{
+	mvec3c tmp;
+	vec2_to_mvec3c(&tmp,V);
+	return toString_mvec3c(&tmp, str, maxLength, fp);
+}
+const char *toString_bvec2(const bvec2 *V, char *str, int maxLength, const char *fp)
+{
+	mvec3c tmp;
+	bvec2_to_mvec3c(&tmp,V);
+	return toString_mvec3c(&tmp, str, maxLength, fp);
+}
+
+#ifdef WIN32
+#define snprintf _snprintf
+#pragma warning(disable:4996) /* quit your whining already */
+#endif /* WIN32 */
+const char *toString_mvec3c(const mvec3c *V, char *str, int maxLength, const char *fp) 
+{
+	int dummyArg = 0; /* prevents compiler warning on some platforms */
+	int stdIdx = 0, l;
+	char tmpBuf[256], tmpFloatBuf[256]; 
+	int i, j, k = 0, bei, ia = 0, s = ga_mvSize[V->gu], p = 0, cnt = 0;
+
+	/* set up the floating point precision */
+	if (fp == NULL) fp = ga_string_fp;
+
+	/* start the string */
+	l = snprintf(tmpBuf, 256, "%s", ga_string_start);
+	if (stdIdx + l <= maxLength) {
+		strcpy(str + stdIdx, tmpBuf);
+		stdIdx += l;
+	}
+	else {
+		snprintf(str, maxLength, "toString_mvec3c: buffer too small");
+		return str;
+	}
+
+	/* print all coordinates */
+	for (i = 0; i <= 6; i++) {
+		if (V->gu & (1 << i)) {
+			for (j = 0; j < ga_groupSize[i]; j++) {
+				double coord = (double)ga_basisElementSignByIndex[ia] * V->c[k];
+				/* goal: print [+|-]V->c[k][* basisVector1 ^ ... ^ basisVectorN] */			
+				snprintf(tmpFloatBuf, 256, fp, fabs(coord));
+				if (atof(tmpFloatBuf) != 0.0) {
+					l = 0;
+
+					/* print [+|-] */
+					l += snprintf(tmpBuf + l, 256-l, "%s", (coord >= 0.0) 
+						? (cnt ? ga_string_plus : "")
+						: ga_string_minus);
+						
+					/* print obj.m_c[k] */
+					l += snprintf(tmpBuf + l, 256-l, tmpFloatBuf, dummyArg);
+
+					if (i) { /* if not grade 0, print [* basisVector1 ^ ... ^ basisVectorN] */
+						l += snprintf(tmpBuf + l, 256-l, "%s", ga_string_mul);
+
+						/* print all basis vectors */
+						bei = 0;
+						while (ga_basisElements[ia][bei] >= 0) {
+							l += snprintf(tmpBuf + l, 256-l, "%s%s", (bei) ? ga_string_wedge : "", 
+							 ga_basisVectorNames[ga_basisElements[ia][bei]]);
+							 bei++;
+						}
+					}
+
+					/* copy all to 'str' */
+					if (stdIdx + l <= maxLength) {
+						strcpy(str + stdIdx, tmpBuf);
+						stdIdx += l;
+					}
+					else {
+						snprintf(str, maxLength, "toString_mvec3c: buffer too small");
+						return str;
+					}
+					cnt++;
+				}
+				k++; ia++;
+			}
+		}
+		else ia += ga_groupSize[i];
+	}
+
+    /* if no coordinates printed: 0 */
+	l = 0;
+	if (cnt == 0) {
+		l += snprintf(tmpBuf + l, 256-l, "0");
+	}
+
+    /* end the string */
+	l += snprintf(tmpBuf + l, 256-l, "%s", ga_string_end);
+	if (stdIdx + l <= maxLength) {
+		strcpy(str + stdIdx, tmpBuf);
+		stdIdx += l;
+		return str;
+	}
+	else {
+		snprintf(str, maxLength, "toString_mvec3c: buffer too small\n");
+		return str;
+	}
+}
+
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 0 -> group 0)
+ */
+void gp_default_0_0_0(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 1 -> group 1)
+ */
+void gp_default_0_1_1(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 2 -> group 2)
+ */
+void gp_default_0_2_2(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 3 -> group 3)
+ */
+void gp_default_0_3_3(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 4 -> group 4)
+ */
+void gp_default_0_4_4(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 0  x  group 5 -> group 5)
+ */
+void gp_default_0_5_5(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 0 -> group 1)
+ */
+void gp_default_1_0_1(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 1 -> group 0)
+ */
+void gp_default_1_1_0(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 1 -> group 2)
+ */
+void gp_default_1_1_2(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 2 -> group 1)
+ */
+void gp_default_1_2_1(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 2 -> group 3)
+ */
+void gp_default_1_2_3(const double *A, const double *B, double *C);
+/**
+ * Computes the partial geometric product of two multivectors (group 1  x  group 3 -> group 2)