[PATCH] Splines for parametric plot

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

This patch cleans up my previous patch (which added splines support for
plotting), and allows splines with parametric plotting,

It should be applied *after* my previous plotting patches.  Try it with

  ploth(x=0,6,[x*sin(x),x*cos(x)],128+256,9)
  ploth(x=0,6,[x*sin(x),x*cos(x)],128+256+1,9)

Enjoy,
Ilya

P.S.  If the lines in the highlvl.c patch are too long, all it does is
      an addition of ", 256 use cubic splines".

--- ./src/graph/plotport.c.with_splines	Sun Nov  8 21:26:24 1998
+++ ./src/graph/plotport.c	Sat Dec  5 03:16:52 1998
@@ -67,6 +67,11 @@ dsprintf9(double d, char *buf)
     return buf;				/* Should not happen? */
 }
 
+static const char *  const	quark = "quark"; /* Used as a special-case */
+static GEN			quark_gen;
+
+#define READ_EXPR(s)	((s)==quark? quark_gen : lisexpr(s))
+
 void
 plot(entree *ep, GEN a, GEN b, char *ch)
 {
@@ -93,7 +98,7 @@ plot(entree *ep, GEN a, GEN b, char *ch)
   av2=avma; limite=(av2+bot)>>1;
   for (i=1; i<=ISCR; i++)
   {
-    gaffect(lisexpr(ch),y[i]);
+    gaffect(READ_EXPR(ch),y[i]);
     if (gcmp(y[i],ysml)<0) ysml=y[i];
     if (gcmp(y[i],ybig)>0) ybig=y[i];
     x = addrr(x,dx);
@@ -826,7 +831,7 @@ single_recursion(dblPointList *pl,char *
   if (depth==RECUR_MAXDEPTH) return;
 
   yy=cgetr(3); xx=gmul2n(gadd(xleft,xright),-1);
-  ep->value = (void*) xx; gaffect(lisexpr(ch), yy);
+  ep->value = (void*) xx; gaffect(READ_EXPR(ch), yy);
 
   if (dy)
   {
@@ -855,7 +860,7 @@ param_recursion(dblPointList *pl,char *c
   if (depth==PARAMR_MAXDEPTH) return;
 
   xx=cgetr(3); yy=cgetr(3); tt=gmul2n(gadd(tleft,tright),-1);
-  ep->value = (void*)tt; p1=lisexpr(ch);
+  ep->value = (void*)tt; p1=READ_EXPR(ch);
   gaffect((GEN)p1[1],xx); gaffect((GEN)p1[2],yy);
 
   if (dx && dy)
@@ -907,7 +912,7 @@ rectplothin(entree *ep, GEN a, GEN b, ch
   dx=gdivgs(gsub(b,a),testpoints-1);
 
   x = cgetr(prec); gaffect(a,x); push_val(ep, x);
-  av2=avma; p1=lisexpr(ch); tx=typ(p1);
+  av2=avma; p1=READ_EXPR(ch); tx=typ(p1);
   if (!is_matvec_t(tx))
   {
     xsml = gtodouble(a); xbig=gtodouble(b);
@@ -957,14 +962,14 @@ rectplothin(entree *ep, GEN a, GEN b, ch
     {
       tleft=cgetr(prec); tright=cgetr(prec);
       av2=avma;
-      gaffect(a,tleft); ep->value = (void*)tleft; p1=lisexpr(ch);
+      gaffect(a,tleft); ep->value = (void*)tleft; p1=READ_EXPR(ch);
       gaffect((GEN)p1[1],xleft); gaffect((GEN)p1[2],yleft);
       for (i=0; i<testpoints-1; i++)
       {
 	if (i) {
 	  gaffect(tright,tleft); gaffect(xright,xleft); gaffect(yright,yleft);
 	 }
-	gaddz(tleft,dx,tright); ep->value = (void*)tright; p1=lisexpr(ch);
+	gaddz(tleft,dx,tright); ep->value = (void*)tright; p1=READ_EXPR(ch);
         gaffect((GEN)p1[1],xright); gaffect((GEN)p1[2],yright);
 
 	Appendx(&pl[0],&pl[0],rtodbl(xleft));
@@ -980,11 +985,11 @@ rectplothin(entree *ep, GEN a, GEN b, ch
     {
       av2=avma;
       gaffect(a,xleft); ep->value = (void*) xleft;
-      gaffect(lisexpr(ch),yleft);
+      gaffect(READ_EXPR(ch),yleft);
       for (i=0; i<testpoints-1; i++)
       {
         gaddz(xleft,dx,xright); ep->value = (void*) xright;
-	gaffect(lisexpr(ch),yright);
+	gaffect(READ_EXPR(ch),yright);
 
 	Appendx(&pl[0],&pl[0],rtodbl(xleft));
 	Appendy(&pl[0],&pl[1],rtodbl(yleft));
@@ -1002,7 +1007,7 @@ rectplothin(entree *ep, GEN a, GEN b, ch
     if (single_c)
       for (i=0; i<testpoints; i++)
       {
-	p1 = lisexpr(ch);
+	p1 = READ_EXPR(ch);
 	pl[0].d[i]=gtodouble(x);
 	Appendy(&pl[0],&pl[1],gtodouble(p1));
 	gaddz(x,dx,x); avma=av2;
@@ -1014,7 +1019,7 @@ rectplothin(entree *ep, GEN a, GEN b, ch
 
       for (i=0; i<testpoints; i++)
       {
-	p1 = lisexpr(ch);
+	p1 = READ_EXPR(ch);
 	for (j=0; j<nl; j=k)
 	{
 	  k=j+1; z=gtodouble((GEN)p1[k]);
@@ -1029,7 +1034,7 @@ rectplothin(entree *ep, GEN a, GEN b, ch
     else /* plothmult */
       for (i=0; i<testpoints; i++)
       {
-	p1 = lisexpr(ch);
+	p1 = READ_EXPR(ch);
 	pl[0].d[i]=gtodouble(x);
 	for (j=1; j<nl; j++) { Appendy(&pl[0],&pl[j],gtodouble((GEN)p1[j])); }
 	gaddz(x,dx,x); avma=av2;
@@ -1039,65 +1044,59 @@ rectplothin(entree *ep, GEN a, GEN b, ch
   return pl;
 }
 
-/* Uses highlevel functions to implement splines in low-level functions.
-
+/* Uses highlevel plotting functions to implement splines as 
+   a low-level plotting function.
    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;
+  long i, j, oldavma = avma;
   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__";
+  GEN ttas = cgetg(4 + 1, t_VEC);
+  GEN param = cgetg(2 + 1, t_VEC);
+  GEN tas;
+  entree *var0 = varentries[ordvar[0]];
   
-  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]);      
   }
+  if (flag & PLOT_PARAMETRIC) {
+      for (j = 1; j <= 4; j++)
+	  ttas[j] = (long)stoi(j);
+      tas = ttas;
+      quark_gen = param;
+  } else
+      tas = xas;
   for (i = 0; i <= lx - 4; i++) {
       long oavma = avma;
-      GEN poly;
 
-      for (j = 1; j <= 4; j++) {
+      for (j = 1; j <= 4; j++) {	/* Need to copy to have
+					   correct headers of xas and yas. */
 	  xas[j] = xa[i + j];
-	  yas[j] = ya[i + j];      
+	  yas[j] = ya[i + j];
       }
-      poly = polint(xas, yas, polx[0], NULL);
-      var1->value = (void*)poly;	/* Temporary make it into expr */
+      if (flag & PLOT_PARAMETRIC) {
+	  param[1] = (long)polint(ttas, xas, polx[0], NULL);
+	  param[2] = (long)polint(ttas, yas, polx[0], NULL);
+      } else
+	  quark_gen = polint(xas, yas, polx[0], NULL);
+      
       rectploth( ne, var0, 
-		 (GEN)(i==0 ? xa[1] : xa[i+2]),
-		 (GEN)(i==lx-4 ? xa[lx] : xa[i+3]),
-		 name, 4,		/* XXXX precision */
+		 (GEN)(i==0 ? tas[1] : tas[2]),
+		 (GEN)(i==lx-4 ? tas[4] : tas[3]),
+		 (char*)quark,		/* Remove const */
+		 4,			/* XXXX precision */
 		 (PLOT_RECURSIVE | PLOT_NO_RESCALE | PLOT_NO_FRAME
-		  | PLOT_NO_AXE_Y | PLOT_NO_AXE_X),
+		  | PLOT_NO_AXE_Y | PLOT_NO_AXE_X)
+		 | (flag & PLOT_PARAMETRIC),
 		 2);			/* Start with 3 points */
       avma = oavma;
   }
-  var1->value = old1;			/* Restore the old value */
-
-#if 0					/* XXXX Does not work??? */
-  if (dokill) {
-      kill0(var1);
-      manage_var(1,NULL);
-  }
-#endif  
-
   avma = oldavma;
 }
 
@@ -1216,11 +1215,17 @@ rectplothrawin(long stringrect, long dra
 	rectpoints0(drawrect,x.d,y.d,nbpoints);
     }
     if ((flags & PLOT_POINTS_LINES) || !(flags & PLOT_POINTS)) {
-	rectlinetype(drawrect, rectline_itype + ltype); 	/* Graphs. */
-	if ((flags & PLOT_PARAMETRIC) || !(flags & PLOT_SPLINES))
-	    rectlines0(drawrect,x.d,y.d,nbpoints,0);
-	else
-	    rectsplines(drawrect,x.d,y.d,nbpoints,0);
+	if (flags & PLOT_SPLINES) {
+	    /* rectsplines will call us back with ltype == 0 */
+	    int old = rectline_itype;
+
+	    rectline_itype = rectline_itype + ltype;
+	    rectsplines(drawrect,x.d,y.d,nbpoints,flags);	    
+	    rectline_itype = old;
+	} else {
+	    rectlinetype(drawrect, rectline_itype + ltype); 	/* Graphs. */
+	    rectlines0(drawrect,x.d,y.d,nbpoints,0);	    
+	}
     }
     ltype++;				/* Graphs. */
   }
--- ./src/language/highlvl.c~	Sun Nov  8 21:26:30 1998
+++ ./src/language/highlvl.c	Sat Dec  5 03:22:56 1998
@@ -177,7 +177,7 @@ char *helpmessages_highlevel[]={
   "plotcursor(w): current position of cursor in rectwindow w",
   "plotdraw(list): draw vector of rectwindows list at indicated x,y positions; list is a vector w1,x1,y1,w2,x2,y2,etc. . ",
   "plotfile(filename): set the output file for plotting output. \"-\" redirects to the same place as PARI output",
-  "ploth(X=a,b,expr,{flags=0},{n=0}): plot of expression expr, X goes from a to b in high resolution. Both flags and n are optional. Binary digits of flags mean : 1 parametric plot, 2 recursive plot, 8 omit x-axis, 16 omit y-axis, 32 omit frame, 64 do not join points, 128 plot both lines and points. n specifies number of reference points on the graph (0=use default value). Returns a vector for the bounding box",
+  "ploth(X=a,b,expr,{flags=0},{n=0}): plot of expression expr, X goes from a to b in high resolution. Both flags and n are optional. Binary digits of flags mean : 1 parametric plot, 2 recursive plot, 8 omit x-axis, 16 omit y-axis, 32 omit frame, 64 do not join points, 128 plot both lines and points, 256 use cubic splines. n specifies number of reference points on the graph (0=use default value). Returns a vector for the bounding box",
   "plothraw(listx,listy,{flag=0}): plot in high resolution points whose x (resp. y) coordinates are in listx (resp. listy). If flag is non zero, join points",
   "plothsizes(): returns array of 6 elements: terminal width and height, sizes for ticks in horizontal and vertical directions (in pixels), width and height of characters",
   "plotinit(w,x,y): initialize rectwindow w to size x,y",