[BACK]Return to plotport.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari / src / graph

Annotation of OpenXM_contrib/pari/src/graph/plotport.c, Revision 1.1.1.1

1.1       maekawa     1: /*******************************************************************/
                      2: /*                                                                 */
                      3: /*                         PLOT ROUTINES                           */
                      4: /*                                                                 */
                      5: /*******************************************************************/
                      6: /* $Id: plotport.c,v 1.1.1.1 1999/09/16 13:47:44 karim Exp $ */
                      7: #include "pari.h"
                      8: #include "rect.h"
                      9:
                     10: void push_val(entree *ep, GEN a);
                     11: void pop_val(entree *ep);
                     12: void postdraw0(long *w, long *x, long *y, long lw);
                     13: void rectdraw0(long *w, long *x, long *y, long lw, long do_free);
                     14: static void PARI_get_psplot();
                     15:
                     16: static long current_color[NUMRECT];
                     17: PariRect **rectgraph = NULL;
                     18: PARI_plot pari_plot, pari_psplot;
                     19: long  rectpoint_itype = 0;
                     20: long  rectline_itype  = 0;
                     21:
                     22: #define STRINGRECT (NUMRECT-2)
                     23: #define DRAWRECT (NUMRECT-1)
                     24:
                     25: #define PLOTH_NUMPOINTS 1000
                     26: #define PARAM_NUMPOINTS 1500
                     27: #define RECUR_NUMPOINTS 8
                     28:
                     29: #define RECUR_MAXDEPTH 10
                     30: #define RECUR_PREC 0.001
                     31: #define PARAMR_MAXDEPTH 10
                     32:
                     33: /********************************************************************/
                     34: /**                                                                **/
                     35: /**                         LOW-RES PLOT                           **/
                     36: /**                                                                **/
                     37: /********************************************************************/
                     38: #define ISCR 64
                     39: #define JSCR 22
                     40: #define BLANK ' '
                     41: #define ZERO1 ','
                     42: #define ZERO2 '-'
                     43: #define ZERO3 '`'
                     44: #define PICTZERO(j) ((j) % 3 ? ((j) % 3 == 2 ? ZERO3 : ZERO2) : ZERO1)
                     45: #define YY '|'
                     46: #define XX_UPPER '\''
                     47: #define XX_LOWER '.'
                     48: #define FF1 '_'
                     49: #define FF2 'x'
                     50: #define FF3 '"'
                     51: #define PICT(j) ((j) % 3 ? ((j) % 3 == 2 ? FF3 : FF2) : FF1)
                     52:
                     53: static char *
                     54: dsprintf9(double d, char *buf)
                     55: {
                     56:   int i = 10;
                     57:
                     58:   while (--i >= 0) {
                     59:     sprintf(buf, "%9.*g", i, d);
                     60:     if (strlen(buf) <= 9) return buf;
                     61:   }
                     62:   return buf; /* Should not happen? */
                     63: }
                     64:
                     65: #define QUARK  ((char*)NULL) /* Used as a special-case */
                     66: static GEN quark_gen;
                     67:
                     68: #define READ_EXPR(s)   ((s)==QUARK? quark_gen : lisexpr(s))
                     69:
                     70: void
                     71: plot(entree *ep, GEN a, GEN b, char *ch, long prec)
                     72: {
                     73:   long av = avma, av2,limite,jz,j,i,sig;
                     74:   int jnew, jpre = 0; /* for lint */
                     75:   GEN p1,p2,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
                     76:   unsigned char scr[ISCR+1][JSCR+1];
                     77:   char buf[80], z;
                     78:
                     79:   sig=gcmp(b,a); if (!sig) return;
                     80:   if (sig<0) { x=a; a=b; b=x; }
                     81:   x = cgetr(prec); gaffect(a,x); push_val(ep, x);
                     82:   for (i=1; i<=ISCR; i++) y[i]=cgetr(3);
                     83:   p1 = gdivgs(gsub(b,a), ISCR-1);
                     84:   dx = cgetr(prec); gaffect(p1, dx);
                     85:   ysml=gzero; ybig=gzero;
                     86:   for (j=1; j<=JSCR; j++) scr[1][j]=scr[ISCR][j]=YY;
                     87:   for (i=2; i<ISCR; i++)
                     88:   {
                     89:     scr[i][1] = XX_LOWER;
                     90:     scr[i][JSCR]=XX_UPPER;
                     91:     for (j=2; j<JSCR; j++) scr[i][j]=BLANK;
                     92:   }
                     93:   av2=avma; limite=stack_lim(av2,1);
                     94:   for (i=1; i<=ISCR; i++)
                     95:   {
                     96:     gaffect(READ_EXPR(ch),y[i]);
                     97:     if (gcmp(y[i],ysml)<0) ysml=y[i];
                     98:     if (gcmp(y[i],ybig)>0) ybig=y[i];
                     99:     x = addrr(x,dx);
                    100:     if (low_stack(limite, stack_lim(av2,1)))
                    101:     {
                    102:       long tetpil=avma;
                    103:       if (DEBUGMEM>1) err(warnmem,"plot");
                    104:       x = gerepile(av2,tetpil,rcopy(x));
                    105:     }
                    106:     ep->value = (void*)x;
                    107:   }
                    108:   avma=av2; diff=gsub(ybig,ysml);
                    109:   if (gcmp0(diff)) { ybig=gaddsg(1,ybig); diff=gun; }
                    110:   dyj = gdivsg((JSCR-1)*3+2,diff);
                    111:   jz = 3-gtolong(gmul(ysml,dyj));
                    112:   av2=avma; z = PICTZERO(jz); jz = jz/3;
                    113:   for (i=1; i<=ISCR; i++)
                    114:   {
                    115:     scr[i][jz]=z; j=3+gtolong(gmul(gsub(y[i],ysml),dyj));
                    116:     jnew = j/3;
                    117:     if (i > 1)
                    118:     {
                    119:       int i_up, i_lo,  mid = (jpre+jnew)/2, up, lo;
                    120:       int domid = 0;
                    121:
                    122:       /* If the gap is 1, leave it as it is. */
                    123:       if (jpre < jnew - 2) {
                    124:         up = jnew - 1; i_up = i;
                    125:         lo = jpre + 1; i_lo = i - 1;
                    126:         domid = 1;
                    127:       } else if (jnew < jpre - 2) {
                    128:         up = jpre - 1; i_up = i - 1;
                    129:         lo = jnew + 1; i_lo = i;
                    130:         domid = 1;
                    131:       }
                    132:       if (domid)
                    133:       {
                    134:         while (lo <= mid) scr[i_lo][lo++] = ':';
                    135:         while (up > mid)  scr[i_up][up--] = ':';
                    136:       }
                    137:     }
                    138:     scr[i][jnew] = PICT(j); avma=av2;
                    139:     jpre = jnew;
                    140:   }
                    141:   p1=cgetr(3); gaffect(ybig,p1); pariputc('\n');
                    142:   pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
                    143:   for (i=1; i<=ISCR; i++) pariputc(scr[i][JSCR]);
                    144:   pariputc('\n');
                    145:   for (j=(JSCR-1); j>=2; j--)
                    146:   {
                    147:     pariputs("          ");
                    148:     for (i=1; i<=ISCR; i++) pariputc(scr[i][j]);
                    149:     pariputc('\n');
                    150:   }
                    151:   p1=cgetr(3); gaffect(ysml,p1);
                    152:   pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
                    153:   for (i=1; i<=ISCR; i++)  pariputc(scr[i][1]);
                    154:   pariputc('\n');
                    155:   p1=cgetr(3); gaffect(a,p1);
                    156:   p2=cgetr(3); gaffect(b,p2);
                    157:   pariputsf("%10s%-9.7g%*.7g\n"," ",rtodbl(p1),ISCR-9,rtodbl(p2));
                    158:   pop_val(ep); avma=av;
                    159: }
                    160:
                    161: /********************************************************************/
                    162: /**                                                                **/
                    163: /**                      RECTPLOT FUNCTIONS                        **/
                    164: /**                                                                **/
                    165: /********************************************************************/
                    166: void
                    167: init_graph()
                    168: {
                    169:   int n;
                    170:
                    171:   rectgraph = (PariRect**) gpmalloc(sizeof(PariRect*)*NUMRECT);
                    172:   for (n=0; n<NUMRECT; n++)
                    173:   {
                    174:     PariRect *e = (PariRect*) gpmalloc(sizeof(PariRect));
                    175:
                    176:     e->head = e->tail = NULL;
                    177:     e->sizex = e->sizey = 0;
                    178:     current_color[n] = DEFAULT_COLOR;
                    179:     rectgraph[n] = e;
                    180:   }
                    181: }
                    182:
                    183: void
                    184: free_graph()
                    185: {
                    186:   int i;
                    187:
                    188:   for (i=0; i<NUMRECT; i++)
                    189:   {
                    190:     PariRect *e=rectgraph[i];
                    191:
                    192:     if (RHead(e)) killrect(i);
                    193:     free((void *)e);
                    194:   }
                    195:   free((void *)rectgraph);
                    196: }
                    197:
                    198: static PariRect *
                    199: check_rect(long ne)
                    200: {
                    201:   if (!GOODRECT(ne))
                    202:     err(talker,"not an rplot vector type in graphic function");
                    203:   return rectgraph[ne];
                    204: }
                    205:
                    206: static PariRect *
                    207: check_rect_init(long ne)
                    208: {
                    209:   PariRect *e = check_rect(ne);
                    210:   if (!RHead(e))
                    211:     err(talker,"you must initialize the rectwindow first");
                    212:   return e;
                    213: }
                    214:
                    215: void
                    216: initrect(long ne, long x, long y)
                    217: {
                    218:   PariRect *e;
                    219:   RectObj *z;
                    220:
                    221:   if (x<=1 || y<=1) err(talker,"incorrect dimensions in initrect");
                    222:   e = check_rect(ne); if (RHead(e)) killrect(ne);
                    223:
                    224:   z = (RectObj*) gpmalloc(sizeof(RectObj));
                    225:   RoNext(z) = NULL;
                    226:   RoType(z) = ROt_NULL;
                    227:   RHead(e)=RTail(e)=z;
                    228:   RXsize(e)=x; RYsize(e)=y;
                    229:   RXcursor(e)=0; RYcursor(e)=0;
                    230:   RXscale(e)=1; RXshift(e)=0;
                    231:   RYscale(e)=1; RYshift(e)=0;
                    232: }
                    233:
                    234: GEN
                    235: rectcursor(long ne)
                    236: {
                    237:   PariRect *e = check_rect_init(ne);
                    238:   GEN z=cgetg(3,t_VEC);
                    239:
                    240:   z[1] = lstoi((long)RXcursor(e)); z[2] = lstoi((long)RYcursor(e));
                    241:   return z;
                    242: }
                    243:
                    244: static void
                    245: rectscale0(long ne, double x1, double x2, double y1, double y2)
                    246: {
                    247:   PariRect *e = check_rect_init(ne);
                    248:   double p2,p3;
                    249:
                    250:   p2 = RXshift(e) + RXscale(e) * RXcursor(e);
                    251:   p3 = RYshift(e) + RYscale(e) * RYcursor(e);
                    252:   RXscale(e) = RXsize(e)/(x2-x1); RXshift(e) = -x1*RXscale(e);
                    253:   RYscale(e) = RYsize(e)/(y1-y2); RYshift(e) = -y2*RYscale(e);
                    254:   RXcursor(e) = (p2 - RXshift(e)) / RXscale(e);
                    255:   RYcursor(e) = (p3 - RYshift(e)) / RYscale(e);
                    256: }
                    257:
                    258: void
                    259: rectscale(long ne, GEN x1, GEN x2, GEN y1, GEN y2)
                    260: {
                    261:   rectscale0(ne, gtodouble(x1), gtodouble(x2), gtodouble(y1), gtodouble(y2));
                    262: }
                    263:
                    264: /*
                    265:  * #ifdef LONG_IS_64BIT
                    266:  * #  define MAXDTOL 9223372036854775807.0
                    267:  * #else
                    268:  * #  define MAXDTOL 2147483647.0
                    269:  * #endif
                    270:  *
                    271:  *
                    272:  * long
                    273:  * DTOL(double t)
                    274:  * {
                    275:  *   if (fabs(t) > MAXDTOL) err(affer2);
                    276:  *   return (long)t;
                    277:  * }
                    278:  */
                    279:
                    280: #define DTOL(t) ((long)(t))
                    281:
                    282: static void
                    283: rectmove0(long ne, double x, double y, long relative)
                    284: {
                    285:   PariRect *e = check_rect_init(ne);
                    286:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj1P));
                    287:
                    288:   if (relative)
                    289:    { RXcursor(e) += x; RYcursor(e) += y; }
                    290:   else
                    291:    { RXcursor(e) = x; RYcursor(e) = y; }
                    292:   RoNext(z) = 0; RoType(z) = ROt_MV;
                    293:   RoMVx(z) = DTOL(RXcursor(e) * RXscale(e) + RXshift(e));
                    294:   RoMVy(z) = DTOL(RYcursor(e) * RYscale(e) + RYshift(e));
                    295:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    296:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    297: }
                    298:
                    299: void
                    300: rectmove(long ne, GEN x, GEN y) /* code = 0 */
                    301: {
                    302:   rectmove0(ne,gtodouble(x),gtodouble(y),0);
                    303: }
                    304:
                    305: void
                    306: rectrmove(long ne, GEN x, GEN y) /* code = 0 */
                    307: {
                    308:   rectmove0(ne,gtodouble(x),gtodouble(y),1);
                    309: }
                    310:
                    311: void
                    312: rectpoint0(long ne, double x, double y,long relative) /* code = 1 */
                    313: {
                    314:   PariRect *e = check_rect_init(ne);
                    315:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj1P));
                    316:
                    317:   if (relative)
                    318:    { RXcursor(e) += x; RYcursor(e) += y; }
                    319:   else
                    320:    { RXcursor(e) = x; RYcursor(e) = y; }
                    321:   RoNext(z)=0;
                    322:   RoPTx(z) = DTOL(RXcursor(e)*RXscale(e) + RXshift(e));
                    323:   RoPTy(z) = DTOL(RYcursor(e)*RYscale(e) + RYshift(e));
                    324:   RoType(z) = ((RoPTx(z)<0)||(RoPTy(z)<0)||(RoPTx(z)>RXsize(e))
                    325:               ||(RoPTy(z)>RYsize(e))) ? ROt_MV : ROt_PT;
                    326:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    327:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    328:   RoCol(z)=current_color[ne];
                    329: }
                    330:
                    331: void
                    332: rectpoint(long ne, GEN x, GEN y) /* code = 1 */
                    333: {
                    334:   rectpoint0(ne,gtodouble(x),gtodouble(y),0);
                    335: }
                    336:
                    337: void
                    338: rectrpoint(long ne, GEN x, GEN y) /* code = 1 */
                    339: {
                    340:   rectpoint0(ne,gtodouble(x),gtodouble(y),1);
                    341: }
                    342:
                    343: void
                    344: rectcolor(long ne, long color)
                    345: {
                    346:   check_rect(ne);
                    347:   if (!GOODCOLOR(color)) err(talker,"This is not a valid color");
                    348:   current_color[ne]=color;
                    349: }
                    350:
                    351: void
                    352: rectline0(long ne, double gx2, double gy2, long relative) /* code = 2 */
                    353: {
                    354:   long dx,dy,dxy,xmin,xmax,ymin,ymax,x1,y1,x2,y2;
                    355:   PariRect *e = check_rect_init(ne);
                    356:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
                    357:
                    358:   x1 = DTOL(RXcursor(e)*RXscale(e) + RXshift(e));
                    359:   y1 = DTOL(RYcursor(e)*RYscale(e) + RYshift(e));
                    360:   if (relative)
                    361:     { RXcursor(e)+=gx2; RYcursor(e)+=gy2; }
                    362:   else
                    363:     { RXcursor(e)=gx2; RYcursor(e)=gy2; }
                    364:   x2 = DTOL(RXcursor(e)*RXscale(e) + RXshift(e));
                    365:   y2 = DTOL(RYcursor(e)*RYscale(e) + RYshift(e));
                    366:   xmin = max(min(x1,x2),0); xmax = min(max(x1,x2),RXsize(e));
                    367:   ymin = max(min(y1,y2),0); ymax = min(max(y1,y2),RYsize(e));
                    368:   dxy = x1*y2 - y1*x2; dx = x2-x1; dy = y2-y1;
                    369:   if (dy)
                    370:   {
                    371:     if (dx*dy<0)
                    372:       { xmin = max(xmin,(dxy+RYsize(e)*dx)/dy); xmax=min(xmax,dxy/dy); }
                    373:     else
                    374:       { xmin=max(xmin,dxy/dy); xmax=min(xmax,(dxy+RYsize(e)*dx)/dy); }
                    375:   }
                    376:   if (dx)
                    377:   {
                    378:     if (dx*dy<0)
                    379:       { ymin=max(ymin,(RXsize(e)*dy-dxy)/dx); ymax=min(ymax,-dxy/dx); }
                    380:     else
                    381:       { ymin=max(ymin,-dxy/dx); ymax=min(ymax,(RXsize(e)*dy-dxy)/dx); }
                    382:   }
                    383:   RoNext(z)=0;
                    384:   RoLNx1(z) = xmin; RoLNx2(z) = xmax;
                    385:   if (dx*dy<0) { RoLNy1(z) = ymax; RoLNy2(z) = ymin; }
                    386:   else { RoLNy1(z) = ymin; RoLNy2(z) = ymax; }
                    387:   RoType(z) = ((xmin>xmax)||(ymin>ymax)) ? ROt_MV : ROt_LN;
                    388:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    389:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    390:   RoCol(z)=current_color[ne];
                    391: }
                    392:
                    393: void
                    394: rectline(long ne, GEN gx2, GEN gy2) /* code = 2 */
                    395: {
                    396:   rectline0(ne, gtodouble(gx2), gtodouble(gy2),0);
                    397: }
                    398:
                    399: void
                    400: rectrline(long ne, GEN gx2, GEN gy2) /* code = 2 */
                    401: {
                    402:   rectline0(ne, gtodouble(gx2), gtodouble(gy2),1);
                    403: }
                    404:
                    405: void
                    406: rectbox0(long ne, double gx2, double gy2, long relative)
                    407: {
                    408:   long x1,y1,x2,y2,xmin,ymin,xmax,ymax;
                    409:   double xx,yy;
                    410:   PariRect *e = check_rect_init(ne);
                    411:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
                    412:
                    413:   x1 = DTOL(RXcursor(e)*RXscale(e) + RXshift(e));
                    414:   y1 = DTOL(RYcursor(e)*RYscale(e) + RYshift(e));
                    415:   if (relative)
                    416:   { xx = RXcursor(e)+gx2; yy = RYcursor(e)+gy2; }
                    417:   else
                    418:   {  xx = gx2; yy = gy2; }
                    419:   x2=DTOL(xx*RXscale(e) + RXshift(e));
                    420:   y2=DTOL(yy*RYscale(e) + RYshift(e));
                    421:   xmin = max(min(x1,x2),0); xmax = min(max(x1,x2),RXsize(e));
                    422:   ymin = max(min(y1,y2),0); ymax = min(max(y1,y2),RYsize(e));
                    423:
                    424:   RoNext(z)=0; RoType(z) = ROt_BX;
                    425:   RoBXx1(z) = xmin; RoBXy1(z) = ymin;
                    426:   RoBXx2(z) = xmax; RoBXy2(z) = ymax;
                    427:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    428:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    429:   RoCol(z)=current_color[ne];
                    430: }
                    431:
                    432: void
                    433: rectbox(long ne, GEN gx2, GEN gy2) /* code = 3 */
                    434: {
                    435:   rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 0);
                    436: }
                    437:
                    438: void
                    439: rectrbox(long ne, GEN gx2, GEN gy2) /* code = 3 */
                    440: {
                    441:   rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 1);
                    442: }
                    443:
                    444: void
                    445: killrect(long ne)
                    446: {
                    447:   RectObj *p1,*p2;
                    448:   PariRect *e = check_rect_init(ne);
                    449:
                    450:   current_color[ne]=DEFAULT_COLOR;
                    451:   p1=RHead(e);
                    452:   RHead(e) = RTail(e) = NULL;
                    453:   RXsize(e) = RYsize(e) = 0;
                    454:   RXcursor(e) = RYcursor(e) = 0;
                    455:   RXscale(e) = RYscale(e) = 1;
                    456:   RXshift(e) = RYshift(e) = 0;
                    457:   while (p1)
                    458:   {
                    459:     if (RoType(p1)==ROt_MP || RoType(p1)==ROt_ML)
                    460:     {
                    461:       free(RoMPxs(p1)); free(RoMPys(p1));
                    462:     }
                    463:     if (RoType(p1)==ROt_ST) free(RoSTs(p1));
                    464:     p2=RoNext(p1); free(p1); p1=p2;
                    465:   }
                    466: }
                    467:
                    468: void
                    469: rectpoints0(long ne, double *listx, double *listy, long lx) /* code = 4 */
                    470: {
                    471:   long *ptx,*pty,x,y,i,cp=0;
                    472:   PariRect *e = check_rect_init(ne);
                    473:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjMP));
                    474:
                    475:   ptx=(long*) gpmalloc(lx*sizeof(long));
                    476:   pty=(long*) gpmalloc(lx*sizeof(long));
                    477:   for (i=0; i<lx; i++)
                    478:   {
                    479:     x=DTOL(RXscale(e)*listx[i] + RXshift(e));
                    480:     y=DTOL(RYscale(e)*listy[i] + RYshift(e));
                    481:     if ((x>=0)&&(y>=0)&&(x<=RXsize(e))&&(y<=RYsize(e)))
                    482:     {
                    483:       ptx[cp]=x; pty[cp]=y; cp++;
                    484:     }
                    485:   }
                    486:   RoNext(z)=0; RoType(z) = ROt_MP;
                    487:   RoMPcnt(z) = cp; RoMPxs(z) = ptx; RoMPys(z) = pty;
                    488:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    489:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    490:   RoCol(z)=current_color[ne];
                    491: }
                    492:
                    493: void
                    494: rectpoints(long ne, GEN listx, GEN listy) /* code = 4 */
                    495: {
                    496:   long i,lx, tx=typ(listx), ty=typ(listy);
                    497:   double *px,*py;
                    498:
                    499:   if (!is_matvec_t(tx) || !is_matvec_t(ty))
                    500:   {
                    501:     rectpoint(ne, listx, listy); return;
                    502:   }
                    503:   if (tx == t_MAT || ty == t_MAT) err(ploter4);
                    504:   lx=lg(listx); if (lg(listy)!=lx) err(ploter5);
                    505:   lx--; if (!lx) return;
                    506:
                    507:   px = (double*) gpmalloc(lx*sizeof(double));
                    508:   py = (double*) gpmalloc(lx*sizeof(double));
                    509:   for (i=0; i<lx; i++)
                    510:   {
                    511:     px[i]=gtodouble((GEN)listx[i+1]); py[i]=gtodouble((GEN)listy[i+1]);
                    512:   }
                    513:   rectpoints0(ne,px,py,lx);
                    514:   free(px); free(py);
                    515: }
                    516:
                    517: void
                    518: rectlines0(long ne, double *x, double *y, long lx, long flag) /* code = 5 */
                    519: {
                    520:   long i,I,*ptx,*pty;
                    521:   PariRect *e = check_rect_init(ne);
                    522:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
                    523:
                    524:   I = flag ? lx+1 : lx;
                    525:   ptx = (long*) gpmalloc(I*sizeof(long));
                    526:   pty = (long*) gpmalloc(I*sizeof(long));
                    527:   for (i=0; i<lx; i++)
                    528:   {
                    529:     ptx[i]=DTOL(RXscale(e)*x[i] + RXshift(e));
                    530:     pty[i]=DTOL(RYscale(e)*y[i] + RYshift(e));
                    531:   }
                    532:   if (flag)
                    533:   {
                    534:     ptx[i]=DTOL(RXscale(e)*x[0] + RXshift(e));
                    535:     pty[i]=DTOL(RYscale(e)*y[0] + RYshift(e));
                    536:   }
                    537:   RoNext(z)=0; RoType(z)=ROt_ML;
                    538:   RoMLcnt(z)=lx; RoMLxs(z)=ptx; RoMLys(z)=pty;
                    539:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    540:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    541:   RoCol(z) = current_color[ne];
                    542: }
                    543:
                    544: void
                    545: rectlines(long ne, GEN listx, GEN listy, long flag) /* code = 5 */
                    546: {
                    547:   long tx=typ(listx), ty=typ(listy), lx=lg(listx), i;
                    548:   double *x, *y;
                    549:
                    550:   if (!is_matvec_t(tx) || !is_matvec_t(ty))
                    551:   {
                    552:     rectline(ne, listx, listy); return;
                    553:   }
                    554:   if (tx == t_MAT || ty == t_MAT) err(ploter4);
                    555:   if (lg(listy)!=lx) err(ploter5);
                    556:   lx--; if (!lx) return;
                    557:
                    558:   x = (double*) gpmalloc(lx*sizeof(double));
                    559:   y = (double*) gpmalloc(lx*sizeof(double));
                    560:   for (i=0; i<lx; i++)
                    561:   {
                    562:     x[i] = gtodouble((GEN)listx[i+1]); y[i] = gtodouble((GEN)listy[i+1]);
                    563:   }
                    564:   rectlines0(ne,x,y,lx,flag);
                    565:   free(x); free(y);
                    566: }
                    567:
                    568: static void
                    569: put_string(long win, long x, long y, char *str)
                    570: {
                    571:   rectmove0(win,(double)x,(double)y,0); rectstring(win,str);
                    572: }
                    573:
                    574: /* Allocate memory, then put string */
                    575: void
                    576: rectstring(long ne, char *str) /* code = 6 */
                    577: {
                    578:   PariRect *e = check_rect_init(ne);
                    579:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjST));
                    580:   long l = strlen(str);
                    581:   char *s = (char *) gpmalloc(l+1);
                    582:
                    583:   strcpy(s,str);
                    584:   RoNext(z)=0; RoType(z) = ROt_ST;
                    585:   RoSTl(z) = l; RoSTs(z) = s;
                    586:   RoSTx(z) = DTOL(RXscale(e)*RXcursor(e)+RXshift(e));
                    587:   RoSTy(z) = DTOL(RYscale(e)*RYcursor(e)+RYshift(e));
                    588:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    589:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    590:   RoCol(z)=current_color[ne];
                    591: }
                    592:
                    593: void
                    594: rectpointtype(long ne, long type) /* code = 0 */
                    595: {
                    596:  if (ne == -1) {
                    597:      rectpoint_itype = type;
                    598:  } else {
                    599:   PariRect *e = check_rect_init(ne);
                    600:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPN));
                    601:
                    602:   RoNext(z) = 0; RoType(z) = ROt_PTT;
                    603:   RoPTTpen(z) = type;
                    604:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    605:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    606:  }
                    607: }
                    608:
                    609: void
                    610: rectpointsize(long ne, GEN size) /* code = 0 */
                    611: {
                    612:  if (ne == -1) {
                    613:      set_pointsize(gtodouble(size));   /* Immediate set */
                    614:  } else {
                    615:      PariRect *e = check_rect_init(ne);
                    616:      RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPS));
                    617:
                    618:      RoNext(z) = 0; RoType(z) = ROt_PTS;
                    619:      RoPTSsize(z) = gtodouble(size);
                    620:      if (!RHead(e)) RHead(e)=RTail(e)=z;
                    621:      else { RoNext(RTail(e))=z; RTail(e)=z; }
                    622:  }
                    623: }
                    624:
                    625: void
                    626: rectlinetype(long ne, long type)
                    627: {
                    628:  if (ne == -1) {
                    629:      rectline_itype = type;
                    630:  } else {
                    631:   PariRect *e = check_rect_init(ne);
                    632:   RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPN));
                    633:
                    634:   RoNext(z) = 0; RoType(z) = ROt_LNT;
                    635:   RoLNTpen(z) = type;
                    636:   if (!RHead(e)) RHead(e)=RTail(e)=z;
                    637:   else { RoNext(RTail(e))=z; RTail(e)=z; }
                    638:  }
                    639: }
                    640:
                    641: void
                    642: rectcopy(long source, long dest, long xoff, long yoff)
                    643: {
                    644:   PariRect *s = check_rect_init(source), *d = check_rect_init(dest);
                    645:   RectObj *p1 = RHead(s);
                    646:   RectObj *tail = RTail(d);
                    647:   RectObj *next;
                    648:   int i;
                    649:
                    650:   while (p1)
                    651:   {
                    652:     switch(RoType(p1))
                    653:     {
                    654:       case ROt_PT:
                    655:        next = (RectObj*) gpmalloc(sizeof(RectObj1P));
                    656:        memcpy(next,p1,sizeof(RectObj1P));
                    657:        RoPTx(next) += xoff; RoPTy(next) += yoff;
                    658:        RoNext(tail) = next; tail = next;
                    659:        break;
                    660:       case ROt_LN: case ROt_BX:
                    661:        next = (RectObj*) gpmalloc(sizeof(RectObj2P));
                    662:        memcpy(next,p1,sizeof(RectObj2P));
                    663:        RoLNx1(next) += xoff; RoLNy1(next) += yoff;
                    664:        RoLNx2(next) += xoff; RoLNy2(next) += yoff;
                    665:        RoNext(tail) = next; tail = next;
                    666:        break;
                    667:       case ROt_MP: case ROt_ML:
                    668:        next = (RectObj*) gpmalloc(sizeof(RectObjMP));
                    669:        memcpy(next,p1,sizeof(RectObjMP));
                    670:        RoMPxs(next) = (long*) gpmalloc(sizeof(long)*RoMPcnt(next));
                    671:        RoMPys(next) = (long*) gpmalloc(sizeof(long)*RoMPcnt(next));
                    672:        memcpy(RoMPxs(next),RoMPxs(p1),sizeof(long)*RoMPcnt(next));
                    673:        memcpy(RoMPys(next),RoMPys(p1),sizeof(long)*RoMPcnt(next));
                    674:        for (i=0; i<RoMPcnt(next); i++)
                    675:        {
                    676:          RoMPxs(next)[i] += xoff; RoMPys(next)[i] += yoff;
                    677:         }
                    678:        RoNext(tail) = next; tail = next;
                    679:        break;
                    680:       case ROt_ST:
                    681:        next = (RectObj*) gpmalloc(sizeof(RectObjMP));
                    682:        memcpy(next,p1,sizeof(RectObjMP));
                    683:        RoSTs(next) = (char*) gpmalloc(RoSTl(p1)+1);
                    684:        memcpy(RoSTs(next),RoSTs(p1),RoSTl(p1)+1);
                    685:        RoSTx(next) += xoff; RoSTy(next) += yoff;
                    686:        RoNext(tail) = next; tail = next;
                    687:        break;
                    688:       case ROt_PTT: case ROt_LNT: case ROt_PTS:
                    689:        next = (RectObj*) gpmalloc(sizeof(RectObjPN));
                    690:        memcpy(next,p1,sizeof(RectObjPN));
                    691:        RoNext(tail) = next; tail = next;
                    692:        break;
                    693:     }
                    694:     p1=RoNext(p1);
                    695:   }
                    696:   RoNext(tail) = NULL; RTail(d) = tail;
                    697: }
                    698:
                    699: #define CLIPLINE_NONEMPTY      1
                    700: #define CLIPLINE_CLIP_1                2
                    701: #define CLIPLINE_CLIP_2                4
                    702:
                    703: /* A simpler way is to clip by 4 half-planes */
                    704: static int
                    705: clipline(long xmin, long xmax, long ymin, long ymax, long *x1p, long *y1p, long *x2p, long *y2p)
                    706: {
                    707:     int xy_exch = 0, rc = CLIPLINE_NONEMPTY;
                    708:     double t;
                    709:     double xi, yi;
                    710:     double sl;
                    711:     double xmn, xmx;
                    712:     double ymn, ymx;
                    713:     int x1_is_ymn, x1_is_xmn;
                    714:     double x1 = *x1p, x2 = *x2p, y1 = *y1p, y2 = *y2p;
                    715:
                    716:     if ((x1 < xmin &&  x2 < xmin) || (x1 > xmax && x2 > xmax))
                    717:        return 0;
                    718:     if (fabs(x1 - x2) < fabs(y1 - y2)) { /* Exchange x and y */
                    719:        xy_exch = 1;
                    720:        t = xmin, xmin = ymin, ymin = t;
                    721:        t = xmax, xmax = ymax, ymax = t;
                    722:        t = x1, x1 = y1, y1 = t;
                    723:        t = x2, x2 = y2, y2 = t;
                    724:     }
                    725:
                    726:     /* Build y as a function of x */
                    727:     xi = x1, yi = y1;
                    728:     sl = ( (x1==x2) ? 0 : (y2 - yi)/(x2 - xi) );
                    729:
                    730:     xmn = x1, xmx = x2;
                    731:     if (x1 > x2) {
                    732:        x1_is_xmn = 0;
                    733:        xmn = x2, xmx = x1;
                    734:     } else
                    735:        x1_is_xmn = 1;
                    736:
                    737:     if (xmn < xmin) {
                    738:        xmn = xmin;
                    739:        rc |= (x1_is_xmn ? CLIPLINE_CLIP_1 : CLIPLINE_CLIP_2);
                    740:     }
                    741:
                    742:     if (xmx > xmax) {
                    743:        xmx = xmax;
                    744:        rc |= (x1_is_xmn ? CLIPLINE_CLIP_2 : CLIPLINE_CLIP_1);
                    745:     }
                    746:
                    747:     if (xmn > xmx)
                    748:        return 0;
                    749:
                    750:     ymn = yi + (xmn - xi)*sl;
                    751:     ymx = yi + (xmx - xi)*sl;
                    752:
                    753:     if (sl < 0)
                    754:        t = ymn, ymn = ymx, ymx = t;
                    755:     if (ymn > ymax || ymx < ymin)
                    756:        return 0;
                    757:
                    758:     if (rc & CLIPLINE_CLIP_1)
                    759:        x1 = (x1_is_xmn ? xmn : xmx);
                    760:     if (rc & CLIPLINE_CLIP_2)
                    761:        x2 = (x1_is_xmn ? xmx : xmn);
                    762:
                    763:     /* Now we know there is an intersection, need to move x1 and x2 */
                    764:     x1_is_ymn = ((sl >= 0) == (x1 < x2));
                    765:     if (ymn < ymin) {
                    766:        double x = (ymin - yi)/sl + xi; /* slope != 0  ! */
                    767:
                    768:        if (x1_is_ymn)
                    769:            x1 = x, rc |= CLIPLINE_CLIP_1;
                    770:        else
                    771:            x2 = x, rc |= CLIPLINE_CLIP_2;
                    772:     }
                    773:     if (ymx > ymax) {
                    774:        double x = (ymax - yi)/sl + xi; /* slope != 0  ! */
                    775:
                    776:        if (x1_is_ymn)
                    777:            x2 = x, rc |= CLIPLINE_CLIP_2;
                    778:        else
                    779:            x1 = x, rc |= CLIPLINE_CLIP_1;
                    780:     }
                    781:     if (rc & CLIPLINE_CLIP_1)
                    782:        y1 = yi + (x1 - xi)*sl;
                    783:     if (rc & CLIPLINE_CLIP_2)
                    784:        y2 = yi + (x2 - xi)*sl;
                    785:     if (xy_exch)                       /* Exchange x and y */
                    786:        *x1p = y1, *x2p = y2, *y1p = x1, *y2p = x2;
                    787:     else
                    788:        *x1p = x1, *x2p = x2, *y1p = y1, *y2p = y2;
                    789:     return rc;
                    790: }
                    791:
                    792: void
                    793: rectclip(long rect)
                    794: {
                    795:   PariRect *s = check_rect_init(rect);
                    796:   RectObj *p1 = RHead(s);
                    797:   RectObj **prevp = &RHead(s);
                    798:   RectObj *next;
                    799:   double xmin = 0;
                    800:   double xmax = RXsize(s);
                    801:   double ymin = 0;
                    802:   double ymax = RYsize(s);
                    803:
                    804:   while (p1) {
                    805:       int did_clip = 0;
                    806:
                    807:       next = RoNext(p1);
                    808:       switch(RoType(p1)) {
                    809:       case ROt_PT:
                    810:          if ( RoPTx(p1) < xmin || RoPTx(p1) > xmax
                    811:               || RoPTy(p1) < ymin || RoPTy(p1) > ymax) {
                    812:                 remove:
                    813:              *prevp = next;
                    814:              free(p1);
                    815:              break;
                    816:          }
                    817:          goto do_next;
                    818:       case ROt_BX:
                    819:          if (RoLNx1(p1) < xmin)
                    820:              RoLNx1(p1) = xmin, did_clip = 1;
                    821:          if (RoLNx2(p1) < xmin)
                    822:              RoLNx2(p1) = xmin, did_clip = 1;
                    823:          if (RoLNy1(p1) < ymin)
                    824:              RoLNy1(p1) = ymin, did_clip = 1;
                    825:          if (RoLNy2(p1) < ymin)
                    826:              RoLNy2(p1) = ymin, did_clip = 1;
                    827:          if (RoLNx1(p1) > xmax)
                    828:              RoLNx1(p1) = xmax, did_clip = 1;
                    829:          if (RoLNx2(p1) > xmax)
                    830:              RoLNx2(p1) = xmax, did_clip = 1;
                    831:          if (RoLNy1(p1) > ymax)
                    832:              RoLNy1(p1) = ymax, did_clip = 1;
                    833:          if (RoLNy2(p1) > ymax)
                    834:              RoLNy2(p1) = ymax, did_clip = 1;
                    835:          /* Remove zero-size clipped boxes */
                    836:          if ( did_clip
                    837:               && RoLNx1(p1) == RoLNx2(p1) && RoLNy1(p1) == RoLNy2(p1) )
                    838:              goto remove;
                    839:          goto do_next;
                    840:       case ROt_LN:
                    841:          if (!clipline(xmin, xmax, ymin, ymax,
                    842:                        &RoLNx1(p1), &RoLNy1(p1), &RoLNx2(p1), &RoLNy2(p1)))
                    843:              goto remove;
                    844:          goto do_next;
                    845:       case ROt_MP:
                    846:       {
                    847:          int c = RoMPcnt(p1);
                    848:          int f = 0, t = 0;
                    849:
                    850:          while (f < c) {
                    851:              if ( RoMPxs(p1)[f] >= xmin && RoMPxs(p1)[f] <= xmax
                    852:                   && RoMPys(p1)[f] >= ymin && RoMPys(p1)[f] <= ymax) {
                    853:                  if (t != f) {
                    854:                      RoMPxs(p1)[t] = RoMPxs(p1)[f];
                    855:                      RoMPys(p1)[t] = RoMPys(p1)[f];
                    856:                  }
                    857:                  t++;
                    858:              }
                    859:              f++;
                    860:          }
                    861:          if (t == 0)
                    862:              goto remove;
                    863:          RoMPcnt(p1) = t;
                    864:          goto do_next;
                    865:       }
                    866:       case ROt_ML:
                    867:       {
                    868:          /* Hard case.  Here we need to break a multiline into
                    869:             several pieces if some part is clipped. */
                    870:          int c = RoMPcnt(p1) - 1;
                    871:          int f = 0, t = 0, had_lines = 0, had_hole = 0, rc;
                    872:          long ox = RoMLxs(p1)[0], oy = RoMLys(p1)[0], oxn, oyn;
                    873:
                    874:          while (f < c) {
                    875:              /* Endpoint of this segment is the startpoint of the
                    876:                 next one, so we need to preserve it if it is clipped. */
                    877:              oxn = RoMLxs(p1)[f + 1], oyn = RoMLys(p1)[f + 1];
                    878:              rc = clipline(xmin, xmax, ymin, ymax,
                    879:                            /* &RoMLxs(p1)[f], &RoMLys(p1)[f], */
                    880:                            &ox, &oy,
                    881:                            &RoMLxs(p1)[f+1], &RoMLys(p1)[f+1]);
                    882:              RoMLxs(p1)[f] = ox; RoMLys(p1)[f] = oy;
                    883:              ox = oxn; oy = oyn;
                    884:              if (!rc) {
                    885:                  if (had_lines)
                    886:                      had_hole = 1;
                    887:                  f++;
                    888:                  continue;
                    889:              }
                    890:
                    891:              if ( !had_lines || (!(rc & CLIPLINE_CLIP_1) && !had_hole) ) {
                    892:                  /* Continuous */
                    893:                  had_lines = 1;
                    894:                  if (t != f) {
                    895:                      if (t == 0) {
                    896:                          RoMPxs(p1)[t] = RoMPxs(p1)[f];
                    897:                          RoMPys(p1)[t] = RoMPys(p1)[f];
                    898:                      }
                    899:                      RoMPxs(p1)[t+1] = RoMPxs(p1)[f+1];
                    900:                      RoMPys(p1)[t+1] = RoMPys(p1)[f+1];
                    901:                  }
                    902:                  t++;
                    903:                  f++;
                    904:                  if ( rc & CLIPLINE_CLIP_2)
                    905:                      had_hole = 1, RoMLcnt(p1) = t + 1;
                    906:                  continue;
                    907:              }
                    908:              /* Is not continuous, automatically p1 is not free()ed.  */
                    909:              RoMLcnt(p1) = t + 1;
                    910:              if ( rc & CLIPLINE_CLIP_2) { /* Needs separate entry */
                    911:                  RectObj *n = (RectObj*) gpmalloc(sizeof(RectObj2P));
                    912:
                    913:                  RoType(n) = ROt_LN;
                    914:                  RoCol(n) = RoCol(p1);
                    915:                  RoLNx1(n) = RoMLxs(p1)[f];    RoLNy1(n) = RoMLys(p1)[f];
                    916:                  RoLNx2(n) = RoMLxs(p1)[f+1];  RoLNy2(n) = RoMLys(p1)[f+1];
                    917:                  RoNext(n) = next;
                    918:                  RoNext(p1) = n;
                    919:                  /* Restore the unclipped value: */
                    920:                  RoMLxs(p1)[f + 1] = oxn;      RoMLys(p1)[f + 1] = oyn;
                    921:                  f++;
                    922:                  prevp = &RoNext(n);
                    923:              }
                    924:              if (f + 1 < c) {          /* Are other lines */
                    925:                  RectObj *n = (RectObj*) gpmalloc(sizeof(RectObjMP));
                    926:                  RoType(n) = ROt_ML;
                    927:                  RoCol(n) = RoCol(p1);
                    928:                  RoMLcnt(n) = c - f;
                    929:                  RoMLxs(n) = (long*) gpmalloc(sizeof(long)*(c - f));
                    930:                  RoMLys(n) = (long*) gpmalloc(sizeof(long)*(c - f));
                    931:                  memcpy(RoMPxs(n),RoMPxs(p1) + f, sizeof(long)*(c - f));
                    932:                  memcpy(RoMPys(n),RoMPys(p1) + f, sizeof(long)*(c - f));
                    933:                  RoMPxs(n)[0] = oxn; RoMPys(n)[0] = oyn;
                    934:                  RoNext(n) = next;
                    935:                  RoNext(p1) = n;
                    936:                  next = n;
                    937:              }
                    938:              goto do_next;
                    939:          }
                    940:          if (t == 0)
                    941:              goto remove;
                    942:          goto do_next;
                    943:       }
                    944: #if 0
                    945:       case ROt_ST:
                    946:          next = (RectObj*) gpmalloc(sizeof(RectObjMP));
                    947:          memcpy(next,p1,sizeof(RectObjMP));
                    948:          RoSTs(next) = (char*) gpmalloc(RoSTl(p1)+1);
                    949:          memcpy(RoSTs(next),RoSTs(p1),RoSTl(p1)+1);
                    950:          RoSTx(next) += xoff; RoSTy(next) += yoff;
                    951:          RoNext(tail) = next; tail = next;
                    952:          break;
                    953: #endif
                    954:       default: {
                    955:        do_next:
                    956:              prevp = &RoNext(p1);
                    957:              break;
                    958:          }
                    959:       }
                    960:       p1 = next;
                    961:   }
                    962: }
                    963:
                    964: /********************************************************************/
                    965: /**                                                                **/
                    966: /**                        HI-RES PLOT                             **/
                    967: /**                                                                **/
                    968: /********************************************************************/
                    969:
                    970: static void
                    971: Appendx(dblPointList *f, dblPointList *l,double x)
                    972: {
                    973:   (l->d)[l->nb++]=x;
                    974:   if (x < f->xsml) f->xsml=x;
                    975:   else if (x > f->xbig) f->xbig=x;
                    976: }
                    977:
                    978: static void
                    979: Appendy(dblPointList *f, dblPointList *l,double y)
                    980: {
                    981:   (l->d)[l->nb++]=y;
                    982:   if (y < f->ysml) f->ysml=y;
                    983:   else if (y > f->ybig) f->ybig=y;
                    984: }
                    985:
                    986: /* Convert data from GEN to double before we call rectplothrawin */
                    987: static dblPointList*
                    988: gtodblList(GEN data, long flags)
                    989: {
                    990:   dblPointList *l;
                    991:   double xsml,xbig,ysml,ybig;
                    992:   long tx=typ(data), ty, nl=lg(data)-1, lx1,lx;
                    993:   register long i,j,u,v;
                    994:   long param=(flags & PLOT_PARAMETRIC);
                    995:   GEN x,y;
                    996:
                    997:   if (! is_vec_t(tx)) err(talker,"not a vector in gtodblList");
                    998:   if (!nl) return NULL;
                    999:   lx1=lg(data[1]);
                   1000:
                   1001:   /* Allocate memory, then convert coord. to double */
                   1002:   l=(dblPointList*) gpmalloc(nl*sizeof(dblPointList));
                   1003:   for (i=0; i<nl-1; i+=2)
                   1004:   {
                   1005:     u=i+1;
                   1006:     x=(GEN)data[u]; tx = typ(x);
                   1007:     y=(GEN)data[u+1]; ty = typ(y);
                   1008:     if (!is_vec_t(tx) || !is_vec_t(ty)) err(ploter4);
                   1009:     lx=lg(x); if (lg(y)!=lx) err(ploter5);
                   1010:     if (!param && lx != lx1) err(ploter5);
                   1011:     lx--;
                   1012:
                   1013:     if (lx)
                   1014:     {
                   1015:       l[i].d = (double*) gpmalloc(lx*sizeof(double));
                   1016:       l[u].d = (double*) gpmalloc(lx*sizeof(double));
                   1017:
                   1018:       for (j=0; j<lx; j=v)
                   1019:       {
                   1020:        v=j+1;
                   1021:        l[i].d[j]=gtodouble((GEN)x[v]);
                   1022:        l[u].d[j]=gtodouble((GEN)y[v]);
                   1023:       }
                   1024:     }
                   1025:     l[i].nb=l[u].nb=lx;
                   1026:   }
                   1027:
                   1028:   xsml=xbig=l[0].d[0]; ysml=ybig=l[1].d[0];
                   1029:
                   1030:   /* Now compute extremas */
                   1031:   if (param)
                   1032:   {
                   1033:     l[0].nb=nl/2;
                   1034:     for (i=0; i<l[0].nb; i+=2)
                   1035:     {
                   1036:       u=i+1;
                   1037:       for (j=0; j<l[u].nb; j++)
                   1038:       {
                   1039:        if (l[i].d[j]<xsml)
                   1040:          xsml=l[i].d[j];
                   1041:        else
                   1042:          if (l[i].d[j]>xbig) xbig=l[i].d[j];
                   1043:
                   1044:        if (l[u].d[j]<ysml)
                   1045:          ysml=l[u].d[j];
                   1046:        else
                   1047:          if (l[u].d[j]>ybig) ybig=l[u].d[j];
                   1048:       }
                   1049:     }
                   1050:   }
                   1051:   else
                   1052:   {
                   1053:     l[0].nb=nl-1;
                   1054:     for (j=0; j<l[1].nb; j++)
                   1055:     {
                   1056:       if (l[0].d[j]<xsml)
                   1057:        xsml=l[0].d[j];
                   1058:       else
                   1059:        if (l[0].d[j]>xbig) xbig=l[0].d[j];
                   1060:     }
                   1061:     for (i=1; i<=l[0].nb; i++)
                   1062:     {
                   1063:       for (j=0; j<l[i].nb; j++)
                   1064:       {
                   1065:        if (l[i].d[j]<ysml)
                   1066:          ysml=l[i].d[j];
                   1067:        else
                   1068:          if (l[i].d[j]>ybig) ybig=l[i].d[j];
                   1069:       }
                   1070:     }
                   1071:   }
                   1072:   l[0].xsml=xsml; l[0].xbig=xbig;
                   1073:   l[0].ysml=ysml; l[0].ybig=ybig;
                   1074:   return l;
                   1075: }
                   1076:
                   1077: static void
                   1078: single_recursion(dblPointList *pl,char *ch,entree *ep,GEN xleft,GEN yleft,
                   1079:   GEN xright,GEN yright,long depth)
                   1080: {
                   1081:   GEN xx,yy;
                   1082:   long av=avma;
                   1083:   double dy=pl[0].ybig - pl[0].ysml;
                   1084:
                   1085:   if (depth==RECUR_MAXDEPTH) return;
                   1086:
                   1087:   yy=cgetr(3); xx=gmul2n(gadd(xleft,xright),-1);
                   1088:   ep->value = (void*) xx; gaffect(READ_EXPR(ch), yy);
                   1089:
                   1090:   if (dy)
                   1091:   {
                   1092:     if (fabs(rtodbl(yleft)+rtodbl(yright)-2*rtodbl(yy))/dy < RECUR_PREC)
                   1093:       return;
                   1094:   }
                   1095:   single_recursion(pl,ch,ep, xleft,yleft, xx,yy, depth+1);
                   1096:
                   1097:   Appendx(&pl[0],&pl[0],rtodbl(xx));
                   1098:   Appendy(&pl[0],&pl[1],rtodbl(yy));
                   1099:
                   1100:   single_recursion(pl,ch,ep, xx,yy, xright,yright, depth+1);
                   1101:
                   1102:   avma=av;
                   1103: }
                   1104:
                   1105: static void
                   1106: param_recursion(dblPointList *pl,char *ch,entree *ep, GEN tleft,GEN xleft,
                   1107:   GEN yleft, GEN tright,GEN xright,GEN yright, long depth)
                   1108: {
                   1109:   GEN tt,xx,yy, p1;
                   1110:   long av=avma;
                   1111:   double dy=pl[0].ybig - pl[0].ysml;
                   1112:   double dx=pl[0].xbig - pl[0].xsml;
                   1113:
                   1114:   if (depth==PARAMR_MAXDEPTH) return;
                   1115:
                   1116:   xx=cgetr(3); yy=cgetr(3); tt=gmul2n(gadd(tleft,tright),-1);
                   1117:   ep->value = (void*)tt; p1=READ_EXPR(ch);
                   1118:   gaffect((GEN)p1[1],xx); gaffect((GEN)p1[2],yy);
                   1119:
                   1120:   if (dx && dy)
                   1121:   {
                   1122:     if (fabs(rtodbl(xleft)+rtodbl(xright)-2*rtodbl(xx))/dx < RECUR_PREC &&
                   1123:        fabs(rtodbl(yleft)+rtodbl(yright)-2*rtodbl(yy))/dy < RECUR_PREC)
                   1124:         return;
                   1125:   }
                   1126:   param_recursion(pl,ch,ep, tleft,xleft,yleft, tt,xx,yy, depth+1);
                   1127:
                   1128:   Appendx(&pl[0],&pl[0],rtodbl(xx));
                   1129:   Appendy(&pl[0],&pl[1],rtodbl(yy));
                   1130:
                   1131:   param_recursion(pl,ch,ep, tt,xx,yy, tright,xright,yright, depth+1);
                   1132:
                   1133:   avma=av;
                   1134: }
                   1135:
                   1136: /*
                   1137:  *  Pure graphing. If testpoints is 0, it is set to the default.
                   1138:  *  Returns a dblPointList of (absolute) coordinates.
                   1139:  */
                   1140: static dblPointList *
                   1141: rectplothin(entree *ep, GEN a, GEN b, char *ch, long prec, ulong flags,
                   1142:             long testpoints)
                   1143: {
                   1144:   long single_c;
                   1145:   long param=flags & PLOT_PARAMETRIC;
                   1146:   long recur=flags & PLOT_RECURSIVE;
                   1147:   GEN p1,dx,x,xleft,xright,yleft,yright,tleft,tright;
                   1148:   dblPointList *pl;
                   1149:   long tx,av = avma,av2,i,j,sig,nc,nl,nbpoints;
                   1150:   double xsml,xbig,ysml,ybig,fx,fy;
                   1151:
                   1152:   if (!testpoints)
                   1153:   {
                   1154:     if (recur)
                   1155:       testpoints = RECUR_NUMPOINTS;
                   1156:     else
                   1157:       testpoints = param? PARAM_NUMPOINTS : PLOTH_NUMPOINTS;
                   1158:   }
                   1159:   if (recur)
                   1160:     nbpoints = testpoints << (param? PARAMR_MAXDEPTH : RECUR_MAXDEPTH);
                   1161:   else
                   1162:     nbpoints = testpoints;
                   1163:
                   1164:   sig=gcmp(b,a); if (!sig) return 0;
                   1165:   if (sig<0) { x=a; a=b; b=x; }
                   1166:   dx=gdivgs(gsub(b,a),testpoints-1);
                   1167:
                   1168:   x = cgetr(prec); gaffect(a,x); push_val(ep, x);
                   1169:   av2=avma; p1=READ_EXPR(ch); tx=typ(p1);
                   1170:   if (!is_matvec_t(tx))
                   1171:   {
                   1172:     xsml = gtodouble(a);
                   1173:     xbig = gtodouble(b);
                   1174:     ysml = ybig = gtodouble(p1);
                   1175:     nc=1; nl=2; /* nc = nb of curves; nl = nb of coord. lists */
                   1176:     if (param)
                   1177:       err(warner,"flag PLOT_PARAMETRIC ignored");
                   1178:     single_c=1; param=0;
                   1179:   }
                   1180:   else
                   1181:   {
                   1182:     if (tx != t_VEC) err(talker,"not a row vector in ploth");
                   1183:     nl=lg(p1)-1; if (!nl) { avma=av; return 0; }
                   1184:     single_c=0;
                   1185:     if (param) nc=nl/2; else { nc=nl; nl++; }
                   1186:     if (recur && nc > 1) err(talker,"multi-curves cannot be plot recursively");
                   1187:
                   1188:     if (param)
                   1189:     {
                   1190:       xbig=xsml=gtodouble((GEN)p1[1]);
                   1191:       ybig=ysml=gtodouble((GEN)p1[2]);
                   1192:       for (i=3; i<=nl; i++)
                   1193:       {
                   1194:        fx=gtodouble((GEN)p1[i]); i++;
                   1195:         fy=gtodouble((GEN)p1[i]);
                   1196:        if (xsml>fx) xsml=fx; else if (xbig<fx) xbig=fx;
                   1197:        if (ysml>fy) ysml=fy; else if (ybig<fy) ybig=fy;
                   1198:       }
                   1199:     }
                   1200:     else
                   1201:     {
                   1202:       xsml=gtodouble(a); xbig=gtodouble(b);
                   1203:       ysml=gtodouble(vecmin(p1)); ybig=gtodouble(vecmax(p1));
                   1204:     }
                   1205:   }
                   1206:
                   1207:   pl=(dblPointList*) gpmalloc(nl*sizeof(dblPointList));
                   1208:   for (i = 0; i < nl; i++)
                   1209:   {
                   1210:     pl[i].d = (double*) gpmalloc((nbpoints+1)*sizeof(dblPointList));
                   1211:     pl[i].nb=0;
                   1212:   }
                   1213:   pl[0].xsml=xsml; pl[0].ysml=ysml; pl[0].xbig=xbig; pl[0].ybig=ybig;
                   1214:
                   1215:   if (recur) /* recursive plot */
                   1216:   {
                   1217:     xleft=cgetr(3); xright=cgetr(3); yleft=cgetr(3); yright=cgetr(3);
                   1218:     if (param)
                   1219:     {
                   1220:       tleft=cgetr(prec); tright=cgetr(prec);
                   1221:       av2=avma;
                   1222:       gaffect(a,tleft); ep->value = (void*)tleft; p1=READ_EXPR(ch);
                   1223:       gaffect((GEN)p1[1],xleft); gaffect((GEN)p1[2],yleft);
                   1224:       for (i=0; i<testpoints-1; i++)
                   1225:       {
                   1226:        if (i) {
                   1227:          gaffect(tright,tleft); gaffect(xright,xleft); gaffect(yright,yleft);
                   1228:         }
                   1229:        gaddz(tleft,dx,tright); ep->value = (void*)tright;
                   1230:         p1 = READ_EXPR(ch);
                   1231:         if (lg(p1) != 3) err(talker,"inconsistent data in rectplothin");
                   1232:         gaffect((GEN)p1[1],xright); gaffect((GEN)p1[2],yright);
                   1233:
                   1234:        Appendx(&pl[0],&pl[0],rtodbl(xleft));
                   1235:        Appendy(&pl[0],&pl[1],rtodbl(yleft));
                   1236:
                   1237:        param_recursion(pl,ch,ep, tleft,xleft,yleft, tright,xright,yright, 0);
                   1238:        avma=av2;
                   1239:       }
                   1240:       Appendx(&pl[0],&pl[0],rtodbl(xright));
                   1241:       Appendy(&pl[0],&pl[1],rtodbl(yright));
                   1242:     }
                   1243:     else /* single_c */
                   1244:     {
                   1245:       av2=avma;
                   1246:       gaffect(a,xleft); ep->value = (void*) xleft;
                   1247:       gaffect(READ_EXPR(ch),yleft);
                   1248:       for (i=0; i<testpoints-1; i++)
                   1249:       {
                   1250:         gaddz(xleft,dx,xright); ep->value = (void*) xright;
                   1251:        gaffect(READ_EXPR(ch),yright);
                   1252:
                   1253:        Appendx(&pl[0],&pl[0],rtodbl(xleft));
                   1254:        Appendy(&pl[0],&pl[1],rtodbl(yleft));
                   1255:
                   1256:         single_recursion(pl,ch,ep,xleft,yleft,xright,yright,0);
                   1257:         avma=av2;
                   1258:         gaffect(xright,xleft); gaffect(yright,yleft);
                   1259:       }
                   1260:       Appendx(&pl[0],&pl[0],rtodbl(xright));
                   1261:       Appendy(&pl[0],&pl[1],rtodbl(yright));
                   1262:     }
                   1263:   }
                   1264:   else /* non-recursive plot */
                   1265:   {
                   1266:     if (single_c)
                   1267:       for (i=0; i<testpoints; i++)
                   1268:       {
                   1269:        p1 = READ_EXPR(ch);
                   1270:        pl[0].d[i]=gtodouble(x);
                   1271:        Appendy(&pl[0],&pl[1],gtodouble(p1));
                   1272:        gaddz(x,dx,x); avma=av2;
                   1273:       }
                   1274:     else if (param)
                   1275:     {
                   1276:       long k;
                   1277:       double z;
                   1278:
                   1279:       for (i=0; i<testpoints; i++)
                   1280:       {
                   1281:        p1 = READ_EXPR(ch);
                   1282:         if (lg(p1) != nl+1) err(talker,"inconsistent data in rectplothin");
                   1283:        for (j=0; j<nl; j=k)
                   1284:        {
                   1285:          k=j+1; z=gtodouble((GEN)p1[k]);
                   1286:          Appendx(&pl[0], &pl[j],z);
                   1287:
                   1288:          j=k; k++; z=gtodouble((GEN)p1[k]);
                   1289:          Appendy(&pl[0], &pl[j],z);
                   1290:         }
                   1291:        gaddz(x,dx,x); avma=av2;
                   1292:       }
                   1293:     }
                   1294:     else /* plothmult */
                   1295:       for (i=0; i<testpoints; i++)
                   1296:       {
                   1297:        p1 = READ_EXPR(ch);
                   1298:         if (lg(p1) != nl) err(talker,"inconsistent data in rectplothin");
                   1299:        pl[0].d[i]=gtodouble(x);
                   1300:        for (j=1; j<nl; j++) { Appendy(&pl[0],&pl[j],gtodouble((GEN)p1[j])); }
                   1301:        gaddz(x,dx,x); avma=av2;
                   1302:       }
                   1303:   }
                   1304:   pl[0].nb=nc; pop_val(ep); avma = av;
                   1305:   return pl;
                   1306: }
                   1307:
                   1308: GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
                   1309:
                   1310: /* Uses highlevel plotting functions to implement splines as
                   1311:    a low-level plotting function.
                   1312:    Most probably one does not need to make varn==0 into pure variable,
                   1313:    since plotting functions should take care of this. */
                   1314: static void
                   1315: rectsplines(long ne, double *x, double *y, long lx, long flag)
                   1316: {
                   1317:   long i, j, oldavma = avma;
                   1318:   GEN tas, xa = cgetg(lx+1, t_VEC), ya = cgetg(lx+1, t_VEC);
                   1319:   entree *var0 = varentries[ordvar[0]];
                   1320:
                   1321:   if (lx < 4)
                   1322:       err(talker, "Too few points (%ld) for spline plot", lx);
                   1323:   for (i = 1; i <= lx; i++) {
                   1324:       xa[i] = (long) dbltor(x[i-1]);
                   1325:       ya[i] = (long) dbltor(y[i-1]);
                   1326:   }
                   1327:   if (flag & PLOT_PARAMETRIC) {
                   1328:       tas = new_chunk(4);
                   1329:       for (j = 1; j <= 4; j++)
                   1330:          tas[j-1] = (long)stoi(j);
                   1331:       quark_gen = cgetg(2 + 1, t_VEC);
                   1332:   }
                   1333:   for (i = 0; i <= lx - 4; i++) {
                   1334:       long oavma = avma;
                   1335:
                   1336:       xa++; ya++;
                   1337:       if (flag & PLOT_PARAMETRIC) {
                   1338:          quark_gen[1] = (long)polint_i(tas, xa, polx[0], 4, NULL);
                   1339:          quark_gen[2] = (long)polint_i(tas, ya, polx[0], 4, NULL);
                   1340:       } else {
                   1341:          quark_gen = polint_i(xa, ya, polx[0], 4, NULL);
                   1342:           tas = xa;
                   1343:       }
                   1344:       rectploth(ne, var0,
                   1345:                  (GEN)(i==0 ? tas[0] : tas[1]),
                   1346:                  (GEN)(i==lx-4 ? tas[3] : tas[2]),
                   1347:                  QUARK,
                   1348:                 DEFAULTPREC,           /* XXXX precision */
                   1349:                 PLOT_RECURSIVE
                   1350:                    | PLOT_NO_RESCALE
                   1351:                    | PLOT_NO_FRAME
                   1352:                   | PLOT_NO_AXE_Y
                   1353:                    | PLOT_NO_AXE_X
                   1354:                    | (flag & PLOT_PARAMETRIC),
                   1355:                 2);                    /* Start with 3 points */
                   1356:       avma = oavma;
                   1357:   }
                   1358:   avma = oldavma;
                   1359: }
                   1360:
                   1361: /*
                   1362:  * Plot a dblPointList. Complete with axes, bounding box, etc.
                   1363:  * We use two drawing rectangles: one for strings, another
                   1364:  * for graphs.
                   1365:  *
                   1366:  * data is an array of structs. Its meaning depends on flags :
                   1367:  *
                   1368:  * + data[0] contains global extremas, the number of curves to plot
                   1369:  *   (data[0].nb) and a list of doubles (first set of x-coordinates).
                   1370:  *
                   1371:  * + data[i].nb (i>0) contains the number of points in the list
                   1372:  *   data[i].d (hopefully, data[2i].nb=data[2i+1].nb when i>0...)
                   1373:  *
                   1374:  * + If flags contain PLOT_PARAMETRIC, the array length should be
                   1375:  *   even, and successive pairs (data[2i].d, data[2i+1].d) represent
                   1376:  *   curves to plot.
                   1377:  *
                   1378:  * + If there is no such flag, the first element is an array with
                   1379:  *   x-coordinates and the following ones contain y-coordinates.
                   1380:  *
                   1381:  * Additional flags: PLOT_NO_AXE_X, PLOT_NO_AXE_Y, PLOT_NO_FRAME.
                   1382:  */
                   1383:
                   1384: static GEN
                   1385: rectplothrawin(long stringrect, long drawrect, dblPointList *data,
                   1386:                long flags, PARI_plot *WW)
                   1387: {
                   1388:   PARI_plot W;
                   1389:   GEN res;
                   1390:   dblPointList y,x;
                   1391:   double xsml,xbig,ysml,ybig,tmp;
                   1392:   long ltype=0, ltop=avma;
                   1393:   long is,js,i,lm,rm,tm,bm,nc,nbpoints, w[2], wx[2], wy[2];
                   1394:
                   1395:   w[0]=stringrect; w[1]=drawrect;
                   1396:   if (!data) return cgetg(1,t_VEC);
                   1397:   x=data[0]; nc=x.nb; xsml=x.xsml; xbig=x.xbig; ysml=x.ysml; ybig=x.ybig;
                   1398:   if (xbig-xsml < 1.e-9)
                   1399:   {
                   1400:     tmp=fabs(xsml)/10; if (!tmp) tmp=0.1;
                   1401:     xbig+=tmp; xsml-=tmp;
                   1402:   }
                   1403:   if (ybig-ysml < 1.e-9)
                   1404:   {
                   1405:     tmp=fabs(ysml)/10; if (!tmp) tmp=0.1;
                   1406:     ybig+=tmp; ysml-=tmp;
                   1407:   }
                   1408:
                   1409:   if (WW) /* no rectwindow has been supplied ==> ps or screen output */
                   1410:   {
                   1411:     W = *WW;
                   1412:     lm = W.fwidth*10; /* left margin   */
                   1413:     rm = W.hunit-1; /* right margin  */
                   1414:     tm = W.vunit-1; /* top margin    */
                   1415:     bm = W.vunit+W.fheight-1; /* bottom margin */
                   1416:
                   1417:     is = W.width - (lm+rm); js = W.height - (tm+bm);
                   1418:     wx[0]=wy[0]=0; wx[1]=lm; wy[1]=tm;
                   1419:   /*
                   1420:    * Window size (W.width x W.height) is given in pixels, and
                   1421:    * correct pixels are 0..w_width-1.
                   1422:    *
                   1423:    * On the other hand, rect functions work with windows whose pixel
                   1424:    * range is [0,width].
                   1425:    */
                   1426:
                   1427:     initrect(stringrect, W.width-1, W.height-1);
                   1428:     if (drawrect != stringrect) initrect(drawrect, is-1, js-1);
                   1429:   }
                   1430:
                   1431:   if (!(flags & PLOT_NO_RESCALE))
                   1432:     rectscale0(drawrect, xsml, xbig, ysml, ybig);
                   1433:
                   1434:   if (!(flags & PLOT_NO_FRAME))
                   1435:   {
                   1436:     rectlinetype(drawrect, -2);                /* Frame. */
                   1437:     current_color[drawrect]=BLACK;
                   1438:     rectmove0(drawrect,xsml,ysml,0);
                   1439:     rectbox0(drawrect,xbig,ybig,0);
                   1440:   }
                   1441:
                   1442:   if (!(flags & PLOT_NO_AXE_Y) && (xsml<=0 && xbig >=0))
                   1443:   {
                   1444:     rectlinetype(drawrect, -1);                /* Axes. */
                   1445:     current_color[drawrect]=BLUE;
                   1446:     rectmove0(drawrect,0.0,ysml,0);
                   1447:     rectline0(drawrect,0.0,ybig,0);
                   1448:   }
                   1449:
                   1450:   if (!(flags & PLOT_NO_AXE_X) && (ysml<=0 && ybig >=0))
                   1451:   {
                   1452:     rectlinetype(drawrect, -1);                /* Axes. */
                   1453:     current_color[drawrect]=BLUE;
                   1454:     rectmove0(drawrect,xsml,0.0,0);
                   1455:     rectline0(drawrect,xbig,0.0,0);
                   1456:   }
                   1457:
                   1458:   if (flags & PLOT_PARAMETRIC) i=0; else i=1;
                   1459:   current_color[drawrect]=RED;
                   1460:   for (; ltype < nc; )
                   1461:   {
                   1462:     if (nc>1)
                   1463:     {
                   1464:       if (ltype & 1) current_color[drawrect]=RED;
                   1465:       else current_color[drawrect]=SIENNA;
                   1466:     }
                   1467:     if (flags & PLOT_PARAMETRIC) x=data[i++];
                   1468:
                   1469:     y=data[i++]; nbpoints=y.nb;
                   1470:     if ((flags & PLOT_POINTS_LINES) || (flags & PLOT_POINTS)) {
                   1471:        rectlinetype(drawrect, rectpoint_itype + ltype);        /* Graphs. */
                   1472:        rectpointtype(drawrect, rectpoint_itype + ltype);       /* Graphs. */
                   1473:        rectpoints0(drawrect,x.d,y.d,nbpoints);
                   1474:     }
                   1475:     if ((flags & PLOT_POINTS_LINES) || !(flags & PLOT_POINTS)) {
                   1476:        if (flags & PLOT_SPLINES) {
                   1477:            /* rectsplines will call us back with ltype == 0 */
                   1478:            int old = rectline_itype;
                   1479:
                   1480:            rectline_itype = rectline_itype + ltype;
                   1481:            rectsplines(drawrect,x.d,y.d,nbpoints,flags);
                   1482:            rectline_itype = old;
                   1483:        } else {
                   1484:            rectlinetype(drawrect, rectline_itype + ltype);     /* Graphs. */
                   1485:            rectlines0(drawrect,x.d,y.d,nbpoints,0);
                   1486:        }
                   1487:     }
                   1488:     ltype++;                           /* Graphs. */
                   1489:   }
                   1490:   for (i--; i>=0; i--) free(data[i].d);
                   1491:   free(data);
                   1492:
                   1493:   if (WW)
                   1494:   {
                   1495:     char c1[16],c2[16],c3[16],c4[16];
                   1496:
                   1497:     sprintf(c1,"%9.3f",ybig); sprintf(c2,"%9.3f",ysml);
                   1498:     sprintf(c3,"%9.3f",xsml); sprintf(c4,"%9.3f",xbig);
                   1499:
                   1500:     rectlinetype(stringrect,-2); /* Frame */
                   1501:     current_color[stringrect]=BLACK;
                   1502:     put_string(stringrect, 0, W.fheight - 1, c1);
                   1503:     put_string(stringrect, 0, W.height - (bm+ 2 * W.vunit), c2);
                   1504:     put_string(stringrect, lm-(W.fwidth*2), W.height-bm+W.fheight-1, c3);
                   1505:     put_string(stringrect, W.width-(W.fwidth*10), W.height-bm+W.fheight-1,c4);
                   1506:
                   1507:     if (flags & PLOT_POSTSCRIPT)
                   1508:       postdraw0(w,wx,wy,2);
                   1509:     else
                   1510:       rectdraw0(w,wx,wy,2, 0);
                   1511:
                   1512:     killrect(drawrect); if (stringrect != drawrect) killrect(stringrect);
                   1513:   }
                   1514:
                   1515:   avma=ltop;
                   1516:   res = cgetg(5,t_VEC);
                   1517:   res[1]=(long)dbltor(xsml); res[2]=(long)dbltor(xbig);
                   1518:   res[3]=(long)dbltor(ysml); res[4]=(long)dbltor(ybig);
                   1519:   return res;
                   1520: }
                   1521:
                   1522: /*************************************************************************/
                   1523: /*                                                                       */
                   1524: /*                          HI-RES FUNCTIONS                             */
                   1525: /*                                                                       */
                   1526: /*************************************************************************/
                   1527:
                   1528: GEN
                   1529: rectploth(long drawrect,entree *ep,GEN a,GEN b,char *ch,
                   1530:           long prec,ulong flags,long testpoints)
                   1531: {
                   1532:   dblPointList *pl=rectplothin(ep, a,b, ch, prec, flags,testpoints);
                   1533:   return rectplothrawin(0,drawrect, pl, flags,NULL);
                   1534: }
                   1535:
                   1536: GEN
                   1537: rectplothraw(long drawrect, GEN data, long flags)
                   1538: {
                   1539:   dblPointList *pl=gtodblList(data,flags);
                   1540:   return rectplothrawin(0,drawrect,pl,flags,NULL);
                   1541: }
                   1542:
                   1543: static PARI_plot*
                   1544: init_output(long flags)
                   1545: {
                   1546:   if (flags & PLOT_POSTSCRIPT)
                   1547:     { PARI_get_psplot(); return &pari_psplot; }
                   1548:   else
                   1549:     { PARI_get_plot(0); return &pari_plot; }
                   1550: }
                   1551:
                   1552: static GEN
                   1553: ploth0(long stringrect,long drawrect,entree *ep,GEN a,GEN b,char *ch,
                   1554:              long prec,ulong flags,long testpoints)
                   1555: {
                   1556:   PARI_plot *output = init_output(flags);
                   1557:   dblPointList *pl=rectplothin(ep, a,b, ch, prec, flags,testpoints);
                   1558:   return rectplothrawin(stringrect,drawrect, pl, flags, output);
                   1559: }
                   1560:
                   1561: static GEN
                   1562: plothraw0(long stringrect, long drawrect, GEN listx, GEN listy, long flags)
                   1563: {
                   1564:   PARI_plot *output = init_output(flags);
                   1565:   long data[] = {evaltyp(t_VEC) | m_evallg(3), 0, 0};
                   1566:   dblPointList *pl;
                   1567:
                   1568:   data[1] = (long) listx;
                   1569:   data[2] = (long) listy;
                   1570:   pl=gtodblList(data,PLOT_PARAMETRIC);
                   1571:   if (!pl) return cgetg(1,t_VEC);
                   1572:   return rectplothrawin(stringrect,drawrect,pl,flags | PLOT_PARAMETRIC,output);
                   1573: }
                   1574:
                   1575: GEN
                   1576: plothraw(GEN listx, GEN listy, long flags)
                   1577: {
                   1578:   flags=(flags)? 0: PLOT_POINTS;
                   1579:   return plothraw0(STRINGRECT, DRAWRECT, listx, listy, flags);
                   1580: }
                   1581:
                   1582: GEN
                   1583: ploth(entree *ep, GEN a, GEN b, char *ch, long prec,long flags,long numpoints)
                   1584: {
                   1585:   return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,flags,numpoints);
                   1586: }
                   1587:
                   1588: GEN
                   1589: ploth2(entree *ep, GEN a, GEN b, char *ch, long prec)
                   1590: {
                   1591:   return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,PLOT_PARAMETRIC,0);
                   1592: }
                   1593:
                   1594: GEN
                   1595: plothmult(entree *ep, GEN a, GEN b, char *ch, long prec)
                   1596: {
                   1597:   return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,0,0);
                   1598: }
                   1599:
                   1600: GEN
                   1601: postplothraw(GEN listx, GEN listy, long flags)
                   1602: {
                   1603:   flags=(flags)? 0: PLOT_POINTS;
                   1604:   return plothraw0(STRINGRECT, DRAWRECT, listx, listy, flags|PLOT_POSTSCRIPT);
                   1605: }
                   1606:
                   1607: GEN
                   1608: postploth(entree *ep, GEN a, GEN b, char *ch, long prec,long flags,
                   1609:            long numpoints)
                   1610: {
                   1611:   return ploth0(STRINGRECT,DRAWRECT,ep,a,b,ch,prec,flags|PLOT_POSTSCRIPT,
                   1612:                 numpoints);
                   1613: }
                   1614:
                   1615: GEN
                   1616: postploth2(entree *ep, GEN a, GEN b, char *ch, long prec,
                   1617:            long numpoints)
                   1618: {
                   1619:   return ploth0(STRINGRECT,DRAWRECT,ep,a,b,ch,prec,
                   1620:                PLOT_PARAMETRIC|PLOT_POSTSCRIPT,numpoints);
                   1621: }
                   1622:
                   1623: GEN
                   1624: plothsizes()
                   1625: {
                   1626:   GEN vect = cgetg(1+6,t_VEC);
                   1627:   int i;
                   1628:
                   1629:   PARI_get_plot(0);
                   1630:   for (i=1; i<7; i++) vect[i]=lgeti(3);
                   1631:   affsi(w_width,(GEN)vect[1]); affsi(w_height,(GEN)vect[2]);
                   1632:   affsi(h_unit, (GEN)vect[3]); affsi(v_unit,  (GEN)vect[4]);
                   1633:   affsi(f_width,(GEN)vect[5]); affsi(f_height,(GEN)vect[6]);
                   1634:   return vect;
                   1635: }
                   1636:
                   1637: /*************************************************************************/
                   1638: /*                                                                       */
                   1639: /*                         POSTSCRIPT OUTPUT                             */
                   1640: /*                                                                       */
                   1641: /*************************************************************************/
                   1642:
                   1643: typedef struct spoint {
                   1644:   int x,y; } SPoint;
                   1645: typedef struct ssegment {
                   1646:   int x1,y1,x2,y2; } SSegment;
                   1647: typedef struct srectangle {
                   1648:   int x,y,width,height; } SRectangle;
                   1649:
                   1650: static void ps_point(FILE *psfile, int x, int y);
                   1651: static void ps_line(FILE *psfile, int x1, int y1, int x2, int y2);
                   1652: static void ps_rect(FILE *psfile, int x1, int y1, int x2, int y2);
                   1653: static void ps_string(FILE *psfile, int x, int y, char *c);
                   1654:
                   1655: #undef ISCR
                   1656: #undef JSCR
                   1657: #define ISCR 1120 /* 1400 en haute resolution */
                   1658: #define JSCR 800  /* 1120 en haute resolution */
                   1659:
                   1660: static void
                   1661: PARI_get_psplot()
                   1662: {
                   1663:   pari_psplot.height = JSCR - 60;
                   1664:   pari_psplot.width  = ISCR - 40;
                   1665:   pari_psplot.fheight = 15;
                   1666:   pari_psplot.fwidth = 6;
                   1667:   pari_psplot.hunit = 5;
                   1668:   pari_psplot.vunit = 5;
                   1669: }
                   1670:
                   1671: static void
                   1672: gendraw(GEN list,long ps)
                   1673: {
                   1674:   long i,n,ne,*w,*x,*y;
                   1675:   GEN x0,y0,win;
                   1676:
                   1677:   if (typ(list) != t_VEC) err(talker,"not a vector in rectdraw");
                   1678:   n = lg(list)-1;
                   1679:   if (n%3) err(talker,"incorrect number of components in rectdraw");
                   1680:   n = n/3; if (!n) return;
                   1681:   w = (long*)gpmalloc(n*sizeof(long));
                   1682:   x = (long*)gpmalloc(n*sizeof(long));
                   1683:   y = (long*)gpmalloc(n*sizeof(long));
                   1684:   for (i=0; i<n; i++)
                   1685:   {
                   1686:     win=(GEN)list[3*i+1]; x0=(GEN)list[3*i+2]; y0=(GEN)list[3*i+3];
                   1687:     if (typ(win)!=t_INT || typ(x0)!=t_INT || typ(y0)!= t_INT)
                   1688:       err(talker, "not an integer type in rectdraw");
                   1689:     ne=itos(win); check_rect(ne);
                   1690:     x[i]=itos(x0); y[i]=itos(y0); w[i]=ne;
                   1691:   }
                   1692:   if (ps) postdraw0(w,x,y,n); else rectdraw0(w,x,y,n, 1);
                   1693:   free(x); free(y); free(w);
                   1694: }
                   1695:
                   1696: void
                   1697: postdraw(GEN list) { gendraw(list,1); }
                   1698:
                   1699: void
                   1700: rectdraw(GEN list) { gendraw(list,0); }
                   1701:
                   1702: static char*
                   1703: zmalloc(size_t x)
                   1704: {
                   1705:   return x? gpmalloc(x): NULL;
                   1706: }
                   1707:
                   1708: void
                   1709: postdraw0(long *w, long *x, long *y, long lw)
                   1710: {
                   1711:   long *ptx,*pty,*numpoints,*numtexts,*xtexts,*ytexts;
                   1712:   RectObj *p1;
                   1713:   PariRect *e;
                   1714:   long i,j,x0,y0;
                   1715:   long a,b,c,d,nd[ROt_MAX+1];
                   1716:   char **texts;
                   1717:   FILE *psfile;
                   1718:
                   1719:   SPoint *points, **lines, *SLine;
                   1720:   SSegment *seg;
                   1721:   SRectangle *rect, SRec;
                   1722:
                   1723:   psfile = fopen(current_psfile, "a");
                   1724:   if (!psfile)
                   1725:     err(openfiler,"postscript",current_psfile);
                   1726:
                   1727:   for (i=0; i<=ROt_MAX; i++) nd[i]=0;
                   1728:
                   1729:   for (i=0; i<lw; i++)
                   1730:   {
                   1731:     e=rectgraph[w[i]]; p1=RHead(e);
                   1732:     while (p1)
                   1733:     {
                   1734:       if (RoType(p1) != ROt_MP) nd[RoType(p1)]++;
                   1735:       else nd[ROt_PT]+=RoMPcnt(p1);
                   1736:       p1=RoNext(p1);
                   1737:     }
                   1738:   }
                   1739:   points=(SPoint*) zmalloc(nd[ROt_PT]*sizeof(SPoint));
                   1740:   seg=(SSegment*) zmalloc(nd[ROt_LN]*sizeof(SSegment));
                   1741:   rect=(SRectangle*) zmalloc(nd[ROt_BX]*sizeof(SRectangle));
                   1742:   lines=(SPoint**) zmalloc(nd[ROt_ML]*sizeof(SPoint*));
                   1743:   numpoints=(long*) zmalloc(nd[ROt_ML]*sizeof(long));
                   1744:   texts=(char**) zmalloc(nd[ROt_ST]*sizeof(char*));
                   1745:   numtexts=(long*) zmalloc(nd[ROt_ST]*sizeof(long));
                   1746:   xtexts = (long*) zmalloc(nd[ROt_ST]*sizeof(long));
                   1747:   ytexts = (long*) zmalloc(nd[ROt_ST]*sizeof(long));
                   1748:   for (i=0; i<=ROt_MAX; i++) nd[i]=0;
                   1749:
                   1750:   for (i=0; i<lw; i++)
                   1751:   {
                   1752:     e=rectgraph[w[i]]; p1=RHead(e); x0=x[i]; y0=y[i];
                   1753:     while (p1)
                   1754:     {
                   1755:       switch(RoType(p1))
                   1756:       {
                   1757:        case ROt_PT:
                   1758:          points[nd[ROt_PT]].x=RoPTx(p1)+x0;
                   1759:          points[nd[ROt_PT]].y=RoPTy(p1)+y0;
                   1760:          nd[ROt_PT]++; break;
                   1761:        case ROt_LN:
                   1762:          seg[nd[ROt_LN]].x1=RoLNx1(p1)+x0;
                   1763:          seg[nd[ROt_LN]].y1=RoLNy1(p1)+y0;
                   1764:          seg[nd[ROt_LN]].x2=RoLNx2(p1)+x0;
                   1765:          seg[nd[ROt_LN]].y2=RoLNy2(p1)+y0;
                   1766:          nd[ROt_LN]++; break;
                   1767:        case ROt_BX:
                   1768:          a=rect[nd[ROt_BX]].x=RoBXx1(p1)+x0;
                   1769:          b=rect[nd[ROt_BX]].y=RoBXy1(p1)+y0;
                   1770:          rect[nd[ROt_BX]].width =RoBXx2(p1)+x0-a;
                   1771:          rect[nd[ROt_BX]].height=RoBXy2(p1)+y0-b;
                   1772:          nd[ROt_BX]++; break;
                   1773:        case ROt_MP:
                   1774:          ptx=RoMPxs(p1); pty=RoMPys(p1);
                   1775:          for (j=0; j<RoMPcnt(p1); j++)
                   1776:          {
                   1777:            points[nd[ROt_PT]+j].x=ptx[j]+x0;
                   1778:            points[nd[ROt_PT]+j].y=pty[j]+y0;
                   1779:          }
                   1780:          nd[ROt_PT]+=RoMPcnt(p1); break;
                   1781:        case ROt_ML:
                   1782:          ptx=RoMLxs(p1); pty=RoMLys(p1);
                   1783:          numpoints[nd[ROt_ML]]=RoMLcnt(p1);
                   1784:          lines[nd[ROt_ML]]=(SPoint*)gpmalloc(RoMLcnt(p1)*sizeof(SPoint));
                   1785:          for (j=0; j<RoMLcnt(p1); j++)
                   1786:          {
                   1787:            lines[nd[ROt_ML]][j].x=ptx[j]+x0;
                   1788:            lines[nd[ROt_ML]][j].y=pty[j]+y0;
                   1789:          }
                   1790:          nd[ROt_ML]++; break;
                   1791:        case ROt_ST:
                   1792:          texts[nd[ROt_ST]]=RoSTs(p1);
                   1793:          numtexts[nd[ROt_ST]]=RoSTl(p1);
                   1794:          xtexts[nd[ROt_ST]]=RoSTx(p1)+x0;
                   1795:          ytexts[nd[ROt_ST]]=RoSTy(p1)+y0;
                   1796:          nd[ROt_ST]++; break;
                   1797:        default: break;
                   1798:       }
                   1799:       p1=RoNext(p1);
                   1800:     }
                   1801:   }
                   1802:   fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont 16 scalefont setfont\n0.65 0.65 scale\n");
                   1803:   for (i=0; i<nd[ROt_PT]; i++)
                   1804:     ps_point(psfile,points[i].x,points[i].y);
                   1805:   for (i=0; i<nd[ROt_LN]; i++)
                   1806:     ps_line(psfile,seg[i].x1,seg[i].y1,seg[i].x2,seg[i].y2);
                   1807:   for (i=0; i<nd[ROt_BX]; i++)
                   1808:   {
                   1809:     SRec=rect[i]; a=SRec.x; b=SRec.y; c=a+SRec.width;
                   1810:     d=b+SRec.height; ps_rect(psfile,a,b,c,d);
                   1811:   }
                   1812:   for (i=0; i<nd[ROt_ML]; i++)
                   1813:   {
                   1814:     SLine=lines[i];
                   1815:     for (j=0; j<numpoints[i]; j++)
                   1816:     {
                   1817:       if (!j) fprintf(psfile,"%d %d moveto\n",SLine[0].y,SLine[0].x);
                   1818:       else fprintf(psfile,"%d %d lineto\n",SLine[j].y,SLine[j].x);
                   1819:     }
                   1820:   }
                   1821:   for (i=0; i<nd[ROt_ST]; i++)
                   1822:     ps_string(psfile,xtexts[i],ytexts[i],texts[i]);
                   1823:   fprintf(psfile,"stroke showpage\n"); fclose(psfile);
                   1824: #define xfree(pointer) if (pointer) free(pointer)
                   1825:   xfree(points); xfree(seg); xfree(rect); xfree(numpoints);
                   1826:   for (i=0; i<nd[ROt_ML]; i++) xfree(lines[i]);
                   1827:   xfree(lines); xfree(texts); xfree(numtexts); xfree(xtexts); xfree(ytexts);
                   1828: #undef xfree
                   1829: }
                   1830:
                   1831: static void
                   1832: ps_point(FILE *psfile, int x, int y)
                   1833: {
                   1834:   fprintf(psfile,"%d %d moveto\n0 2 rlineto 2 0 rlineto 0 -2 rlineto closepath fill\n",y,x);
                   1835: }
                   1836:
                   1837: static void
                   1838: ps_line(FILE *psfile, int x1, int y1, int x2, int y2)
                   1839: {
                   1840:   fprintf(psfile,"%d %d moveto\n%d %d lineto\n",y1,x1,y2,x2);
                   1841: }
                   1842:
                   1843: static void
                   1844: ps_rect(FILE *psfile, int x1, int y1, int x2, int y2)
                   1845: {
                   1846:   fprintf(psfile,"%d %d moveto\n%d %d lineto\n%d %d lineto\n%d %d lineto\nclosepath\n",y1,x1,y1,x2,y2,x2,y2,x1);
                   1847: }
                   1848:
                   1849: static void
                   1850: ps_string(FILE *psfile, int x, int y, char *c)
                   1851: {
                   1852:   fprintf(psfile,"%d %d moveto 90 rotate\n(%s) show -90 rotate\n",y,x,c);
                   1853: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>