2.0.12: half-baked support of splines in plotting

[Ilya Zakharevich <ilya@math.ohio-state.edu>]

This patch:

  a) when gaffect(non_constant, constant_type), calls geval() on non_constant
	if this will not lead to an obvious loop;
  b) Implements (in some situations) new plotting flag 256 to use splines.
	Try
	    ploth(x=0,3,x^4,128+256,5)

(It should be applied after my previous patch, to have PLOT_SPLINES defined.)

However, it contains a bug.  Look for #ifdef I_KNOW_HOW_TO_DO_THIS

I want to have a named variable, then delete it without a trace.
I know how to do creating part (up to some extend - I think it will fail if
the var#1 was previously deleled), but the part in #ifdef fails: apparently
delete_var() is a NOP.

Can somebody fix this?  I need a named variable to pass its name to
highlevel plotting functions (due to "a" above, now
    y = x^17
    ploth(x=2,5,y)
works).

Ilya

P.S. A year ago I had some Perl code to make Math::Pari to
     autosubdivide so that splines will give necessary precision
     (similar to current recursive linear approximation).  Now, when
     PARI proper is not that much behind Math::Pari module in plotting
     support, maybe I could reimplement it in PARI.

--- ./src/graph/plotport.c~	Sun Nov  8 04:22:06 1998
+++ ./src/graph/plotport.c	Sun Nov  8 17:56:08 1998
@@ -1037,6 +1037,68 @@ rectplothin(entree *ep, GEN a, GEN b, ch
   return pl;
 }
 
+/* Uses highlevel functions to implement splines in low-level functions.
+
+   Most probably one does not need to make varn==0 into pure variable,
+   since plotting functions should take care of this. */
+static void
+rectsplines(long ne, double *x, double *y, long lx, long flag)
+{
+  long i, j, oldavma = avma, dokill = 0;
+  GEN xa = cgetg(lx + 1, t_VEC), ya = cgetg(lx + 1, t_VEC);
+  GEN xas = cgetg(4 + 1, t_VEC), yas = cgetg(4 + 1, t_VEC);
+  void *old1;
+  entree *var0 = varentries[ordvar[0]], *var1 = varentries[ordvar[1]];
+  char *name = 0;
+  static const char * const fakename = "splines__";
+  
+  if (!var1) {
+      GEN res = lisexpr((char*)fakename);
+      dokill = 1;
+      var1 = varentries[ordvar[1]];
+      if (!var1 || !(name = var1->name) && strcmp(name, fakename))
+	  err(talker, "panic: no new var in splines");
+  }
+  if (lx < 4)
+      err(talker, "Too few points (%ld) for spline plot", lx);
+  old1 = var1->value;
+  if (!name)
+      name = var1->name;
+  for (i = 1; i <= lx; i++) {
+      xa[i] = (long) dbltor(x[i-1]);
+      ya[i] = (long) dbltor(y[i-1]);      
+  }
+  for (i = 0; i <= lx - 4; i++) {
+      long oavma = avma;
+      GEN poly;
+
+      for (j = 1; j <= 4; j++) {
+	  xas[j] = xa[i + j];
+	  yas[j] = ya[i + j];      
+      }
+      poly = polint(xas, yas, polx[0], NULL);
+      var1->value = (void*)poly;	/* Temporary make it into expr */
+      rectploth( ne, var0, 
+		 (GEN)(i==0 ? xa[1] : xa[i+2]),
+		 (GEN)(i==lx-4 ? xa[lx] : xa[i+3]),
+		 name, 4,		/* XXXX precision */
+		 (PLOT_RECURSIVE | PLOT_NO_RESCALE | PLOT_NO_FRAME
+		  | PLOT_NO_AXE_Y | PLOT_NO_AXE_X),
+		 2);			/* Start with 3 points */
+      avma = oavma;
+  }
+  var1->value = old1;			/* Restore the old value */
+
+#ifdef I_KNOW_HOW_TO_DO_THIS			/* XXXX Does not work??? */
+  if (dokill) {
+      kill0(var1);
+      manage_var(1,NULL);
+  }
+#endif  
+
+  avma = oldavma;
+}
+
 /*
  * Plot a dblPointList. Complete with axes, bounding box, etc.
  * We use two drawing rectangles: one for strings, another
@@ -1153,7 +1215,10 @@ rectplothrawin(long stringrect, long dra
     }
     if ((flags & PLOT_POINTS_LINES) || !(flags & PLOT_POINTS)) {
 	rectlinetype(drawrect, rectline_itype + ltype); 	/* Graphs. */
-	rectlines0(drawrect,x.d,y.d,nbpoints,0);
+	if ((flags & PLOT_PARAMETRIC) || !(flags & PLOT_SPLINES))
+	    rectlines0(drawrect,x.d,y.d,nbpoints,0);
+	else
+	    rectsplines(drawrect,x.d,y.d,nbpoints,0);
     }
     ltype++;				/* Graphs. */
   }
--- ./src/basemath/gen2.c~	Fri Nov  6 10:08:14 1998
+++ ./src/basemath/gen2.c	Sun Nov  8 16:47:22 1998
@@ -1453,7 +1453,33 @@ gaffect(GEN x, GEN y)
     return;
   }
 
-  if (is_const_t(ty)) err(assignerf,tx,ty);
+  if (is_const_t(ty)) {
+      /* initial_value macro is not defined now... */
+#define initial_value(ep) ((ep)+1)
+      if (tx == t_POL) {
+	  entree *var = varentries[ordvar[varn(x)]];
+	  /* Is a bound expression, thus should not loop: */
+	  if (var && var->value != (void*)initial_value(var)) {
+	      gaffect(geval(x),y);
+	      return;
+	  }	  
+      } else if (is_rfrac_t(tx)) {
+	  GEN num = (GEN)x[1], den = (GEN)x[2];
+	  entree *varnum = varentries[ordvar[varn(num)]];
+	  entree *varden = varentries[ordvar[varn(den)]];
+
+	  /* Are bound expressions, thus should not loop: */
+	  if (varnum && varnum->value != (void*)initial_value(varnum)
+	      && varden && varden->value !=
+	      (void*)initial_value(varden)) {
+	      gaffect(geval(x),y);
+	      return;
+	  }
+      }
+#undef initial_value
+      err(assignerf,tx,ty);
+  }
+  
   switch(tx)
   {
     case t_POL: