xref: /dports/math/pari/pari-2.13.3/src/graph/plotport.c

/* Copyright (C) 2000  The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */

/*******************************************************************/
/*                                                                 */
/*                         PLOT ROUTINES                           */
/*                                                                 */
/*******************************************************************/
#include "pari.h"
#include "paripriv.h"
#include "rect.h"

static void (*pari_get_plot)(PARI_plot *);

/* no need for THREAD: OK to share this */
static hashtable *rgb_colors = NULL;

THREAD PariRect rectgraph[18]; /*NUMRECT*/
static THREAD long current_color[18]; /*NUMRECT*/

static long plotpoint_itype = 0, rectline_itype  = 0;

const long NUMRECT = 18;
const long RECUR_MAXDEPTH = 10;
const double RECUR_PREC = 0.001;
const long DEFAULT_COLOR = 1, AXIS_COLOR = 2;

enum {
  ROt_MV = 1, /* Move */
  ROt_PT,  /* Point */
  ROt_LN,  /* Line */
  ROt_BX,  /* Box */
  ROt_FBX, /* Filled Box */
  ROt_MP,  /* Multiple point */
  ROt_ML,  /* Multiple lines */
  ROt_ST,  /* String */
  ROt_PTT,  /* Point type change */
  ROt_LNT,  /* Line type change */
  ROt_PTS,  /* Point size change */
  ROt_NULL, /* To be the start of the chain */
};

/* string justification */
#define RoSTdirLEFT       0x00
#define RoSTdirRIGHT      0x02
#define RoSTdirHPOS_mask  0x03
#define RoSTdirBOTTOM     0x00
#define RoSTdirTOP        0x08
#define RoSTdirVPOS_mask  0x0c
#define RoSTdirHGAP       0x10
#define RoSTdirVGAP       0x20

/* ploth flags */
#define PLOT_PARAMETRIC   0x00001
#define PLOT_RECURSIVE    0x00002
#define PLOT_NO_RESCALE   0x00004
#define PLOT_NO_AXE_X     0x00008
#define PLOT_NO_AXE_Y     0x00010
#define PLOT_NO_FRAME     0x00020
#define PLOT_POINTS       0x00040
#define PLOT_POINTS_LINES 0x00080
#define PLOT_SPLINES      0x00100
#define PLOT_NO_TICK_X    0x00200
#define PLOT_NO_TICK_Y    0x00400
#define PLOT_NODOUBLETICK 0x00800
#define PLOT_COMPLEX      0x01000
#define PLOT_PARA         0x02000

INLINE long
DTOL(double t) { return (long)(t + 0.5); }

static const long PS_WIDTH = 1120 - 60; /* 1400 - 60 for hi-res */
static const long PS_HEIGH = 800 - 40; /* 1120 - 60 for hi-res */
static const long PS_SCALE = 1000; /* Allowing 64x zoom on 500ppi */

static void
_psdraw_scale(PARI_plot *T, GEN w, GEN x, GEN y)
{
  pari_sp av = avma;
  FILE *F = fopen(current_psfile, "a");
  if (!F) pari_err_FILE("postscript file",current_psfile);
  fputs(rect2ps(w,x,y,T), F);
  fclose(F); set_avma(av);
}
static void
_psdraw(PARI_plot *T, GEN w, GEN x, GEN y)
{ (void)T; _psdraw_scale(NULL,w,x,y); }
static void
pari_get_psplot(PARI_plot *T)
{
  T->width  = PS_WIDTH;
  T->height = PS_HEIGH;
  T->fheight= 15;
  T->fwidth = 6;
  T->hunit  = 5;
  T->vunit  = 5;
  T->dwidth = 0;
  T->dheight= 0;
  T->draw = NULL;
}
static void
pari_get_svgplot(PARI_plot *T)
{
  T->width   = 480;
  T->height  = 320;
  T->fheight = 12;
  T->fwidth  = 6;
  T->hunit   = 3;
  T->vunit   = 3;
  T->dwidth  = 0;
  T->dheight = 0;
  T->draw = NULL;
}

/********************************************************************/
/**                                                                **/
/**                      RECTPLOT FUNCTIONS                        **/
/**                                                                **/
/********************************************************************/
void
pari_init_graphics(void) { pari_get_plot = &pari_get_svgplot; }

void
pari_set_plot_engine(void (*plot)(PARI_plot *))
{
  long n;
  pari_get_plot = plot;
  for (n = 0; n < NUMRECT; n++)
  {
    PariRect *e = &rectgraph[n];
    RHead(e) = RTail(e) = NULL;
    RXsize(e) = RYsize(e) = 0;
  }
}

void
pari_kill_plot_engine(void)
{
  int i;
  for (i=0; i<NUMRECT; i++)
  {
    PariRect *e = &rectgraph[i];
    if (RHead(e)) plotkill(i);
  }
  if (rgb_colors)
  {
    pari_free((void*)rgb_colors->table);
    pari_free((void*)rgb_colors);
  }
}

static PariRect *
check_rect(long ne)
{
  const char *f = "graphic function";
  const long m = NUMRECT-1;
  if (ne < 0) pari_err_DOMAIN(f, "rectwindow", "<", gen_0, stoi(ne));
  if (ne > m) pari_err_DOMAIN(f, "rectwindow", ">", stoi(m), stoi(ne));
  return &rectgraph[ne];
}

static PariRect *
check_rect_init(long ne)
{
  PariRect *e = check_rect(ne);
  if (!RHead(e)) pari_err_TYPE("graphic function [use plotinit()]", stoi(ne));
  return e;
}
static void
Rchain(PariRect *e, RectObj *z)
{
  if (!RHead(e)) RHead(e) = z; else RoNext(RTail(e)) = z;
  RTail(e) = z;
  RoNext(z) = NULL;
}

static long
rgb_to_long(long r, long g, long b)
{ return (r << 16) | (g << 8) | b; }
/* c from graphcolormap */
static long
colormap_to_color(long i)
{
  GEN c = GP_DATA->colormap;
  long k = i+1, l = lg(c)-1;
  int r, g, b;
  if (k > l)
    pari_err_COMPONENT("graphcolormap",">", stoi(l), stoi(k));
  color_to_rgb(gel(c, k), &r,&g,&b);
  return rgb_to_long(r, g, b);
}

static void
initrect_i(long ne, long x, long y)
{
  PariRect *e;
  RectObj *z;

  if (x <= 1) pari_err_DOMAIN("plotinit", "x", "<=", gen_1, stoi(x));
  if (y <= 1) pari_err_DOMAIN("plotinit", "y", "<=", gen_1, stoi(y));
  e = check_rect(ne); if (RHead(e)) plotkill(ne);

  current_color[ne] = colormap_to_color(DEFAULT_COLOR);
  z = (RectObj*) pari_malloc(sizeof(RectObj));
  RoType(z) = ROt_NULL;
  Rchain(e, z);
  RXsize(e) = x; RXcursor(e) = 0; RXscale(e) = 1; RXshift(e) = 0;
  RYsize(e) = y; RYcursor(e) = 0; RYscale(e) = 1; RYshift(e) = 0;
}
static long
initrect_get_arg(GEN x, long dft)
{
  if (!x) return dft;
  if (typ(x) != t_INT) pari_err_TYPE("plotinit",x);
  return itos(x);
}
void
plotinit(long ne, GEN x, GEN y, long flag)
{
  const long m = NUMRECT-3;
  long xi, yi;
  PARI_plot T;

  if (flag)
  {
    pari_get_plot(&T);
    xi = T.width -1; if (x) xi = DTOL(xi * gtodouble(x));
    yi = T.height-1; if (y) yi = DTOL(yi * gtodouble(y));
  }
  else
  {
    if (!x || !y) pari_get_plot(&T);
    xi = initrect_get_arg(x, T.width -1);
    yi = initrect_get_arg(y, T.height-1);
  }
  if (ne > m) pari_err_DOMAIN("plotinit", "rectwindow", ">", stoi(m), stoi(ne));
  initrect_i(ne, xi, yi);
}

GEN
plotcursor(long ne)
{
  PariRect *e = check_rect_init(ne);
  return mkvec2s((long)RXcursor(e), (long)RYcursor(e));
}

static void
plotscale0(long ne, double x1, double x2, double y1, double y2)
{
  PariRect *e = check_rect_init(ne);
  double x, y;

  x = RXshift(e) + RXscale(e) * RXcursor(e);
  y = RYshift(e) + RYscale(e) * RYcursor(e);
  RXscale(e) = RXsize(e)/(x2-x1); RXshift(e) = -x1*RXscale(e);
  RYscale(e) = RYsize(e)/(y1-y2); RYshift(e) = -y2*RYscale(e);
  RXcursor(e) = (x - RXshift(e)) / RXscale(e);
  RYcursor(e) = (y - RYshift(e)) / RYscale(e);
}
void
plotscale(long ne, GEN x1, GEN x2, GEN y1, GEN y2)
{ plotscale0(ne, gtodouble(x1), gtodouble(x2), gtodouble(y1), gtodouble(y2)); }

static void
plotmove0(long ne, double x, double y, long relative)
{
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj1P));

  if (relative) { RXcursor(e) += x; RYcursor(e) += y; }
  else          { RXcursor(e) = x; RYcursor(e) = y; }
  RoType(z) = ROt_MV;
  RoMVx(z) = RXcursor(e) * RXscale(e) + RXshift(e);
  RoMVy(z) = RYcursor(e) * RYscale(e) + RYshift(e);
  Rchain(e, z);
}
static void
_move(long ne, double x, double y)
{ plotmove0(ne,x,y,0); }
void
plotmove(long ne, GEN x, GEN y)
{ plotmove0(ne,gtodouble(x),gtodouble(y),0); }
void
plotrmove(long ne, GEN x, GEN y)
{ plotmove0(ne,gtodouble(x),gtodouble(y),1); }

/* ROt_MV/ROt_PT */
static void
plotpoint0(long ne, double x, double y,long relative)
{
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj1P));

  if (relative) { RXcursor(e) += x; RYcursor(e) += y; }
  else          { RXcursor(e) = x; RYcursor(e) = y; }
  RoPTx(z) = RXcursor(e)*RXscale(e) + RXshift(e);
  RoPTy(z) = RYcursor(e)*RYscale(e) + RYshift(e);
  RoType(z) = ( DTOL(RoPTx(z)) < 0
                || DTOL(RoPTy(z)) < 0 || DTOL(RoPTx(z)) > RXsize(e)
                || DTOL(RoPTy(z)) > RYsize(e) ) ? ROt_MV : ROt_PT;
  Rchain(e, z);
  RoCol(z) = current_color[ne];
}
static void
plotpoint(long ne, GEN x, GEN y)
{ plotpoint0(ne,gtodouble(x),gtodouble(y),0); }
void
plotrpoint(long ne, GEN x, GEN y)
{ plotpoint0(ne,gtodouble(x),gtodouble(y),1); }

GEN
plotcolor(long ne, GEN c)
{
  long t = typ(c), n = lg(GP_DATA->colormap)-2;
  int r, g, b;
  check_rect(ne);
  if (t == t_INT)
  {
    long i = itos(c);
    if (i < 0) pari_err_DOMAIN("plotcolor", "color", "<", gen_0, c);
    if (i > n) pari_err_DOMAIN("plotcolor", "color", ">", stoi(n), c);
    c = gel(GP_DATA->colormap,i+1);
  }
  else
  {
    if (t == t_VEC) { c = ZV_to_zv(c); t = typ(c); }
    if (t != t_VECSMALL && t != t_STR) pari_err_TYPE("plotcolor",c);
  }
  color_to_rgb(c, &r,&g,&b);
  current_color[ne] = rgb_to_long(r,g,b);
  return mkvec3s(r, g, b);
}

/* ROt_MV/ROt_LN */
static void
rectline0(long ne, double gx2, double gy2, long relative)
{
  double dx, dy, dxy, xmin, xmax, ymin, ymax, x1, y1, x2, y2;
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj2P));
  const double c = 1 + 1e-10;

  x1 = RXcursor(e)*RXscale(e) + RXshift(e);
  y1 = RYcursor(e)*RYscale(e) + RYshift(e);
  if (relative)
    { RXcursor(e)+=gx2; RYcursor(e)+=gy2; }
  else
    { RXcursor(e)=gx2; RYcursor(e)=gy2; }
  x2 = RXcursor(e)*RXscale(e) + RXshift(e);
  y2 = RYcursor(e)*RYscale(e) + RYshift(e);
  xmin = maxdd(mindd(x1,x2),0); xmax = mindd(maxdd(x1,x2),RXsize(e));
  ymin = maxdd(mindd(y1,y2),0); ymax = mindd(maxdd(y1,y2),RYsize(e));
  dxy = x1*y2 - y1*x2; dx = x2-x1; dy = y2-y1;
  if (dy)
  {
    double a = (dxy + RYsize(e)*dx) / dy, b = dxy / dy;
    if (dx*dy < 0)
    { xmin=maxdd(xmin,a); xmax=mindd(xmax,b); }
    else
    { xmin=maxdd(xmin,b); xmax=mindd(xmax,a); }
  }
  if (dx)
  {
    double a = (RXsize(e)*dy - dxy) / dx, b = -dxy / dx;
    if (dx*dy < 0)
    { ymin=maxdd(ymin,a); ymax=mindd(ymax,b); }
    else
    { ymin=maxdd(ymin,b); ymax=mindd(ymax,a); }
  }
  RoLNx1(z) = xmin;
  RoLNx2(z) = xmax;
  if (dx*dy < 0) { RoLNy1(z) = ymax; RoLNy2(z) = ymin; }
  else         { RoLNy1(z) = ymin; RoLNy2(z) = ymax; }
  RoType(z) = (xmin>xmax*c || ymin>ymax*c) ? ROt_MV : ROt_LN;
  Rchain(e, z);
  RoCol(z) = current_color[ne];
}
static void
_line(long ne, double x, double y)
{ rectline0(ne, x, y, 0); }
void
plotline(long ne, GEN gx2, GEN gy2)
{ rectline0(ne, gtodouble(gx2), gtodouble(gy2),0); }
void
plotrline(long ne, GEN gx2, GEN gy2)
{ rectline0(ne, gtodouble(gx2), gtodouble(gy2),1); }

enum {
  TICKS_CLOCKW   = 1, /* Draw in clockwise direction */
  TICKS_ACLOCKW  = 2, /* Draw in anticlockwise direction */
  TICKS_ENDSTOO  = 4, /* Draw at endspoints if needed */
  TICKS_NODOUBLE = 8  /* Do not draw double-length ticks */
};

/* Given coordinates of ends of a line, and labels l1 l2 attached to the
 * ends, plot ticks where the label coordinate takes "round" values */
static void
rectticks(PARI_plot *WW, long ne, double dx1, double dy1, double dx2,
          double dy2, double l1, double l2, long flags)
{
  long dx, dy, dxy, dxy1, x1, y1, x2, y2, nticks, n, n1, dn;
  double minstep, maxstep, step, l_min, l_max, minl, maxl, dl, dtx, dty, x, y;
  double ddx, ddy;
  const double mult[3] = { 2./1., 5./2., 10./5. };
  PariRect *e = check_rect_init(ne);
  int do_double = !(flags & TICKS_NODOUBLE);

  x1 = DTOL(dx1*RXscale(e) + RXshift(e));
  y1 = DTOL(dy1*RYscale(e) + RYshift(e));
  x2 = DTOL(dx2*RXscale(e) + RXshift(e));
  y2 = DTOL(dy2*RYscale(e) + RYshift(e));
  dx = x2 - x1; if (dx < 0) dx = -dx;
  dy = y2 - y1; if (dy < 0) dy = -dy;
  dxy1 = maxss(dx, dy);
  dx /= WW->hunit;
  dy /= WW->vunit;
  if (dx > 1000 || dy > 1000)
    dxy = 1000; /* avoid overflow */
  else
    dxy = usqrt(dx*dx + dy*dy);
  nticks = (long) ((dxy + 2.5)/4);
  if (!nticks) return;

  /* Find nticks (or less) "round" numbers between l1 and l2. For our purpose
   * round numbers have "last significant" decimal digit either
   *    - any;
   *    - even;
   *    - divisible by 5.
   * We need to choose which alternative is better. */
  if (l1 < l2)
    l_min = l1, l_max = l2;
  else
    l_min = l2, l_max = l1;
  minstep = (l_max - l_min)/(nticks + 1);
  maxstep = 2.5*(l_max - l_min);
  step = exp(log(10.) * floor(log10(minstep)));
  if (!(flags & TICKS_ENDSTOO)) {
    double d = 2*(l_max - l_min)/dxy1; /* Two pixels off */
    l_min += d;
    l_max -= d;
  }
  for (n = 0; ; n++)
  {
    if (step >= maxstep) return;

    if (step >= minstep) {
      minl = ceil(l_min/step);
      maxl = floor(l_max/step);
      if (minl <= maxl && maxl - minl + 1 <= nticks) {
        nticks = (long) (maxl - minl + 1);
        l_min = minl * step;
        l_max = maxl * step; break;
      }
    }
    step *= mult[ n % 3 ];
  }
  /* Where to position doubleticks. Variants:
   * small: each 5, double: each 10  ; n=2 mod 3
   * small: each 2, double: each 10  ; n=1 mod 3
   * small: each 1, double: each  5 */
  dn = (n % 3 == 2)? 2: 5;
  n1 = ((long)minl) % dn; /* unused if do_double = FALSE */

  /* now l_min and l_max keep min/max values of l with ticks, and nticks is
     the number of ticks to draw. */
  if (nticks == 1) ddx = ddy = 0; /* -Wall */
  else {
    dl = (l_max - l_min)/(nticks - 1);
    ddx = (dx2 - dx1) * dl / (l2 - l1);
    ddy = (dy2 - dy1) * dl / (l2 - l1);
  }
  x = dx1 + (dx2 - dx1) * (l_min - l1) / (l2 - l1);
  y = dy1 + (dy2 - dy1) * (l_min - l1) / (l2 - l1);
  /* assume hunit and vunit form a square. For clockwise ticks: */
  dtx = WW->hunit * dy/dxy * (y2 > y1 ? 1 : -1); /* y-coord runs down */
  dty = WW->vunit * dx/dxy * (x2 > x1 ? 1 : -1);
  for (n = 0; n < nticks; n++, x += ddx, y += ddy) {
    RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj2P));
    double lunit = WW->hunit > 1 ? 1.5 : 2;
    double l = (do_double && (n + n1) % dn == 0) ? lunit: 1;
    double x1, x2, y1, y2;
    x1 = x2 = x*RXscale(e) + RXshift(e);
    y1 = y2 = y*RYscale(e) + RYshift(e);
    if (flags & TICKS_CLOCKW)  { x1 += dtx*l; y1 -= dty*l; }
    if (flags & TICKS_ACLOCKW) { x2 -= dtx*l; y2 += dty*l; }
    RoLNx1(z) = x1; RoLNy1(z) = y1;
    RoLNx2(z) = x2; RoLNy2(z) = y2;
    RoType(z) = ROt_LN;
    Rchain(e, z);
    RoCol(z) = current_color[ne];
  }
}

static void
rectbox0(long ne, double gx2, double gy2, long relative, long filled)
{
  double xx, yy, x1, y1, x2, y2, xmin, ymin, xmax, ymax;
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj2P));

  x1 = RXcursor(e)*RXscale(e) + RXshift(e);
  y1 = RYcursor(e)*RYscale(e) + RYshift(e);
  if (relative)
  { xx = RXcursor(e)+gx2; yy = RYcursor(e)+gy2; }
  else
  {  xx = gx2; yy = gy2; }
  x2 = xx*RXscale(e) + RXshift(e);
  y2 = yy*RYscale(e) + RYshift(e);
  xmin = maxdd(mindd(x1,x2),0); xmax = mindd(maxdd(x1,x2),RXsize(e));
  ymin = maxdd(mindd(y1,y2),0); ymax = mindd(maxdd(y1,y2),RYsize(e));

  RoType(z) = filled ? ROt_FBX: ROt_BX;
  RoBXx1(z) = xmin; RoBXy1(z) = ymin;
  RoBXx2(z) = xmax; RoBXy2(z) = ymax;
  Rchain(e, z);
  RoCol(z) = current_color[ne];
}
static void
_box(long ne, double x, double y)
{ rectbox0(ne, x, y, 0, 0); }
void
plotbox(long ne, GEN gx2, GEN gy2, long f)
{ rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 0, f); }
void
plotrbox(long ne, GEN gx2, GEN gy2, long f)
{ rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 1, f); }

static void
freeobj(RectObj *z) {
  switch(RoType(z)) {
    case ROt_MP: case ROt_ML:
      pari_free(RoMPxs(z));
      pari_free(RoMPys(z)); break;
    case ROt_ST:
      pari_free(RoSTs(z)); break;
  }
  pari_free(z);
}

void
plotkill(long ne)
{
  RectObj *z, *t;
  PariRect *e = check_rect_init(ne);

  z = RHead(e);
  RHead(e) = RTail(e) = NULL;
  RXsize(e) = RYsize(e) = 0;
  RXcursor(e) = RYcursor(e) = 0;
  RXscale(e) = RYscale(e) = 1;
  RXshift(e) = RYshift(e) = 0;
  while (z) { t = RoNext(z); freeobj(z); z = t; }
}

/* ROt_MP */
static void
plotpoints0(long ne, double *X, double *Y, long lx)
{
  double *px, *py;
  long i, cp=0;
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObjMP));

  RoMPxs(z) = px = (double*) pari_malloc(lx*sizeof(double));
  RoMPys(z) = py = (double*) pari_malloc(lx*sizeof(double));
  for (i=0; i<lx; i++)
  {
    double x = RXscale(e)*X[i] + RXshift(e);
    double y = RYscale(e)*Y[i] + RYshift(e);
    if (x >= 0 && y >= 0 && x <= RXsize(e) && y <= RYsize(e))
    {
      px[cp] = x; py[cp] = y; cp++;
    }
  }
  RoType(z) = ROt_MP;
  RoMPcnt(z) = cp;
  Rchain(e, z);
  RoCol(z) = current_color[ne];
}
void
plotpoints(long ne, GEN X, GEN Y)
{
  pari_sp av = avma;
  double *px, *py;
  long i, lx;

  if (!is_vec_t(typ(X)) || !is_vec_t(typ(Y))) { plotpoint(ne, X, Y); return; }
  lx = lg(X); if (lg(Y) != lx) pari_err_DIM("plotpoints");
  lx--; if (!lx) return;

  px = (double*)stack_malloc_align(lx*sizeof(double), sizeof(double)); X++;
  py = (double*)stack_malloc_align(lx*sizeof(double), sizeof(double)); Y++;
  for (i=0; i<lx; i++)
  {
    px[i] = gtodouble(gel(X,i));
    py[i] = gtodouble(gel(Y,i));
  }
  plotpoints0(ne,px,py,lx); set_avma(av);
}

/* ROt_ML */
static void
rectlines0(long ne, double *x, double *y, long lx, long flag)
{
  long i,I;
  double *ptx,*pty;
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObj2P));

  I = flag ? lx+1 : lx;
  ptx = (double*) pari_malloc(I*sizeof(double));
  pty = (double*) pari_malloc(I*sizeof(double));
  for (i=0; i<lx; i++)
  {
    ptx[i] = RXscale(e)*x[i] + RXshift(e);
    pty[i] = RYscale(e)*y[i] + RYshift(e);
  }
  if (flag)
  {
    ptx[i] = RXscale(e)*x[0] + RXshift(e);
    pty[i] = RYscale(e)*y[0] + RYshift(e);
  }
  Rchain(e, z);
  RoType(z) = ROt_ML;
  RoMLcnt(z) = I;
  RoMLxs(z) = ptx;
  RoMLys(z) = pty;
  RoCol(z) = current_color[ne];
}
void
plotlines(long ne, GEN X, GEN Y, long flag)
{
  pari_sp av = avma;
  double *x, *y;
  long i, lx;

  if (!is_vec_t(typ(X)) || !is_vec_t(typ(Y))) { plotline(ne, X, Y); return; }
  lx = lg(X); if (lg(Y) != lx) pari_err_DIM("plotlines");
  lx--; if (!lx) return;

  x = (double*)stack_malloc_align(lx*sizeof(double), sizeof(double)); X++;
  y = (double*)stack_malloc_align(lx*sizeof(double), sizeof(double)); Y++;
  for (i=0; i<lx; i++)
  {
    x[i] = gtodouble(gel(X,i));
    y[i] = gtodouble(gel(Y,i));
  }
  rectlines0(ne,x,y,lx,flag); set_avma(av);
}

/* ROt_ST */
void
plotstring(long ne, char *str, long dir)
{
  PariRect *e = check_rect_init(ne);
  RectObj *z = (RectObj*) pari_malloc(sizeof(RectObjST));
  long l = strlen(str);
  char *s = (char *) pari_malloc(l+1);

  memcpy(s,str,l+1);
  RoType(z) = ROt_ST;
  RoSTl(z) = l;
  RoSTs(z) = s;
  RoSTx(z) = RXscale(e)*RXcursor(e)+RXshift(e);
  RoSTy(z) = RYscale(e)*RYcursor(e)+RYshift(e);
  RoSTdir(z) = dir;
  Rchain(e, z);
  RoCol(z) = current_color[ne];
}

/* ROt_PTT */
void
plotpointtype(long ne, long type)
{
 if (ne == -1) plotpoint_itype = type;
 else {
   PariRect *e = check_rect_init(ne);
   RectObj *z = (RectObj*) pari_malloc(sizeof(RectObjPN));
   RoType(z) = ROt_PTT;
   RoPTTpen(z) = type;
   Rchain(e, z);
 }
}

/* ROt_PTS. FIXME: this function is a noop, since no graphic driver implement
 * this code. ne == -1 (change globally). */
void
plotpointsize(long ne, GEN size)
{
 if (ne == -1) { /*do nothing*/ }
 else {
   PariRect *e = check_rect_init(ne);
   RectObj *z = (RectObj*) pari_malloc(sizeof(RectObjPS));
   RoType(z) = ROt_PTS;
   RoPTSsize(z) = gtodouble(size);
   Rchain(e, z);
 }
}

void
plotlinetype(long ne, long type)
{
 if (ne == -1) rectline_itype = type;
 else {
   PariRect *e = check_rect_init(ne);
   RectObj *z = (RectObj*) pari_malloc(sizeof(RectObjPN));
   RoType(z) = ROt_LNT;
   RoLNTpen(z) = type;
   Rchain(e, z);
 }
}

#define RECT_CP_RELATIVE  0x1
#define RECT_CP_NW        0x0
#define RECT_CP_SW        0x2
#define RECT_CP_SE        0x4
#define RECT_CP_NE        0x6

static double*
cpd(double* R, size_t t)
{ void *o = pari_malloc(t * sizeof(double)); memcpy(o,R,t); return (double*)o; }
static void*
cp(void* R, size_t t)
{ void *o = pari_malloc(t); memcpy(o,R,t); return o; }
void
plotcopy(long source, long dest, GEN xoff, GEN yoff, long flag)
{
  PariRect *s = check_rect_init(source), *d = check_rect_init(dest);
  RectObj *R, *tail = RTail(d);
  long i, x, y;
  if (flag & RECT_CP_RELATIVE) {
    double xd = gtodouble(xoff), yd = gtodouble(yoff);
    PARI_plot T;
    if (xd > 1) pari_err_DOMAIN("plotcopy","dx",">",gen_1,xoff);
    if (xd < 0) pari_err_DOMAIN("plotcopy","dx","<",gen_0,xoff);
    if (yd > 1) pari_err_DOMAIN("plotcopy","dy",">",gen_1,yoff);
    if (yd < 0) pari_err_DOMAIN("plotcopy","dy","<",gen_0,yoff);
    pari_get_plot(&T);
    x = DTOL(xd * (T.width-1));
    y = DTOL(yd * (T.height-1));
  } else {
    if (typ(xoff) != t_INT) pari_err_TYPE("plotcopy",xoff);
    if (typ(yoff) != t_INT) pari_err_TYPE("plotcopy",yoff);
    x = itos(xoff);
    y = itos(yoff);
  }
  switch (flag & ~RECT_CP_RELATIVE)
  {
    case RECT_CP_NW: break;
    case RECT_CP_SW: y = RYsize(d) - RYsize(s) - y; break;
    case RECT_CP_SE: y = RYsize(d) - RYsize(s) - y; /* fall through */
    case RECT_CP_NE: x = RXsize(d) - RXsize(s) - x; break;
  }
  for (R = RHead(s); R; R = RoNext(R))
  {
    RectObj *o;
    switch(RoType(R))
    {
      case ROt_PT:
        o = (RectObj*)cp(R, sizeof(RectObj1P));
        RoPTx(o) += x; RoPTy(o) += y;
        break;
      case ROt_LN: case ROt_BX: case ROt_FBX:
        o = (RectObj*)cp(R, sizeof(RectObj2P));
        RoLNx1(o) += x; RoLNy1(o) += y;
        RoLNx2(o) += x; RoLNy2(o) += y;
        break;
      case ROt_MP: case ROt_ML:
        o = (RectObj*)cp(R, sizeof(RectObjMP));
        RoMPxs(o) = (double*)cp(RoMPxs(R), sizeof(double)*RoMPcnt(o));
        RoMPys(o) = (double*)cp(RoMPys(R), sizeof(double)*RoMPcnt(o));
        for (i=0; i<RoMPcnt(o); i++) { RoMPxs(o)[i] += x; RoMPys(o)[i] += y; }
        break;
      case ROt_ST:
        o = (RectObj*)cp(R, sizeof(RectObjST));
        RoSTs(o) = (char*)cp(RoSTs(R),RoSTl(R)+1);
        RoSTx(o) += x; RoSTy(o) += y;
        break;
      default: /* ROt_PTT, ROt_LNT, ROt_PTS */
        o = (RectObj*)cp(R, sizeof(RectObjPN));
        break;
    }
    RoNext(tail) = o; tail = o;
  }
  RoNext(tail) = NULL; RTail(d) = tail;
}

enum {CLIPLINE_NONEMPTY = 1, CLIPLINE_CLIP_1 = 2, CLIPLINE_CLIP_2 = 4};
/* A simpler way is to clip by 4 half-planes */
static int
clipline(double xmin, double xmax, double ymin, double ymax,
         double *x1p, double *y1p, double *x2p, double *y2p)
{
  int xy_exch = 0, rc = CLIPLINE_NONEMPTY;
  double t, sl;
  double xi, xmn, xmx;
  double yi, ymn, ymx;
  int x1_is_ymn, x1_is_xmn;
  double x1 = *x1p, x2 = *x2p, y1 = *y1p, y2 = *y2p;

  if ((x1 < xmin &&  x2 < xmin) || (x1 > xmax && x2 > xmax))
    return 0;
  if (fabs(x1 - x2) < fabs(y1 - y2)) { /* Exchange x and y */
    xy_exch = 1;
    dswap(xmin, ymin); dswap(x1, y1);
    dswap(xmax, ymax); dswap(x2, y2);
  }

  /* Build y as a function of x */
  xi = x1;
  yi = y1;
  sl = x1==x2? 0: (y2 - yi)/(x2 - xi);

  if (x1 > x2) {
    x1_is_xmn = 0;
    xmn = x2;
    xmx = x1;
  } else {
    x1_is_xmn = 1;
    xmn = x1;
    xmx = x2;
  }

  if (xmn < xmin) {
    xmn = xmin;
    rc |= x1_is_xmn? CLIPLINE_CLIP_1: CLIPLINE_CLIP_2;
  }
  if (xmx > xmax) {
    xmx = xmax;
    rc |= x1_is_xmn? CLIPLINE_CLIP_2: CLIPLINE_CLIP_1;
  }
  if (xmn > xmx) return 0;

  ymn = yi + (xmn - xi)*sl;
  ymx = yi + (xmx - xi)*sl;

  if (sl < 0) t = ymn, ymn = ymx, ymx = t;
  if (ymn > ymax || ymx < ymin) return 0;

  if (rc & CLIPLINE_CLIP_1) x1 = x1_is_xmn? xmn: xmx;
  if (rc & CLIPLINE_CLIP_2) x2 = x1_is_xmn? xmx: xmn;

  /* Now we know there is an intersection, need to move x1 and x2 */
  x1_is_ymn = ((sl >= 0) == (x1 < x2));
  if (ymn < ymin) {
    double x = (ymin - yi)/sl + xi; /* slope != 0  ! */
    if (x1_is_ymn) x1 = x, rc |= CLIPLINE_CLIP_1;
    else           x2 = x, rc |= CLIPLINE_CLIP_2;
  }
  if (ymx > ymax) {
    double x = (ymax - yi)/sl + xi; /* slope != 0  ! */
    if (x1_is_ymn) x2 = x, rc |= CLIPLINE_CLIP_2;
    else           x1 = x, rc |= CLIPLINE_CLIP_1;
  }
  if (rc & CLIPLINE_CLIP_1) y1 = yi + (x1 - xi)*sl;
  if (rc & CLIPLINE_CLIP_2) y2 = yi + (x2 - xi)*sl;
  if (xy_exch) /* Exchange x and y */
    *x1p = y1, *x2p = y2, *y1p = x1, *y2p = x2;
  else
    *x1p = x1, *x2p = x2, *y1p = y1, *y2p = y2;
  return rc;
}

void
plotclip(long rect)
{
  PariRect *s = check_rect_init(rect);
  RectObj *next, *R = RHead(s), **prevp = &RHead(s);
  double xmin = 0, xmax = RXsize(s);
  double ymin = 0, ymax = RYsize(s);

  for (; R; R = next) {
    int did_clip = 0;
#define REMOVE() { *prevp = next; freeobj(R); break; }
#define NEXT() { prevp = &RoNext(R); break; }

    next = RoNext(R);
    switch(RoType(R)) {
      case ROt_PT:
        if ( DTOL(RoPTx(R)) < xmin || DTOL(RoPTx(R)) > xmax
          || DTOL(RoPTy(R)) < ymin || DTOL(RoPTy(R)) > ymax) REMOVE();
        NEXT();
      case ROt_BX: case ROt_FBX:
        if (RoLNx1(R) < xmin) RoLNx1(R) = xmin, did_clip = 1;
        if (RoLNx2(R) < xmin) RoLNx2(R) = xmin, did_clip = 1;
        if (RoLNy1(R) < ymin) RoLNy1(R) = ymin, did_clip = 1;
        if (RoLNy2(R) < ymin) RoLNy2(R) = ymin, did_clip = 1;
        if (RoLNx1(R) > xmax) RoLNx1(R) = xmax, did_clip = 1;
        if (RoLNx2(R) > xmax) RoLNx2(R) = xmax, did_clip = 1;
        if (RoLNy1(R) > ymax) RoLNy1(R) = ymax, did_clip = 1;
        if (RoLNy2(R) > ymax) RoLNy2(R) = ymax, did_clip = 1;
        /* Remove zero-size clipped boxes */
        if (did_clip && RoLNx1(R) == RoLNx2(R)
                     && RoLNy1(R) == RoLNy2(R)) REMOVE();
        NEXT();
      case ROt_LN:
        if (!clipline(xmin, xmax, ymin, ymax,
                      &RoLNx1(R), &RoLNy1(R),
                      &RoLNx2(R), &RoLNy2(R))) REMOVE();
        NEXT();
      case ROt_MP: {
        int c = RoMPcnt(R), f = 0, t = 0;

        while (f < c) {
          if ( DTOL(RoMPxs(R)[f]) >= xmin && DTOL(RoMPxs(R)[f]) <= xmax
            && DTOL(RoMPys(R)[f]) >= ymin && DTOL(RoMPys(R)[f]) <= ymax) {
            if (t != f) {
              RoMPxs(R)[t] = RoMPxs(R)[f];
              RoMPys(R)[t] = RoMPys(R)[f];
            }
            t++;
          }
          f++;
        }
        if (t == 0) REMOVE();
        RoMPcnt(R) = t;
        NEXT();
      }
      case ROt_ML: {
        /* Hard case. Break a multiline into several pieces
         * if some part is clipped. */
        int c = RoMPcnt(R) - 1;
        int f = 0, t = 0, had_lines = 0, had_hole = 0, rc;
        double ox = RoMLxs(R)[0], oy = RoMLys(R)[0], oxn, oyn;

        while (f < c) {
        /* Endpoint of this segment is startpoint of next one: need to
         * preserve it if it is clipped. */
          oxn = RoMLxs(R)[f+1];
          oyn = RoMLys(R)[f+1];
          rc = clipline(xmin, xmax, ymin, ymax,
                  &ox, &oy, /* &RoMLxs(R)[f], &RoMLys(R)[f], */
                  &RoMLxs(R)[f+1], &RoMLys(R)[f+1]);
          RoMLxs(R)[f] = ox; ox = oxn;
          RoMLys(R)[f] = oy; oy = oyn;
          if (!rc) {
            if (had_lines) had_hole = 1;
            f++; continue;
          }

          if (!had_lines || (!(rc & CLIPLINE_CLIP_1) && !had_hole) ) {
            /* Continuous */
            had_lines = 1;
            if (t != f) {
              if (t == 0) {
                RoMPxs(R)[t] = RoMPxs(R)[f];
                RoMPys(R)[t] = RoMPys(R)[f];
              }
              RoMPxs(R)[t+1] = RoMPxs(R)[f+1];
              RoMPys(R)[t+1] = RoMPys(R)[f+1];
            }
            t++;
            f++;
            if (rc & CLIPLINE_CLIP_2) had_hole = 1, RoMLcnt(R) = t+1;
            continue;
          }
          /* Is not continuous, automatically R is not pari_free()ed.  */
          t++;
          RoMLcnt(R) = t;
          if (rc & CLIPLINE_CLIP_2) { /* Needs separate entry */
            RectObj *n = (RectObj*) pari_malloc(sizeof(RectObj2P));
            RoType(n) = ROt_LN;
            RoCol(n) = RoCol(R);
            RoLNx1(n) = RoMLxs(R)[f];   RoLNy1(n) = RoMLys(R)[f];
            RoLNx2(n) = RoMLxs(R)[f+1]; RoLNy2(n) = RoMLys(R)[f+1];
            RoNext(n) = next;
            RoNext(R) = n;
            /* Restore the unclipped value: */
            RoMLxs(R)[f+1] = oxn; RoMLys(R)[f+1] = oyn;
            f++;
            prevp = &RoNext(n);
          }
          if (f + 1 < c) { /* Are other lines */
            RectObj *n = (RectObj*) pari_malloc(sizeof(RectObjMP));
            RoType(n) = ROt_ML;
            RoCol(n) = RoCol(R);
            RoMLcnt(n) = c - f;
            RoMLxs(n) = cpd(RoMPxs(R) + f, c-f);
            RoMLys(n) = cpd(RoMPys(R) + f, c-f);
            RoMPxs(n)[0] = oxn;
            RoMPys(n)[0] = oyn;
            RoNext(n) = next;
            RoNext(R) = n;
            next = n;
          }
          break;
        }
        if (t == 0) REMOVE();
        NEXT();
      }
    }
#undef REMOVE
#undef NEXT
  }
}

/********************************************************************/
/**                                                                **/
/**                        HI-RES PLOT                             **/
/**                                                                **/
/********************************************************************/
static void
set_xrange(dblPointList *f, double x)
{ if (x < f->xsml) f->xsml = x;
  if (x > f->xbig) f->xbig = x; }
static void
Appendx(dblPointList *f, dblPointList *l, double x)
{ (l->d)[l->nb++] = x; set_xrange(f,x); }
static void
set_yrange(dblPointList *f, double y)
{ if (y < f->ysml) f->ysml = y;
  if (y > f->ybig) f->ybig = y; }
static void
Appendy(dblPointList *f, dblPointList *l, double y)
{ (l->d)[l->nb++] = y; set_yrange(f,y); }

static void
get_xy(long cplx, GEN t, double *x, double *y)
{
  GEN a, b;
  if (cplx)
  {
    if (typ(t) == t_VEC)
    {
      if (lg(t) != 2) pari_err_DIM("get_xy");
      t = gel(t,1);
    }
    a = real_i(t); b = imag_i(t);
  }
  else
  {
    if (typ(t) != t_VEC || lg(t) != 3) pari_err_DIM("get_xy");
    a = gel(t,1); b = gel(t,2);
  }
  *x = gtodouble(a);
  *y = gtodouble(b);
}
/* t a t_VEC (possibly a scalar if cplx), get next (x,y) coordinate starting
 * at index *i [update i] */
static void
get_xy_from_vec(long cplx, GEN t, long *i, double *x, double *y)
{
  GEN a, b;
  if (cplx)
  {
    if (typ(t) == t_VEC) t = gel(t,*i);
    a = real_i(t); b = imag_i(t); (*i)++;
  }
  else
  {
    a = gel(t, (*i)++);
    b = gel(t, (*i)++);
  }
  *x = gtodouble(a);
  *y = gtodouble(b);
}

static void
dblV_from_RgV(dblPointList *L, GEN x)
{
  long j, l = lg(x);
  double *X = (double*)pari_malloc(l*sizeof(double));
  for (j = 1; j < l; j++) X[j-1] = gtodouble(gel(x,j));
  L->d = X; L->nb = l-1;
}

/* Convert data from GEN to double before we call plotrecthrawin. */
static dblPointList*
gtodblList(GEN data, long flags)
{
  dblPointList *l, *L;
  double *X, *Y;
  long nl = lg(data)-1, i, j;
  const long param = (flags & (PLOT_PARAMETRIC|PLOT_COMPLEX));
  const long cplx = (flags & PLOT_COMPLEX);

  if (! is_vec_t(typ(data))) pari_err_TYPE("gtodblList",data);
  if (!nl) return NULL;

  if (nl == 1 && !cplx) pari_err_DIM("gtodblList");
  for (i = 1; i <= nl; i++)
  { /* Check input first */
    GEN x = gel(data,i);
    if (!is_vec_t(typ(x))) pari_err_TYPE("gtodblList",x);
  }
  if (flags & PLOT_PARAMETRIC)
  {
    if (odd(nl))
      pari_err_TYPE("gtodbllist [odd #components in parametric plot]", data);
    for (i = 1; i < nl; i += 2)
    {
      GEN x = gel(data,i), y = gel(data,i+1);
      if (lg(y) != lg(x)) pari_err_DIM("gtodblList");
    }
  }
  else if (!cplx)
  {
    long l1 = lg(gel(data,1)); if (l1 == 1) return NULL;
    for (i = 2; i <= nl; i++)
      if (lg(gel(data,i)) != l1) pari_err_DIM("gtodblList");
  }

  /* Now allocate memory and convert coord. to double */
  if (cplx)
  {
    l = (dblPointList*)pari_malloc((2*nl)*sizeof(dblPointList));
    for (i = 0; i < nl; i++)
    {
      pari_sp av = avma;
      GEN x = gel(data,i+1);
      dblV_from_RgV(&l[2*i],   real_i(x));
      dblV_from_RgV(&l[2*i+1], imag_i(x)); set_avma(av);
    }
  }
  else if (param)
  {
    l = (dblPointList*)pari_malloc(nl*sizeof(dblPointList));
    for (i = 1; i < nl; i += 2)
    {
      dblV_from_RgV(&l[i-1], gel(data,i));
      dblV_from_RgV(&l[i], gel(data,i+1));
    }
  }
  else
  {
    l = (dblPointList*)pari_malloc(nl*sizeof(dblPointList));
    dblV_from_RgV(&l[0], gel(data,1));
    for (i = 2; i <= nl; i++) dblV_from_RgV(&l[i-1], gel(data, i));
  }
  /* Compute extremas */
  L = &l[0];
  if (param)
  {
    L->nb = cplx? nl: nl/2;
    for (i=0; i < L->nb; i+=2)
      if (l[i+1].nb) break;
    if (i >= L->nb) { pari_free(l); return NULL; }
    L->xsml = L->xbig = l[i  ].d[0];
    L->ysml = L->ybig = l[i+1].d[0];
    for (; i < L->nb; i+=2)
    {
      long nbi = l[i+1].nb; X = l[i].d; Y = l[i+1].d;
      for (j = 0; j < nbi; j++) { set_xrange(L, X[j]); set_yrange(L, Y[j]); }
    }
  }
  else
  {
    L->nb = nl-1;
    X = L->d;   L->xsml = L->xbig = X[0];
    Y = l[1].d; L->ysml = L->ybig = Y[0];
    for (j=0; j < l[1].nb; j++) set_xrange(L, X[j]);
    for (i=1; i <= L->nb; i++)
    {
      long nbi = l[i].nb; Y = l[i].d;
      for (j = 0; j < nbi; j++) set_yrange(L, Y[j]);
    }
  }
  return l;
}

/* x,y t_REAL; return (x+y)/2,  */
static GEN
rmiddle(GEN x, GEN y) { GEN z = addrr(x,y); shiftr_inplace(z,-1); return z; }

static void
single_recursion(void *E, GEN(*eval)(void*,GEN), dblPointList *pl,
                 GEN xl,double yl, GEN xr,double yr,long depth)
{
  GEN xx;
  pari_sp av = avma;
  double yy, dy=pl[0].ybig - pl[0].ysml;

  if (depth==RECUR_MAXDEPTH) return;

  xx = rmiddle(xl,xr);
  yy = gtodouble(eval(E,xx));

  if (dy && fabs(yl+yr-2*yy) < dy*RECUR_PREC) return;
  single_recursion(E,eval, pl,xl,yl, xx,yy, depth+1);
  Appendx(&pl[0],&pl[0],rtodbl(xx));
  Appendy(&pl[0],&pl[1],yy);
  single_recursion(E,eval, pl,xx,yy, xr,yr, depth+1);
  set_avma(av);
}

static void
param_recursion(void *E,GEN(*eval)(void*,GEN), long cplx, dblPointList *pl,
  GEN tl,double xl, double yl, GEN tr,double xr,double yr, long depth)
{
  GEN t;
  pari_sp av = avma;
  double xx, dy=pl[0].ybig - pl[0].ysml;
  double yy, dx=pl[0].xbig - pl[0].xsml;

  if (depth==RECUR_MAXDEPTH) return;

  t = rmiddle(tl,tr);
  get_xy(cplx, eval(E,t), &xx,&yy);

  if (dx && dy && fabs(xl+xr-2*xx) < dx*RECUR_PREC
               && fabs(yl+yr-2*yy) < dy*RECUR_PREC) return;
  param_recursion(E,eval, cplx, pl, tl,xl,yl, t,xx,yy, depth+1);
  Appendx(&pl[0],&pl[0],xx);
  Appendy(&pl[0],&pl[1],yy);
  param_recursion(E,eval,cplx, pl, t,xx,yy, tr,xr,yr, depth+1);
  set_avma(av);
}

/* Graph 'code' for parameter values in [a,b], using 'N' sample points
 * (0 = use a default value); code is either a t_CLOSURE or a t_POL or a
 * t_VEC of two t_POLs from rectsplines. Returns a dblPointList of
 * (absolute) coordinates. */
static dblPointList *
plotrecthin(void *E, GEN(*eval)(void*, GEN), GEN a, GEN b, ulong flags,
            long N, long prec)
{
  const double INF = 1.0/0.0;
  const long param = flags & (PLOT_PARAMETRIC|PLOT_COMPLEX);
  const long recur = flags & PLOT_RECURSIVE;
  const long cplx = flags & PLOT_COMPLEX;
  GEN t, dx, x;
  dblPointList *pl;
  long tx, i, j, sig, nc, nl, ncoords, nbpoints, non_vec = 0;

  sig = gcmp(b,a); if (!sig) return NULL;
  if (sig < 0) swap(a, b);
  if (N == 1) pari_err_DOMAIN("ploth", "#points", "<", gen_2, stoi(N));
  if (!N) N = recur? 8: (param? 1500: 1000);
  /* compute F(a) to determine nc = #curves; nl = #coord. lists */
  x = gtofp(a, prec);
  t = eval(E, x); tx = typ(t);
  if (cplx) nc = nl = (tx == t_VEC)? lg(t)-1: 1;
  else if (param)
  {
    if (tx != t_VEC) pari_err_TYPE("ploth [not a t_VEC in parametric plot]", t);
    nl = lg(t)-1; nc = nl >> 1;
    if (odd(nl)) pari_err_TYPE("ploth [odd #components in parametric plot]",t);
  }
  else if (!is_matvec_t(tx)) { nl = 2; non_vec = 1; nc = 1; }
  else
  {
    if (tx != t_VEC) pari_err_TYPE("ploth [not a t_VEC]",t);
    nl = lg(t);
    nc = nl-1;
  }
  if (!nc) return NULL;
  if (recur && nc > 1) pari_err_TYPE("ploth [multi-curves + recursive]",t);

  ncoords = cplx? 2*nl: nl;
  nbpoints = recur? N << RECUR_MAXDEPTH: N;
  pl=(dblPointList*) pari_malloc(ncoords*sizeof(dblPointList));
  /* set [xy]sml,[xy]big to default values */
  if (param)
  {
    pl[0].xsml = INF;
    pl[0].xbig =-INF;
  } else {
    pl[0].xsml = gtodouble(a);
    pl[0].xbig = gtodouble(b);
  }
  pl[0].ysml = INF;
  pl[0].ybig =-INF;
  for (i = 0; i < ncoords; i++)
  {
    pl[i].d = (double*)pari_malloc((nbpoints+1)*sizeof(double));
    pl[i].nb=0;
  }
  dx = divru(gtofp(gsub(b,a),prec), N-1);
  if (recur)
  { /* recursive plot */
    double yl, yr = 0;
    if (param)
    {
      GEN tl = cgetr(prec), tr = cgetr(prec);
      double xl, xr = 0;
      pari_sp av2 = avma;
      affgr(a, tl);
      t = eval(E, tl);
      get_xy(cplx,t, &xl,&yl);
      for (i=0; i<N-1; i++, set_avma(av2))
      {
        if (i) { affrr(tr,tl); xl = xr; yl = yr; }
        addrrz(tl,dx,tr);
        t = eval(E, tr);
        get_xy(cplx,t, &xr,&yr);
        Appendx(&pl[0],&pl[0],xl);
        Appendy(&pl[0],&pl[1],yl);
        param_recursion(E,eval, cplx, pl, tl,xl,yl, tr,xr,yr, 0);
      }
      Appendx(&pl[0],&pl[0],xr);
      Appendy(&pl[0],&pl[1],yr);
    }
    else /* single curve */
    {
      GEN xl = cgetr(prec), xr = cgetr(prec);
      pari_sp av2 = avma;
      affgr(a,xl);
      yl = gtodouble(eval(E,xl));
      for (i=0; i<N-1; i++, set_avma(av2))
      {
        addrrz(xl,dx,xr);
        yr = gtodouble(eval(E,xr));
        Appendx(&pl[0],&pl[0],rtodbl(xl));
        Appendy(&pl[0],&pl[1],yl);
        single_recursion(E,eval, pl,xl,yl,xr,yr,0);
        affrr(xr,xl); yl = yr;
      }
      Appendx(&pl[0],&pl[0],rtodbl(xr));
      Appendy(&pl[0],&pl[1],yr);
    }
  }
  else /* nonrecursive plot */
  {
    GEN V, X = cgetg(N+1, t_VEC);
    for (i = 1; i <= N; i++) { gel(X,i) = x; x = addrr(x,dx); }
    if (flags & PLOT_PARA && eval == gp_call)
    {
      GEN worker = snm_closure(is_entry("_parapply_slice_worker"),
                               mkvec((GEN)E));
      V = gen_parapply_slice(worker, X, mt_nbthreads());
    }
    else
    {
      V = cgetg(N+1, t_VEC);
      for (i = 1; i <= N; i++) gel(V,i) = eval(E,gel(X,i));
    }
    if (param)
    {
      for (i = 1; i <= N; i++)
      {
        long nt, k, j;
        t = gel(V,i);
        if (typ(t) != t_VEC)
        {
          if (cplx) nt = 1;
          else nt = 0; /* trigger error */
        }
        else
          nt = lg(t)-1;
        if (nt != nl) pari_err_DIM("plotrecth");
        k = 0; j = 1;
        while (j <= nl)
        {
          double xx, yy;
          get_xy_from_vec(cplx, t, &j, &xx, &yy);
          Appendx(&pl[0], &pl[k++], xx);
          Appendy(&pl[0], &pl[k++], yy);
        }
      }
    }
    else if (non_vec)
      for (i = 1; i <= N; i++)
      {
        Appendy(&pl[0], &pl[1], gtodouble(gel(V,i)));
        pl[0].d[i-1] = gtodouble(gel(X,i));
      }
    else /* vector of nonparametric curves */
      for (i = 1; i <= N; i++)
      {
        t = gel(V,i);
        if (typ(t) != t_VEC || lg(t) != nl) pari_err_DIM("plotrecth");
        for (j = 1; j < nl; j++) Appendy(&pl[0], &pl[j], gtodouble(gel(t,j)));
        pl[0].d[i-1] = gtodouble(gel(X,i));
      }
  }
  pl[0].nb = nc; return pl;
}

static GEN
spline_eval(void* E, GEN x) { return gsubst((GEN)E,0,x); }

/* Uses highlevel plotting functions to implement splines as
   a low-level plotting function. */
static void
rectsplines(long ne, double *x, double *y, long lx, long flag)
{
  long i, j;
  pari_sp av0 = avma;
  GEN X = pol_x(0), xa = cgetg(lx+1, t_VEC), ya = cgetg(lx+1, t_VEC);
  GEN tas, pol3;
  long param = flag & PLOT_PARAMETRIC;
  const long fl = param | PLOT_RECURSIVE | PLOT_NO_RESCALE | PLOT_NO_FRAME
                        | PLOT_NO_AXE_Y | PLOT_NO_AXE_X;

  if (lx < 4) pari_err(e_MISC, "Too few points (%ld) for spline plot", lx);
  for (i = 1; i <= lx; i++) {
    gel(xa,i) = dbltor(x[i-1]);
    gel(ya,i) = dbltor(y[i-1]);
  }
  if (param) {
    tas = new_chunk(4);
    for (j = 1; j <= 4; j++) gel(tas,j-1) = utoipos(j);
    pol3 = cgetg(3, t_VEC);
  }
  else
    tas = pol3 = NULL; /* gcc -Wall */
  for (i = 0; i <= lx - 4; i++) {
    pari_sp av = avma;

    xa++; ya++;
    if (param) {
      gel(pol3,1) = polintspec(tas, xa, X, 4, NULL);
      gel(pol3,2) = polintspec(tas, ya, X, 4, NULL);
    } else {
      pol3 = polintspec(xa, ya, X, 4, NULL);
      tas = xa;
    }
    /* Start with 3 points */
    plotrecth((void*)pol3, &spline_eval, ne,
               i== 0 ? gel(tas,0) : gel(tas,1),
               i==lx-4 ? gel(tas,3) : gel(tas,2),
               fl, 2, DEFAULTPREC);
    set_avma(av);
  }
  set_avma(av0);
}

static void
pari_get_fmtplot(GEN fmt, PARI_plot *T)
{
  char *f = GSTR(fmt);
  if (!strcmp(f, "svg")) pari_get_svgplot(T);
  else if (!strcmp(f, "ps")) pari_get_psplot(T);
  else pari_err_TYPE("plotexport [unknown format]", fmt);
}
static GEN
fmt_convert(GEN fmt, GEN w, GEN x, GEN y, PARI_plot *T)
{
  char *f, *s = NULL;
  if (typ(fmt) != t_STR) pari_err_TYPE("plotexport",fmt);
  f = GSTR(fmt);
  if (!strcmp(f, "svg"))
    s = rect2svg(w,x,y,T);
  else if (!strcmp(f, "ps"))
    s = rect2ps(w,x,y,T);
  else
    pari_err_TYPE("plotexport [unknown format]", fmt);
  return strtoGENstr(s);
}

static void
Draw(PARI_plot *T, GEN w, GEN x, GEN y)
{
  if (!T->draw) pari_err(e_MISC,"high resolution graphics disabled");
  T->draw(T, w,x,y);
}
static void
set_range(double m, double M, double *sml, double *big)
{
  if (M - m < 1.e-9)
  {
    double d = fabs(m)/10; if (!d) d = 0.1;
    M += d; m -= d;
  }
  *sml = m; *big = M;
}
/* Plot a dblPointList. Complete with axes, bounding box, etc.
 *
 * data is an array of structs. Its meaning depends on flags :
 *
 * + data[0] contains global extremas, the number of curves to plot
 *   (data[0].nb) and a list of doubles (first set of x-coordinates).
 *
 * + data[i].nb (i>0) contains the number of points in the list
 *   data[i].d (hopefully, data[2i].nb=data[2i+1].nb when i>0...)
 *
 * + If flags contain PLOT_PARAMETRIC, the array length should be
 *   even, and successive pairs (data[2i].d, data[2i+1].d) represent
 *   curves to plot.
 *
 * + If there is no such flag, the first element is an array with
 *   x-coordinates and the following ones contain y-coordinates.
 * If W != NULL, output wrt this PARI_plot using two drawing rectangles:
 * one for labels, another for graphs. Else draw to rectwindow ne without
 * labels.
 * If fmt != NULL (requires W != NULL), output is a t_STR containing the
 * converted picture, else a bounding box */
static GEN
plotrecthrawin(GEN fmt, PARI_plot *W, long ne, dblPointList *data, long flags)
{
  const long param = flags & (PLOT_PARAMETRIC|PLOT_COMPLEX);
  const long max_graphcolors = lg(GP_DATA->graphcolors)-1;
  const pari_sp av = avma;
  dblPointList x, y;
  double xsml, xbig, ysml, ybig;
  long ltype, i, nc, w[3], wx[3], wy[3];

  if (!data) return cgetg(1,t_VEC);
  x = data[0]; nc = x.nb;
  set_range(x.xsml, x.xbig, &xsml, &xbig);
  set_range(x.ysml, x.ybig, &ysml, &ybig);
  if (W)
  { /* actual output; else output to rectwindow: no labels */
    const long se = NUMRECT-2;
    long lm, rm, tm, bm;
    char YBIG[16], YSML[16], XSML[16], XBIG[16];
    /* left/right/top/bottom margin */
    sprintf(YSML,"%.5g", ysml); sprintf(YBIG,"%.5g", ybig);
    sprintf(XSML,"%.5g", xsml); sprintf(XBIG,"%.5g", xbig);
    /* left margin has y labels with hgap on both sides of text */
    lm = maxss(strlen(YSML),strlen(YBIG))*W->fwidth + 2*W->hunit-1;
    rm = W->hunit-1;
    tm = W->vunit-1;
    bm = W->vunit+W->fheight-1;
    w[0] = wx[0] = wy[0] = evaltyp(t_VECSMALL) | evallg(3);
    w[1] = se; wx[1] = 0;  wy[1] = 0;
    w[2] = ne; wx[2] = lm; wy[2] = tm;
   /* Window (width x height) is given in pixels, correct pixels are 0..n-1,
    * whereas rect functions work with windows whose pixel range is [0,n] */
    initrect_i(se, W->width - 1, W->height - 1);
    initrect_i(ne, W->width - (lm+rm) - 1, W->height - (tm+bm) - 1);
    /* draw labels on se */
    _move(se,lm,0); plotstring(se, YBIG, RoSTdirRIGHT|RoSTdirHGAP|RoSTdirTOP);
    _move(se,lm,W->height-bm); plotstring(se,YSML, RoSTdirRIGHT|RoSTdirHGAP|RoSTdirVGAP);
    _move(se,lm,W->height-bm); plotstring(se, XSML, RoSTdirLEFT|RoSTdirTOP);
    _move(se,W->width-rm-1, W->height-bm); plotstring(se, XBIG, RoSTdirRIGHT|RoSTdirTOP);
  }
  if (!(flags & PLOT_NO_RESCALE)) plotscale0(ne, xsml, xbig, ysml, ybig);
  if (!(flags & PLOT_NO_FRAME))
  {
    long fl = (flags & PLOT_NODOUBLETICK)? TICKS_CLOCKW|TICKS_NODOUBLE
                                         : TICKS_CLOCKW;
    PARI_plot T, *pl;
    if (W) pl = W; else { pl = &T; pari_get_plot(pl); }
    plotlinetype(ne, -2); /* frame */
    current_color[ne] = colormap_to_color(DEFAULT_COLOR);
    _move(ne,xsml,ysml);
    _box(ne,xbig,ybig);
    if (!(flags & PLOT_NO_TICK_X)) {
      rectticks(pl, ne, xsml, ysml, xbig, ysml, xsml, xbig, fl);
      rectticks(pl, ne, xbig, ybig, xsml, ybig, xbig, xsml, fl);
    }
    if (!(flags & PLOT_NO_TICK_Y)) {
      rectticks(pl, ne, xbig, ysml, xbig, ybig, ysml, ybig, fl);
      rectticks(pl, ne, xsml, ybig, xsml, ysml, ybig, ysml, fl);
    }
  }
  if (!(flags & PLOT_NO_AXE_Y) && (xsml<=0 && xbig >=0))
  {
    plotlinetype(ne, -1); /* axes */
    current_color[ne] = colormap_to_color(AXIS_COLOR);
    _move(ne,0.0,ysml);
    _line(ne,0.0,ybig);
  }
  if (!(flags & PLOT_NO_AXE_X) && (ysml<=0 && ybig >=0))
  {
    plotlinetype(ne, -1); /* axes */
    current_color[ne] = colormap_to_color(AXIS_COLOR);
    _move(ne,xsml,0.0);
    _line(ne,xbig,0.0);
  }

  if (param) {
    i = 0;
    flags |= PLOT_PARAMETRIC;
    flags &= (~PLOT_COMPLEX); /* turn COMPLEX to PARAMETRIC*/
  } else i = 1;
  for (ltype = 0; ltype < nc; ltype++)
  {
    long c = GP_DATA->graphcolors[1+(ltype%max_graphcolors)];
    current_color[ne] = colormap_to_color(c);
    if (param) x = data[i++];

    y = data[i++];
    if (flags & (PLOT_POINTS_LINES|PLOT_POINTS)) {
      plotlinetype(ne, plotpoint_itype + ltype); /* Graphs */
      plotpointtype(ne,plotpoint_itype + ltype); /* Graphs */
      plotpoints0(ne, x.d, y.d, y.nb);
      if (!(flags & PLOT_POINTS_LINES)) continue;
    }

    if (flags & PLOT_SPLINES) {
      /* rectsplines will call us back with ltype == 0 */
      int old = rectline_itype;
      rectline_itype = rectline_itype + ltype;
      rectsplines(ne, x.d, y.d, y.nb, flags);
      rectline_itype = old;
    } else {
      plotlinetype(ne, rectline_itype + ltype); /* Graphs */
      rectlines0(ne, x.d, y.d, y.nb, 0);
    }
  }
  for (i--; i>=0; i--) pari_free(data[i].d);
  pari_free(data);

  if (W)
  {
    GEN s = NULL;
    if (fmt) s = fmt_convert(fmt, w, wx, wy, W); else Draw(W, w,wx,wy);
    plotkill(w[1]);
    plotkill(w[2]);
    if (fmt) return s;
  }
  set_avma(av);
  retmkvec4(dbltor(xsml), dbltor(xbig), dbltor(ysml), dbltor(ybig));
}

/*************************************************************************/
/*                                                                       */
/*                          HI-RES FUNCTIONS                             */
/*                                                                       */
/*************************************************************************/
/* If T != NULL, draw using the attached graphic (using rectwindow ne as a temp)
 * Else write to rectwindow 'ne'.
 * Graph y=f(x), x=a..b, use n points */
static GEN
plotrecth_i(GEN fmt, void *E, GEN(*f)(void*,GEN), PARI_plot *T, long ne,
            GEN a,GEN b, ulong flags,long n, long prec)
{
  pari_sp av = avma;
  dblPointList *pl = plotrecthin(E,f, a,b, flags, n, prec);
  set_avma(av); return plotrecthrawin(fmt, T, ne, pl, flags);
}
GEN
plotrecth(void *E, GEN(*f)(void*,GEN), long ne, GEN a,GEN b,
          ulong flags, long n, long prec)
{ return plotrecth_i(NULL, E,f, NULL, ne, a,b, flags&~PLOT_PARA, n, prec); }
GEN
plotrecth0(long ne, GEN a,GEN b,GEN code,ulong flags,long n, long prec)
{ EXPR_WRAP(code, plotrecth(EXPR_ARG, ne, a,b, flags, n, prec)); }
static GEN
_ploth(void *E, GEN(*f)(void*,GEN), GEN a, GEN b,long flags, long n, long prec)
{
  PARI_plot T; pari_get_plot(&T);
  return plotrecth_i(NULL, E,f, &T, NUMRECT-1, a,b, flags,n, prec);
}
GEN
ploth(void *E, GEN(*f)(void*,GEN), GEN a, GEN b,long flags, long n, long prec)
{ return _ploth(E, f, a, b, flags&~PLOT_PARA, n, prec); }
GEN
parploth(GEN a, GEN b, GEN code, long flags, long n, long prec)
{ return _ploth(code, gp_call, a, b, flags|PLOT_PARA, n, prec); }
GEN
ploth0(GEN a, GEN b, GEN code, long flags,long n, long prec)
{ EXPR_WRAP(code, ploth(EXPR_ARG, a,b,flags,n, prec)); }

GEN
psploth(void *E, GEN(*f)(void*,GEN), GEN a,GEN b, long flags, long n, long prec)
{
  PARI_plot T; pari_get_psplot(&T); T.draw = &_psdraw;
  return plotrecth_i(NULL, E,f, &T, NUMRECT-1, a,b, flags&~PLOT_PARA,n, prec);
}
GEN
psploth0(GEN a, GEN b, GEN code, long flags, long n, long prec)
{ EXPR_WRAP(code, psploth(EXPR_ARG, a, b, flags, n, prec)); }

static GEN
_plothexport(GEN fmt, void *E, GEN(*f)(void*,GEN), GEN a,GEN b, long flags,
            long n, long prec)
{
  pari_sp av = avma;
  GEN s;
  PARI_plot T; pari_get_fmtplot(fmt, &T);
  s = plotrecth_i(fmt, E,f, &T, NUMRECT-1, a,b, flags,n, prec);
  return gerepileuptoleaf(av, s);
}
GEN
plothexport(GEN fmt, void *E, GEN(*f)(void*,GEN), GEN a,GEN b, long flags,
            long n, long prec)
{ return _plothexport(fmt, E, f, a, b, flags&~PLOT_PARA, n, prec); }
GEN
plothexport0(GEN fmt, GEN a, GEN b, GEN code, long flags, long n, long prec)
{ EXPR_WRAP(code, plothexport(fmt, EXPR_ARG, a, b, flags, n, prec)); }
GEN
parplothexport(GEN fmt, GEN a, GEN b, GEN code, long flags, long n, long prec)
{ return _plothexport(fmt, code, gp_call, a, b, flags|PLOT_PARA, n, prec); }

/* Draw list of points */
static GEN
plotrecthraw_i(GEN fmt, PARI_plot *T, long ne, GEN data, long flags)
{
  dblPointList *pl = gtodblList(data,flags);
  return plotrecthrawin(fmt, T, ne, pl, flags);
}
static GEN
plothraw_i(GEN fmt, PARI_plot *T, GEN X, GEN Y, long flag)
{
  pari_sp av = avma;
  switch (flag) {
    case 0: flag = PLOT_PARAMETRIC|PLOT_POINTS; break;
    case 1: flag = PLOT_PARAMETRIC; break;
    default: flag |= PLOT_PARAMETRIC; break;
  }
  return gerepileupto(av, plotrecthraw_i(fmt, T, NUMRECT-1, mkvec2(X,Y), flag));
}
GEN
plothraw(GEN X, GEN Y, long flags)
{ PARI_plot T; pari_get_plot(&T); return plothraw_i(NULL,&T,X,Y,flags); }
GEN
psplothraw(GEN X, GEN Y, long flags)
{ PARI_plot T; pari_get_psplot(&T); T.draw = &_psdraw;
  return plothraw_i(NULL,&T,X,Y,flags); }
GEN
plotrecthraw(long ne, GEN data, long flags)
{ return plotrecthraw_i(NULL, NULL, ne, data, flags); }
GEN
plothrawexport(GEN fmt, GEN X, GEN Y, long flags)
{ PARI_plot T; pari_get_fmtplot(fmt,&T); return plothraw_i(fmt,&T,X,Y,flags); }

GEN
plothsizes(long flag)
{
  GEN vect = cgetg(1+8,t_VEC);
  PARI_plot T;

  pari_get_plot(&T);
  gel(vect,1) = stoi(T.width);
  gel(vect,2) = stoi(T.height);
  if (flag) {
    gel(vect,3) = dbltor(T.hunit*1.0/T.width);
    gel(vect,4) = dbltor(T.vunit*1.0/T.height);
    gel(vect,5) = dbltor(T.fwidth*1.0/T.width);
    gel(vect,6) = dbltor(T.fheight*1.0/T.height);
  } else {
    gel(vect,3) = stoi(T.hunit);
    gel(vect,4) = stoi(T.vunit);
    gel(vect,5) = stoi(T.fwidth);
    gel(vect,6) = stoi(T.fheight);
  }
  gel(vect,7) = stoi(T.dwidth);
  gel(vect,8) = stoi(T.dheight);
  return vect;
}

/*************************************************************************/
/*                                                                       */
/*                         POSTSCRIPT OUTPUT                             */
/*                                                                       */
/*************************************************************************/
static long
wxy_n(GEN wxy)
{
  long n;
  switch(typ(wxy))
  {
    case t_INT: return 1;
    case t_VEC:
      n = lg(wxy)-1;
      if (n%3) pari_err_DIM("plotdraw");
      return n/3;
  }
  pari_err_TYPE("plotdraw",wxy);
  return 0;/*LCOV_EXCL_LINE*/
}
static void
wxy_init(GEN wxy, GEN *pW, GEN *pX, GEN *pY, PARI_plot *T)
{
  long i, n = wxy_n(wxy);
  GEN W, X, Y;
  *pW = W = cgetg(n+1, t_VECSMALL); /* win number */
  *pX = X = cgetg(n+1, t_VECSMALL);
  *pY = Y = cgetg(n+1, t_VECSMALL); /* (x,y)-offset */
  if (typ(wxy) == t_INT)
  {
    W[1] = itos(wxy); check_rect_init(W[1]);
    X[1] = 0;
    Y[1] = 0; return;
  }
  for (i = 1; i <= n; i++)
  {
    GEN w = gel(wxy,3*i-2), x = gel(wxy,3*i-1), y = gel(wxy,3*i);
    if (typ(w) != t_INT) pari_err_TYPE("plotdraw",w);
    W[i] = itos(w); check_rect_init(W[i]);
    if (T) {
      X[i] = DTOL(gtodouble(x)*(T->width - 1));
      Y[i] = DTOL(gtodouble(y)*(T->height - 1));
    } else {
      X[i] = gtos(x);
      Y[i] = gtos(y);
    }
  }
}
/* if flag is set, rescale wrt T */
static void
gendraw(PARI_plot *T, GEN wxy, long flag)
{
  GEN w, x, y, W, X, Y;
  long i, l;
  wxy_init(wxy, &w,&x,&y, flag? T: NULL);
  l = lg(w);
  /* malloc mandatory in case draw() forks then pari_close(). Done after
   * wxy_init to avoid leak on error */
  W = cgetalloc(t_VECSMALL, l);
  X = cgetalloc(t_VECSMALL, l);
  Y = cgetalloc(t_VECSMALL, l);
  for (i = 1; i < l; i++) { W[i] = w[i]; X[i] = x[i]; Y[i] = y[i]; }
  Draw(T,W,X,Y);
  pari_free(W); pari_free(X); pari_free(Y);
}
void
psdraw(GEN wxy, long flag)
{ PARI_plot T; pari_get_psplot(&T); T.draw = flag? &_psdraw: &_psdraw_scale;
  gendraw(&T, wxy, flag); }
void
plotdraw(GEN wxy, long flag)
{ PARI_plot T; pari_get_plot(&T); gendraw(&T, wxy, flag); }
GEN
plotexport(GEN fmt, GEN wxy, long flag)
{
  pari_sp av = avma;
  GEN w, x, y;
  PARI_plot _T, *T = flag? &_T: NULL;
  if (T) pari_get_plot(T);
  wxy_init(wxy, &w, &x, &y, T);
  return gerepileuptoleaf(av, fmt_convert(fmt, w, x, y, T));
}

/* may be called after pari_close(): don't use the PARI stack */
void
gen_draw(struct plot_eng *eng, GEN w, GEN x, GEN y, double xs, double ys)
{
  void *data = eng->data;
  long i, j, lw = lg(w);
  long hgapsize = eng->pl->hunit, fheight = eng->pl->fheight;
  long vgapsize = eng->pl->vunit,  fwidth = eng->pl->fwidth;
  for(i = 1; i < lw; i++)
  {
    PariRect *e = &rectgraph[w[i]];
    RectObj *R;
    long x0 = x[i], y0 = y[i];
    for (R = RHead(e); R; R = RoNext(R))
    {
      long col = RoCol(R);
      switch(RoType(R))
      {
      case ROt_PT:
        eng->sc(data,col);
        eng->pt(data, DTOL((RoPTx(R)+x0)*xs), DTOL((RoPTy(R)+y0)*ys));
        break;
      case ROt_LN:
        eng->sc(data,col);
        eng->ln(data, DTOL((RoLNx1(R)+x0)*xs), DTOL((RoLNy1(R)+y0)*ys),
                      DTOL((RoLNx2(R)+x0)*xs), DTOL((RoLNy2(R)+y0)*ys));
        break;
      case ROt_BX:
        eng->sc(data,col);
        eng->bx(data,
                DTOL((RoBXx1(R)+x0)*xs),
                DTOL((RoBXy1(R)+y0)*ys),
                DTOL((RoBXx2(R)-RoBXx1(R))*xs),
                DTOL((RoBXy2(R)-RoBXy1(R))*ys));
        break;
      case ROt_FBX:
        eng->sc(data,col);
        eng->fb(data,
                DTOL((RoBXx1(R)+x0)*xs),
                DTOL((RoBXy1(R)+y0)*ys),
                DTOL((RoBXx2(R)-RoBXx1(R))*xs),
                DTOL((RoBXy2(R)-RoBXy1(R))*ys));
        break;
      case ROt_MP:
        {
          double *ptx = RoMPxs(R);
          double *pty = RoMPys(R);
          long     nb = RoMPcnt(R);
          struct plot_points *points =
            (struct plot_points *) pari_malloc(sizeof(*points)*nb);
          for(j=0;j<nb;j++)
          {
            points[j].x = DTOL((ptx[j]+x0)*xs);
            points[j].y = DTOL((pty[j]+y0)*ys);
          }
          eng->sc(data,col);
          eng->mp(data, nb, points);
          pari_free(points);
          break;
        }
      case ROt_ML:
        {
          double *ptx = RoMLxs(R);
          double *pty = RoMLys(R);
          long     nb = RoMLcnt(R);
          struct plot_points *points =
            (struct plot_points *) pari_malloc(sizeof(*points)*nb);
          for(j=0;j<nb;j++)
          {
            points[j].x = DTOL((ptx[j]+x0)*xs);
            points[j].y = DTOL((pty[j]+y0)*ys);
          }
          eng->sc(data,col);
          eng->ml(data, nb, points);
          pari_free(points);
          break;
        }
      case ROt_ST:
        {
          long dir = RoSTdir(R);
          long h = dir & RoSTdirHPOS_mask, hgap  = 0;
          long v = dir & RoSTdirVPOS_mask, vgap  = 0;
          long x, y, l = RoSTl(R);
          long shift = (h == RoSTdirLEFT ? 0 : (h == RoSTdirRIGHT? 2: 1));
          long vshift= (v == RoSTdirBOTTOM? 0: (v == RoSTdirTOP? 2: 1));
          if (dir & RoSTdirHGAP)
            hgap = (h == RoSTdirLEFT) ? hgapsize : -hgapsize;
          if (dir & RoSTdirVGAP)
            vgap = (v == RoSTdirBOTTOM) ? 2*vgapsize : -2*vgapsize;
          x = DTOL(xs * (RoSTx(R) + x0 + hgap - (l * fwidth * shift)/2));
          y = DTOL(ys * (RoSTy(R) + y0 - (vgap - vshift*(fheight-1))/2));
          eng->sc(data,col);
          eng->st(data, x, y, RoSTs(R), l);
          break;
        }
      default:
        break;
      }
    }
  }
}
/*************************************************************************/
/*                               SVG                                     */
/*************************************************************************/

struct svg_data {
  pari_str str;
  char hexcolor[8];  /* "#rrggbb\0" */
};
#define data_str(d) (&((struct svg_data*)(d))->str)
#define data_hexcolor(d) (((struct svg_data*)(d))->hexcolor)

/* Work with precision 1/scale */
static const float SVG_SCALE = 1024.0;

static float
svg_rescale(float x) { return x / SVG_SCALE; }

static void
svg_point(void *data, long x, long y)
{
  pari_str *S = data_str(data);

  str_printf(S, "<circle cx='%.2f' cy='%.2f' r='0.5' ",
    svg_rescale(x), svg_rescale(y));
  str_printf(S, "style='fill:%s;stroke:none;'/>", data_hexcolor(data));
}

static void
svg_line(void *data, long x1, long y1, long x2, long y2)
{
  pari_str *S = data_str(data);

  str_printf(S, "<line x1='%.2f' y1='%.2f' x2='%.2f' y2='%.2f' ",
    svg_rescale(x1), svg_rescale(y1), svg_rescale(x2), svg_rescale(y2));
  str_printf(S, "style='fill:none;stroke:%s;'/>", data_hexcolor(data));
}

static void
svg_rect(void *data, long x, long y, long w, long h)
{
  pari_str *S = data_str(data);

  str_printf(S, "<rect x='%.2f' y='%.2f' width='%.2f' height='%.2f' ",
    svg_rescale(x), svg_rescale(y), svg_rescale(w), svg_rescale(h));
  str_printf(S, "style='fill:none;stroke:%s;'/>", data_hexcolor(data));
}

static void
svg_fillrect(void *data, long x, long y, long w, long h)
{
  pari_str *S = data_str(data);
  const char * color = data_hexcolor(data);
  str_printf(S, "<rect x='%.2f' y='%.2f' width='%.2f' height='%.2f' ",
    svg_rescale(x), svg_rescale(y), svg_rescale(w), svg_rescale(h));
  str_printf(S, "style='fill:%s;stroke:%s;'/>", color, color);
}

static void
svg_points(void *data, long nb, struct plot_points *p)
{
  long i;
  for (i = 0; i < nb; i++)
    svg_point(data, p[i].x, p[i].y);
}

static void
svg_color(void *data, long col)
{
  static const char hex[] = "0123456789abcdef";
  char *c = data_hexcolor(data);
  int r, g, b;
  long_to_rgb(col, &r, &g, &b);
  c[0] = '#';
  c[1] = hex[r / 16];
  c[2] = hex[r & 15];
  c[3] = hex[g / 16];
  c[4] = hex[g & 15];
  c[5] = hex[b / 16];
  c[6] = hex[b & 15];
  c[7] = '\0';
}

static void
svg_lines(void *data, long nb, struct plot_points *p)
{
  long i;
  pari_str *S = data_str(data);

  str_printf(S, "<polyline points='");
  for (i = 0; i < nb; i++)
  {
    if (i > 0) str_printf(S, " ");
    str_printf(S, "%.2f,%.2f", svg_rescale(p[i].x), svg_rescale(p[i].y));
  }
  str_printf(S, "' style='fill:none;stroke:%s;'/>", data_hexcolor(data));
}

static void
svg_text(void *data, long x, long y, char *text, long numtext)
{
  pari_str *S = data_str(data);
  (void)numtext;
  str_printf(S, "<text x='%.5f' y='%.5f' font-size='%ld' style='fill:%s;'>%s</text>",
    svg_rescale(x),svg_rescale(y), 12, data_hexcolor(data), text);
}

static void
svg_head(PARI_plot *T, pari_str *S)
{
  str_printf(S, "<svg width='%ld' height='%ld' version='1.1' xmlns='http://www.w3.org/2000/svg'>", T->width, T->height);
}

static void
svg_tail(pari_str *S)
{
  str_printf(S, "</svg>");
}

char *
rect2svg(GEN w, GEN x, GEN y, PARI_plot *T)
{
  struct plot_eng pl;
  struct svg_data data;
  PARI_plot U;

  str_init(&data.str, 1);
  svg_color(&data, 0);
  if (!T)
  {
    long i, l = lg(w), xmax = 0, ymax = 0;
    T = &U; pari_get_svgplot(T);
    for (i = 1; i < l; i++)
    {
      PariRect *e = check_rect_init(w[i]);
      xmax = maxss(xmax, RXsize(e) + x[i]);
      ymax = maxss(ymax, RYsize(e) + y[i]);
    }
    T->width = xmax;
    T->height = ymax;
  }
  pl.data = &data;
  pl.sc = &svg_color;
  pl.pt = &svg_point;
  pl.ln = &svg_line;
  pl.bx = &svg_rect;
  pl.fb = &svg_fillrect;
  pl.mp = &svg_points;
  pl.ml = &svg_lines;
  pl.st = &svg_text;
  pl.pl = T;

  svg_head(T, &data.str);
  gen_draw(&pl, w, x, y, SVG_SCALE, SVG_SCALE);
  svg_tail(&data.str);

  return data.str.string;
}

/*************************************************************************/
/*                            POSTSCRIPT                                 */
/*************************************************************************/
static void
ps_sc(void *data, long col)
{
  pari_str *S = (pari_str*)data;
  int r, g, b; long_to_rgb(col, &r, &g, &b);
  if (!r && !g && !b)
    str_puts(S,"c0\n");
  else
    str_printf(S,"%.6f %.6f %.6f c\n", r/255., g/255., b/255.);
}

static void
ps_point(void *data, long x, long y)
{
  pari_str *S = (pari_str*)data;
  str_printf(S,"%ld %ld p\n",y,x);
}

static void
ps_line(void *data, long x1, long y1, long x2, long y2)
{
  pari_str *S = (pari_str*)data;
  str_printf(S,"%ld %ld m %ld %ld l\n",y1,x1,y2,x2);
  str_printf(S,"stroke\n");
}

static void
ps_rect(void *data, long x, long y, long w, long h)
{
  pari_str *S = (pari_str*)data;
  str_printf(S,"%ld %ld m %ld %ld l %ld %ld l %ld %ld l closepath currentlinejoin 0 setlinejoin stroke setlinejoin\n",
             y,x, y,x+w, y+h,x+w, y+h,x);
}

static void
ps_fillrect(void *data, long x, long y, long w, long h)
{
  pari_str *S = (pari_str*)data;
  str_printf(S,"%ld %ld m %ld %ld l %ld %ld l %ld %ld l closepath currentlinejoin 0 setlinejoin fill setlinejoin\n",
             y,x, y,x+w, y+h,x+w, y+h,x);
}

static void
ps_points(void *data, long nb, struct plot_points *p)
{
  long i;
  for (i=0; i<nb; i++) ps_point(data, p[i].x, p[i].y);
}

static void
ps_lines(void *data, long nb, struct plot_points *p)
{
  pari_str *S = (pari_str*)data;
  long i;
  str_printf(S,"%ld %ld m\n",p[0].y,p[0].x);
  for (i=1; i<nb; i++) str_printf(S, "%ld %ld l\n", p[i].y, p[i].x);
  str_printf(S,"stroke\n");
}

static void
ps_string(void *data, long x, long y, char *s, long length)
{
  pari_str *S = (pari_str*)data;
  (void)length;
  if (strpbrk(s, "(\\)")) {
    str_printf(S,"(");
    while (*s) {
      if ( *s=='(' || *s==')' || *s=='\\' ) str_putc(S,'\\');
      str_putc(S, *s);
      s++;
    }
  } else
    str_printf(S,"(%s", s);
  str_printf(S,") %ld %ld m 90 rotate show -90 rotate\n", y, x);
}

char *
rect2ps_i(GEN w, GEN x, GEN y, PARI_plot *T, int plotps)
{
  struct plot_eng pl;
  PARI_plot U;
  pari_str S;
  double xs = 0.65*PS_SCALE, ys = 0.65*PS_SCALE;
  if (T) /* res wrt T dimens */
  {
    if (plotps)
      xs = ys = PS_SCALE;
    else
    {
      xs *= ((double)PS_WIDTH) / T->width;
      ys *= ((double)PS_HEIGH) / T->height;
    }
  }
  else
  {
    T = &U; pari_get_psplot(T);
  }
  str_init(&S, 1);
  /* Definitions taken from post terminal of Gnuplot. */
  str_printf(&S, "%%!\n\
50 50 translate\n\
1 %d div 1 %d div scale\n\
1 setlinejoin\n\
/p {moveto 0 2 rlineto 2 0 rlineto 0 -2 rlineto closepath fill} def\n\
/c0 {0 0 0 setrgbcolor} def\n\
/c {setrgbcolor} def\n\
/l {lineto} def\n\
/m {moveto} def\n"
"/Times-Roman findfont %ld scalefont setfont\n",
PS_SCALE, PS_SCALE, DTOL(T->fheight * xs));

  pl.sc = &ps_sc;
  pl.pt = &ps_point;
  pl.ln = &ps_line;
  pl.bx = &ps_rect;
  pl.fb = &ps_fillrect;
  pl.mp = &ps_points;
  pl.ml = &ps_lines;
  pl.st = &ps_string;
  pl.pl = T;
  pl.data = (void*)&S;

  if (plotps) str_printf(&S,"0 %ld translate -90 rotate\n", (T->height - 50)*PS_SCALE);
  gen_draw(&pl, w, x, y, xs, ys);
  str_puts(&S,"stroke showpage\n");
  *S.cur = 0; return S.string;
}
char *
rect2ps(GEN w, GEN x, GEN y, PARI_plot *T)
{ return rect2ps_i(w,x,y,T,0); }

void
pari_plot_by_file(const char *env, const char *suf, const char *img)
{
  const char *cmd, *s = pari_unique_filename_suffix("plotfile", suf);
  FILE *f = fopen(s, "w");
  if (!f) pari_err_FILE("image file", s);
  fputs(img, f); (void)fclose(f);
  cmd = os_getenv(env);
#ifdef GP_MIME_OPEN
  if (!cmd) cmd = GP_MIME_OPEN;
#else
  if (!cmd) cmd = "open -W";
#endif
  cmd = pari_sprintf("%s \"%s\" 2>/dev/null", cmd, s);
  gpsystem(cmd);
  pari_unlink(s);
  pari_free((char*)s);
}

/*************************************************************************/
/*                                                                       */
/*                           RGB COLORS                                  */
/*                                                                       */
/*************************************************************************/
/* generated from /etc/X11/rgb.txt by the following perl script
#!/usr/bin/perl
while(<>)
{
  ($hex, $name) = split(/\t\t/, $_);
  $hex =~ s/^ +//; chomp($name); $name =~ s, *,,g;
  $hex = sprintf("0x%02x%02x%02x", split(/\s+/, $hex));
  $name = lc($name); next if ($done{$name});
  $done{$name} = 1;
  print "COL(\"$name\", $hex),\n";
}
*/

#define COL(x,y) {(void*)x,(void*)y,0,NULL}
static hashentry col_list[] = {
COL("", 0x000000),
COL("snow", 0xfffafa),
COL("ghostwhite", 0xf8f8ff),
COL("whitesmoke", 0xf5f5f5),
COL("gainsboro", 0xdcdcdc),
COL("floralwhite", 0xfffaf0),
COL("oldlace", 0xfdf5e6),
COL("linen", 0xfaf0e6),
COL("antiquewhite", 0xfaebd7),
COL("papayawhip", 0xffefd5),
COL("blanchedalmond", 0xffebcd),
COL("bisque", 0xffe4c4),
COL("peachpuff", 0xffdab9),
COL("navajowhite", 0xffdead),
COL("moccasin", 0xffe4b5),
COL("cornsilk", 0xfff8dc),
COL("ivory", 0xfffff0),
COL("lemonchiffon", 0xfffacd),
COL("seashell", 0xfff5ee),
COL("honeydew", 0xf0fff0),
COL("mintcream", 0xf5fffa),
COL("azure", 0xf0ffff),
COL("aliceblue", 0xf0f8ff),
COL("lavender", 0xe6e6fa),
COL("lavenderblush", 0xfff0f5),
COL("mistyrose", 0xffe4e1),
COL("white", 0xffffff),
COL("black", 0x000000),
COL("darkslategray", 0x2f4f4f),
COL("darkslategrey", 0x2f4f4f),
COL("dimgray", 0x696969),
COL("dimgrey", 0x696969),
COL("slategray", 0x708090),
COL("slategrey", 0x708090),
COL("lightslategray", 0x778899),
COL("lightslategrey", 0x778899),
COL("gray", 0xbebebe),
COL("grey", 0xbebebe),
COL("lightgrey", 0xd3d3d3),
COL("lightgray", 0xd3d3d3),
COL("midnightblue", 0x191970),
COL("navy", 0x000080),
COL("navyblue", 0x000080),
COL("cornflowerblue", 0x6495ed),
COL("darkslateblue", 0x483d8b),
COL("slateblue", 0x6a5acd),
COL("mediumslateblue", 0x7b68ee),
COL("lightslateblue", 0x8470ff),
COL("mediumblue", 0x0000cd),
COL("royalblue", 0x4169e1),
COL("blue", 0x0000ff),
COL("dodgerblue", 0x1e90ff),
COL("deepskyblue", 0x00bfff),
COL("skyblue", 0x87ceeb),
COL("lightskyblue", 0x87cefa),
COL("steelblue", 0x4682b4),
COL("lightsteelblue", 0xb0c4de),
COL("lightblue", 0xadd8e6),
COL("powderblue", 0xb0e0e6),
COL("paleturquoise", 0xafeeee),
COL("darkturquoise", 0x00ced1),
COL("mediumturquoise", 0x48d1cc),
COL("turquoise", 0x40e0d0),
COL("cyan", 0x00ffff),
COL("lightcyan", 0xe0ffff),
COL("cadetblue", 0x5f9ea0),
COL("mediumaquamarine", 0x66cdaa),
COL("aquamarine", 0x7fffd4),
COL("darkgreen", 0x006400),
COL("darkolivegreen", 0x556b2f),
COL("darkseagreen", 0x8fbc8f),
COL("seagreen", 0x2e8b57),
COL("mediumseagreen", 0x3cb371),
COL("lightseagreen", 0x20b2aa),
COL("palegreen", 0x98fb98),
COL("springgreen", 0x00ff7f),
COL("lawngreen", 0x7cfc00),
COL("green", 0x00ff00),
COL("chartreuse", 0x7fff00),
COL("mediumspringgreen", 0x00fa9a),
COL("greenyellow", 0xadff2f),
COL("limegreen", 0x32cd32),
COL("yellowgreen", 0x9acd32),
COL("forestgreen", 0x228b22),
COL("olivedrab", 0x6b8e23),
COL("darkkhaki", 0xbdb76b),
COL("khaki", 0xf0e68c),
COL("palegoldenrod", 0xeee8aa),
COL("lightgoldenrodyellow", 0xfafad2),
COL("lightyellow", 0xffffe0),
COL("yellow", 0xffff00),
COL("gold", 0xffd700),
COL("lightgoldenrod", 0xeedd82),
COL("goldenrod", 0xdaa520),
COL("darkgoldenrod", 0xb8860b),
COL("rosybrown", 0xbc8f8f),
COL("indianred", 0xcd5c5c),
COL("saddlebrown", 0x8b4513),
COL("sienna", 0xa0522d),
COL("peru", 0xcd853f),
COL("burlywood", 0xdeb887),
COL("beige", 0xf5f5dc),
COL("wheat", 0xf5deb3),
COL("sandybrown", 0xf4a460),
COL("tan", 0xd2b48c),
COL("chocolate", 0xd2691e),
COL("firebrick", 0xb22222),
COL("brown", 0xa52a2a),
COL("darksalmon", 0xe9967a),
COL("salmon", 0xfa8072),
COL("lightsalmon", 0xffa07a),
COL("orange", 0xffa500),
COL("darkorange", 0xff8c00),
COL("coral", 0xff7f50),
COL("lightcoral", 0xf08080),
COL("tomato", 0xff6347),
COL("orangered", 0xff4500),
COL("red", 0xff0000),
COL("hotpink", 0xff69b4),
COL("deeppink", 0xff1493),
COL("pink", 0xffc0cb),
COL("lightpink", 0xffb6c1),
COL("palevioletred", 0xdb7093),
COL("maroon", 0xb03060),
COL("mediumvioletred", 0xc71585),
COL("violetred", 0xd02090),
COL("magenta", 0xff00ff),
COL("violet", 0xee82ee),
COL("plum", 0xdda0dd),
COL("orchid", 0xda70d6),
COL("mediumorchid", 0xba55d3),
COL("darkorchid", 0x9932cc),
COL("darkviolet", 0x9400d3),
COL("blueviolet", 0x8a2be2),
COL("purple", 0xa020f0),
COL("mediumpurple", 0x9370db),
COL("thistle", 0xd8bfd8),
COL("snow1", 0xfffafa),
COL("snow2", 0xeee9e9),
COL("snow3", 0xcdc9c9),
COL("snow4", 0x8b8989),
COL("seashell1", 0xfff5ee),
COL("seashell2", 0xeee5de),
COL("seashell3", 0xcdc5bf),
COL("seashell4", 0x8b8682),
COL("antiquewhite1", 0xffefdb),
COL("antiquewhite2", 0xeedfcc),
COL("antiquewhite3", 0xcdc0b0),
COL("antiquewhite4", 0x8b8378),
COL("bisque1", 0xffe4c4),
COL("bisque2", 0xeed5b7),
COL("bisque3", 0xcdb79e),
COL("bisque4", 0x8b7d6b),
COL("peachpuff1", 0xffdab9),
COL("peachpuff2", 0xeecbad),
COL("peachpuff3", 0xcdaf95),
COL("peachpuff4", 0x8b7765),
COL("navajowhite1", 0xffdead),
COL("navajowhite2", 0xeecfa1),
COL("navajowhite3", 0xcdb38b),
COL("navajowhite4", 0x8b795e),
COL("lemonchiffon1", 0xfffacd),
COL("lemonchiffon2", 0xeee9bf),
COL("lemonchiffon3", 0xcdc9a5),
COL("lemonchiffon4", 0x8b8970),
COL("cornsilk1", 0xfff8dc),
COL("cornsilk2", 0xeee8cd),
COL("cornsilk3", 0xcdc8b1),
COL("cornsilk4", 0x8b8878),
COL("ivory1", 0xfffff0),
COL("ivory2", 0xeeeee0),
COL("ivory3", 0xcdcdc1),
COL("ivory4", 0x8b8b83),
COL("honeydew1", 0xf0fff0),
COL("honeydew2", 0xe0eee0),
COL("honeydew3", 0xc1cdc1),
COL("honeydew4", 0x838b83),
COL("lavenderblush1", 0xfff0f5),
COL("lavenderblush2", 0xeee0e5),
COL("lavenderblush3", 0xcdc1c5),
COL("lavenderblush4", 0x8b8386),
COL("mistyrose1", 0xffe4e1),
COL("mistyrose2", 0xeed5d2),
COL("mistyrose3", 0xcdb7b5),
COL("mistyrose4", 0x8b7d7b),
COL("azure1", 0xf0ffff),
COL("azure2", 0xe0eeee),
COL("azure3", 0xc1cdcd),
COL("azure4", 0x838b8b),
COL("slateblue1", 0x836fff),
COL("slateblue2", 0x7a67ee),
COL("slateblue3", 0x6959cd),
COL("slateblue4", 0x473c8b),
COL("royalblue1", 0x4876ff),
COL("royalblue2", 0x436eee),
COL("royalblue3", 0x3a5fcd),
COL("royalblue4", 0x27408b),
COL("blue1", 0x0000ff),
COL("blue2", 0x0000ee),
COL("blue3", 0x0000cd),
COL("blue4", 0x00008b),
COL("dodgerblue1", 0x1e90ff),
COL("dodgerblue2", 0x1c86ee),
COL("dodgerblue3", 0x1874cd),
COL("dodgerblue4", 0x104e8b),
COL("steelblue1", 0x63b8ff),
COL("steelblue2", 0x5cacee),
COL("steelblue3", 0x4f94cd),
COL("steelblue4", 0x36648b),
COL("deepskyblue1", 0x00bfff),
COL("deepskyblue2", 0x00b2ee),
COL("deepskyblue3", 0x009acd),
COL("deepskyblue4", 0x00688b),
COL("skyblue1", 0x87ceff),
COL("skyblue2", 0x7ec0ee),
COL("skyblue3", 0x6ca6cd),
COL("skyblue4", 0x4a708b),
COL("lightskyblue1", 0xb0e2ff),
COL("lightskyblue2", 0xa4d3ee),
COL("lightskyblue3", 0x8db6cd),
COL("lightskyblue4", 0x607b8b),
COL("slategray1", 0xc6e2ff),
COL("slategray2", 0xb9d3ee),
COL("slategray3", 0x9fb6cd),
COL("slategray4", 0x6c7b8b),
COL("lightsteelblue1", 0xcae1ff),
COL("lightsteelblue2", 0xbcd2ee),
COL("lightsteelblue3", 0xa2b5cd),
COL("lightsteelblue4", 0x6e7b8b),
COL("lightblue1", 0xbfefff),
COL("lightblue2", 0xb2dfee),
COL("lightblue3", 0x9ac0cd),
COL("lightblue4", 0x68838b),
COL("lightcyan1", 0xe0ffff),
COL("lightcyan2", 0xd1eeee),
COL("lightcyan3", 0xb4cdcd),
COL("lightcyan4", 0x7a8b8b),
COL("paleturquoise1", 0xbbffff),
COL("paleturquoise2", 0xaeeeee),
COL("paleturquoise3", 0x96cdcd),
COL("paleturquoise4", 0x668b8b),
COL("cadetblue1", 0x98f5ff),
COL("cadetblue2", 0x8ee5ee),
COL("cadetblue3", 0x7ac5cd),
COL("cadetblue4", 0x53868b),
COL("turquoise1", 0x00f5ff),
COL("turquoise2", 0x00e5ee),
COL("turquoise3", 0x00c5cd),
COL("turquoise4", 0x00868b),
COL("cyan1", 0x00ffff),
COL("cyan2", 0x00eeee),
COL("cyan3", 0x00cdcd),
COL("cyan4", 0x008b8b),
COL("darkslategray1", 0x97ffff),
COL("darkslategray2", 0x8deeee),
COL("darkslategray3", 0x79cdcd),
COL("darkslategray4", 0x528b8b),
COL("aquamarine1", 0x7fffd4),
COL("aquamarine2", 0x76eec6),
COL("aquamarine3", 0x66cdaa),
COL("aquamarine4", 0x458b74),
COL("darkseagreen1", 0xc1ffc1),
COL("darkseagreen2", 0xb4eeb4),
COL("darkseagreen3", 0x9bcd9b),
COL("darkseagreen4", 0x698b69),
COL("seagreen1", 0x54ff9f),
COL("seagreen2", 0x4eee94),
COL("seagreen3", 0x43cd80),
COL("seagreen4", 0x2e8b57),
COL("palegreen1", 0x9aff9a),
COL("palegreen2", 0x90ee90),
COL("palegreen3", 0x7ccd7c),
COL("palegreen4", 0x548b54),
COL("springgreen1", 0x00ff7f),
COL("springgreen2", 0x00ee76),
COL("springgreen3", 0x00cd66),
COL("springgreen4", 0x008b45),
COL("green1", 0x00ff00),
COL("green2", 0x00ee00),
COL("green3", 0x00cd00),
COL("green4", 0x008b00),
COL("chartreuse1", 0x7fff00),
COL("chartreuse2", 0x76ee00),
COL("chartreuse3", 0x66cd00),
COL("chartreuse4", 0x458b00),
COL("olivedrab1", 0xc0ff3e),
COL("olivedrab2", 0xb3ee3a),
COL("olivedrab3", 0x9acd32),
COL("olivedrab4", 0x698b22),
COL("darkolivegreen1", 0xcaff70),
COL("darkolivegreen2", 0xbcee68),
COL("darkolivegreen3", 0xa2cd5a),
COL("darkolivegreen4", 0x6e8b3d),
COL("khaki1", 0xfff68f),
COL("khaki2", 0xeee685),
COL("khaki3", 0xcdc673),
COL("khaki4", 0x8b864e),
COL("lightgoldenrod1", 0xffec8b),
COL("lightgoldenrod2", 0xeedc82),
COL("lightgoldenrod3", 0xcdbe70),
COL("lightgoldenrod4", 0x8b814c),
COL("lightyellow1", 0xffffe0),
COL("lightyellow2", 0xeeeed1),
COL("lightyellow3", 0xcdcdb4),
COL("lightyellow4", 0x8b8b7a),
COL("yellow1", 0xffff00),
COL("yellow2", 0xeeee00),
COL("yellow3", 0xcdcd00),
COL("yellow4", 0x8b8b00),
COL("gold1", 0xffd700),
COL("gold2", 0xeec900),
COL("gold3", 0xcdad00),
COL("gold4", 0x8b7500),
COL("goldenrod1", 0xffc125),
COL("goldenrod2", 0xeeb422),
COL("goldenrod3", 0xcd9b1d),
COL("goldenrod4", 0x8b6914),
COL("darkgoldenrod1", 0xffb90f),
COL("darkgoldenrod2", 0xeead0e),
COL("darkgoldenrod3", 0xcd950c),
COL("darkgoldenrod4", 0x8b6508),
COL("rosybrown1", 0xffc1c1),
COL("rosybrown2", 0xeeb4b4),
COL("rosybrown3", 0xcd9b9b),
COL("rosybrown4", 0x8b6969),
COL("indianred1", 0xff6a6a),
COL("indianred2", 0xee6363),
COL("indianred3", 0xcd5555),
COL("indianred4", 0x8b3a3a),
COL("sienna1", 0xff8247),
COL("sienna2", 0xee7942),
COL("sienna3", 0xcd6839),
COL("sienna4", 0x8b4726),
COL("burlywood1", 0xffd39b),
COL("burlywood2", 0xeec591),
COL("burlywood3", 0xcdaa7d),
COL("burlywood4", 0x8b7355),
COL("wheat1", 0xffe7ba),
COL("wheat2", 0xeed8ae),
COL("wheat3", 0xcdba96),
COL("wheat4", 0x8b7e66),
COL("tan1", 0xffa54f),
COL("tan2", 0xee9a49),
COL("tan3", 0xcd853f),
COL("tan4", 0x8b5a2b),
COL("chocolate1", 0xff7f24),
COL("chocolate2", 0xee7621),
COL("chocolate3", 0xcd661d),
COL("chocolate4", 0x8b4513),
COL("firebrick1", 0xff3030),
COL("firebrick2", 0xee2c2c),
COL("firebrick3", 0xcd2626),
COL("firebrick4", 0x8b1a1a),
COL("brown1", 0xff4040),
COL("brown2", 0xee3b3b),
COL("brown3", 0xcd3333),
COL("brown4", 0x8b2323),
COL("salmon1", 0xff8c69),
COL("salmon2", 0xee8262),
COL("salmon3", 0xcd7054),
COL("salmon4", 0x8b4c39),
COL("lightsalmon1", 0xffa07a),
COL("lightsalmon2", 0xee9572),
COL("lightsalmon3", 0xcd8162),
COL("lightsalmon4", 0x8b5742),
COL("orange1", 0xffa500),
COL("orange2", 0xee9a00),
COL("orange3", 0xcd8500),
COL("orange4", 0x8b5a00),
COL("darkorange1", 0xff7f00),
COL("darkorange2", 0xee7600),
COL("darkorange3", 0xcd6600),
COL("darkorange4", 0x8b4500),
COL("coral1", 0xff7256),
COL("coral2", 0xee6a50),
COL("coral3", 0xcd5b45),
COL("coral4", 0x8b3e2f),
COL("tomato1", 0xff6347),
COL("tomato2", 0xee5c42),
COL("tomato3", 0xcd4f39),
COL("tomato4", 0x8b3626),
COL("orangered1", 0xff4500),
COL("orangered2", 0xee4000),
COL("orangered3", 0xcd3700),
COL("orangered4", 0x8b2500),
COL("red1", 0xff0000),
COL("red2", 0xee0000),
COL("red3", 0xcd0000),
COL("red4", 0x8b0000),
COL("debianred", 0xd70751),
COL("deeppink1", 0xff1493),
COL("deeppink2", 0xee1289),
COL("deeppink3", 0xcd1076),
COL("deeppink4", 0x8b0a50),
COL("hotpink1", 0xff6eb4),
COL("hotpink2", 0xee6aa7),
COL("hotpink3", 0xcd6090),
COL("hotpink4", 0x8b3a62),
COL("pink1", 0xffb5c5),
COL("pink2", 0xeea9b8),
COL("pink3", 0xcd919e),
COL("pink4", 0x8b636c),
COL("lightpink1", 0xffaeb9),
COL("lightpink2", 0xeea2ad),
COL("lightpink3", 0xcd8c95),
COL("lightpink4", 0x8b5f65),
COL("palevioletred1", 0xff82ab),
COL("palevioletred2", 0xee799f),
COL("palevioletred3", 0xcd6889),
COL("palevioletred4", 0x8b475d),
COL("maroon1", 0xff34b3),
COL("maroon2", 0xee30a7),
COL("maroon3", 0xcd2990),
COL("maroon4", 0x8b1c62),
COL("violetred1", 0xff3e96),
COL("violetred2", 0xee3a8c),
COL("violetred3", 0xcd3278),
COL("violetred4", 0x8b2252),
COL("magenta1", 0xff00ff),
COL("magenta2", 0xee00ee),
COL("magenta3", 0xcd00cd),
COL("magenta4", 0x8b008b),
COL("orchid1", 0xff83fa),
COL("orchid2", 0xee7ae9),
COL("orchid3", 0xcd69c9),
COL("orchid4", 0x8b4789),
COL("plum1", 0xffbbff),
COL("plum2", 0xeeaeee),
COL("plum3", 0xcd96cd),
COL("plum4", 0x8b668b),
COL("mediumorchid1", 0xe066ff),
COL("mediumorchid2", 0xd15fee),
COL("mediumorchid3", 0xb452cd),
COL("mediumorchid4", 0x7a378b),
COL("darkorchid1", 0xbf3eff),
COL("darkorchid2", 0xb23aee),
COL("darkorchid3", 0x9a32cd),
COL("darkorchid4", 0x68228b),
COL("purple1", 0x9b30ff),
COL("purple2", 0x912cee),
COL("purple3", 0x7d26cd),
COL("purple4", 0x551a8b),
COL("mediumpurple1", 0xab82ff),
COL("mediumpurple2", 0x9f79ee),
COL("mediumpurple3", 0x8968cd),
COL("mediumpurple4", 0x5d478b),
COL("thistle1", 0xffe1ff),
COL("thistle2", 0xeed2ee),
COL("thistle3", 0xcdb5cd),
COL("thistle4", 0x8b7b8b),
COL("gray0", 0x000000),
COL("grey0", 0x000000),
COL("gray1", 0x030303),
COL("grey1", 0x030303),
COL("gray2", 0x050505),
COL("grey2", 0x050505),
COL("gray3", 0x080808),
COL("grey3", 0x080808),
COL("gray4", 0x0a0a0a),
COL("grey4", 0x0a0a0a),
COL("gray5", 0x0d0d0d),
COL("grey5", 0x0d0d0d),
COL("gray6", 0x0f0f0f),
COL("grey6", 0x0f0f0f),
COL("gray7", 0x121212),
COL("grey7", 0x121212),
COL("gray8", 0x141414),
COL("grey8", 0x141414),
COL("gray9", 0x171717),
COL("grey9", 0x171717),
COL("gray10", 0x1a1a1a),
COL("grey10", 0x1a1a1a),
COL("gray11", 0x1c1c1c),
COL("grey11", 0x1c1c1c),
COL("gray12", 0x1f1f1f),
COL("grey12", 0x1f1f1f),
COL("gray13", 0x212121),
COL("grey13", 0x212121),
COL("gray14", 0x242424),
COL("grey14", 0x242424),
COL("gray15", 0x262626),
COL("grey15", 0x262626),
COL("gray16", 0x292929),
COL("grey16", 0x292929),
COL("gray17", 0x2b2b2b),
COL("grey17", 0x2b2b2b),
COL("gray18", 0x2e2e2e),
COL("grey18", 0x2e2e2e),
COL("gray19", 0x303030),
COL("grey19", 0x303030),
COL("gray20", 0x333333),
COL("grey20", 0x333333),
COL("gray21", 0x363636),
COL("grey21", 0x363636),
COL("gray22", 0x383838),
COL("grey22", 0x383838),
COL("gray23", 0x3b3b3b),
COL("grey23", 0x3b3b3b),
COL("gray24", 0x3d3d3d),
COL("grey24", 0x3d3d3d),
COL("gray25", 0x404040),
COL("grey25", 0x404040),
COL("gray26", 0x424242),
COL("grey26", 0x424242),
COL("gray27", 0x454545),
COL("grey27", 0x454545),
COL("gray28", 0x474747),
COL("grey28", 0x474747),
COL("gray29", 0x4a4a4a),
COL("grey29", 0x4a4a4a),
COL("gray30", 0x4d4d4d),
COL("grey30", 0x4d4d4d),
COL("gray31", 0x4f4f4f),
COL("grey31", 0x4f4f4f),
COL("gray32", 0x525252),
COL("grey32", 0x525252),
COL("gray33", 0x545454),
COL("grey33", 0x545454),
COL("gray34", 0x575757),
COL("grey34", 0x575757),
COL("gray35", 0x595959),
COL("grey35", 0x595959),
COL("gray36", 0x5c5c5c),
COL("grey36", 0x5c5c5c),
COL("gray37", 0x5e5e5e),
COL("grey37", 0x5e5e5e),
COL("gray38", 0x616161),
COL("grey38", 0x616161),
COL("gray39", 0x636363),
COL("grey39", 0x636363),
COL("gray40", 0x666666),
COL("grey40", 0x666666),
COL("gray41", 0x696969),
COL("grey41", 0x696969),
COL("gray42", 0x6b6b6b),
COL("grey42", 0x6b6b6b),
COL("gray43", 0x6e6e6e),
COL("grey43", 0x6e6e6e),
COL("gray44", 0x707070),
COL("grey44", 0x707070),
COL("gray45", 0x737373),
COL("grey45", 0x737373),
COL("gray46", 0x757575),
COL("grey46", 0x757575),
COL("gray47", 0x787878),
COL("grey47", 0x787878),
COL("gray48", 0x7a7a7a),
COL("grey48", 0x7a7a7a),
COL("gray49", 0x7d7d7d),
COL("grey49", 0x7d7d7d),
COL("gray50", 0x7f7f7f),
COL("grey50", 0x7f7f7f),
COL("gray51", 0x828282),
COL("grey51", 0x828282),
COL("gray52", 0x858585),
COL("grey52", 0x858585),
COL("gray53", 0x878787),
COL("grey53", 0x878787),
COL("gray54", 0x8a8a8a),
COL("grey54", 0x8a8a8a),
COL("gray55", 0x8c8c8c),
COL("grey55", 0x8c8c8c),
COL("gray56", 0x8f8f8f),
COL("grey56", 0x8f8f8f),
COL("gray57", 0x919191),
COL("grey57", 0x919191),
COL("gray58", 0x949494),
COL("grey58", 0x949494),
COL("gray59", 0x969696),
COL("grey59", 0x969696),
COL("gray60", 0x999999),
COL("grey60", 0x999999),
COL("gray61", 0x9c9c9c),
COL("grey61", 0x9c9c9c),
COL("gray62", 0x9e9e9e),
COL("grey62", 0x9e9e9e),
COL("gray63", 0xa1a1a1),
COL("grey63", 0xa1a1a1),
COL("gray64", 0xa3a3a3),
COL("grey64", 0xa3a3a3),
COL("gray65", 0xa6a6a6),
COL("grey65", 0xa6a6a6),
COL("gray66", 0xa8a8a8),
COL("grey66", 0xa8a8a8),
COL("gray67", 0xababab),
COL("grey67", 0xababab),
COL("gray68", 0xadadad),
COL("grey68", 0xadadad),
COL("gray69", 0xb0b0b0),
COL("grey69", 0xb0b0b0),
COL("gray70", 0xb3b3b3),
COL("grey70", 0xb3b3b3),
COL("gray71", 0xb5b5b5),
COL("grey71", 0xb5b5b5),
COL("gray72", 0xb8b8b8),
COL("grey72", 0xb8b8b8),
COL("gray73", 0xbababa),
COL("grey73", 0xbababa),
COL("gray74", 0xbdbdbd),
COL("grey74", 0xbdbdbd),
COL("gray75", 0xbfbfbf),
COL("grey75", 0xbfbfbf),
COL("gray76", 0xc2c2c2),
COL("grey76", 0xc2c2c2),
COL("gray77", 0xc4c4c4),
COL("grey77", 0xc4c4c4),
COL("gray78", 0xc7c7c7),
COL("grey78", 0xc7c7c7),
COL("gray79", 0xc9c9c9),
COL("grey79", 0xc9c9c9),
COL("gray80", 0xcccccc),
COL("grey80", 0xcccccc),
COL("gray81", 0xcfcfcf),
COL("grey81", 0xcfcfcf),
COL("gray82", 0xd1d1d1),
COL("grey82", 0xd1d1d1),
COL("gray83", 0xd4d4d4),
COL("grey83", 0xd4d4d4),
COL("gray84", 0xd6d6d6),
COL("grey84", 0xd6d6d6),
COL("gray85", 0xd9d9d9),
COL("grey85", 0xd9d9d9),
COL("gray86", 0xdbdbdb),
COL("grey86", 0xdbdbdb),
COL("gray87", 0xdedede),
COL("grey87", 0xdedede),
COL("gray88", 0xe0e0e0),
COL("grey88", 0xe0e0e0),
COL("gray89", 0xe3e3e3),
COL("grey89", 0xe3e3e3),
COL("gray90", 0xe5e5e5),
COL("grey90", 0xe5e5e5),
COL("gray91", 0xe8e8e8),
COL("grey91", 0xe8e8e8),
COL("gray92", 0xebebeb),
COL("grey92", 0xebebeb),
COL("gray93", 0xededed),
COL("grey93", 0xededed),
COL("gray94", 0xf0f0f0),
COL("grey94", 0xf0f0f0),
COL("gray95", 0xf2f2f2),
COL("grey95", 0xf2f2f2),
COL("gray96", 0xf5f5f5),
COL("grey96", 0xf5f5f5),
COL("gray97", 0xf7f7f7),
COL("grey97", 0xf7f7f7),
COL("gray98", 0xfafafa),
COL("grey98", 0xfafafa),
COL("gray99", 0xfcfcfc),
COL("grey99", 0xfcfcfc),
COL("gray100", 0xffffff),
COL("grey100", 0xffffff),
COL("darkgrey", 0xa9a9a9),
COL("darkgray", 0xa9a9a9),
COL("darkblue", 0x00008b),
COL("darkcyan", 0x008b8b),
COL("darkmagenta", 0x8b008b),
COL("darkred", 0x8b0000),
COL("lightgreen", 0x90ee90),
COL(NULL,0) /* sentinel */
};
#undef COL

void
long_to_rgb(long c, int *r, int *g, int *b)
{
  *b = c & 0xff; c >>= 8;
  *g = c & 0xff; c >>= 8;
  *r = c;
}
static int
hex2(const char *s)
{
  int m = 0, c = 0, i;
  for (i = 0; i < 2; i++, s++)
  {
    if (*s >= '0' && *s <= '9')
      c = *s - '0';
    else if (*s >= 'A' && *s <= 'F')
      c = *s - 'A' + 10;
    else if (*s >= 'a' && *s <= 'f')
      c = *s - 'a' + 10;
    else pari_err(e_MISC,"incorrect hexadecimal number: %s", s);
    m = 16*m + c;
  }
  return m;
}
void
colorname_to_rgb(const char *s, int *r, int *g, int *b)
{
  if (!rgb_colors) rgb_colors = hashstr_import_static(col_list, 1000);
  if (*s == '#' && strlen(s) == 7)
  {
    *r = hex2(s+1);
    *g = hex2(s+3);
    *b = hex2(s+5);
  }
  else
  {
    hashentry *ep = hash_search(rgb_colors, (void*)s);
    if (!ep) pari_err(e_MISC, "unknown color %s", s);
    long_to_rgb((long)ep->val, r,g,b);
  }
}

static void
chk_8bit(int v, GEN c)
{ if (v & ~0xff) pari_err(e_MISC, "invalid RGB code: %Ps", c); }
void
color_to_rgb(GEN c, int *r, int *g, int *b)
{
  switch(typ(c))
  {
    case t_STR:
      colorname_to_rgb(GSTR(c), r,g,b);
      break;
    default: /* t_VECSMALL: */
      *r = c[1]; chk_8bit(*r, c);
      *g = c[2]; chk_8bit(*g, c);
      *b = c[3]; chk_8bit(*b, c);
      break;
  }
}