Commits

Barry Schwartz committed 34d693a

Checking in files generated by g25, because it is an obscure tool that might not be available someday.

  • Participants
  • Parent commits de630a6

Comments (0)

Files changed (3)

 CPPFLAGS += -I"${top_srcdir}/c" -I"${top_builddir}/c"
 
 lib_LTLIBRARIES = libebeno.la
-libebeno_la_SOURCES = c/nan.c ${builddir}/c/e2ga.c					\
+libebeno_la_SOURCES = c/nan.c ${srcdir}/c/e2ga.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
 
-${builddir}/c/e2ga.h: ${srcdir}/c/e2ga.xml
+${srcdir}/c/e2ga.h: ${srcdir}/c/e2ga.xml
 	$(G25) $<
-	-$(MKDIR_P) ${builddir}/c
-	mv e2ga.c e2ga.h c
-${builddir}/c/e2ga.c: ${builddir}/c/e2ga.h ${srcdir}/c/e2ga.xml
+#	-$(MKDIR_P) ${builddir}/c
+	mv e2ga.c e2ga.h ${srcdir}/c
+${srcdir}/c/e2ga.c: ${srcdir}/c/e2ga.h ${srcdir}/c/e2ga.xml
 	@$(force_rebuild)
 
-EXTRA_DIST = ${srcdir}/c/e2ga.xml c/e2ga.c c/e2ga.h
+EXTRA_DIST = ${srcdir}/c/e2ga.xml ${srcdir}/c/e2ga.c ${srcdir}/c/e2ga.h
 
 #--------------------------------------------------------------------------
 
 EXTRA_DIST += COPYING INSTALL
 MOSTLYCLEANFILES = $(foreach f, $(MODFILES), $(call modfile,${f}))	\
 	Doxyfile
-MAINTAINERCLEANFILES = c/e2ga.c c/e2ga.h								\
+MAINTAINERCLEANFILES = ${srcdir}/c/e2ga.c ${srcdir}/c/e2ga.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
+/*
+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);
+	}
+}
+/*
+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_ */