home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / bibli2 < prev    next >
Encoding:
Text File  |  1991-11-28  |  22.7 KB  |  895 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.   affsr(1,y=cgetr(prec));
  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)-32*(prec-2)-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.      
  661. #define JMIT 10000
  662.      
  663. {
  664.   long av,tetpil,j,nterm,jterm;
  665.   GEN  p1,sum,sum0,q1,tmp,dum,*tab,unreel;
  666.   
  667.   newvalue(ep,gun); q1=(GEN)ep->value; gaffect(a,q1);
  668.   tab=(GEN *)newbloc(JMIT);
  669.   av=avma;
  670.   p1=lisexpr(ch);affsr(1,unreel=cgetr(prec));
  671.   nterm=1;tab[1]=p1;sum=gmul2n(p1,-1);
  672.   for(jterm=2;jterm<=JMIT;jterm++)
  673.     {
  674.       tmp=tab[1];a=gaddsg(1,a);gaffect(a,q1);tab[1]=gmul(unreel,lisexpr(ch));
  675.       for(j=1;j<nterm;j++)
  676.     {
  677.       dum=tab[j+1];tab[j+1]=gmul2n(gadd(tab[j],tmp),-1);tmp=dum;
  678.     }
  679.       tab[nterm+1]=gmul2n(gadd(tab[nterm],tmp),-1);
  680.       if(gcmp(gabs(tab[nterm+1],prec),gabs(tab[nterm],prec))<=0)
  681.     p1=gmul2n(tab[++nterm],-1);
  682.       else p1=tab[nterm+1];
  683.       sum0=sum;sum=gadd(sum,p1);
  684.       if((gcmp0(p1)||(gexpo(p1)<-32*(prec-2)+5)||(gegal(sum,sum0)))&&(jterm>=10))
  685.     {
  686.       tetpil=avma;
  687.       killbloc((GEN)tab);
  688.       killvalue(ep);
  689.       return gerepile(av,tetpil,gcopy(sum));
  690.     }
  691.     }
  692.   err(eulsumer1);
  693. }
  694.  
  695. GEN  sumpos(ep,a,ch,prec)
  696.      entree *ep;
  697.      GEN a;
  698.      char *ch;
  699.      
  700. {
  701.   long av,tetpil,k,jterm,dsum;
  702.   GEN  sum,term,q1,p1,*tab,unreel,r;
  703.   
  704.   tab=(GEN *)newbloc(JMIT);
  705.   av = avma; a=gsubgs(a,1);
  706.   affsr(1,unreel=cgetr(prec));term=gzero;r=gun;k=0;
  707.   p1=cgeti(prec+1);p1[1]=0x1000001+prec;
  708.   newvalue(ep,p1); q1=(GEN)ep->value;
  709.   do
  710.     {
  711.       gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
  712.       term=gadd(term,p1);r=shifti(r,1);k++;
  713.     }
  714.   while(gexpo(p1)>=(-32*(prec-2)+5));
  715.   sum=gzero;
  716.   eulsum(&sum,term,1,tab,&dsum,prec);
  717.   for(jterm=2;jterm<=JMIT;jterm++)
  718.     {
  719.       term=gzero;r=stoi(jterm);k=0;
  720.       do
  721.     {
  722.       gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
  723.       term=gadd(term,p1);r=shifti(r,1);k++;
  724.     }
  725.       while(gexpo(p1)>=(-32*(prec-2)+5));
  726.       if(!odd(jterm)) term=gneg(term);
  727.       eulsum(&sum,term,jterm,tab,&dsum,prec);
  728.       if(dsum< -32*(prec-2)+5)
  729.     {
  730.       tetpil=avma;
  731.       killbloc((GEN)tab);
  732.       killvalue(ep);
  733.       return gerepile(av,tetpil,gcopy(sum));
  734.     }
  735.     }
  736.   killbloc((GEN)tab);
  737.   err(eulsumer1);
  738. }
  739.  
  740. /********************************************************************/
  741. /********************************************************************/
  742. /**                                                                **/
  743. /**                        TRACE GROSSIER                          **/
  744. /**                                                                **/
  745. /********************************************************************/
  746. /********************************************************************/
  747.  
  748.  
  749. GEN plot(ep,a,b,ch)
  750.      entree *ep;
  751.      GEN a,b;
  752.      char *ch;
  753.      
  754. #define ISCR 60
  755. #define JSCR 21
  756. #define BLANK ' '
  757. #define ZERO '-'
  758. #define YY '|'
  759. #define XX '-'
  760. #define FF 'x'
  761.      
  762. {
  763.   long av,av2,jz,j,i,sig;
  764.   GEN p1,p2,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
  765.   char scr[ISCR+1][JSCR+1], thestring[100];
  766.   
  767.   sig=gcmp(b,a);if(!sig) return gnil;
  768.   av=avma; pariputc('\n');
  769.   if(sig<0) {x=a;a=b;b=x;}
  770.   newvalue(ep,cgetr(3)); x=(GEN)ep->value;
  771.   for(i=1;i<=ISCR;i++) y[i]=cgetr(3);
  772.   gaffect(a,x);
  773.   dx=gdivgs(gsub(b,a),(ISCR-1));ysml=gzero;ybig=gzero;
  774.   for(j=1;j<=JSCR;j++) scr[1][j]=scr[ISCR][j]=YY;
  775.   for(i=2;i<ISCR;i++)
  776.     {
  777.       scr[i][1]=scr[i][JSCR]=XX;
  778.       for(j=2;j<JSCR;j++) scr[i][j]=BLANK;
  779.     }
  780.   av2=avma;
  781.   for(i=1;i<=ISCR;i++)
  782.     {
  783.       gaffect(lisexpr(ch),y[i]);
  784.       if(gcmp(y[i],ysml)<0) ysml=y[i];
  785.       if(gcmp(y[i],ybig)>0) ybig=y[i];
  786.       gaddz(x,dx,x);avma=av2;
  787.     }
  788.   diff=gsub(ybig,ysml);
  789.   if(gcmp0(diff)) {ybig=gaddsg(1,ybig);diff=gun;}
  790.   dyj=gdivsg(JSCR-1,diff);jz=1-itos(ground(gmul(ysml,dyj)));
  791.   av2=avma;
  792.   for(i=1;i<=ISCR;i++)
  793.     {
  794.       scr[i][jz]=ZERO;j=1+itos(ground(gmul(gsub(y[i],ysml),dyj)));
  795.       scr[i][j]=FF;avma=av2;
  796.     }
  797.   p1=cgetr(4);gaffect(ybig,p1);
  798.   sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
  799.   for(i=1;i<=ISCR;i++)  pariputc(scr[i][JSCR]); pariputc('\n');
  800.   for(j=(JSCR-1);j>=2;j--)
  801.     {
  802.       pariputs("            ");
  803.       for(i=1;i<=ISCR;i++) pariputc(scr[i][j]);
  804.       pariputc('\n');
  805.     }
  806.   p1=cgetr(4);gaffect(ysml,p1);
  807.   sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
  808.   for(i=1;i<=ISCR;i++)  pariputc(scr[i][1]);
  809.   pariputc('\n');
  810.   p1=cgetr(4);p2=cgetr(4);gaffect(a,p1);gaffect(b,p2);
  811.   sprintf(thestring, "%8s %10.3lf %44s %10.3lf\n"," ",rtodbl(p1)," ",rtodbl(p2));
  812.   pariputs(thestring);
  813.   killvalue(ep);
  814.   avma=av; pariputc('\n'); return gnil;
  815. }
  816.  
  817. /********************************************************************/
  818. /********************************************************************/
  819. /**                                                                **/
  820. /**                  RECHERCHE DE ZEROS REELS                      **/
  821. /**                                                                **/
  822. /********************************************************************/
  823. /********************************************************************/
  824.  
  825. GEN zbrent(ep,a,b,ch,prec)
  826.      entree *ep;
  827.      GEN a,b;
  828.      char *ch;
  829.      long prec;
  830.      
  831. {
  832.   long av=avma,tetpil,sig,iter,itmax;
  833.   GEN q1,c,d,e,tol,toli,min1,min2,fa,fb,fc,p,q,r,s,xm;
  834.   
  835.   if(typ(a)!=2) {p=cgetr(prec);gaffect(a,p);a=p;}
  836.   if(typ(b)!=2) {p=cgetr(prec);gaffect(b,p);b=p;}
  837.   sig=cmprr(b,a);if(!sig) {avma = av; return gzero;}
  838.   if(sig<0) {c=a;a=b;b=c;}
  839.   newvalue(ep,a); fa=lisexpr(ch);
  840.   changevalue(ep,b); q1=(GEN)ep->value; fb=lisexpr(ch);
  841.   if(gsigne(fa)*gsigne(fb)>0) err(zbrenter1);
  842.   itmax=32*prec+50;affsr(1,tol=cgetr(3));tol=shiftr(tol,-(32*(prec-2)-5));
  843.   fc=fb;
  844.   for(iter=1;iter<=itmax;iter++)
  845.     {
  846.       if(gsigne(fb)*gsigne(fc)>0)
  847.     {
  848.       c=a;fc=fa;d=subrr(b,a);e=d;
  849.     }
  850.       if(gcmp(gabs(fc),gabs(fb))<0)
  851.     {
  852.       a=b;b=c;c=a;fa=fb;fb=fc;fc=fa;
  853.     }
  854.       toli=mulrr(tol,absr(b));
  855.       xm=shiftr(subrr(c,b),-1);
  856.       if((cmprr(absr(xm),toli)<=0)||gcmp0(fb))
  857.     {
  858.       tetpil=avma;
  859.       killvalue(ep);
  860.       return gerepile(av,tetpil,gcopy(b));
  861.     }
  862.       if((cmprr(absr(e),toli)>=0)&&(gcmp(gabs(fa),gabs(fb))>0))
  863.     {
  864.       s=gdiv(fb,fa);
  865.       if(cmprr(a,b)==0)
  866.         {
  867.           p=gmul2n(gmul(xm,s),1);q=gsubsg(1,s);
  868.         }
  869.       else
  870.         {
  871.           q=gdiv(fa,fc);r=gdiv(fb,fc);
  872.           p=gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
  873.           p=gmul(s,gsub(p,gmul(gsub(b,a),gsubgs(r,1))));
  874.           q=gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
  875.         }
  876.       if(gsigne(p)>0) q=gneg(q);
  877.       p=gabs(p);
  878.       min1=gsub(gmulsg(3,gmul(xm,q)),gabs(gmul(q,toli)));
  879.       min2=gabs(gmul(e,q));
  880.       if(gcmp(gmul2n(p,1),gmin(min1,min2))<0) { e=d;d=gdiv(p,q);}
  881.       else {d=xm;e=d;}
  882.     }
  883.       else {d=xm;e=d;}
  884.       a=b;fa=fb;
  885.       if(gcmp(gabs(d),toli)>0) b=addrr(b,d);
  886.       else
  887.     {
  888.       if(gsigne(xm)>0) b=addrr(b,absr(toli));
  889.       else b=subrr(b,absr(toli));
  890.     }
  891.       gaffect(b,q1); ;fb=lisexpr(ch);
  892.     }
  893.   err(zbrenter2);
  894. }
  895.