home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / bibli2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  46.1 KB  |  1,764 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                  BIBLIOTHEQUE  MATHEMATIQUE                    **/
  5. /**                     (deuxieme partie)                          **/
  6. /**                                                                **/
  7. /**                     Copyright Babe Cool                        **/
  8. /**                                                                **/
  9. /********************************************************************/
  10. /********************************************************************/
  11.  
  12. #include "genpari.h"
  13.  
  14. /********************************************************************/
  15. /********************************************************************/
  16. /**                                                                **/
  17. /**                        ITERATIONS                              **/
  18. /**                                                                **/
  19. /********************************************************************/
  20. /********************************************************************/
  21.  
  22. GEN forpari(ep,a,b,ch)
  23.      entree *ep;
  24.      GEN a,b;
  25.      char *ch;
  26. {
  27.   GEN p1;
  28.   long av=avma,tx=typ(a),la,lx,i,e;
  29.   if(tx>2) err(forer1);
  30.   if(tx==1)
  31.     {
  32.       la=lx=lgef(a);if(lx<3) lx=3;
  33.       if(!gcmp0(b)) {e=(gexpo(b)>>5)+3;if(e>lx) lx=e;}
  34.       p1=newbloc(lx);for(i=0;i<la;i++) p1[i]=a[i];
  35.       setlg(p1,lx);p1[-1]=(long)ep->value;ep->value=(void*)p1;
  36.     }
  37.   else newvalue(ep,a);
  38.   p1=(GEN)ep->value;
  39.   while (gcmp(p1,b)<=0)
  40.     {
  41.       lisseq(ch);gaddgsz(p1,1,p1);avma=av;
  42.     }
  43.   killvalue(ep);return gnil;
  44. }
  45.  
  46. GEN forstep(ep,a,b,s,ch)
  47.      entree *ep;
  48.      GEN a,b,s;
  49.      char *ch;
  50. {
  51.   GEN p1;
  52.   long av=avma,tx=typ(a),lx,la,s1,i,e;
  53.   s1=gsigne(s);if(!s1) err(forer2);
  54.   if(tx>2) err(forer1);
  55.   if(tx==1)
  56.     {
  57.       la=lx=lgef(a);if(lx<3) lx=3;
  58.       if(!gcmp0(b)) {e=(gexpo(b)>>5)+3;if(e>lx) lx=e;}
  59.       if(typ(s)==2) {affir(a,p1=cgetr(lg(s)));newvalue(ep,p1);avma=av;}
  60.       else
  61.     {
  62.       p1=newbloc(lx);for(i=0;i<la;i++) p1[i]=a[i];
  63.       setlg(p1,lx);p1[-1]=(long)ep->value;ep->value=(void*)p1;
  64.     }
  65.     }
  66.   else newvalue(ep,a);
  67.   p1=(GEN)ep->value;
  68.   if(s1>0)
  69.     while (gcmp(p1,b)<= 0)
  70.       {
  71.     lisseq(ch);gaddz(p1,s,p1);avma = av;
  72.       }
  73.   else
  74.     while (gcmp(p1,b)>= 0)
  75.       {
  76.     lisseq(ch);gaddz(p1,s,p1);avma = av;
  77.       }
  78.   killvalue(ep);return gnil;
  79. }
  80.  
  81. GEN forprime(ep,a,b,ch)
  82.      entree *ep;
  83.      GEN a,b;
  84.      char *ch;
  85. {
  86.   GEN p1;
  87.   long av=avma,prime=0;
  88.   byteptr p=diffptr;
  89.   
  90.   newvalue(ep,gun); p1 = (GEN)ep->value;
  91.   while((*p)&&gcmpgs(a,prime)>0) prime += *p++;
  92.   if(!*p) err(recprimer);
  93.   affsi(prime,p1);
  94.   while(gcmp(p1,b)<=0)
  95.     {
  96.       if(!*p) err(recprimer);
  97.       lisseq(ch);addsiz(*p++,p1,p1);avma=av;
  98.     }
  99.   killvalue(ep);return gnil;
  100. }
  101.  
  102. GEN fordiv(ep,a,ch)
  103.      entree *ep;
  104.      GEN a;
  105.      char *ch;
  106. {
  107.   long i,l,av2,av=avma;
  108.   GEN p1,t=divisors(a);
  109.   l=lg(t);
  110.   newvalue(ep,a); p1=(GEN)ep->value;
  111.   av2=avma;
  112.   for(i=1;i<l;i++)
  113.     {
  114.       gaffect(t[i],p1);
  115.       lisseq(ch);
  116.       avma=av2;
  117.     }
  118.   avma=av;
  119.   killvalue(ep);return gnil;
  120. }
  121.  
  122. /********************************************************************/
  123. /********************************************************************/
  124. /**                                                                **/
  125. /**       SOMMES,PRODUITS,VECTEURS,MATRICES ET RECURRENCES         **/
  126. /**                                                                **/
  127. /********************************************************************/
  128. /********************************************************************/
  129.  
  130. GEN     somme(ep,x,a,b,ch)
  131.      
  132.      entree *ep;
  133.      GEN x,a,b;
  134.      char *ch;
  135.      
  136. {
  137.   GEN   y,z,p1;
  138.   long  av = avma,tetpil,limite=(av+bot)/2;
  139.   
  140.   newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a, p1);
  141.   y = x;
  142.   tetpil = avma;
  143.   if(gcmp(a,b)>0)
  144.     {
  145.       if(!gcmp1(gsub(a,b))) err(sommeer1);
  146.       avma = av;
  147.       y=gcopy(x);
  148.     }
  149.   else
  150.     do
  151.       {
  152.     if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
  153.     z=lisexpr(ch);
  154.     tetpil=avma; y=gadd(y,z);
  155.     gaddgsz(p1,1,p1);
  156.       }
  157.   while(gcmp(p1,b)<=0);
  158.   killvalue(ep);
  159.   return gerepile(av,tetpil,y);
  160. }
  161.  
  162. GEN     suminf(ep,a,ch,prec)
  163.      
  164.      entree *ep;
  165.      GEN a;
  166.      char *ch;
  167.      long prec;
  168.      
  169. {
  170.   GEN   y,z,p1;
  171.   long  av=avma,tetpil,limite=(bot+avma)/2,fl=0;
  172.   
  173.   newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a, p1);
  174.   y=gun;
  175.   do
  176.     {
  177.       if (avma < limite) {tetpil = avma; y=gerepile(av,tetpil,gcopy(y));}
  178.       z=lisexpr(ch); y=gadd(y,z);
  179.       gaddgsz(p1,1,p1);
  180.       if((!gcmp0(z))&&(gexpo(z)>gexpo(y)-((prec-2)<<5)-5)) fl=0;
  181.       else fl++;
  182.     }
  183.   while(fl<3);
  184.   killvalue(ep);
  185.   tetpil=avma; return gerepile(av,tetpil,gsubgs(y,1));
  186. }
  187.  
  188. GEN     produit(ep,x,a,b,ch)
  189.      
  190.      entree *ep;
  191.      GEN        x,a,b;
  192.      char       *ch;
  193.      
  194. {
  195.   GEN   y,z,p1;
  196.   long  av=avma,tetpil,limite=(av+bot)/2;
  197.   
  198.   newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a,p1);
  199.   y = x;
  200.   tetpil = avma;
  201.   if(gcmp(a,b)>0)
  202.     {
  203.       if(!gcmp1(gsub(a,b))) err(proder1);
  204.       avma=av;y=gcopy(x);
  205.     }
  206.   else
  207.     do
  208.       {
  209.     if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
  210.     z=lisexpr(ch);
  211.     tetpil=avma;y=gmul(y,z);
  212.     gaddgsz(p1,1,p1);
  213.       }
  214.   while(gcmp(p1,b)<=0);
  215.   killvalue(ep);
  216.   return gerepile(av,tetpil,y);
  217. }
  218.  
  219. GEN     prodinf(ep,a,ch,prec)
  220.      
  221.      entree *ep;
  222.      GEN a;
  223.      char *ch;
  224.      long prec;
  225.      
  226. {
  227.   GEN   y,z,p1;
  228.   long  av=avma,tetpil,limite=(av+bot)/2,fl=0;
  229.   
  230.   newvalue(ep, gun); p1=(GEN)ep->value; gaffect(a,p1);
  231.   affsr(1,y=cgetr(prec));
  232.   do
  233.     {
  234.       if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
  235.       z=lisexpr(ch); y=gmul(y,z);
  236.       gaddgsz(p1,1,p1);
  237.       if(gexpo(gsubgs(z,1))>-32*(prec-2)-5) fl=0;else fl++;
  238.     }
  239.   while(fl<3);
  240.   killvalue(ep);
  241.   tetpil = avma; return gerepile(av,tetpil,gcopy(y));
  242. }
  243.  
  244. GEN     prodinf1(ep,a,ch,prec)
  245.      
  246.      entree *ep;
  247.      GEN a;
  248.      char *ch;
  249.      long prec;
  250.      
  251. {
  252.   GEN   y,z,p1,p2;
  253.   long  av=avma,tetpil,limite=(av+bot)/2,fl=0;
  254.   
  255.   newvalue(ep,gun); p1=(GEN)ep->value; gaffect(a,p1);
  256.   affsr(1,y=cgetr(prec));
  257.   do
  258.     {
  259.       if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
  260.       p2=lisexpr(ch);z=gaddsg(1,p2);
  261.       y=gmul(y,z);
  262.       gaddgsz(p1,1,p1);
  263.       if((!gcmp0(z))&&(gexpo(p2)>-32*(prec-2)-5)) fl=0;else fl++;
  264.     }
  265.   while(fl<3);
  266.   killvalue(ep);
  267.   tetpil = avma; return gerepile(av,tetpil,gcopy(y));
  268. }
  269.  
  270. GEN prodeuler(ep,a,b,ch,prec)
  271.      entree *ep;
  272.      GEN a,b;
  273.      char *ch;
  274.      long prec;
  275.      
  276. {
  277.   GEN y,z,p1;
  278.   long av=avma,tetpil,prime=0;
  279.   byteptr p=diffptr;
  280.   
  281.   newvalue(ep,gun); p1 = (GEN)ep->value;
  282.   while((*p)&&gcmpgs(a,prime)>0) prime += *p++;
  283.   if(!*p) err(recprimer);
  284.   affsi(prime,p1);affsr(1,y=cgetr(prec));
  285.   if(gcmp(p1,b)>0)
  286.     {
  287.       if(!gcmp1(gsub(a,b))) err(recer1);
  288.       tetpil=avma;
  289.       y=gcopy(y);
  290.     }
  291.   else
  292.     do
  293.       {
  294.     if(!*p) err(recprimer);
  295.     z=lisexpr(ch);
  296.     tetpil=avma; y=gmul(y,z);
  297.     addsiz(*p++,p1,p1);
  298.       }
  299.   while(gcmp(p1,b)<=0);
  300.   killvalue(ep);
  301.   return gerepile(av,tetpil,y);
  302. }
  303.  
  304. GEN     vecteur(ep,nmax,ch)
  305.      entree *ep;
  306.      GEN        nmax;
  307.      char       *ch;
  308.      
  309. {
  310.   GEN   y,p1,t;
  311.   long  i,m;
  312.   
  313.   if((typ(nmax)!=1) || (signe(nmax)<0)) err(vecer1);
  314.   m=itos(nmax);
  315.   y=cgetg(m+1 ,17);
  316.   newvalue(ep, gun); p1 = (GEN)ep->value;
  317.   for(i=1;i<=m;i++)
  318.     {
  319.       affsi(i,p1);
  320.       t=lisexpr(ch);
  321.       y[i] = isonstack(t) ? (long)t : lcopy(t);
  322.     }
  323.   killvalue(ep);
  324.   return y;
  325. }
  326.  
  327. GEN     vvecteur(ep,nmax,ch)
  328.      entree *ep;
  329.      GEN        nmax;
  330.      char       *ch;
  331.      
  332. {
  333.   GEN   y=vecteur(ep,nmax,ch);
  334.   settyp(y,18);
  335.   return y;
  336. }
  337.  
  338. GEN     matrice(ep1,ep2,nlig,ncol,ch)
  339.      entree *ep1,*ep2;
  340.      GEN        nlig,ncol;
  341.      char       *ch;
  342.      
  343. {
  344.   GEN   y,z,t,p1,p2;
  345.   long  i,j,m,n;
  346.   
  347.   if((typ(nlig)!=1) || (signe(nlig)<=0)) err(mater1);
  348.   if((typ(ncol)!=1) || (signe(ncol)<=0)) err(mater1);
  349.   n=itos(nlig); m=itos(ncol);
  350.   newvalue(ep1, gun); p1 = (GEN)ep1->value;
  351.   newvalue(ep2, gun); p2 = (GEN)ep2->value;
  352.   y=cgetg(m+1 ,19);
  353.   for(i=1;i<=m;i++)
  354.     {
  355.       affsi(i,p2);
  356.       z=cgetg(n+1 ,18);
  357.       y[i]=(long)z;
  358.       for(j=1;j<=n;j++)
  359.     {
  360.       affsi(j,p1);
  361.       t=lisexpr(ch);
  362.       z[j] = isonstack(t) ? (long)t : lcopy(t);
  363.     }
  364.     }
  365.   killvalue(ep1); killvalue(ep2); 
  366.   return y;
  367. }
  368.  
  369. GEN     divsomme(ep,num,ch)
  370.      entree *ep;
  371.      GEN        num;
  372.      char       *ch;
  373.      
  374. {
  375.   GEN   y,z,p1;
  376.   long  av=avma,tetpil,d,n,d2;
  377.   
  378.   newvalue(ep, gun); p1 = (GEN)ep->value;
  379.   n=itos(num);/* provisoire */
  380.   tetpil=av=avma;y=gzero;
  381.   for(d = d2 = 1; d2 < n; d++, d2 += d+d-1)
  382.     if (!(n%d))
  383.       {
  384.     affsi(d,p1);y=gadd(y, lisexpr(ch));
  385.     affsi(n/d,p1); z = lisexpr(ch);
  386.     tetpil=avma; y=gadd(y,z);
  387.       }
  388.   if (d2 == n)
  389.     {
  390.       affsi(d,p1); z = lisexpr(ch);
  391.       tetpil=avma; y=gadd(y,z);
  392.     }
  393.   killvalue(ep);
  394.   return gerepile(av,tetpil,y);
  395. }
  396.  
  397. /********************************************************************/
  398. /********************************************************************/
  399. /**                                                                **/
  400. /**                   CALCUL D'UNE INTEGRALE                       **/
  401. /**                   ( Methode de ROMBERG )                       **/
  402. /**                                                                **/
  403. /********************************************************************/
  404. /********************************************************************/
  405.  
  406. GEN     qromb(ep,a,b,ch,prec)
  407.      entree *ep;
  408.      GEN a,b;
  409.      char *ch;
  410.      long prec;
  411.      
  412. #define JMAX 25
  413. #define JMAXP JMAX+3
  414. #define KLOC 5
  415.      
  416. {
  417.   GEN     ss,dss,s,h,q1,p1,p2,p3,p4,p,qlint,del,sz,x,sum;
  418.   long    av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
  419.   
  420.   av=avma;l=prec;
  421.   if(typ(a)!=2) { p=cgetr(prec);gaffect(a,p);a=p;}
  422.   if(typ(b)!=2) { p=cgetr(prec);gaffect(b,p);b=p;}
  423.   qlint=subrr(b,a);sig=signe(qlint);
  424.   if(!sig) {avma=av;return gzero;}
  425.   newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
  426.   if(sig<0) {setsigne(qlint,1);s=a;a=b;b=s;}
  427.   p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
  428.   s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
  429.   affsr(1,h[1]=lgetr(l));
  430.   gaffect(a,q1); p1=lisexpr(ch); if (!isonstack(p1)) p1=gcopy(p1);
  431.   gaffect(b,q1); p2=lisexpr(ch);
  432.   s[1]=lmul2n(gmul(qlint,gadd(p1,p2)),-1);it=1;sz=gmul(gzero,s[1]);
  433.   s[2]=s[1];h[2]=lshiftr(h[1],-2);
  434.   for(j=2;j<=JMAX;j++)
  435.     {
  436.       tnm=it;del=divrs(qlint,tnm);x=addrr(a,shiftr(del,-1));
  437.       for(sum=sz,j1=1;j1<=it;j1++,x=addrr(x,del))
  438.     {
  439.       affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);
  440.     }
  441.       it*=2;
  442.       s[j]=lmul2n(gadd(s[j-1],gmul(sum,del)),-1);
  443.       if(j>=KLOC)
  444.     {
  445.       for(j1=1;j1<=KLOC;j1++)
  446.         {
  447.           p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
  448.         }
  449.       tetpil=avma;ss=polint(p4,p3,gzero,&dss);
  450.       j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-j-5;
  451.       if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
  452.         {
  453.           if(gcmp0(gimag(ss))) ss=greal(ss);
  454.           tetpil=avma;
  455.           killvalue(ep);
  456.           return gerepile(av,tetpil,gmulsg(sig,ss));
  457.         }
  458.     }
  459.       s[j+1]=s[j];h[j+1]=lshiftr(h[j],-2);
  460.     }
  461.   err(intger2);
  462. }         
  463.  
  464. GEN     qromo(ep,a,b,ch,prec)
  465.      entree *ep;
  466.      GEN a,b;
  467.      char *ch;
  468.      long prec;
  469.      
  470. #undef JMAX
  471. #define JMAX 16
  472. #define JMAXP JMAX+3
  473. #define KLOC 5
  474.      
  475. {
  476.   GEN     ss,dss,s,h,q1,sz,p1,p3,p4,p,qlint,del,ddel,x,sum;
  477.   long    av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
  478.   
  479.   av=avma;l=prec;
  480.   if(typ(a)!=2) { p=cgetr(prec);gaffect(a,p);a=p;}
  481.   if(typ(b)!=2) { p=cgetr(prec);gaffect(b,p);b=p;}
  482.   qlint=subrr(b,a);sig=signe(qlint);
  483.   if(!sig) {avma=av;return gzero;}
  484.   if(sig<0) {setsigne(qlint,1);s=a;a=b;b=s;}
  485.   newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
  486.   p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
  487.   s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
  488.   affsr(1,h[1]=lgetr(l));affrr(shiftr(addrr(a,b),-1),q1);p1=lisexpr(ch);
  489.   s[1]=lmul(qlint,p1);it=1;sz=gmul(gzero,s[1]);
  490.   s[2]=s[1];h[2]=ldivrs(h[1],9);
  491.   for(j=2;j<=JMAX;j++)
  492.     {
  493.       tnm=it;del=divrs(qlint,3*tnm);ddel=shiftr(del,1);
  494.       x=addrr(a,shiftr(del,-1));sum=sz;
  495.       for(j1=1;j1<=it;j1++)
  496.     {
  497.       affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);x=addrr(x,ddel);
  498.       affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);x=addrr(x,del);
  499.     }
  500.       it*=3;
  501.       s[j]=ladd(gdivgs(s[j-1],3),gmul(sum,del));
  502.       if(j>=KLOC)
  503.     {
  504.       for(j1=1;j1<=KLOC;j1++)
  505.         {
  506.           p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
  507.         }
  508.       tetpil=avma;ss=polint(p4,p3,gzero,&dss);
  509.       j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-(3*j/2)-5;
  510.       if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
  511.         {
  512.           if(gcmp0(gimag(ss))) ss=greal(ss);
  513.           tetpil=avma; killvalue(ep);
  514.           return gerepile(av,tetpil,gmulsg(sig,ss));
  515.         }
  516.     }
  517.       s[j+1]=s[j];h[j+1]=ldivrs(h[j],9);
  518.     }
  519.   err(intger2);
  520. }         
  521.  
  522. GEN     qromi(ep,a,b,ch,prec)
  523.      entree *ep;
  524.      GEN a,b;
  525.      char *ch;
  526.      long prec;
  527.      
  528. #undef JMAX
  529. #define JMAX 16
  530. #define JMAXP JMAX+3
  531. #define KLOC 5
  532.      
  533. {
  534.   GEN     ss,dss,s,h,q1,sz,p1,p3,p4,p,qlint,del,ddel,x,sum;
  535.   long    av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
  536.   
  537.   av=avma;l=prec;
  538.   p=cgetr(prec);gaffect(ginv(a),p);a=p;
  539.   p=cgetr(prec);gaffect(ginv(b),p);b=p;
  540.   qlint=subrr(b,a);sig= -signe(qlint);
  541.   if(!sig) {avma=av;return gzero;}
  542.   if(sig>0) {setsigne(qlint,1);s=a;a=b;b=s;}
  543.   newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
  544.   p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
  545.   s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
  546.   affsr(1,h[1]=lgetr(l));x=divsr(2,addrr(a,b));
  547.   affrr(x,q1);
  548.   p1=gmul(lisexpr(ch),mulrr(x,x));s[1]=lmul(qlint,p1);it=1;
  549.   sz=gmul(gzero,s[1]);s[2]=s[1];h[2]=ldivrs(h[1],9);
  550.   for(j=2;j<=JMAX;j++)
  551.     {
  552.       tnm=it;del=divrs(qlint,3*tnm);ddel=shiftr(del,1);
  553.       x=addrr(a,shiftr(del,-1));sum=sz;
  554.       for(j1=1;j1<=it;j1++)
  555.     {
  556.       divsrz(1,x,q1);p1=gmul(lisexpr(ch),mulrr(q1,q1));
  557.       sum=gadd(sum,p1);x=addrr(x,ddel);
  558.       divsrz(1,x,q1);p1=gmul(lisexpr(ch),mulrr(q1,q1));
  559.       sum=gadd(sum,p1);x=addrr(x,del);
  560.     }
  561.       it*=3;
  562.       s[j]=ladd(gdivgs(s[j-1],3),gmul(sum,del));
  563.       if(j>=KLOC)
  564.     {
  565.       for(j1=1;j1<=KLOC;j1++)
  566.         {
  567.           p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
  568.         }
  569.       tetpil=avma;ss=polint(p4,p3,gzero,&dss);
  570.       j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-(3*j/2)-5;
  571.       if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
  572.         {
  573.           if(gcmp0(gimag(ss))) ss=greal(ss);
  574.           tetpil=avma; killvalue(ep);
  575.           return gerepile(av,tetpil,gmulsg(sig,ss));
  576.         }
  577.     }
  578.       s[j+1]=s[j];h[j+1]=ldivrs(h[j],9);
  579.     }
  580.   err(intger2);
  581. }         
  582.  
  583. GEN     rombint(ep,a,b,ch,prec)
  584.      entree *ep;
  585.      GEN a,b;
  586.      char *ch;
  587.      long prec;
  588.      
  589. {
  590.   GEN     s,p1,p2,p3;
  591.   static  long p4[3]={0x1ff0003,0xff000003,1};
  592.   long    l,av,tetpil;
  593.   
  594.   l=gcmp(b,a);
  595.   if(!l) return gzero;
  596.   if(l<0) {s=a;a=b;b=s;}
  597.   av=avma;
  598.   if(gcmpgs(b,100)>=0)
  599.     {
  600.       if(gcmpgs(a,1)>=0) return qromi(ep,a,b,ch,prec);
  601.       p1=qromi(ep,gun,b,ch,prec);
  602.       if(gcmpgs(a,-100)>=0)
  603.     {
  604.       p2=qromo(ep,a,gun,ch,prec);tetpil=avma;
  605.       return gerepile(av,tetpil,gadd(p1,p2));
  606.     }
  607.       p2=qromo(ep,p4,gun,ch,prec);p3=gadd(p2,qromi(ep,a,p4,ch,prec));
  608.       tetpil=avma;return gerepile(av,tetpil,gadd(p1,p3));
  609.     }
  610.   if(gcmpgs(a,-100)>=0) return qromo(ep,a,b,ch,prec);
  611.   if(gcmpgs(b,-1)>=0)
  612.     {
  613.       p1=qromi(ep,a,p4,ch,prec);p2=qromo(ep,p4,b,ch,prec);tetpil=avma;
  614.       return gerepile(av,tetpil,gadd(p1,p2));
  615.     }
  616.   return qromi(ep,a,b,ch,prec);
  617. }         
  618.  
  619. /********************************************************************/
  620. /********************************************************************/
  621. /**                                                                **/
  622. /**                    SOMMATION DE SERIES                         **/
  623. /**                                                                **/
  624. /********************************************************************/
  625. /********************************************************************/
  626.  
  627. void eulsum(sum,term,jterm,tab,dsum,prec)
  628.      GEN *sum,term,tab[];
  629.      long jterm,prec,*dsum;
  630.      
  631. {
  632.   long j;
  633.   static long nterm;
  634.   GEN  tmp,dum,p1;
  635.   static GEN unreel;
  636.   
  637.   if(jterm==1)
  638.     {
  639.       nterm=1;tab[1]=term;*sum=gmul2n(term,-1);affsr(1,unreel=cgetr(prec));
  640.     }
  641.   else
  642.     {
  643.       tmp=tab[1];tab[1]=gmul(unreel,term);
  644.       for(j=1;j<nterm;j++)
  645.     {
  646.       dum=tab[j+1];tab[j+1]=gmul2n(gadd(tab[j],tmp),-1);tmp=dum;
  647.     }
  648.       tab[nterm+1]=gmul2n(gadd(tab[nterm],tmp),-1);
  649.       if(gcmp(gabs(tab[nterm+1],prec),gabs(tab[nterm],prec))<=0)
  650.     p1=gmul2n(tab[++nterm],-1);
  651.       else p1=tab[nterm+1];
  652.       *sum=gadd(*sum,p1);*dsum=gexpo(p1);
  653.     }
  654. }
  655.  
  656. GEN  sumalt(ep,a,ch,prec)
  657.      entree *ep;
  658.      GEN a;
  659.      char *ch;
  660.      long prec;
  661.  
  662. #define JMIT 10000
  663.      
  664. {
  665.   long av,tetpil,j,nterm,jterm;
  666.   GEN  p1,sum,sum0,q1,tmp,dum,*tab,unreel;
  667.   
  668.   newvalue(ep,gun); q1=(GEN)ep->value; gaffect(a,q1);
  669.   tab=(GEN *)newbloc(JMIT);
  670.   av=avma;
  671.   p1=lisexpr(ch);affsr(1,unreel=cgetr(prec));
  672.   nterm=1;tab[1]=p1;sum=gmul2n(p1,-1);
  673.   for(jterm=2;jterm<=JMIT;jterm++)
  674.     {
  675.       tmp=tab[1];a=gaddsg(1,a);gaffect(a,q1);tab[1]=gmul(unreel,lisexpr(ch));
  676.       for(j=1;j<nterm;j++)
  677.     {
  678.       dum=tab[j+1];tab[j+1]=gmul2n(gadd(tab[j],tmp),-1);tmp=dum;
  679.     }
  680.       tab[nterm+1]=gmul2n(gadd(tab[nterm],tmp),-1);
  681.       if(gcmp(gabs(tab[nterm+1],prec),gabs(tab[nterm],prec))<=0)
  682.     p1=gmul2n(tab[++nterm],-1);
  683.       else p1=tab[nterm+1];
  684.       sum0=sum;sum=gadd(sum,p1);
  685.       if((gcmp0(p1)||(gexpo(p1)<-32*(prec-2)+5)||(gegal(sum,sum0)))&&(jterm>=10))
  686.     {
  687.       tetpil=avma;
  688.       killbloc((GEN)tab);
  689.       killvalue(ep);
  690.       return gerepile(av,tetpil,gcopy(sum));
  691.     }
  692.     }
  693.   err(eulsumer1);
  694. }
  695.  
  696. GEN  sumpos(ep,a,ch,prec)
  697.      entree *ep;
  698.      GEN a;
  699.      char *ch;
  700.      long prec;
  701.  
  702. {
  703.   long av,tetpil,k,jterm,dsum;
  704.   GEN  sum,term,q1,p1,*tab,unreel,r;
  705.   
  706.   tab=(GEN *)newbloc(JMIT);
  707.   av = avma; a=gsubgs(a,1);
  708.   affsr(1,unreel=cgetr(prec));term=gzero;r=gun;k=0;
  709.   p1=cgeti(prec+1);p1[1]=0x1000001+prec;
  710.   newvalue(ep,p1); q1=(GEN)ep->value;
  711.   do
  712.     {
  713.       gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
  714.       term=gadd(term,p1);r=shifti(r,1);k++;
  715.     }
  716.   while(gexpo(p1)>=(-32*(prec-2)+5));
  717.   sum=gzero;
  718.   eulsum(&sum,term,1,tab,&dsum,prec);
  719.   for(jterm=2;jterm<=JMIT;jterm++)
  720.     {
  721.       term=gzero;r=stoi(jterm);k=0;
  722.       do
  723.     {
  724.       gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
  725.       term=gadd(term,p1);r=shifti(r,1);k++;
  726.     }
  727.       while(gexpo(p1)>=(-32*(prec-2)+5));
  728.       if(!odd(jterm)) term=gneg(term);
  729.       eulsum(&sum,term,jterm,tab,&dsum,prec);
  730.       if(dsum< -32*(prec-2)+5)
  731.     {
  732.       tetpil=avma;
  733.       killbloc((GEN)tab);
  734.       killvalue(ep);
  735.       return gerepile(av,tetpil,gcopy(sum));
  736.     }
  737.     }
  738.   killbloc((GEN)tab);
  739.   err(eulsumer1);
  740. }
  741.  
  742. /********************************************************************/
  743. /********************************************************************/
  744. /**                                                                **/
  745. /**                        TRACE GROSSIER                          **/
  746. /**                                                                **/
  747. /********************************************************************/
  748. /********************************************************************/
  749.  
  750.  
  751. GEN plot(ep,a,b,ch)
  752.      entree *ep;
  753.      GEN a,b;
  754.      char *ch;
  755.      
  756. #define ISCR 60
  757. #define JSCR 21
  758. #define BLANK ' '
  759. #define ZERO '-'
  760. #define YY '|'
  761. #define XX '-'
  762. #define FF 'x'
  763.      
  764. {
  765.   long av,av2,jz,j,i,sig;
  766.   GEN p1,p2,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
  767.   char scr[ISCR+1][JSCR+1], thestring[100];
  768.   
  769.   sig=gcmp(b,a);if(!sig) return gnil;
  770.   av=avma; pariputc('\n');
  771.   if(sig<0) {x=a;a=b;b=x;}
  772.   newvalue(ep,cgetr(3)); x=(GEN)ep->value;
  773.   for(i=1;i<=ISCR;i++) y[i]=cgetr(3);
  774.   gaffect(a,x);
  775.   dx=gdivgs(gsub(b,a),(ISCR-1));ysml=gzero;ybig=gzero;
  776.   for(j=1;j<=JSCR;j++) scr[1][j]=scr[ISCR][j]=YY;
  777.   for(i=2;i<ISCR;i++)
  778.     {
  779.       scr[i][1]=scr[i][JSCR]=XX;
  780.       for(j=2;j<JSCR;j++) scr[i][j]=BLANK;
  781.     }
  782.   av2=avma;
  783.   for(i=1;i<=ISCR;i++)
  784.     {
  785.       gaffect(lisexpr(ch),y[i]);
  786.       if(gcmp(y[i],ysml)<0) ysml=y[i];
  787.       if(gcmp(y[i],ybig)>0) ybig=y[i];
  788.       gaddz(x,dx,x);avma=av2;
  789.     }
  790.   diff=gsub(ybig,ysml);
  791.   if(gcmp0(diff)) {ybig=gaddsg(1,ybig);diff=gun;}
  792.   dyj=gdivsg(JSCR-1,diff);jz=1-gtolong(gmul(ysml,dyj));
  793.   av2=avma;
  794.   for(i=1;i<=ISCR;i++)
  795.     {
  796.       scr[i][jz]=ZERO;j=1+gtolong(gmul(gsub(y[i],ysml),dyj));
  797.       scr[i][j]=FF;avma=av2;
  798.     }
  799.   p1=cgetr(4);gaffect(ybig,p1);
  800.   sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
  801.   for(i=1;i<=ISCR;i++)  pariputc(scr[i][JSCR]); pariputc('\n');
  802.   for(j=(JSCR-1);j>=2;j--)
  803.     {
  804.       pariputs("            ");
  805.       for(i=1;i<=ISCR;i++) pariputc(scr[i][j]);
  806.       pariputc('\n');
  807.     }
  808.   p1=cgetr(4);gaffect(ysml,p1);
  809.   sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
  810.   for(i=1;i<=ISCR;i++)  pariputc(scr[i][1]);
  811.   pariputc('\n');
  812.   p1=cgetr(4);p2=cgetr(4);gaffect(a,p1);gaffect(b,p2);
  813.   sprintf(thestring, "%8s %10.3lf %44s %10.3lf\n"," ",rtodbl(p1)," ",rtodbl(p2));
  814.   pariputs(thestring);
  815.   killvalue(ep);
  816.   avma=av; pariputc('\n'); return gnil;
  817. }
  818.  
  819. /********************************************************************/
  820. /********************************************************************/
  821. /**                                                                **/
  822. /**                  RECHERCHE DE ZEROS REELS                      **/
  823. /**                                                                **/
  824. /********************************************************************/
  825. /********************************************************************/
  826.  
  827. GEN zbrent(ep,a,b,ch,prec)
  828.      entree *ep;
  829.      GEN a,b;
  830.      char *ch;
  831.      long prec;
  832.      
  833. {
  834.   long av=avma,tetpil,sig,iter,itmax;
  835.   GEN q1,c,d,e,tol,toli,min1,min2,fa,fb,fc,p,q,r,s,xm;
  836.   
  837.   if(typ(a)!=2) {p=cgetr(prec);gaffect(a,p);a=p;}
  838.   if(typ(b)!=2) {p=cgetr(prec);gaffect(b,p);b=p;}
  839.   sig=cmprr(b,a);if(!sig) {avma = av; return gzero;}
  840.   if(sig<0) {c=a;a=b;b=c;}
  841.   newvalue(ep,a); fa=lisexpr(ch);
  842.   changevalue(ep,b); q1=(GEN)ep->value; fb=lisexpr(ch);
  843.   if(gsigne(fa)*gsigne(fb)>0) err(zbrenter1);
  844.   itmax=32*prec+50;affsr(1,tol=cgetr(3));tol=shiftr(tol,-(32*(prec-2)-5));
  845.   fc=fb;
  846.   for(iter=1;iter<=itmax;iter++)
  847.     {
  848.       if(gsigne(fb)*gsigne(fc)>0)
  849.     {
  850.       c=a;fc=fa;d=subrr(b,a);e=d;
  851.     }
  852.       if(gcmp(gabs(fc),gabs(fb))<0)
  853.     {
  854.       a=b;b=c;c=a;fa=fb;fb=fc;fc=fa;
  855.     }
  856.       toli=mulrr(tol,absr(b));
  857.       xm=shiftr(subrr(c,b),-1);
  858.       if((cmprr(absr(xm),toli)<=0)||gcmp0(fb))
  859.     {
  860.       tetpil=avma;
  861.       killvalue(ep);
  862.       return gerepile(av,tetpil,gcopy(b));
  863.     }
  864.       if((cmprr(absr(e),toli)>=0)&&(gcmp(gabs(fa),gabs(fb))>0))
  865.     {
  866.       s=gdiv(fb,fa);
  867.       if(cmprr(a,b)==0)
  868.         {
  869.           p=gmul2n(gmul(xm,s),1);q=gsubsg(1,s);
  870.         }
  871.       else
  872.         {
  873.           q=gdiv(fa,fc);r=gdiv(fb,fc);
  874.           p=gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
  875.           p=gmul(s,gsub(p,gmul(gsub(b,a),gsubgs(r,1))));
  876.           q=gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
  877.         }
  878.       if(gsigne(p)>0) q=gneg(q);
  879.       p=gabs(p);
  880.       min1=gsub(gmulsg(3,gmul(xm,q)),gabs(gmul(q,toli)));
  881.       min2=gabs(gmul(e,q));
  882.       if(gcmp(gmul2n(p,1),gmin(min1,min2))<0) { e=d;d=gdiv(p,q);}
  883.       else {d=xm;e=d;}
  884.     }
  885.       else {d=xm;e=d;}
  886.       a=b;fa=fb;
  887.       if(gcmp(gabs(d),toli)>0) b=addrr(b,d);
  888.       else
  889.     {
  890.       if(gsigne(xm)>0) b=addrr(b,absr(toli));
  891.       else b=subrr(b,absr(toli));
  892.     }
  893.       gaffect(b,q1); ;fb=lisexpr(ch);
  894.     }
  895.   err(zbrenter2);
  896. }
  897.  
  898. GEN    lllgen(x)
  899.      GEN   x;
  900.  
  901. {
  902.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  903.   GEN g;
  904.  
  905.   if(tx!=19) err(lller1);
  906.   av=avma;
  907.   g=cgetg(lx,19);
  908.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  909.   for(i=1;i<lx;i++)
  910.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  911.   av1=avma;return gerepile(av,av1,lllgramallgen(g,2));
  912. }
  913.  
  914. GEN lllkerimgen(x)
  915.      GEN x;
  916. {
  917.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  918.   GEN g;
  919.  
  920.   if(tx!=19) err(lller1);
  921.   av=avma;
  922.   g=cgetg(lx,19);
  923.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  924.   for(i=1;i<lx;i++)
  925.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  926.   av1=avma;return gerepile(av,av1,lllgramallgen(g,0));
  927. }
  928.  
  929. GEN lllgramgen(x)
  930.      GEN x;
  931. {
  932.   return lllgramallgen(x,2);
  933. }
  934.  
  935. GEN lllgramkerimgen(x)
  936.      GEN x;
  937. {
  938.   return lllgramallgen(x,0);
  939. }
  940.  
  941. int pslg(x)
  942.      GEN x;
  943. {
  944.   if(gcmp0(x)) return 2;
  945.   return (typ(x)<10)?3:lgef(x);
  946. }
  947.  
  948. GEN lllgramallgen(x,all)
  949.      GEN x;
  950.      long all;
  951. {
  952.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec;
  953.   GEN u,B,lam,q,cq,h,la,bb,p2,p3,p4,y,fl;
  954.   int ps1,ps2,flc;
  955.  
  956.   if(tx!=19) err(lllger1);
  957.   n=lx-1;if(n<=1) return idmat(n);
  958.   if(lg(x[1])!=lx) err(lllger2);
  959.   av=avma;lim=(avma+bot)>>1;
  960.   B=cgetg(lx+1,18);
  961.   fl=cgetg(lx,17);
  962.   B[1]=un;lam=cgetg(lx,19);
  963.   for(j=1;j<lx;j++) lam[j]=lgetg(lx,18);
  964.   for(i=1;i<lx;i++) 
  965.     for(j=1;j<=i;j++)
  966.       {
  967.     if((j<i)&&(!signe(fl[j]))) coeff(lam,i,j)=coeff(lam,j,i)=zero;
  968.     else
  969.       {
  970.         u=(GEN)coeff(x,i,j);
  971.         if(typ(u)>10) err(lllger4);
  972.         for(k=1;k<j;k++)
  973.           if(signe(fl[k])) u=gdiv(gsub(gmul(B[k+1],u),gmul(coeff(lam,i,k),coeff(lam,j,k))),B[k]);
  974.         if(j<i) {coeff(lam,i,j)=(long)u;coeff(lam,j,i)=zero;}
  975.         else 
  976.           {
  977.         if(!gcmp0(u)) {B[i+1]=(long)u;coeff(lam,i,i)=fl[i]=un;}
  978.         else {B[i+1]=B[i];coeff(lam,i,i)=fl[i]=zero;}
  979.           }
  980.       }
  981.       }
  982.   k=2;h=idmat(n);flc=0;
  983.   for(;;)
  984.     {
  985.       u=(GEN)coeff(lam,k,k-1);
  986.       if(pslg(u)>=pslg(B[k]))
  987.     {
  988.       q=gdeuc(u,B[k]);cq=gdivsg(1,content(q));q=gmul(q,cq);flc=1;
  989.       h[k]=lsub(gmul(cq,h[k]),gmul(q,h[k-1]));
  990.       coeff(lam,k,k-1)=lsub(gmul(cq,coeff(lam,k,k-1)),gmul(q,B[k]));
  991.       for(i=1;i<k-1;i++) 
  992.         coeff(lam,k,i)=lsub(gmul(cq,coeff(lam,k,i)),gmul(q,coeff(lam,k-1,i)));
  993.     }
  994.       ps1=pslg(gmul(B[k],B[k]));
  995.       ps2=pslg(gadd(p3=gmul(B[k-1],B[k+1]),p4=gmul(la=(GEN)coeff(lam,k,k-1),coeff(lam,k,k-1))));
  996.       if(signe(fl[k-1])&&((ps1>ps2)||((ps1==ps2)&&flc)||(!signe(fl[k]))))
  997.     {
  998.       p1=h[k-1];h[k-1]=h[k];h[k]=p1;if((ps1==ps2)&&flc) flc=0;else flc=1;
  999.       for(j=1;j<=k-2;j++) 
  1000.         {
  1001.           p1=coeff(lam,k-1,j);coeff(lam,k-1,j)=coeff(lam,k,j);
  1002.           coeff(lam,k,j)=p1;
  1003.         }
  1004.       if(signe(fl[k]))
  1005.         {
  1006.           for(i=k+1;i<=n;i++)
  1007.         {
  1008.           bb=(GEN)coeff(lam,i,k);
  1009.           coeff(lam,i,k)=ldiv(gsub(gmul(B[k+1],coeff(lam,i,k-1)),gmul(la,bb)),B[k]);
  1010.           coeff(lam,i,k-1)=ldiv(gadd(gmul(la,coeff(lam,i,k-1)),gmul(B[k-1],bb)),B[k]);
  1011.         }
  1012.           B[k]=ldiv(gadd(p3,p4),B[k]);
  1013.         }
  1014.       else
  1015.         {
  1016.           if(!gcmp0(la))
  1017.         {
  1018.           p2=(GEN)B[k];p1=ldiv(p4,p2);
  1019.           for(i=k+1;i<lx;i++)
  1020.             coeff(lam,i,k-1)=ldiv(gmul(la,coeff(lam,i,k-1)),p2);
  1021.           for(j=k+1;j<lx-1;j++)
  1022.             for(i=j+1;i<lx;i++)
  1023.               coeff(lam,i,j)=ldiv(gmul(p1,coeff(lam,i,j)),p2);
  1024.           B[k+1]=B[k]=p1;
  1025.           for(i=k+2;i<=lx;i++)
  1026.             B[i]=ldiv(gmul(p1,B[i]),p2);
  1027.         }
  1028.           else
  1029.         {
  1030.           coeff(lam,k,k-1)=zero;
  1031.           for(i=k+1;i<lx;i++)
  1032.             {coeff(lam,i,k)=coeff(lam,i,k-1);coeff(lam,i,k-1)=zero;}
  1033.           B[k]=B[k-1];fl[k]=un;fl[k-1]=zero;
  1034.         }
  1035.         }
  1036.       if(k>2) k--;
  1037.     }
  1038.       else
  1039.     {
  1040.       for(l=k-2;l>=1;l--)
  1041.         {
  1042.           u=(GEN)coeff(lam,k,l);
  1043.           if(pslg(u)>=pslg(B[l+1]))
  1044.         {
  1045.           q=gdeuc(u,B[l+1]);cq=gdivsg(1,content(q));q=gmul(q,cq);flc=1;
  1046.           h[k]=lsub(gmul(cq,h[k]),gmul(q,h[l]));
  1047.           coeff(lam,k,l)=lsub(gmul(cq,coeff(lam,k,l)),gmul(q,B[l+1]));
  1048.           for(i=1;i<l;i++) coeff(lam,k,i)=lsub(gmul(cq,coeff(lam,k,i)),gmul(q,coeff(lam,l,i)));
  1049.         }
  1050.         }
  1051.       k++;
  1052.       if(k>n) 
  1053.         {
  1054.           for(k=1;(k<=n)&&(!signe(fl[k]));k++);
  1055.           tetpil=avma;
  1056.           if(!all)
  1057.         {
  1058.           y=cgetg(3,17);
  1059.           p2=cgetg(k,19);for(i=1;i<k;i++) p2[i]=lcopy(h[i]);
  1060.           y[1]=(long)p2;p2=cgetg(n-k+2,19);y[2]=(long)p2;
  1061.           for(i=k;i<=n;i++) p2[i-k+1]=lcopy(h[i]);
  1062.         }
  1063.           else
  1064.         {
  1065.           if(all==1)
  1066.             {
  1067.               y=cgetg(k,19);for(i=1;i<k;i++) y[i]=lcopy(h[i]);
  1068.             }
  1069.           else
  1070.             {
  1071.               y=cgetg(n-k+2,19);
  1072.               for(i=k;i<=n;i++) y[i-k+1]=lcopy(h[i]);
  1073.             }
  1074.         }
  1075.           return gerepile(av,tetpil,y);
  1076.         }
  1077.     }
  1078.       if(avma<lim)
  1079.     {
  1080.       tetpil=avma;
  1081.       B=gcopy(B);h=gcopy(h);lam=gcopy(lam);fl=gcopy(fl);
  1082.       dec=lpile(av,tetpil,0)>>2;
  1083.       B+=dec;h+=dec;lam+=dec;fl+=dec;
  1084.     }
  1085.     }
  1086. }
  1087.  
  1088. /********************************************************************/
  1089. /********************************************************************/
  1090. /**                                                                **/
  1091. /**                   FONCTIONS DE RECTPLOT                        **/
  1092. /**                                                                **/
  1093. /********************************************************************/
  1094. /********************************************************************/
  1095.  
  1096. GEN initrect(ne,x,y)
  1097.      long ne,x,y;
  1098.  
  1099. {
  1100.   long *e;
  1101.  
  1102.   if((ne<0)||(ne>15)) err(rploter2);
  1103.   if((x<=1)||(y<=1)) err(rploter1);
  1104.   e=rectgraph[ne];e[0]=e[1]=e[4]=e[5]=0;e[2]=x;e[3]=y;
  1105.   return gnil;
  1106. }
  1107.  
  1108. GEN rectcursor(ne)
  1109.      long ne;
  1110. {
  1111.   GEN z;long *e;
  1112.  
  1113.   if((ne<0)||(ne>15)) err(rploter2);
  1114.   e=rectgraph[ne];
  1115.   z=cgetg(3,17);z[1]=lstoi(e[4]);z[2]=lstoi(e[5]);
  1116.   return z;
  1117. }
  1118.  
  1119. GEN rectmove(ne,x,y) /* code = 0 */
  1120.      long ne,x,y;
  1121. {
  1122.   long *e,*z;
  1123.  
  1124.   if((ne<0)||(ne>15)) err(rploter2);
  1125.   e=rectgraph[ne];
  1126.   if(!(z=(long*)malloc(16))) err(memer);
  1127.   z[0]=z[1]=0;e[4]=z[2]=x;e[5]=z[3]=y;
  1128.   if(!e[0]) e[0]=e[1]=(long)z;
  1129.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1130.   return gnil;
  1131. }
  1132.  
  1133. GEN rectrmove(ne,x,y) /* code = 0 */
  1134.      long ne,x,y;
  1135. {
  1136.   long *e,*z;
  1137.  
  1138.   if((ne<0)||(ne>15)) err(rploter2);
  1139.   e=rectgraph[ne];
  1140.   if(!(z=(long*)malloc(16))) err(memer);
  1141.   z[0]=z[1]=0;e[4]=z[2]=x+e[4];e[5]=z[3]=y+e[5];
  1142.   if(!e[0]) e[0]=e[1]=(long)z;
  1143.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1144.   return gnil;
  1145. }
  1146.  
  1147. GEN rectpoint(ne,x,y) /* code = 1 */
  1148.      long ne,x,y;
  1149. {
  1150.   long *e,*z;
  1151.  
  1152.   if((ne<0)||(ne>15)) err(rploter2);
  1153.   e=rectgraph[ne];
  1154.   if(!(z=(long*)malloc(16))) err(memer);
  1155.   z[0]=0;e[4]=z[2]=x;e[5]=z[3]=y;
  1156.   z[1]=((x<0)||(y<0)||(x>e[2])||(y>e[3]))?0:1;
  1157.   if(!e[0]) e[0]=e[1]=(long)z;
  1158.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1159.   return gnil;
  1160. }
  1161.  
  1162. GEN rectrpoint(ne,x,y) /* code = 1 */
  1163.      long ne,x,y;
  1164. {
  1165.   long *e,*z,a,b;
  1166.  
  1167.   if((ne<0)||(ne>15)) err(rploter2);
  1168.   e=rectgraph[ne];
  1169.   if(!(z=(long*)malloc(16))) err(memer);
  1170.   z[0]=0;a=e[4]=z[2]=x+e[4];b=e[5]=z[3]=y+e[5];
  1171.   z[1]=((a<0)||(b<0)||(a>e[2])||(b>e[3]))?0:1;
  1172.   if(!e[0]) e[0]=e[1]=(long)z;
  1173.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1174.   return gnil;
  1175. }
  1176.  
  1177. GEN rectline(ne,x1,y1,x2,y2) /* code = 2 */
  1178.      long ne,x1,y1,x2,y2;
  1179. {
  1180.   long *e,*z,dx,dy,dxy,xmin,xmax,ymin,ymax;
  1181.  
  1182.   if((ne<0)||(ne>15)) err(rploter2);
  1183.   e=rectgraph[ne];
  1184.   if(!(z=(long*)malloc(24))) err(memer);
  1185.   xmin=max(min(x1,x2),0);xmax=min(max(x1,x2),e[2]);
  1186.   ymin=max(min(y1,y2),0);ymax=min(max(y1,y2),e[3]);
  1187.   dxy=x1*y2-y1*x2;dx=x2-x1;dy=y2-y1;
  1188.   if(dy)
  1189.     {
  1190.       if(dx*dy<0) {xmin=max(xmin,(dxy+e[3]*dx)/dy);xmax=min(xmax,dxy/dy);}
  1191.       else {xmin=max(xmin,dxy/dy);xmax=min(xmax,(dxy+e[3]*dx)/dy);}
  1192.     }
  1193.   if(dx)
  1194.     {
  1195.       if(dx*dy<0) {ymin=max(ymin,(e[2]*dy-dxy)/dx);ymax=min(ymax,-dxy/dx);}
  1196.       else {ymin=max(ymin,-dxy/dx);ymax=min(ymax,(e[2]*dy-dxy)/dx);}
  1197.     }
  1198.   z[0]=0;z[2]=xmin;z[4]=xmax;e[4]=x2;e[5]=y2;
  1199.   if(dx*dy<0) {z[3]=ymax;z[5]=ymin;}
  1200.   else {z[3]=ymin;z[5]=ymax;}
  1201.   z[1]=((xmin>xmax)||(ymin>ymax))?0:2;
  1202.   if(!e[0]) e[0]=e[1]=(long)z;
  1203.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1204.   return gnil;
  1205. }
  1206.  
  1207. GEN rectrline(ne,x2,y2) /* code = 2 */
  1208.      long ne,x2,y2;
  1209. {
  1210.   long *e,*z,x1,y1,dx,dy,dxy,xmin,xmax,ymin,ymax;
  1211.  
  1212.   if((ne<0)||(ne>15)) err(rploter2);
  1213.   e=rectgraph[ne];
  1214.   if(!(z=(long*)malloc(24))) err(memer);
  1215.   x1=e[4];y1=e[5];
  1216.   xmin=max(min(x1,x2),0);xmax=min(max(x1,x2),e[2]);
  1217.   ymin=max(min(y1,y2),0);ymax=min(max(y1,y2),e[3]);
  1218.   dxy=x1*y2-y1*x2;dx=x2-x1;dy=y2-y1;
  1219.   if(dy)
  1220.     {
  1221.       if(dx*dy<0) {xmin=max(xmin,(dxy+e[3]*dx)/dy);xmax=min(xmax,dxy/dy);}
  1222.       else {xmin=max(xmin,dxy/dy);xmax=min(xmax,(dxy+e[3]*dx)/dy);}
  1223.     }
  1224.   if(dx)
  1225.     {
  1226.       if(dx*dy<0) {ymin=max(ymin,(e[2]*dy-dxy)/dx);ymax=min(ymax,-dxy/dx);}
  1227.       else {ymin=max(ymin,-dxy/dx);ymax=min(ymax,(e[2]*dy-dxy)/dx);}
  1228.     }
  1229.   z[0]=0;z[2]=xmin;z[4]=xmax;e[4]=x2;e[5]=y2;
  1230.   if(dx*dy<0) {z[3]=ymax;z[5]=ymin;}
  1231.   else {z[3]=ymin;z[5]=ymax;}
  1232.   z[1]=((xmin>xmax)||(ymin>ymax))?0:2;
  1233.   if(!e[0]) e[0]=e[1]=(long)z;
  1234.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1235.   return gnil;
  1236. }
  1237.  
  1238.  
  1239. GEN rectbox(ne,x1,y1,x2,y2) /* code = 3 */
  1240.      long ne,x1,y1,x2,y2;
  1241. {
  1242.   long *e,*z;
  1243.   long xmin,ymin,xmax,ymax;
  1244.  
  1245.   if((ne<0)||(ne>15)) err(rploter2);
  1246.   e=rectgraph[ne];
  1247.   if(!(z=(long*)malloc(24))) err(memer);
  1248.   xmin=max(min(x1,x2),0);xmax=min(max(x1,x2),e[2]);
  1249.   ymin=max(min(y1,y2),0);ymax=min(max(y1,y2),e[3]);
  1250.   z[0]=0;z[1]=3;z[2]=xmin;z[3]=ymin;z[4]=xmax;z[5]=ymax;
  1251.   if(!e[0]) e[0]=e[1]=(long)z;
  1252.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1253.   return gnil;
  1254. }
  1255.  
  1256. GEN rectrbox(ne,x2,y2) /* code = 3 */
  1257.      long ne,x2,y2;
  1258. {
  1259.   long *e,*z;
  1260.   long x1,y1,xmin,ymin,xmax,ymax;
  1261.  
  1262.   if((ne<0)||(ne>15)) err(rploter2);
  1263.   e=rectgraph[ne];x1=e[4];y1=e[5];
  1264.   if(!(z=(long*)malloc(24))) err(memer);
  1265.   xmin=max(min(x1,x2),0);xmax=min(max(x1,x2),e[2]);
  1266.   ymin=max(min(y1,y2),0);ymax=min(max(y1,y2),e[3]);
  1267.   z[0]=0;z[1]=3;z[2]=xmin;z[3]=ymin;z[4]=xmax;z[5]=ymax;
  1268.   if(!e[0]) e[0]=e[1]=(long)z;
  1269.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1270.   return gnil;
  1271. }
  1272.  
  1273. GEN killrect(ne)
  1274.      long ne;
  1275. {
  1276.   long *e,*p1,*p2;
  1277.  
  1278.   if((ne<0)||(ne>15)) err(rploter2);
  1279.   e=rectgraph[ne];
  1280.   p1=(long*)e[0];e[0]=e[1]=e[2]=e[3]=e[4]=e[5]=0;
  1281.   while((long)p1) 
  1282.     {
  1283.       if((p1[1]==4)||(p1[1]==5)) {free((long *)p1[3]);free((long *)p1[4]);}
  1284.       if((p1[1]==6)) free((long *)p1[3]);
  1285.       p2=(long*)p1[0];free(p1);p1=p2;
  1286.     }
  1287.   return gnil;
  1288. }
  1289.  
  1290. GEN rectpoints(ne,listx,listy) /* code = 4 */
  1291.      long ne;
  1292.      GEN listx,listy;
  1293. {
  1294.   long *e,*z,lx,*px,*py,*ptx,*pty,x,y,i,cp;
  1295.  
  1296.   if((ne<0)||(ne>15)) err(rploter2);
  1297.   e=rectgraph[ne];
  1298.   if((typ(listx)<17)||(typ(listx)>18)||(typ(listy)<17)||(typ(listy)>18))
  1299.     err(ploter4);
  1300.   lx=lg(listx);if(lg(listy)!=lx) err(ploter5);
  1301.   lx--;if(!lx) return gnil;
  1302.   if(!(px=(long*)malloc(lx<<2))) err(memer);
  1303.   if(!(py=(long*)malloc(lx<<2))) err(memer);
  1304.   cp=0;
  1305.   for(i=0;i<lx;i++) 
  1306.     {
  1307.       x=px[i]=gtolong(listx[i+1]);y=py[i]=gtolong(listy[i+1]);
  1308.       if((x>=0)&&(y>=0)&&(x<=e[2])&&(y<=e[3])) cp++;
  1309.     }
  1310.   if(!cp) {free(px);free(py);return gnil;}
  1311.   if(!(ptx=(long*)malloc(cp<<2))) err(memer);
  1312.   if(!(pty=(long*)malloc(cp<<2))) err(memer);
  1313.   cp=0;
  1314.   for(i=0;i<lx;i++) 
  1315.     {
  1316.       x=px[i];y=py[i];
  1317.       if((x>=0)&&(y>=0)&&(x<=e[2])&&(y<=e[3])) {ptx[cp]=x;pty[cp]=y;}
  1318.       cp++;
  1319.     }
  1320.   free(px);free(py);
  1321.   if(!(z=(long*)malloc(20))) err(memer);
  1322.   z[0]=0;z[1]=4;z[2]=cp;z[3]=(long)ptx;z[4]=(long)pty;
  1323.   if(!e[0]) e[0]=e[1]=(long)z;
  1324.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1325.   return gnil;
  1326. }
  1327.  
  1328. GEN rectlines(ne,listx,listy) /* code = 5 */
  1329.      long ne;
  1330.      GEN listx,listy;
  1331. {
  1332.   long *e,*z,lx,*ptx,*pty,i;
  1333.  
  1334.   if((ne<0)||(ne>15)) err(rploter2);
  1335.   e=rectgraph[ne];
  1336.   if((typ(listx)<17)||(typ(listx)>18)||(typ(listy)<17)||(typ(listy)>18))
  1337.     err(ploter4);
  1338.   lx=lg(listx);if(lg(listy)!=lx) err(ploter5);
  1339.   lx--;if(!lx) return gnil;
  1340.   if(!(ptx=(long*)malloc(lx<<2))) err(memer);
  1341.   if(!(pty=(long*)malloc(lx<<2))) err(memer);
  1342.   for(i=0;i<lx;i++) {ptx[i]=gtolong(listx[i+1]);pty[i]=gtolong(listy[i+1]);}
  1343.   if(!(z=(long*)malloc(20))) err(memer);
  1344.   z[0]=0;z[1]=5;z[2]=lx;z[3]=(long)ptx;z[4]=(long)pty;
  1345.   if(!e[0]) e[0]=e[1]=(long)z;
  1346.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1347.   return gnil;
  1348. }
  1349.  
  1350. GEN rectstring(ne,x) /* code = 6 */
  1351.      long ne;
  1352.      GEN x;
  1353. {
  1354.   long *e,*z,i,a,lx;
  1355.   char *c;
  1356.   GEN p1;
  1357.  
  1358.   if((ne<0)||(ne>15)) err(rploter2);
  1359.   e=rectgraph[ne];
  1360.   if(!(z=(long*)malloc(24))) err(memer);
  1361.   if(typ(x)<17)
  1362.     {
  1363.       if(!(c=(char*)malloc(20))) err(memer);
  1364.       sprintf(c,"%9.3lf",gtodouble(x));
  1365.     }
  1366.   else
  1367.     {
  1368.       lx=lg(x);if(!(c=(char*)malloc(lx))) err(memer);
  1369.       for(i=1;i<lx;i++)
  1370.     {
  1371.       p1=(GEN)x[i];if(typ(p1)!=1) err(rploter6);
  1372.       a=itos(p1);if((a<0)||(a>255)) err(rploter6);
  1373.       c[i-1]=(char)a;
  1374.     }
  1375.       c[lx-1]=(char)0;
  1376.     }
  1377.   z[0]=0;z[1]=6;z[2]=strlen(c);z[3]=(long)c;z[4]=e[4];z[5]=e[5];
  1378.   if(!e[0]) e[0]=e[1]=(long)z;
  1379.   else {((long*)e[1])[0]=(long)z;e[1]=(long)z;}
  1380.   return gnil;
  1381. }
  1382.  
  1383. /*************************************************************************/
  1384. /*                                                                       */
  1385. /*                                                                       */
  1386. /*                         POSTSCRIPT OUTPUT                             */
  1387. /*                                                                       */
  1388. /*                                                                       */
  1389. /*************************************************************************/
  1390.  
  1391. typedef struct spoint {
  1392.   int x,y;} SPoint; 
  1393. typedef struct ssegment {
  1394.   int x1,y1,x2,y2;} SSegment;
  1395. typedef struct srectangle {
  1396.   int x,y,width,height;} SRectangle;
  1397.  
  1398. void ps_point(),ps_line(),ps_rect(),ps_string();
  1399.  
  1400. #undef ISCR
  1401. #undef JSCR
  1402. #define ISCR 1120 /* 1400 en haute resolution */     
  1403. #define JSCR 800  /* 1120 en haute resolution */     
  1404. #define DECI 100  /* 140 en haute resolution  */
  1405. #define DECJ  50  /* 70 en haute resolution   */
  1406. #define IDEC 90
  1407. #define JDEC 5
  1408.  
  1409. GEN  postploth(ep,a,b,ch)
  1410.      entree *ep;
  1411.      GEN a,b;
  1412.      char *ch;   
  1413.  
  1414. {
  1415.   long av,av2,jz,j,j1,i,sig,is,is2,js,js2;
  1416.   GEN p1,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
  1417.   char c1[20];
  1418.   char *c2;
  1419.   FILE *psfile;
  1420.  
  1421.   is=ISCR-DECI;js=JSCR-DECJ;is2=is-DECI;js2=js-DECJ;
  1422.   sig=gcmp(b,a); if(!sig) return gnil;
  1423.   psfile = fopen("pari.ps", "a");
  1424.   if (!psfile) err(talker,"cannot open psfile");
  1425.   fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont 16 scalefont setfont\n0.65 0.65 scale\n");
  1426.   ps_line(psfile,DECI,DECJ,DECI,js);
  1427.   ps_line(psfile,DECI,DECJ,is,DECJ);
  1428.   ps_line(psfile,is,DECJ,is,js);
  1429.   ps_line(psfile,DECI,js,is,js);
  1430.  
  1431.   av=avma;
  1432.   if(sig<0) {x=a;a=b;b=x;}
  1433.   for(i=1;i<=is2;i++) y[i]=cgetr(3);
  1434.   newvalue(ep,cgetr(3)); x=(GEN)ep->value; gaffect(a,x);
  1435.   dx=gdivgs(gsub(b,a),is2-1);ysml=gzero;ybig=gzero;
  1436.   av2=avma;
  1437.   for(i=1;i<=is2;i++)
  1438.   {
  1439.     gaffect(lisexpr(ch),y[i]);
  1440.     if(gcmp(y[i],ysml)<0) ysml=y[i];
  1441.     if(gcmp(y[i],ybig)>0) ybig=y[i];
  1442.     gaddz(x,dx,x);avma=av2;
  1443.   }
  1444.   diff=gsub(ybig,ysml);
  1445.   if(gcmp0(diff)) {ybig=gaddsg(1,ybig);diff=gun;}
  1446.   dyj=gdivsg(js2-1,diff);jz=js+itos(ground(gmul(ysml,dyj)));
  1447.   ps_line(psfile,DECI,jz,is,jz);
  1448.   if(gsigne(a)*gsigne(b)<0)
  1449.   {
  1450.     jz=1-itos(ground(gdiv(a,dx)))+DECI;
  1451.     ps_line(psfile,jz,DECJ,jz,js);
  1452.   }
  1453.   av2=avma;
  1454.   for(i=1;i<=is2;i++)
  1455.   {
  1456.     j1=js-itos(ground(gmul(gsub(y[i],ysml),dyj)));
  1457.     if(i==1) fprintf(psfile,"%d %d moveto\n",j1,DECI);
  1458.     else fprintf(psfile,"%d %d lineto\n",j1,i-1+DECI);
  1459.     avma=av2;
  1460.   }
  1461.   p1=cgetr(4);gaffect(ysml,p1);c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1462.   ps_string(psfile,5,js,c2);gaffect(ybig,p1);
  1463.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1464.   ps_string(psfile,5,DECJ,c2);gaffect(a,p1);
  1465.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1466.   ps_string(psfile,DECI-36,js+20,c2);gaffect(b,p1);
  1467.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1468.   ps_string(psfile,is-36,js+20,c2);
  1469.   fprintf(psfile,"stroke showpage\n");fclose(psfile);
  1470.   avma = av;killvalue(ep);return gnil;
  1471. }
  1472.  
  1473. GEN  postploth2(ep,a,b,ch)
  1474.      entree *ep;
  1475.      GEN a,b;
  1476.      char *ch;   
  1477.  
  1478. {
  1479.   long av,av2,jz,iz,k1,k,j,j1,i,sig,is,is2,js,js2;
  1480.   GEN p1,ysml,ybig,xsml,xbig,diffx,diffy,dxj,t,dyj,dt,y[ISCR+1],x[ISCR+1];
  1481.   char c1[20];
  1482.   char *c2;
  1483.   FILE *psfile;
  1484.  
  1485.   is=ISCR-DECI;js=JSCR-DECJ;is2=is-DECI;js2=js-DECJ;
  1486.   sig=gcmp(b,a); if(!sig) return gnil;
  1487.   psfile = fopen("pari.ps", "a");
  1488.   if (!psfile) err(talker,"cannot open psfile");
  1489.   fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont 16 scalefont setfont\n0.65 0.65 scale\n");
  1490.   ps_line(psfile,DECI,DECJ,DECI,js);
  1491.   ps_line(psfile,DECI,DECJ,is,DECJ);
  1492.   ps_line(psfile,is,DECJ,is,js);
  1493.   ps_line(psfile,DECI,js,is,js);
  1494.  
  1495.   av=avma;
  1496.   if(sig<0) {p1=a;a=b;b=p1;}
  1497.   for(i=1;i<=is2;i++) {x[i]=cgetr(3);y[i]=cgetr(3);}
  1498.   newvalue(ep,cgetr(3)); t=(GEN)ep->value; gaffect(a,t);
  1499.   dt=gdivgs(gsub(b,a),is2-1);ysml=ybig=xsml=xbig=gzero;
  1500.   av2=avma;
  1501.   for(i=1;i<=is2;i++)
  1502.   {
  1503.     p1=lisexpr(ch);gaffect(p1[1],x[i]);gaffect(p1[2],y[i]);
  1504.     if(gcmp(y[i],ysml)<0) ysml=y[i];
  1505.     if(gcmp(y[i],ybig)>0) ybig=y[i];
  1506.     if(gcmp(x[i],xsml)<0) xsml=x[i];
  1507.     if(gcmp(x[i],xbig)>0) xbig=x[i];
  1508.     gaddz(t,dt,t);avma=av2;
  1509.   }
  1510.   diffy=gsub(ybig,ysml);
  1511.   if(gcmp0(diffy)) {ybig=gaddsg(1,ybig);diffy=gun;}
  1512.   diffx=gsub(xbig,xsml);
  1513.   if(gcmp0(diffx)) {xbig=gaddsg(1,xbig);diffx=gun;}
  1514.   dyj=gdivsg(js2-1,diffy);jz=js+itos(ground(gmul(ysml,dyj)));
  1515.   dxj=gdivsg(is2-1,diffx);iz=DECI-itos(ground(gmul(xsml,dxj)));
  1516.   if(gsigne(ysml)*gsigne(ybig)<0)
  1517.     ps_line(psfile,DECI,jz,is,jz);
  1518.   if(gsigne(xsml)*gsigne(xbig)<0)
  1519.     ps_line(psfile,iz,DECJ,iz,js);
  1520.   av2=avma;
  1521.   for(i=1;i<=is2;i++)
  1522.   {
  1523.     j1=js-itos(ground(gmul(gsub(y[i],ysml),dyj)));
  1524.     k1=DECI+itos(ground(gmul(gsub(x[i],xsml),dxj)));
  1525.     if(i==1) fprintf(psfile,"%d %d moveto\n",j1,k1);
  1526.     else fprintf(psfile,"%d %d lineto\n",j1,k1);
  1527.     avma=av2;
  1528.   }
  1529.   p1=cgetr(4);gaffect(ysml,p1);c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1530.   ps_string(psfile,5,js,c2);gaffect(ybig,p1);
  1531.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1532.   ps_string(psfile,5,DECJ,c2);gaffect(xsml,p1);
  1533.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1534.   ps_string(psfile,DECI-36,js+20,c2);gaffect(xbig,p1);
  1535.   c2=(char *)sprintf(c1," %9.3lf ",rtodbl(p1));
  1536.   ps_string(psfile,is-36,js+20,c2);
  1537.   fprintf(psfile,"stroke showpage\n");fclose(psfile);
  1538.   avma = av;killvalue(ep);return gnil;
  1539. }
  1540.  
  1541. GEN  postplothraw(listx,listy)
  1542.      GEN listx,listy;
  1543. {
  1544.   long av = avma,av2,i,lx,is,js,is2,js2;
  1545.   SPoint *points;
  1546.   GEN p1,xsml,xbig,ysml,ybig,dx,dy,scal,scaly;
  1547.   FILE *psfile;
  1548.  
  1549.   if((typ(listx)<17)||(typ(listx)>18)||(typ(listy)<17)||(typ(listy)>18))
  1550.     err(ploter4);
  1551.   lx=lg(listx);
  1552.   if(lg(listy)!=lx) err(ploter5);
  1553.   if(lx==1) return gnil;
  1554.   points = (SPoint*)malloc(lx*sizeof(SPoint));
  1555.   if(!points) err(ploter6);
  1556.   is=ISCR-IDEC-45;js=JSCR-85;is2=is+IDEC;js2=js+JDEC;
  1557.   av=avma;xsml=xbig=(GEN)listx[1];ysml=ybig=(GEN)listy[1];
  1558.   for(i = 0; i < lx-1; i++)
  1559.     {
  1560.       p1=(GEN)listx[i+1];
  1561.       if(gcmp(p1,xsml)<0) xsml=p1;if(gcmp(p1,xbig)>0) xbig=p1;
  1562.       p1=(GEN)listy[i+1];
  1563.       if(gcmp(p1,ysml)<0) ysml=p1;if(gcmp(p1,ybig)>0) ybig=p1;
  1564.     }
  1565.   dx=gsub(xbig,xsml);dy=gsub(ybig,ysml);
  1566.   if(gcmp0(dx))
  1567.     {
  1568.       if(gcmp0(dy)) 
  1569.     {
  1570.       scal=gun;dx=gsubsg(is>>1,xsml);
  1571.       dy=gsubsg(js>>1,ysml);
  1572.     }
  1573.       else
  1574.     {
  1575.       scal=gdivsg(js,dy);
  1576.       dx=gneg(gmul(scal,xsml));dy=gneg(gmul(scal,ysml));
  1577.     }
  1578.     }
  1579.   else
  1580.     {
  1581.       scal=gdivsg(is,dx);
  1582.       if(!gcmp0(dy))
  1583.     {
  1584.       scaly=gdivsg(js,dy);if(gcmp(scaly,scal)<0) scal=scaly;
  1585.     }
  1586.       dx=gneg(gmul(scal,xsml));dy=gneg(gmul(scal,ysml));
  1587.     }
  1588.   for(i = 0; i < lx-1; i++)
  1589.     {
  1590.       av2=avma;
  1591.       points[i].x = IDEC + itos(ground(gadd(gmul(listx[i+1],scal),dx)));
  1592.       points[i].y = JDEC + itos(ground(gadd(gmul(listy[i+1],scal),dy)));
  1593.       avma=av2;
  1594.     }
  1595.   psfile = fopen("pari.ps", "a");
  1596.   if (!psfile) err(talker,"cannot open psfile");
  1597.   fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont 16 scalefont setfont\n0.65 0.65 scale\n");
  1598.   ps_line(psfile,IDEC,JDEC,IDEC,js2);
  1599.   ps_line(psfile,IDEC,JDEC,is2,JDEC);
  1600.   ps_line(psfile,is2,JDEC,is2,js2);
  1601.   ps_line(psfile,IDEC,js2,is2,js2);
  1602.   for(i = 0; i < lx-1; i++)
  1603.     ps_point(psfile,points[i].x,points[i].y,1);
  1604.   free(points);avma = av;
  1605.   fprintf(psfile,"stroke showpage\n");fclose(psfile);
  1606.   return gnil;
  1607. }
  1608.  
  1609. GEN postdraw(list)
  1610.      GEN list;
  1611. {
  1612.   long *e,*p1,*ptx,*pty,*numpoints,*numtexts,*xtexts,*ytexts;
  1613.   long n,i,j,x0,y0,av=avma;
  1614.   long a,b,c,d,nd[10],ne;
  1615.   char **texts;
  1616.   FILE *psfile;
  1617.  
  1618.   SPoint *points, **lines, *SLine;
  1619.   SSegment *segments; 
  1620.   SRectangle *rectangles, SRec;
  1621.  
  1622.   if(typ(list)!=17) err(rploter3);
  1623.   n=lg(list)-1;if(n%3) err(rploter4);
  1624.   n=n/3;if(!n) return gnil;
  1625.   nd[0]=nd[1]=nd[2]=nd[3]=nd[4]=nd[5]=nd[6]=0;
  1626.   for(i=0;i<n;i++)
  1627.     {
  1628.       if(typ(list[3*i+1])!=1) err(rploter5);
  1629.       ne=itos(list[3*i+1]);if((ne<0)||(ne>15)) err(rploter2);
  1630.       e=rectgraph[ne];
  1631.       p1=(long*)e[0];while((long)p1) 
  1632.     {
  1633.       if(p1[1]!=4) nd[p1[1]]++;
  1634.       else nd[1]+=p1[2];
  1635.       p1=(long*)p1[0];
  1636.     }
  1637.     }
  1638.   points=(SPoint*)malloc(nd[1]*sizeof(SPoint));
  1639.   segments=(SSegment*)malloc(nd[2]*sizeof(SSegment));
  1640.   rectangles=(SRectangle*)malloc(nd[3]*sizeof(SRectangle));
  1641.   lines=(SPoint**)malloc(nd[5]*sizeof(SPoint*));
  1642.   numpoints=(long*)malloc(nd[5]*sizeof(long));
  1643.   texts=(char**)malloc(nd[6]*sizeof(char*));
  1644.   numtexts=(long*)malloc(nd[6]*sizeof(long));
  1645.   xtexts=(long*)malloc(nd[6]*sizeof(long));
  1646.   ytexts=(long*)malloc(nd[6]*sizeof(long));
  1647.   nd[1]=nd[2]=nd[3]=nd[5]=nd[6]=0;
  1648.   for(i=0;i<n;i++)
  1649.     {
  1650.       e=rectgraph[itos(list[3*i+1])];x0=list[3*i+2];y0=list[3*i+3];
  1651.       if((typ(x0)!=1)||(typ(y0)!=1)) err(rploter5);
  1652.       x0=itos(x0);y0=itos(y0);
  1653.       p1=(long*)e[0];
  1654.       while((long)p1)
  1655.     {
  1656.       switch(p1[1])
  1657.         {
  1658.         case 1: 
  1659.           points[nd[1]].x=p1[2]+x0;
  1660.           points[nd[1]].y=p1[3]+y0;
  1661.           nd[1]++;break;
  1662.         case 2:
  1663.           segments[nd[2]].x1=p1[2]+x0;
  1664.           segments[nd[2]].y1=p1[3]+y0;
  1665.           segments[nd[2]].x2=p1[4]+x0;
  1666.           segments[nd[2]].y2=p1[5]+y0;
  1667.           nd[2]++;break;
  1668.         case 3:
  1669.           a=rectangles[nd[3]].x=p1[2]+x0;
  1670.           b=rectangles[nd[3]].y=p1[3]+y0;
  1671.           rectangles[nd[3]].width=p1[4]+x0-a;
  1672.           rectangles[nd[3]].height=p1[5]+y0-b;
  1673.           nd[3]++;break;
  1674.         case 4:
  1675.           ptx=(long*)p1[3];pty=(long*)p1[4];
  1676.           for(j=0;j<p1[2];j++)
  1677.         {
  1678.           points[nd[1]+j].x=ptx[j]+x0;
  1679.           points[nd[1]+j].y=pty[j]+y0;
  1680.         }
  1681.           nd[1]+=p1[2];break;
  1682.         case 5:
  1683.           ptx=(long*)p1[3];pty=(long*)p1[4];
  1684.           numpoints[nd[5]]=p1[2];
  1685.           lines[nd[5]]=(SPoint*)malloc(p1[2]*sizeof(SPoint));
  1686.           for(j=0;j<p1[2];j++)
  1687.         {
  1688.           lines[nd[5]][j].x=ptx[j]+x0;
  1689.           lines[nd[5]][j].y=pty[j]+y0;
  1690.         }
  1691.           nd[5]++;break;
  1692.         case 6: 
  1693.           texts[nd[6]]=(char*)p1[3];numtexts[nd[6]]=p1[2];
  1694.           xtexts[nd[6]]=p1[4]+x0;ytexts[nd[6]]=p1[5]+y0;
  1695.           nd[6]++;break;
  1696.         default: break;
  1697.         }
  1698.       p1=(long*)p1[0];
  1699.     }
  1700.     }
  1701.   psfile = fopen("pari.ps", "a");
  1702.   if (!psfile) err(talker,"cannot open psfile");
  1703.   fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont 16 scalefont setfont\n0.65 0.65 scale\n");
  1704.   for(i=0;i<nd[1];i++) ps_point(psfile,points[i].x,points[i].y);
  1705.   for(i=0;i<nd[2];i++) ps_line(psfile,segments[i].x1,segments[i].y1,segments[i].x2,segments[i].y2);
  1706.   for(i=0;i<nd[3];i++) 
  1707.     {
  1708.       SRec=rectangles[i];a=SRec.x;b=SRec.y;c=a+SRec.width;
  1709.       d=b+SRec.height;ps_rect(psfile,a,b,c,d);
  1710.     }
  1711.   for(i=0;i<nd[5];i++) 
  1712.     {
  1713.       SLine=lines[i];
  1714.       for(j=0;j<numpoints[i];j++)
  1715.     {
  1716.       if(!j) fprintf(psfile,"%d %d moveto\n",SLine[0].y,SLine[0].x);
  1717.       else fprintf(psfile,"%d %d lineto\n",SLine[j].y,SLine[j].x);
  1718.     }
  1719.     }
  1720.   for(i=0;i<nd[6];i++) 
  1721.       ps_string(psfile,xtexts[i],ytexts[i],texts[i]);
  1722.   fprintf(psfile,"stroke showpage\n");fclose(psfile);
  1723.   free(points);free(segments);free(rectangles);
  1724.   free(numpoints);for(i=0;i<nd[5];i++) free(lines[i]);
  1725.   free(lines);free(texts);free(numtexts);free(xtexts);free(ytexts);
  1726.   avma = av;return gnil;
  1727. }
  1728.  
  1729. void ps_point(psfile,x,y)
  1730.      int x,y;
  1731.      FILE *psfile;
  1732. {
  1733.   fprintf(psfile,"%d %d moveto\n0 2 rlineto 2 0 rlineto 0 -2 rlineto closepath fill\n",y,x);
  1734.   return;
  1735. }
  1736.  
  1737. void ps_line(psfile,x1,y1,x2,y2)
  1738.      int x1,y1,x2,y2;
  1739.      FILE *psfile;
  1740. {
  1741.   fprintf(psfile,"%d %d moveto\n%d %d lineto\n",y1,x1,y2,x2);
  1742.   return;
  1743. }
  1744.  
  1745. void ps_rect(psfile,x1,y1,x2,y2)
  1746.      int x1,y1,x2,y2;
  1747.      FILE *psfile;
  1748. {
  1749.   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);
  1750.   return;
  1751. }
  1752.  
  1753. void ps_string(psfile,x,y,c)
  1754.      int x,y;
  1755.      char *c;
  1756.      FILE *psfile;
  1757.   
  1758. {
  1759.   fprintf(psfile,"%d %d moveto 90 rotate\n(",y,x);
  1760.   fputs(c,psfile);fprintf(psfile,") show -90 rotate\n");
  1761.   return;
  1762. }
  1763.  
  1764.