curveopt: symbolic: wip, before abandon provision for fdf

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

diff --git a/findcurve.c b/findcurve.c
index 0fa5c10a7691f085d7d1d53b6e2302a94a12e2a2..eec17d4da6ac715b5606819310c2ae220cb31534 100644 (file)
--- a/findcurve.c
+++ b/findcurve.c
@@ -25,15 +25,16 @@ static inline double sinc(double x) {
 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) {
+  const double *X = xp;
+
+  double e = 0;
+  for (P=0; P<N-3; P++) {
+    double d;
+    // A is point #p, B #p+1, C #p+2, etc.
+    if (P == 0) {
+    } else if (P==N-2) {
     } else {
-      E_CALCULATE_ED;
+      E_CALCULATE_MID;
 
   double F[N], e;
   int i;

diff --git a/symbolic.py b/symbolic.py
index bda742ace9fd496d8f29deb2a3a8ac07adaf1205..f12a96d18d9496645394195604ef286abb1d7f61 100644 (file)
--- a/symbolic.py
+++ b/symbolic.py
@@ -34,32 +34,31 @@ def vector_symbols(vnames):
     out.append(v)
   return out
 
-abcd = vector_symbols('A B C D')
+A, B, C, D = vector_symbols('A B C D')
 p = vector_symbols('p')
 
-E, F = vector_symbols('E F')
-En, Fn = vector_symbols('En Fn')
+E, H = vector_symbols('E H')
+En, Hn = vector_symbols('En Hn')
+
+EFl, GHl = symbols('EFl GHl')
 
 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
+#    N x 3    coordinates of points 0..N-3
+#    1        EFl = length parameter |EF| for point 1
+#    1        GHl = length parameter |GH| for point N-2
+
+# fixed array in numerical algorithm has:
+#    4 x 3    E, En, H, Hn
 
-def point(p):
-  return Piecewise((p == 0, E + En * 
+#def subst_vect():
 
 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
@@ -69,12 +68,13 @@ Piecewise(( P
   dbg('cost_ABCD')
   dprint(A)
 
-  cost_ABCD = subst_vect(
-  
-
-  A_end = E + 
+  F = E + En * EFl
+  G = H + Hn * GHl
 
-  cost_EBCD = subst_vect(cost_ABCD, A, A_end)
+  #cost_EFCD = subst_vect(cost_ABCD, (A,E), (B,F))
+  #cost_FBCD = subst_vect(cost_ABCD,        (A,F))
+  #cost_ABCG = subst_vect(cost_ABCD,        (D,G))
+  #cost_ABGH = subst_vect(cost_ABCD, (C,G), (D,H))
 
   # diff_A_i = cost_ABCD.diff(vector_component(A, 0))
   # global diff_A_i
@@ -137,19 +137,29 @@ def gen_diff(current, smalls):
 def gen_misc():
   cprintraw('// AUTOGENERATED - DO NOT EDIT\n')
 
+def gen_point_index_macro(macro_basename, c_array_name, base_index):
+  cprintraw('#define %s (&%s[%s])'
+            % (macro_basename, c_array_name, base_index))
+  ijk = 'i j k'.split(' ')
+  for ii in range(0, len(ijk)):
+    cprintraw('#define %s_%s (%s[%d])'
+              % (macro_basename, ijk[ii], macro_basename, ii))
+
 def gen_abcd_ijk_macros():
   abcd = 'A B C D'.split(' ')
-  ijk = 'i j k'.split(' ')
+  eh = 'E En H Hn'.split(' ')
+  for ehi in range(0, len(eh)):
+    gen_point_index_macro(eh[ehi], 'FIXED', ehi * 3)
+
   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)))
+    gen_point_index_macro(abcd[ai], 'X', '(P%+d) * 3' % (ai - 2))
 
   cprint('#define E_CALCULATE_MID')
   cassign(cost_ABCD,'d','dtmp_mid')
 
 def gen_e_calculate():
   cprint('#define E_CALCULATE_MID')
+  cprint('if (P==0) '
   cassign(cost_ABCD,'d','dtmp_mid')
 
 def gen_f_populate():