curveopt: symbolic: wip, before go back to conditional in C

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>

diff --git a/findcurve.c b/findcurve.c
index f6df6e6fd28f15ba856b2ee5d0f598f393464d40..0fa5c10a7691f085d7d1d53b6e2302a94a12e2a2 100644 (file)
--- a/findcurve.c
+++ b/findcurve.c
@@ -26,6 +26,15 @@ static double target[N];
 
 static double cb_Efunc(void *xp) {
   const double *x = xp;
+
+  e = 0;
+  for (i=0; i<N-3; i++) {
+    // A is point i, B i+1, C i+2, etc.
+    if (i == 0) {
+    } else if (i==N-2) {
+    } else {
+      E_CALCULATE_ED;
+
   double F[N], e;
   int i;
   X_EXTRACT;
@@ -40,6 +49,11 @@ static double cb_Efunc(void *xp) {
   return e;
 }
 
+static int cb_fdf(const gsl_vector *x, void *params,
+                 gsl_vector *f, gsl_matrix *J) {
+  
+}
+
 static void cb_step(const gsl_rng *rng, void *xp, double step_size) {
   double *x = xp;
   int i;

diff --git a/symbolic.py b/symbolic.py
index 38ef5e6ffeb463337a83950a21f7257f187b1925..bda742ace9fd496d8f29deb2a3a8ac07adaf1205 100644 (file)
--- a/symbolic.py
+++ b/symbolic.py
@@ -34,12 +34,32 @@ def vector_symbols(vnames):
     out.append(v)
   return out
 
-A, B, C, D = vector_symbols('A B C D')
+abcd = vector_symbols('A B C D')
+p = vector_symbols('p')
+
+E, F = vector_symbols('E F')
+En, Fn = vector_symbols('En Fn')
+
+def vector_component(v, ix):
+  return v.components[N.base_vectors()[ix]]
+
+# x array in numerical algorithm has:
+#    N x 3    coordinates of points 1..N-2
+#    1        length parameter |EA| for point 0
+#    1        length parameter |DE| for point N-1
+
+def point(p):
+  return Piecewise((p == 0, E + En * 
 
 def calculate():
   global calculated
   if calculated: return
 
+  for i in range(0..len(abcd)):
+    abcd[i] = 
+
+Piecewise(( P 
+
   Q = B + 0.5 * (B - A)
   R = C + 0.5 * (C - D)
   QR = R - Q
@@ -49,6 +69,17 @@ def calculate():
   dbg('cost_ABCD')
   dprint(A)
 
+  cost_ABCD = subst_vect(
+  
+
+  A_end = E + 
+
+  cost_EBCD = subst_vect(cost_ABCD, A, A_end)
+
+  # diff_A_i = cost_ABCD.diff(vector_component(A, 0))
+  # global diff_A_i
+  # dbg('diff_A_i')
+
   calculated = True
 
 def ourccode(*a, **kw):
@@ -106,8 +137,20 @@ def gen_diff(current, smalls):
 def gen_misc():
   cprintraw('// AUTOGENERATED - DO NOT EDIT\n')
 
-def gen_x_extract():
-  pass
+def gen_abcd_ijk_macros():
+  abcd = 'A B C D'.split(' ')
+  ijk = 'i j k'.split(' ')
+  for ai in range(0, len(abcd)):
+    for ii in range(0, len(ijk)):
+      cprintraw(('#define %s_%s' % (abcd[ai], ijk[ii]) +
+                 ' X[ ((P) - 1 + %d) * 3 + %d ]' % (ai, ii)))
+
+  cprint('#define E_CALCULATE_MID')
+  cassign(cost_ABCD,'d','dtmp_mid')
+
+def gen_e_calculate():
+  cprint('#define E_CALCULATE_MID')
+  cassign(cost_ABCD,'d','dtmp_mid')
 
 def gen_f_populate():
   cprint('#define F_POPULATE')
@@ -121,9 +164,8 @@ def gen_j_populate():
 
 def gen_C():
   gen_misc()
-  gen_x_extract()
-  gen_f_populate()
-  gen_j_populate()
+  gen_abcd_ijk_macros()
+  gen_e_calculate()
 
 def get_python():
   # https://github.com/sympy/sympy/issues/13642