Commits

Barry Schwartz  committed 940d6a2

Geometric algebra work.

  • Participants
  • Parent commits fc4d0eb

Comments (0)

Files changed (5)

 	_dst->c[1] = (a->c[0]*b->c[1]-a->c[1]*b->c[0]);
 
 }
+void gp_rotor2_vec2(vec2 *_dst, const rotor2 *a, const vec2 *b)
+{
+	_dst->c[0] = (a->c[0]*b->c[0]+a->c[1]*b->c[1]);
+	_dst->c[1] = (a->c[0]*b->c[1]-a->c[1]*b->c[0]);
+
+}
+void gp_vec2_rotor2(vec2 *_dst, const vec2 *a, const rotor2 *b)
+{
+	_dst->c[0] = (a->c[0]*b->c[0]-a->c[1]*b->c[1]);
+	_dst->c[1] = (a->c[0]*b->c[1]+a->c[1]*b->c[0]);
+
+}
+void gp_rotor2_rotor2(rotor2 *_dst, const rotor2 *a, const rotor2 *b)
+{
+	_dst->c[0] = (a->c[0]*b->c[0]-a->c[1]*b->c[1]);
+	_dst->c[1] = (a->c[0]*b->c[1]+a->c[1]*b->c[0]);
+
+}
+double gp_bvec2_bvec2(const bvec2 *a, const bvec2 *b)
+{
+	return -a->c[0]*b->c[0];
+
+}
+void gp_mvec3c_double(mvec3c *_dst, const mvec3c *a, const double b)
+{
+	double c[32];
+	const double* _a[6];
+	const double* _b[1] = {&b};
+	expand(_a, a);
+	ga_double_zero_N(c, 32);
+	if (a->gu & 1) {
+			gp_default_0_0_0(_a[0], _b[0], c + 0);
+	}
+	if (a->gu & 2) {
+			gp_default_1_0_1(_a[1], _b[0], c + 1);
+	}
+	if (a->gu & 4) {
+			gp_default_2_0_2(_a[2], _b[0], c + 6);
+	}
+	if (a->gu & 8) {
+			gp_default_3_0_3(_a[3], _b[0], c + 16);
+	}
+	if (a->gu & 16) {
+			gp_default_4_0_4(_a[4], _b[0], c + 26);
+	}
+	if (a->gu & 32) {
+			gp_default_5_0_5(_a[5], _b[0], c + 31);
+	}
+	compress(c, _dst->c, &(_dst->gu), 0.0, 63);
+}
+void gp_mvec2_double(mvec2 *_dst, const mvec2 *a, const double b)
+{
+	_dst->c[0] = a->c[0]*b;
+	_dst->c[1] = a->c[1]*b;
+	_dst->c[2] = a->c[2]*b;
+	_dst->c[3] = a->c[3]*b;
+
+}
+void gp_vec2_double(vec2 *_dst, const vec2 *a, const double b)
+{
+	_dst->c[0] = a->c[0]*b;
+	_dst->c[1] = a->c[1]*b;
+
+}
+void gp_bvec2_double(bvec2 *_dst, const bvec2 *a, const double b)
+{
+	_dst->c[0] = a->c[0]*b;
+
+}
+void gp_rotor2_double(rotor2 *_dst, const rotor2 *a, const double b)
+{
+	_dst->c[0] = a->c[0]*b;
+	_dst->c[1] = a->c[1]*b;
+
+}
+void gp_double_mvec3c(mvec3c *_dst, const double a, const mvec3c *b)
+{
+	double c[32];
+	const double* _a[1] = {&a};
+	const double* _b[6];
+	expand(_b, b);
+	ga_double_zero_N(c, 32);
+		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 + 6);
+		}
+		if (b->gu & 8) {
+			gp_default_0_3_3(_a[0], _b[3], c + 16);
+		}
+		if (b->gu & 16) {
+			gp_default_0_4_4(_a[0], _b[4], c + 26);
+		}
+		if (b->gu & 32) {
+			gp_default_0_5_5(_a[0], _b[5], c + 31);
+		}
+	compress(c, _dst->c, &(_dst->gu), 0.0, 63);
+}
+void gp_double_mvec2(mvec2 *_dst, const double a, const mvec2 *b)
+{
+	_dst->c[0] = a*b->c[0];
+	_dst->c[1] = a*b->c[1];
+	_dst->c[2] = a*b->c[2];
+	_dst->c[3] = a*b->c[3];
+
+}
+void gp_double_vec2(vec2 *_dst, const double a, const vec2 *b)
+{
+	_dst->c[0] = a*b->c[0];
+	_dst->c[1] = a*b->c[1];
+
+}
+void gp_double_bvec2(bvec2 *_dst, const double a, const bvec2 *b)
+{
+	_dst->c[0] = a*b->c[0];
+
+}
+void gp_double_rotor2(rotor2 *_dst, const double a, const rotor2 *b)
+{
+	_dst->c[0] = a*b->c[0];
+	_dst->c[1] = a*b->c[1];
+
+}
 void div_mvec2_double(mvec2 *_dst, const mvec2 *a, const double b)
 {
 	_dst->c[0] = a->c[0]/((b));
  */
 void gp_vec2_vec2(rotor2 *_dst, const vec2 *a, const vec2 *b);
 /**
+ * Returns geometric product of rotor2 and vec2.
+ */
+void gp_rotor2_vec2(vec2 *_dst, const rotor2 *a, const vec2 *b);
+/**
+ * Returns geometric product of vec2 and rotor2.
+ */
+void gp_vec2_rotor2(vec2 *_dst, const vec2 *a, const rotor2 *b);
+/**
+ * Returns geometric product of rotor2 and rotor2.
+ */
+void gp_rotor2_rotor2(rotor2 *_dst, const rotor2 *a, const rotor2 *b);
+/**
+ * Returns geometric product of bvec2 and bvec2.
+ */
+double gp_bvec2_bvec2(const bvec2 *a, const bvec2 *b);
+/**
+ * Returns geometric product of mvec3c and double.
+ */
+void gp_mvec3c_double(mvec3c *_dst, const mvec3c *a, const double b);
+/**
+ * Returns geometric product of mvec2 and double.
+ */
+void gp_mvec2_double(mvec2 *_dst, const mvec2 *a, const double b);
+/**
+ * Returns geometric product of vec2 and double.
+ */
+void gp_vec2_double(vec2 *_dst, const vec2 *a, const double b);
+/**
+ * Returns geometric product of bvec2 and double.
+ */
+void gp_bvec2_double(bvec2 *_dst, const bvec2 *a, const double b);
+/**
+ * Returns geometric product of rotor2 and double.
+ */
+void gp_rotor2_double(rotor2 *_dst, const rotor2 *a, const double b);
+/**
+ * Returns geometric product of double and mvec3c.
+ */
+void gp_double_mvec3c(mvec3c *_dst, const double a, const mvec3c *b);
+/**
+ * Returns geometric product of double and mvec2.
+ */
+void gp_double_mvec2(mvec2 *_dst, const double a, const mvec2 *b);
+/**
+ * Returns geometric product of double and vec2.
+ */
+void gp_double_vec2(vec2 *_dst, const double a, const vec2 *b);
+/**
+ * Returns geometric product of double and bvec2.
+ */
+void gp_double_bvec2(bvec2 *_dst, const double a, const bvec2 *b);
+/**
+ * Returns geometric product of double and rotor2.
+ */
+void gp_double_rotor2(rotor2 *_dst, const double a, const rotor2 *b);
+/**
  * Returns a / b
  */
 void div_mvec2_double(mvec2 *_dst, const mvec2 *a, const double b);
   <function name="gp" arg1="mvec3c" arg2="mvec3c"/>
   <function name="gp" arg1="mvec2" arg2="mvec2"/>
   <function name="gp" arg1="vec2" arg2="vec2"/>
+  <function name="gp" arg1="rotor2" arg2="vec2"/>
+  <function name="gp" arg1="vec2" arg2="rotor2"/>
+  <function name="gp" arg1="rotor2" arg2="rotor2"/>
+  <function name="gp" arg1="bvec2" arg2="bvec2"/>
+  <function name="gp" arg1="mvec3c" arg2="double"/>
+  <function name="gp" arg1="mvec2" arg2="double"/>
+  <function name="gp" arg1="vec2" arg2="double"/>
+  <function name="gp" arg1="bvec2" arg2="double"/>
+  <function name="gp" arg1="rotor2" arg2="double"/>
+  <function name="gp" arg1="double" arg2="mvec3c"/>
+  <function name="gp" arg1="double" arg2="mvec2"/>
+  <function name="gp" arg1="double" arg2="vec2"/>
+  <function name="gp" arg1="double" arg2="bvec2"/>
+  <function name="gp" arg1="double" arg2="rotor2"/>
 
   <function name="div" arg1="mvec2" arg2="double"/>
 

File pure/ebeno/geometric_algebra.pure

     _ = throw (bad_multivector_grade grade mv);
   end;
 
+//-------------------------------------------------------------------------
+
 // Geometric products.
+
 mvec2 c1@{a0,a1,a2,a3} ::* mvec2 c2@{b0,b1,b2,b3} =
   if dmatrixp c1 && dmatrixp c2 then
     (gp_mvec2_mvec2 (pointer c, pointer cp1, pointer cp2) $$
            a0*b1 + a1*b0 - a2*b3 + a3*b2,
            a0*b2 + a1*b3 + a2*b0 - a3*b1,
            a0*b3 + a1*b2 - a2*b1 + a3*b0};
+
 vec2 c1@{a0,a1} ::* vec2 c2@{b0,b1} =
   if dmatrixp c1 && dmatrixp c2 then
     (gp_vec2_vec2 (pointer c, pointer cp1, pointer cp2) $$
     mvec2 {a0*b0 + a1*b1, 0, 0, a0*b1 - a1*b0};
 
 //-------------------------------------------------------------------------
-
-/*
-infix (*) ⋀ :^: ; // Outer product.
-infix (*) ⋅ :.: ; // Dot product.
-(:^:) = (⋀);
-(:.:) = (⋅);
-*/

File pure/generate_ga_code.pure

   "       cp2::dmatrix = pack c2;\n" +
   "     end)\n" +
   "  else\n" +
-  "    " + str type_out + " {" + elemlist_out + "}\n" +
-  "  if ~cmatrixp c1 && ~cmatrixp c2;\n"
+  "    " + str type_out + " {" + elemlist_out + "};\n"
   when
     ebeno_func = "add_" + str type1 + "_" + str type2;
     a_elem = ["a" + str i | i = indices1];
     b_elem = ["b" + str i | i = indices2];
     a_elemlist = foldl (\el s -> el + ("," + s)) (a_elem!0) (drop 1 a_elem);
-    a_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (a_elem!0 + "::double") (drop 1 a_elem);
     b_elemlist = foldl (\el s -> el + ("," + s)) (b_elem!0) (drop 1 b_elem);
-    b_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (b_elem!0 + "::double") (drop 1 b_elem);
     m::int = foldl max (-1) (list indices1 + list indices2);
     elem_out = map elem_sum (0..m)
       with
   end;
 
 make_multivector_addition type1 type2 indices1::list indices2 type_out =
-  str type1 + " c1@{" + a_elemlist + "} ::+ r::real |\n" +
-  "r::real ::+ " + str type1 + " c1@{" + a_elemlist + "} =\n" +
-  "  if dmatrixp c1 then\n" +
+  str type1 + " c1@{" + a_elemlist + "} ::+ r |\n" +
+  "r ::+ " + str type1 + " c1@{" + a_elemlist + "} =\n" +
+  "  if doublep r && dmatrixp c1 then\n" +
   "    (" + ebeno_func + " (pointer c, pointer cp1, double r) $$ " + str type_out + " c\n" +
   "     when\n" +
   "       c::dmatrix = rowvector $ repeatn (#c1) 0.0;\n" +
   "     end)\n" +
   "  else\n" +
   "    " + str type_out + " {" + elemlist_out + "}\n" +
-  "  if ~cmatrixp c1;\n"
+  "  if ~multivectorp r;\n"
   when
     ebeno_func = "add_" + str type1 + "_double";
     a_elem = ["a" + str i | i = indices1];
     a_elemlist = foldl (\el s -> el + ("," + s)) (a_elem!0) (drop 1 a_elem);
-    a_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (a_elem!0 + "::double") (drop 1 a_elem);
     elemlist_out = foldl (\el s -> el + ("," + s)) (a_elem!0 + " + r") (drop 1 a_elem);
   end;
 
 //-------------------------------------------------------------------------
 
-// FIXME: Consider calling the C code instead of compiling ‘double’
-// versions.
-
 // A function that creates multivector subtraction functions.
 
 make_multivector_subtraction type1 type2 indices1::list indices2::list type_out =
-  sprintf "%s {%s} ::- %s {%s} |\n"
-  (str type1, a_elemlist_double, str type2, b_elemlist_double) +
-  sprintf "%s {%s} ::- %s {%s} = %s {%s};\n"
-  (str type1, a_elemlist, str type2, b_elemlist, str type_out, elemlist_out)
+  (if type1 === type2 then
+     str type1 + " c1@{" + a_elemlist + "} ::- " +
+     str type2 + " c2@{" + b_elemlist + "} =\n"
+   else
+     str type1 + " c1@{" + a_elemlist + "} ::- " +
+     str type2 + " c2@{" + b_elemlist + "} =\n") +
+  "  if dmatrixp c1 && dmatrixp c2 then\n" +
+  "    (" + ebeno_func + " (pointer c, pointer cp1, pointer cp2) $$ " + str type_out + " c\n" +
+  "     when\n" +
+  "       c::dmatrix = rowvector $ repeatn " + str (m + 1) + " 0.0;\n" +
+  "       cp1::dmatrix = pack c1;\n" +
+  "       cp2::dmatrix = pack c2;\n" +
+  "     end)\n" +
+  "  else\n" +
+  "    " + str type_out + " {" + elemlist_out + "};\n\n" +
+  (if type1 === type2 then
+     ""
+   else
+     str type2 + " c2@{" + b_elemlist + "} ::- " +
+     str type1 + " c1@{" + a_elemlist + "} =\n" +
+     "  if dmatrixp c1 && dmatrixp c2 then\n" +
+     "    (" + ebeno_func_inverse + " (pointer c, pointer cp2, pointer cp1) $$ " + str type_out + " c\n" +
+     "     when\n" +
+     "       c::dmatrix = rowvector $ repeatn " + str (m + 1) + " 0.0;\n" +
+     "       cp1::dmatrix = pack c1;\n" +
+     "       cp2::dmatrix = pack c2;\n" +
+     "     end)\n" +
+     "  else\n" +
+     "    " + str type_out + " {" + elemlist_out_inverse + "};\n")
   when
+    ebeno_func = "subtract_" + str type1 + "_" + str type2;
+    ebeno_func_inverse = "subtract_" + str type2 + "_" + str type1;
     a_elem = ["a" + str i | i = indices1];
     b_elem = ["b" + str i | i = indices2];
     a_elemlist = foldl (\el s -> el + ("," + s)) (a_elem!0) (drop 1 a_elem);
-    a_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (a_elem!0 + "::double") (drop 1 a_elem);
     b_elemlist = foldl (\el s -> el + ("," + s)) (b_elem!0) (drop 1 b_elem);
-    b_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (b_elem!0 + "::double") (drop 1 b_elem);
     m::int = foldl max (-1) (list indices1 + list indices2);
     elem_out = map elem_sum (0..m)
       with
           end;
       end;
     elemlist_out = foldl (\el s -> el + ("," + s)) (elem_out!0) (drop 1 elem_out);
+    elem_out_inverse = map elem_sum (0..m)
+      with
+        elem_sum i =
+          case index indices2 i, index indices1 i of
+            -1,-1 = nan;
+            j,-1 = b_elem!j;
+            -1,j = "-" + a_elem!j;
+            j,k = b_elem!j + " - " + a_elem!k;
+          end;
+      end;
+    elemlist_out_inverse = foldl (\el s -> el + ("," + s)) (elem_out_inverse!0) (drop 1 elem_out_inverse);
   end;
+
 make_multivector_subtraction type1 type2 indices1::list indices2 type_out =
-  sprintf "%s {%s} ::- r::double |\n" (str type1, a_elemlist_double) +
-  sprintf "%s {%s} ::- r::real = %s {%s};\n" (str type1, a_elemlist, str type_out, elemlist_out)
+  str type1 + " c1@{" + a_elemlist + "} ::- r =\n" +
+  "  if doublep r && dmatrixp c1 then\n" +
+  "    (" + ebeno_func + " (pointer c, pointer cp1, double r) $$ " + str type_out + " c\n" +
+  "     when\n" +
+  "       c::dmatrix = rowvector $ repeatn (#c1) 0.0;\n" +
+  "       cp1::dmatrix = pack c1;\n" +
+  "     end)\n" +
+  "  else\n" +
+  "    " + str type_out + " {" + elemlist_out + "}\n" +
+  "  if ~multivectorp r;\n\n" +
+  "r ::- " + str type1 + " c1@{" + a_elemlist + "} =\n" +
+  "  if doublep r && dmatrixp c1 then\n" +
+  "    (" + ebeno_func_inverse + " (pointer c, double r, pointer cp1) $$ " + str type_out + " c\n" +
+  "     when\n" +
+  "       c::dmatrix = rowvector $ repeatn (#c1) 0.0;\n" +
+  "       cp1::dmatrix = pack c1;\n" +
+  "     end)\n" +
+  "  else\n" +
+  "    " + str type_out + " {" + elemlist_out_inverse + "}\n" +
+  "  if ~multivectorp r;\n"
   when
+    ebeno_func = "subtract_" + str type1 + "_double";
+    ebeno_func_inverse = "subtract_double_" + str type1;
     a_elem = ["a" + str i | i = indices1];
     a_elemlist = foldl (\el s -> el + ("," + s)) (a_elem!0) (drop 1 a_elem);
-    a_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (a_elem!0 + "::double") (drop 1 a_elem);
     elemlist_out = foldl (\el s -> el + ("," + s)) (a_elem!0 + " - r") (drop 1 a_elem);
+    elemlist_out_inverse = foldl (\el s -> el + (",-" + s)) ("r - " + a_elem!0) (drop 1 a_elem);
   end;
-make_multivector_subtraction type1 type2 indices1 indices2::list type_out =
-  sprintf "r::double ::- %s {%s} |\n" (str type2, b_elemlist_double) +
-  sprintf "r::real ::- %s {%s} = %s {%s};\n" (str type2, b_elemlist, str type_out, elemlist_out)
+
+//-------------------------------------------------------------------------
+
+// A function that makes geometric product functions.
+
+make_multivector_geometric_product type1 type2 indices1::list indices2::list type_out =
+  throw "NOT IMPLEMENTED";
+
+make_multivector_geometric_product type1 type2 indices1::list indices2 type_out =
+  str type1 + " c1@{" + a_elemlist + "} ::* r |\n" +
+  "r ::* " + str type1 + " c1@{" + a_elemlist + "} =\n" +
+  "  if doublep r && dmatrixp c1 then\n" +
+  "    (" + ebeno_func + " (pointer c, pointer cp1, double r) $$ " + str type_out + " c\n" +
+  "     when\n" +
+  "       c::dmatrix = rowvector $ repeatn (#c1) 0.0;\n" +
+  "       cp1::dmatrix = pack c1;\n" +
+  "     end)\n" +
+  "  else\n" +
+  "    " + str type_out + " {" + elemlist_out + "}\n" +
+  "  if ~multivectorp r;\n"
   when
-    b_elem = ["b" + str i | i = indices2];
-    b_elemlist = foldl (\el s -> el + ("," + s)) (b_elem!0) (drop 1 b_elem);
-    b_elemlist_double = foldl (\el s -> el + ("," + s + "::double")) (b_elem!0 + "::double") (drop 1 b_elem);
-    elemlist_out = foldl (\el s -> el + (",-" + s)) ("r - " + b_elem!0) (drop 1 b_elem);
+    ebeno_func = "gp_" + str type1 + "_double";
+    a_elem = ["a" + str i | i = indices1];
+    a_elemlist = foldl (\el s -> el + ("," + s)) (a_elem!0) (drop 1 a_elem);
+    elemlist_out = foldl (\el s -> el + (",r*" + s)) ("r*" + a_elem!0) (drop 1 a_elem);
   end;
 
 //-------------------------------------------------------------------------
 
+// A function that makes outer product functions.
+
+make_multivector_outer_product type1 type2 indices1::list indices2::list type_out =
+  throw "NOT IMPLEMENTED";
+
+make_multivector_outer_product type1 type2 indices1::list indices2 type_out =
+  str type1 + " c1@{" + a_elemlist + "} ⋀ r |\n" +
+  "r ⋀ " + str type1 + " c1@{" + a_elemlist + "} =\n" +
+  "  if doublep r && dmatrixp c1 then 0.0 else 0;\n" +
+  "  if ~multivectorp r;\n";
+
+//-------------------------------------------------------------------------
+
 puts "using ebeno::geometric_algebra_h;";
 puts "using math;";
-puts "using system;";
 
 puts "";
 puts "namespace ebeno;";
 puts "nonfix mvec2 vec2 bvec2 rotor2;";
 puts "public bad_multivector_grade;";
 
+puts "// Wedge, and an ASCII equivalent. These represesent the outer";
+puts "// product.";
+puts "infix (*) ⋀ :^: ;";
+
+// puts "// Math dot, middle dot, ASCII equivalent. These represent the";
+// puts "// scalar product.";
+// puts "infix (*) ⋅ · :.: ;";
+
 puts "";
-puts "type mvec2 (mvec2 c::matrix) = dim c == (1,4);";
-puts "type vec2 (vec2 c::matrix) = dim c == (1,2);";
-puts "type bvec2 (bvec2 c::matrix) = dim c == (1,1);";
-puts "type rotor2 (rotor2 c::matrix) = dim c == (1,2);";
+puts "(:^:) = (⋀);"; //  :^: is the same as wedge.
+// puts "(·) = (⋅);";    //  Middle dot is the same as math dot.
+// puts "(:.:) = (⋅);";  //  :.: is the same as math dot.
+
+puts "";
+puts "type mvec2 (mvec2 {_,_,_,_});";
+puts "type vec2 (vec2 {_,_});";
+puts "type bvec2 (bvec2 {_,_});";
+puts "type rotor2 (rotor2 {_});";
+puts "type multivector (mvec2 {_,_,_,_});";
+puts "type multivector (vec2 {_,_});";
+puts "type multivector (bvec2 {_,_});";
+puts "type multivector (rotor2 {_});";
+
+puts "";
+puts "mvec2p = typep mvec2;";
+puts "vec2p = typep vec2;";
+puts "bvec2p = typep bvec2;";
+puts "rotor2p = typep rotor2;";
+puts "multivectorp = typep multivector;";
 
 //-------------------------------------------------------------------------
 
 puts "";
 
-puts $ make_multivector_addition mvec2 mvec2 [0,1,2,3] [0,1,2,3] mvec2;
-puts $ make_multivector_addition mvec2 vec2 [0,1,2,3] [1,2] mvec2;
-puts $ make_multivector_addition mvec2 bvec2 [0,1,2,3] [3] mvec2;
-puts $ make_multivector_addition mvec2 rotor2 [0,1,2,3] [0,3] mvec2;
-puts $ make_multivector_addition vec2 vec2 [0,1] [0,1] vec2;
-puts $ make_multivector_addition bvec2 bvec2 [0] [0] bvec2;
-puts $ make_multivector_addition bvec2 rotor2 [1] [0,1] rotor2;
-puts $ make_multivector_addition rotor2 rotor2 [0,1] [0,1] rotor2;
-
-puts $ make_multivector_addition mvec2 real [0,1,2,3] () mvec2;
-puts $ make_multivector_addition rotor2 real [0,1] () rotor2;
-
-puts $ make_multivector_addition real mvec2 () [0,1,2,3] mvec2;
-puts $ make_multivector_addition real rotor2 () [0,1] rotor2;
-
-//-------------------------------------------------------------------------
+flip do [make_multivector_addition, make_multivector_subtraction]
+(\func ->
+ puts $ func mvec2 mvec2 [0,1,2,3] [0,1,2,3] mvec2 $$
+ puts $ func mvec2 vec2 [0,1,2,3] [1,2] mvec2 $$
+ puts $ func mvec2 bvec2 [0,1,2,3] [3] mvec2 $$
+ puts $ func mvec2 rotor2 [0,1,2,3] [0,3] mvec2 $$
+ puts $ func vec2 vec2 [0,1] [0,1] vec2 $$
+ puts $ func bvec2 bvec2 [0] [0] bvec2 $$
+ puts $ func rotor2 rotor2 [0,1] [0,1] rotor2 $$
+ puts $ func rotor2 bvec2 [0,1] [1] rotor2 $$
+ puts $ func mvec2 () [0,1,2,3] () mvec2 $$
+ puts $ func rotor2 () [0,1] () rotor2);
 
 puts "";
 
-puts $ make_multivector_subtraction mvec2 mvec2 [0,1,2,3] [0,1,2,3] mvec2;
-puts $ make_multivector_subtraction mvec2 vec2 [0,1,2,3] [1,2] mvec2;
-puts $ make_multivector_subtraction mvec2 bvec2 [0,1,2,3] [3] mvec2;
-puts $ make_multivector_subtraction mvec2 rotor2 [0,1,2,3] [0,3] mvec2;
-
-puts $ make_multivector_subtraction vec2 mvec2 [1,2] [0,1,2,3] mvec2;
-puts $ make_multivector_subtraction vec2 vec2 [0,1] [0,1] vec2;
-
-puts $ make_multivector_subtraction bvec2 mvec2 [3] [0,1,2,3] mvec2;
-puts $ make_multivector_subtraction bvec2 bvec2 [0] [0] bvec2;
-puts $ make_multivector_subtraction bvec2 rotor2 [1] [0,1] rotor2;
-
-puts $ make_multivector_subtraction rotor2 mvec2 [0,3] [0,1,2,3] mvec2;
-puts $ make_multivector_subtraction rotor2 bvec2 [0,1] [1] rotor2;
-puts $ make_multivector_subtraction rotor2 rotor2 [0,1] [0,1] rotor2;
-
-puts $ make_multivector_subtraction mvec2 real [0,1,2,3] () mvec2;
-puts $ make_multivector_subtraction rotor2 real [0,1] () rotor2;
-
-puts $ make_multivector_subtraction real mvec2 () [0,1,2,3] mvec2;
-puts $ make_multivector_subtraction real rotor2 () [0,1] rotor2;
+puts $ make_multivector_geometric_product mvec2 () [0,1,2,3] () mvec2;
+puts $ make_multivector_geometric_product vec2 () [0,1] () vec2;
+puts $ make_multivector_geometric_product bvec2 () [0,1] () bvec2;
+puts $ make_multivector_geometric_product rotor2 () [0] () rotor2;
 
 //-------------------------------------------------------------------------