400% speedup by using maxima's optimize()

4 files changed, 54 insertions, 35 deletions

diff --git a/src/Makefile.lsq_funcs b/src/Makefile.lsq_funcs
index 8f92113..c33fe5f 100644
--- a/src/Makefile.lsq_funcs
+++ b/src/Makefile.lsq_funcs
@@ -1,7 +1,7 @@
 all: ProjectionRectilinear_funcs.cxx ProjectionCylindrical_funcs.cxx
 
-ProjectionRectilinear_funcs.cxx: lsq_rectilinear.mac
-	maxima -b lsq_rectilinear.mac | grep "^double" > ProjectionRectilinear_funcs.cxx
+ProjectionRectilinear_funcs.cxx: lsq_rectilinear.mac expr2c.mac
+	maxima -b lsq_rectilinear.mac | grep "^double" | sed s/%/_/g > ProjectionRectilinear_funcs.cxx
 
-ProjectionCylindrical_funcs.cxx: lsq_cylindrical.mac
-	maxima -b lsq_cylindrical.mac | grep "^double" > ProjectionCylindrical_funcs.cxx
+ProjectionCylindrical_funcs.cxx: lsq_cylindrical.mac expr2c.mac
+	maxima -b lsq_cylindrical.mac | grep "^double" | sed s/%/_/g > ProjectionCylindrical_funcs.cxx
diff --git a/src/expr2c.mac b/src/expr2c.mac
new file mode 100644
index 0000000..dba5987
--- /dev/null
+++ b/src/expr2c.mac
@@ -0,0 +1,20 @@
+printdecl(vars) := (
+	sprint("double", first(vars)),
+	for v in rest(vars) do sprint(",", v),
+	sprint(";"))$
+
+printdef(expr) := sprint(first(expr), "=", second(expr), ";")$
+
+block2c(expr) := 
+	if ?equal(op(expr), block) then (
+		printdecl(first(expr)),
+		for d in reverse(rest(reverse(rest(expr)))) do printdef(d),
+		sprint("return", string(last(expr)), ";}", "
+"))
+	else
+		sprint("return", string(expr), ";}", "
+")$
+
+expr2c(funcname, args, expr) := (
+	sprint("double", funcname, "(", args, ") {"),
+	block2c(subst(pow, "^", optimize(expr))))$
diff --git a/src/lsq_cylindrical.mac b/src/lsq_cylindrical.mac
index 541e4ad..0bf9499 100644
--- a/src/lsq_cylindrical.mac
+++ b/src/lsq_cylindrical.mac
@@ -13,22 +13,22 @@ y_scale : y_rot * scale$
  * Some mangling for C code generation
  */
 
+load("expr2c.mac")$
+
 x_expand : trigexpand(x_scale)$
 y_expand : trigexpand(y_scale)$
 
-args: "(double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick)"$
+args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick"$
 
-printfunc(name, expression) := sprint("double ProjectionCylindrical::", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
-")$
+expr2c("ProjectionCylindrical::mac_x", args, x_expand)$
+expr2c("ProjectionCylindrical::mac_y", args, y_expand)$
 
-printfunc("mac_x", x_expand)$
-printfunc("mac_y", y_expand)$
+expr2c("ProjectionCylindrical::mac_x_dc_view", args, diff(x_expand, c_view))$
+expr2c("ProjectionCylindrical::mac_x_dc_nick", args, diff(x_expand, c_nick))$
+expr2c("ProjectionCylindrical::mac_x_dc_tilt", args, diff(x_expand, c_tilt))$
+expr2c("ProjectionCylindrical::mac_x_dscale", args, diff(x_expand, scale))$
+expr2c("ProjectionCylindrical::mac_y_dc_view", args, diff(y_expand, c_view))$
+expr2c("ProjectionCylindrical::mac_y_dc_nick", args, diff(y_expand, c_nick))$
+expr2c("ProjectionCylindrical::mac_y_dc_tilt", args, diff(y_expand, c_tilt))$
+expr2c("ProjectionCylindrical::mac_y_dscale", args, diff(y_expand, scale))$
 
-printfunc("mac_x_dc_view", diff(x_expand, c_view))$
-printfunc("mac_x_dc_nick", diff(x_expand, c_nick))$
-printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt))$
-printfunc("mac_x_dscale", diff(x_expand, scale))$
-printfunc("mac_y_dc_view", diff(y_expand, c_view))$
-printfunc("mac_y_dc_nick", diff(y_expand, c_nick))$
-printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt))$
-printfunc("mac_y_dscale", diff(y_expand, scale))$
diff --git a/src/lsq_rectilinear.mac b/src/lsq_rectilinear.mac
index 69a9642..69980c0 100644
--- a/src/lsq_rectilinear.mac
+++ b/src/lsq_rectilinear.mac
@@ -15,26 +15,25 @@ y_dist : y_rot * (1 + dist_fact) * scale$
  * Some mangling for C code generation
  */
 
+load("expr2c.mac")$
+
 x_expand : trigexpand(x_dist)$
 y_expand : trigexpand(y_dist)$
 
-args: "(double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick)"$
-
-printfunc(name, expression) := sprint("double ProjectionRectilinear::", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
-")$
+args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick"$
 
-printfunc("mac_x", x_expand)$
-printfunc("mac_y", y_expand)$
+expr2c("ProjectionRectilinear::mac_x", args, x_expand)$
+expr2c("ProjectionRectilinear::mac_y", args, y_expand)$
 
-printfunc("mac_x_dc_view", diff(x_expand, c_view))$
-printfunc("mac_x_dc_nick", diff(x_expand, c_nick))$
-printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt))$
-printfunc("mac_x_dscale", diff(x_expand, scale))$
-printfunc("mac_x_dk0", diff(x_expand, k0))$
-printfunc("mac_x_dk1", diff(x_expand, k1))$
-printfunc("mac_y_dc_view", diff(y_expand, c_view))$
-printfunc("mac_y_dc_nick", diff(y_expand, c_nick))$
-printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt))$
-printfunc("mac_y_dscale", diff(y_expand, scale))$
-printfunc("mac_y_dk0", diff(y_expand, k0))$
-printfunc("mac_y_dk1", diff(y_expand, k1))$
+expr2c("ProjectionRectilinear::mac_x_dc_view", args, diff(x_expand, c_view))$
+expr2c("ProjectionRectilinear::mac_x_dc_nick", args, diff(x_expand, c_nick))$
+expr2c("ProjectionRectilinear::mac_x_dc_tilt", args, diff(x_expand, c_tilt))$
+expr2c("ProjectionRectilinear::mac_x_dscale", args, diff(x_expand, scale))$
+expr2c("ProjectionRectilinear::mac_x_dk0", args, diff(x_expand, k0))$
+expr2c("ProjectionRectilinear::mac_x_dk1", args, diff(x_expand, k1))$
+expr2c("ProjectionRectilinear::mac_y_dc_view", args, diff(y_expand, c_view))$
+expr2c("ProjectionRectilinear::mac_y_dc_nick", args, diff(y_expand, c_nick))$
+expr2c("ProjectionRectilinear::mac_y_dc_tilt", args, diff(y_expand, c_tilt))$
+expr2c("ProjectionRectilinear::mac_y_dscale", args, diff(y_expand, scale))$
+expr2c("ProjectionRectilinear::mac_y_dk0", args, diff(y_expand, k0))$
+expr2c("ProjectionRectilinear::mac_y_dk1", args, diff(y_expand, k1))$