home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / trans3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  31.8 KB  |  1,266 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                +++++++++++++++++++++++++++++++                 **/
  5. /**                +                             +                 **/
  6. /**                +  FONCTIONS TRANSCENDANTES   +                 **/
  7. /**                +     (troisieme partie)      +                 **/
  8. /**                +                             +                 **/
  9. /**                +     copyright Babe Cool     +                 **/
  10. /**                +                             +                 **/
  11. /**                +++++++++++++++++++++++++++++++                 **/
  12. /**                                                                **/
  13. /**                                                                **/
  14. /********************************************************************/
  15. /********************************************************************/
  16.  
  17. # include "genpari.h"
  18.  
  19. GEN qq(),inteta();
  20.  
  21. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  22. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  23. /*~                                                                     ~*/
  24. /*~                     FONCTION K DE BESSEL                            ~*/
  25. /*~                                                                     ~*/
  26. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  27. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  28.  
  29. GEN     kbessel(nu,gx,prec)
  30.     
  31.      GEN nu,gx;
  32.      long prec;
  33.     
  34. {
  35.   GEN  x,y,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp;
  36.   GEN   nu2,w;
  37.   long  l,lbin,av,av1,av2,k,k2,l1,n2,n,ex;
  38.  
  39.   if(typ(nu)==6) return kbessel2(nu,gx,prec);
  40.   if(typ(gx)!=2) {l=prec;k=1;}
  41.   else {l=lg(gx);k=0;x=gx;} 
  42.   y=cgetr(l);l1=l+1;
  43.   av=avma;if(k) gaffect(gx,x=cgetr(l));
  44.   u=cgetr(l1);v=cgetr(l1);c=cgetr(l1);d=cgetr(l1);
  45.   e=cgetr(l1);f=cgetr(l1);
  46.   nu2=gmulgs(gmul(nu,nu),-4);
  47.   n=16*(l-2)*LOG2+PI*sqrt(gtodouble(gnorm(nu)))/2;
  48.   n2=(n<<1);pitemp=mppi(l1);
  49.   lbin=74-(l<<5);av1=avma;
  50.   if (gcmpgs(x,n)<0)
  51.     {
  52.       zf=gsqrt(gdivgs(pitemp,n2));
  53.       zz=cgetr(l1);gdivgsz(un,(n2<<2),zz);
  54.       s=gun;t=gzero;k2=2*n2+1;
  55.       for (k=n2;k>0;--k)
  56.         {
  57.           k2-=2;
  58.           if(k2>32767) p1=gadd(mulss(k2,k2),nu2);
  59.           else p1=gaddsg(k2*k2,nu2);
  60.           ak=gdivgs(gmul(p1,zz),-k);s=gaddsg(1,gmul(ak,s));
  61.           t=gaddsg(k2,gmul(ak,t));
  62.         }
  63.       gmulz(s,zf,u);t=gmul2n(t,-1);
  64.       gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v);avma=av1;
  65.       affsr(n2,q=cgetr(l1));
  66.       r=gmul2n(x,1);av1=avma;
  67.       do
  68.         {
  69.           p1=divsr(5,q);
  70.           if (expo(p1)>= -1) p1=ghalf;
  71.           p2=subsr(1,divrr(r,q));
  72.           if (gcmp(p1,p2)>0) p1=p2;
  73.           gnegz(p1,c);avma=av1;
  74.           k=1;gaffsg(1,d);
  75.           affrr(u,e);affrr(v,f);av2=avma;
  76.           do
  77.             {
  78.               w=gadd(gmul(gsubsg(k,ghalf),u),gmul(gsubgs(q,k),v));
  79.               w=gadd(w,gmul(nu,gsub(u,gmul2n(v,1))));
  80.               gdivgsz(gmul(q,v),k,u);
  81.               gdivgsz(w,k,v);
  82.               gmulz(d,c,d);
  83.               gaddz(e,gmul(d,u),e);
  84.               gaddz(f,p1=gmul(d,v),f);
  85.               k++;ex=gexpo(p1)-gexpo(f);
  86.               avma=av2;
  87.             }
  88.           while(ex>lbin);
  89.           p1=u;u=e;e=p1;p1=v;v=f;f=p1;
  90.           gmulz(q,gaddsg(1,c),q);
  91.           p1=subrr(q,r);ex=expo(p1);avma=av1;
  92.         }
  93.       while(ex>lbin);
  94.       gmulz(u,gpui(gdivgs(x,n),nu),y);
  95.     }
  96.   else
  97.     {
  98.       p2=gmul2n(x,1);
  99.       zf=gsqrt(gdiv(pitemp,p2));
  100.       zz=gdiv(un,gmul2n(p2,2));
  101.       s=gun;k2=2*n2+1;
  102.       for (k=n2;k>0;--k)
  103.         {
  104.           k2-=2;if(k2>32767) p1=gadd(mulss(k2,k2),nu2);
  105.           else p1=gaddsg(k2*k2,nu2);
  106.           ak=gdivgs(gmul(p1,zz),k);
  107.           s=gsubsg(1,gmul(ak,s));
  108.         }
  109.       gmulz(s,zf,y);
  110.     }
  111.   gdivz(y,mpexp(x),y);avma=av;
  112.   return y;
  113. }
  114.  
  115.  
  116. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  117. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  118. /*~                                                                     ~*/
  119. /*~                     FONCTION U(a,b,z) GENERALE                      ~*/
  120. /*~                                                                     ~*/
  121. /*~                         ET CAS PARTICULIERS                         ~*/
  122. /*~                                                                     ~*/
  123. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  124. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  125.  
  126. GEN     hyperu(a,b,gx,prec)
  127.     
  128.      GEN  a,b,gx;
  129.      long prec;
  130.     
  131.      /* On suppose gx>0 et a,b reels   pour l'instant   */
  132. {
  133.   GEN  x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d;
  134.   GEN  w,a1,gn;
  135.   long  l,lbin,av,av1,av2,k,l1,n,ex;
  136.  
  137.   if(typ(gx)!=2) {l=prec;k=1;}
  138.   else {l=lg(gx);k=0;x=gx;} 
  139.   ex=(iscomplex(a)||iscomplex(b));
  140.   if(ex) {y=cgetg(3,6);y[1]=lgetr(l);y[2]=lgetr(l);}
  141.   else y=cgetr(l);
  142.   l1=l+1;av=avma;if(k) gaffect(gx,x=cgetr(l));
  143.   a1=gaddsg(1,gsub(a,b));
  144.   if(ex)
  145.   {
  146.     u=cgetg(3,6);u[1]=lgetr(l1);u[2]=lgetr(l1);
  147.     v=cgetg(3,6);v[1]=lgetr(l1);v[2]=lgetr(l1);
  148.     c=cgetg(3,6);c[1]=lgetr(l1);c[2]=lgetr(l1);
  149.     d=cgetg(3,6);d[1]=lgetr(l1);d[2]=lgetr(l1);
  150.     e=cgetg(3,6);e[1]=lgetr(l1);e[2]=lgetr(l1);
  151.     f=cgetg(3,6);f[1]=lgetr(l1);f[2]=lgetr(l1);
  152.   }
  153.   else
  154.   {
  155.     u=cgetr(l1);v=cgetr(l1);c=cgetr(l1);
  156.     d=cgetr(l1);e=cgetr(l1);f=cgetr(l1);
  157.   }
  158.   n=32*(l-2)*LOG2+PI*sqrt(gtodouble(gabs(gmul(a,a1),l1)));
  159.   lbin=74-(l<<5);av1=avma;
  160.   if (gcmpgs(x,n)<0)
  161.     {
  162.       zf=gpui(gn=stoi(n),gneg(a),l1);
  163.       zz=gdivsg(-1,gn);
  164.       s=gun;t=gzero;
  165.       for (k=n-1;k>=0;--k)
  166.         {
  167.           p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
  168.           s=gaddsg(1,gmul(p1,s));
  169.           t=gadd(gaddsg(k,a),gmul(p1,t));
  170.         }
  171.       gmulz(s,zf,u);t=gmul(t,zz);gmulz(t,zf,v);avma=av1;
  172.       affsr(n,q=cgetr(l1));
  173.       r=x;av1=avma;
  174.       do
  175.         {
  176.           p1=divsr(5,q);
  177.           if (expo(p1)>= -1) p1=ghalf;
  178.           p2=subsr(1,divrr(r,q));
  179.           if (gcmp(p1,p2)>0) p1=p2;
  180.           gnegz(p1,c);avma=av1;
  181.           k=0;gaffsg(1,d);
  182.           gaffect(u,e);gaffect(v,f);
  183.           p3=gsub(q,b);av2=avma;
  184.           do
  185.             {
  186.               w=gadd(gmul(gaddsg(k,a),u),gmul(gaddsg(-k,p3),v));
  187.               k++;gdivgsz(gmul(q,v),k,u);
  188.               gdivgsz(w,k,v);
  189.               gmulz(d,c,d);
  190.               gaddz(e,gmul(d,u),e);
  191.               gaddz(f,p1=gmul(d,v),f);
  192.               ex=gexpo(p1)-gexpo(f);
  193.               avma=av2;
  194.             }
  195.           while(ex>lbin);
  196.           p1=u;u=e;e=p1;p1=v;v=f;f=p1;
  197.           gmulz(q,gaddsg(1,c),q);
  198.           p1=subrr(q,r);ex=expo(p1);avma=av1;
  199.         }
  200.       while(ex>lbin);
  201.       gaffect(u,y);
  202.     }
  203.   else
  204.     {
  205.       zf=gpui(x,gneg(a));
  206.       zz=gdivsg(-1,x);
  207.       s=gun;
  208.       for (k=n-1;k>=0;--k)
  209.         {
  210.           p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
  211.           s=gaddsg(1,gmul(p1,s));
  212.         }
  213.       gmulz(s,zf,y);
  214.     }
  215.   avma=av;
  216.   return y;
  217. }
  218.  
  219. GEN   kbessel2(nu,x,prec)
  220.      GEN nu,x;
  221.      long prec;
  222. {
  223.   long av,tetpil,l;
  224.   GEN  p1,p2,x2,a,pitemp;
  225.  
  226.   av=avma;x2=gshift(x,1);
  227.   if(typ(x)==2) l=lg(x);else l=prec;
  228.   pitemp=mppi(l);
  229.   if(gcmp0(gimag(nu))) a=cgetr(l);
  230.   else
  231.     {a=cgetg(3,6);a[1]=lgetr(l);a[2]=lgetr(l);}
  232.   gaddsgz(1,gshift(nu,1),a);
  233.   p1=hyperu(gshift(a,-1),a,x2,prec);
  234.   p1=gmul(gmul(p1,gpui(x2,nu,prec)),mpsqrt(pitemp));p2=gexp(x,prec);
  235.   tetpil=avma;
  236.   return gerepile(av,tetpil,gdiv(p1,p2));
  237. }
  238.  
  239.  
  240. GEN eint1(x,prec)
  241.      GEN x;
  242.      long prec;
  243. {
  244.   long av,tetpil,l,n,i;
  245.   GEN p1,p2,p3,p4,y;
  246.  
  247.   av=avma;
  248.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  249.   if(expo(x)>=4)
  250.     {
  251.       tetpil=avma;y=gerepile(av,tetpil,incgam2(gzero,x,prec));
  252.     }
  253.   else
  254.     {
  255.       l=lg(x);
  256.       n= -32*(l-2)-1;
  257.       p1=cgetr(l);affsr(1,p1);p2=p1;p3=p1;p4=p1;i=1;
  258.       while(expo(p2)>=n)
  259.         {
  260.           i++;p1=gadd(p1,gdivgs(gun,i));p4=gdivgs(gmul(x,p4),i);
  261.           p2=gmul(p4,p1);p3=gadd(p2,p3);
  262.         }
  263.       p3=gmul(x,gmul(mpexp(negr(x)),p3));
  264.       consteuler(l);p1=gadd(mplog(x),geuler);tetpil=avma;
  265.       y=gerepile(av,tetpil,gsub(p3,p1));
  266.     }
  267.   return y;
  268. }
  269.  
  270. GEN   gerfc(x,prec)
  271.      GEN x;
  272.      long prec;
  273. {
  274.   long av,tetpil,l;
  275.   GEN p1,p2,x2,pitemp;
  276.  
  277.   av=avma;x2=gmul(x,x);p1=incgam(ghalf,x2,prec);
  278.   if(typ(x)==2) l=lg(x);else l=prec;
  279.   pitemp=mppi(l);
  280.   p2=mpsqrt(pitemp);tetpil=avma;
  281.   return gerepile(av,tetpil,gdiv(p1,p2));
  282. }
  283.  
  284. GEN   incgam(a,x,prec)
  285.      GEN a,x;
  286.      long prec;
  287. {
  288.   long av,tetpil;
  289.   GEN p1,p2,y;
  290.  
  291.   av=avma;
  292.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  293.   if((gcmp(subrs(x,1),a)>0)||(gsigne(greal(a))<=0))
  294.     {
  295.       tetpil=avma;y=gerepile(av,tetpil,incgam2(a,x,prec));
  296.     }
  297.   else
  298.     {
  299.       p2=ggamma(a,prec);p1=incgam3(a,x,prec);tetpil=avma;
  300.       y=gerepile(av,tetpil,gsub(p2,p1));
  301.     }
  302.   return y;
  303. }
  304.  
  305. GEN  incgam1(a,x,prec)
  306.      GEN a,x;
  307.      long prec;
  308. {
  309.   long av=avma,tetpil,l,n,i;
  310.   double m,mx;
  311.   GEN p1,p2,p3,y;
  312.  
  313.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  314.   l=lg(x);mx=rtodbl(x);
  315.   m=32*(l-2)*LOG2;n=m/(log(m)-(1+log(mx)));
  316.   gaffect(gaddsg(1,gsub(x,a)),p2=cgetr(l));p3=subrs(p2,n+1);
  317.   for(i=n;i>=1;i--) p3=addrr(subrs(p2,i),divrr(mulsr(i,x),p3));
  318.   y=gmul(mpexp(negr(x),prec),gpui(x,a,prec));tetpil=avma;
  319.   return gerepile(av,tetpil,divrr(y,p3));
  320. }
  321.  
  322. GEN  incgam2(a,x,prec)
  323.      GEN a,x;
  324.      long prec;
  325. {
  326.   long av=avma,tetpil,l,n,i;
  327.   double m,mx;
  328.   GEN p1,p2,p3,y;
  329.  
  330.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  331.   l=lg(x);mx=rtodbl(x);
  332.   m=(8*(l-2)*LOG2+mx/4);n=1+m*m/mx;
  333.   gaffect(gsub(x,a),p2=cgetr(l));
  334.   p3=gdiv(gsubgs(a,n),addrs(p2,(n<<1)));
  335.   for(i=n-1;i>=1;i--) p3=gdiv(gsubgs(a,i),addrr(addrs(p2,(i<<1)),gmulsg(i,p3)));
  336.   p3=gaddsg(1,p3);y=gmul(mpexp(negr(x),prec),gpui(x,gsubgs(a,1),prec));
  337.   tetpil=avma; return gerepile(av,tetpil,gmul(y,p3));
  338. }
  339.  
  340. GEN  incgam3(a,x,prec)
  341.      GEN a,x;
  342.      long prec;
  343.  
  344. {
  345.   long av=avma,tetpil,l,n,i;
  346.   GEN p1,p2,p3,y;
  347.  
  348.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  349.   l=lg(x);
  350.   n= -32*(l-2)-1;
  351.   p3=cgetr(l);affsr(1,p3);p2=gcopy(p3);i=0;
  352.   while(expo(p2)>=n)
  353.     {
  354.       i++;p2=gdiv(gmul(x,p2),gaddgs(a,i));
  355.       p3=gadd(p2,p3);
  356.     }
  357.   y=gdiv(gmul(mpexp(negr(x),prec),gpui(x,a,prec)),a);tetpil=avma;
  358.   return gerepile(av,tetpil,gmul(y,p3));
  359. }
  360.  
  361. GEN   incgam4(a,x,z,prec)
  362.      GEN a,x,z;
  363.      long prec;
  364.  
  365. /* One assumes that z=gamma(a,prec) but no check */
  366.  
  367. {
  368.   long av,tetpil;
  369.   GEN p1,p2,y;
  370.  
  371.   av=avma;
  372.   if(typ(x)!=2) {gaffect(x,p1=cgetr(prec));x=p1;}
  373.   if((gcmp(subrs(x,1),a)>0)||(gsigne(greal(a))<=0))
  374.     {
  375.       tetpil=avma;y=gerepile(av,tetpil,incgam2(a,x,prec));
  376.     }
  377.   else
  378.     {
  379.       p2=incgam3(a,x,prec);tetpil=avma;
  380.       y=gerepile(av,tetpil,gsub(z,p2));
  381.     }
  382.   return y;
  383. }
  384.  
  385. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  386. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  387. /*~                                    ~*/
  388. /*~                FONCTION ZETA DE RIEMANN            ~*/
  389. /*~                                    ~*/
  390. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  391. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  392.  
  393. GEN     czeta(s,prec)
  394.  
  395.      GEN s;
  396.      long prec;
  397.  
  398. {
  399.   long av,n,p,n1,l,l1,flag1,flag2,flag3,flag4,i,i2;
  400.   double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf,x00,x11,y00,epsil,fepsil;
  401.   GEN y,z,res,res2,sig,ms,t,z1,p1,p2,p3,p31,pitemp;
  402.  
  403.   if(gcmp1(s)) err(zetaer1);
  404.   if(typ(s)==6)
  405.     {
  406.       l=16383;
  407.       if(typ(s[1])==2) l=precision(s[1]);
  408.       if(typ(s[2])==2) {l1=precision(s[2]);if(l1<l) l=l1;}
  409.       if(l==16383) l=prec;
  410.       res=cgetg(3,6);res[1]=lgetr(l);res[2]=lgetr(l);
  411.       av=avma;
  412.       res2=cgetg(3,6);res2[1]=lgetr(l+1);res2[2]=lgetr(l+1);
  413.       gaffect(s,res2);s=res2;
  414.       sig=(GEN)s[1];
  415.     }
  416.   else
  417.     {
  418.       if(typ(s)!=2) err(zetaer2);
  419.       res=(signe(s)) ? cgetr(lg(s)) : cgetr(((-expo(s))>>5)+2);
  420.       av=avma;
  421.       res2=(signe(s)) ? cgetr(lg(s)+1) : cgetr(((-expo(s))>>5)+3);
  422.       affrr(s,res2);sig=s=res2;
  423.     }
  424.   if((signe(sig)<=0)||(expo(sig)<-1))
  425.     {
  426.       if(gcmp0(gimag(s))&&gcmp0(gfrac(gmul2n(sig,-1))))
  427.     {
  428.       if(gcmp0(sig)) {gaffect(gneg(ghalf),res);avma=av;return res;}
  429.       else {gaffsg(0,res);avma=av;return res;}
  430.     }
  431.       else {flag1=1;s=gsubsg(1,s);sig=greal(s);}
  432.     }
  433.   else flag1=0;
  434.   t=gimag(s);
  435.   ssig=rtodbl(sig);st=fabs(rtodbl(t));maxbeta=pow(3.0,-2.5);
  436.   if(st)
  437.     {
  438.       ns=ssig*ssig+st*st;
  439.       alpha=C2*(prec-2)-0.39-0.5*(ssig-1.0)*log(ns)-log(ssig)+ssig*2*C1;
  440.       beta=(alpha+ssig)/st-atan(ssig/st);
  441.       if(beta<=0)
  442.     {
  443.       if(ssig>=1.0)
  444.         {
  445.           p=0;sn=sqrt(ns);
  446.           n=ceil(exp(C2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
  447.         }
  448.       else
  449.         {
  450.           p=1;sn=ssig+1;n=ceil(sqrt(sn*sn+st*st)/(2*PI));
  451.         }
  452.     }
  453.       else
  454.     {
  455.       if(beta<maxbeta) xinf=beta+pow(3*beta,1.0/3.0);
  456.       else
  457.         {
  458.           epsil=0.0001;fepsil=0.0087;flag4=1;
  459.           x00=beta+PI/2.0;
  460.           while(flag4)
  461.         {
  462.           y00=x00*x00;x11=(beta+atan(x00))*(1+y00)/y00-1/x00;
  463.           if((x00-x11)<fepsil) flag4=0;
  464.           else x00=x11;
  465.         }
  466.           xinf=x11;
  467.         }
  468.       sp=1.0-ssig+st*xinf;
  469.       if(sp>0)
  470.         {
  471.           p=ceil(sp/2.0);sn=ssig+2*p-1;
  472.           n=ceil(sqrt(sn*sn+st*st)/(2*PI));
  473.         }
  474.       else
  475.         {
  476.           p=0;sn=sqrt(ns);
  477.           n=ceil(exp(C2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
  478.         }
  479.     }
  480.     }
  481.   else
  482.     {
  483.       beta=C2*(prec-2)+0.61+ssig*2*C1-ssig*log(ssig);
  484.       if(beta>0)
  485.     {
  486.       p=ceil(beta/2.0);sn=ssig+2*p-1;
  487.       n=ceil(sqrt(sn*sn+st*st)/(2*PI));
  488.     }
  489.       else
  490.     {
  491.       p=0;sn=sqrt(ssig*ssig+st*st);
  492.       n=ceil(exp(C2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
  493.     }
  494.     }
  495.   if(n<46340) {flag2=1;n1=n*n;} else flag2=0;
  496.   y=gun;ms=gneg(s);
  497.   for(i=2;i<=n;i++) y=gadd(y,p2=gexp(gmul(ms,glog(stoi(i),prec+1))));
  498.   flag3=((2*p)<46340);
  499.   mpbern(p,prec+1);p31=cgetr(prec+1);
  500.   z=gzero;
  501.   for(i=p;i>=1;i--)
  502.     {
  503.       i2=i<<1;
  504.       p1=gmul(gaddsg(i2-1,s),gaddsg(i2,s));
  505.       p1=flag3 ? gdivgs(p1,i2*(i2+1)) : gdivgs(gdivgs(p1,i2),i2+1);
  506.       p1=flag2 ? gdivgs(p1,n1) : gdivgs(gdivgs(p1,n),n);
  507.       if(bernzone[2]>prec+1) {affrr(bern(i),p31);p3=p31;} else p3=(GEN)bern(i);
  508.       z=gadd(divrs(p3,i),gmul(p1,z));
  509.     }
  510.   z1=gsub(gdivsg(n,gsubgs(s,1)),ghalf);
  511.   z=gmul(gadd(z1,gmul(s,gdivgs(z,n<<1))),p2);
  512.   if(!flag1) {gaffect(gadd(y,z),res);avma=av;}
  513.   else
  514.     {
  515.       pitemp=mppi(prec+1);setexpo(pitemp,2);
  516.       y=gmul(gmul(gadd(y,z),ggamma(s,prec+1)),gpui(pitemp,gneg(s)));
  517.       setexpo(pitemp,0);
  518.       y=gmul2n(gmul(y,gcos(gmul(pitemp,s),prec+1)),1);
  519.       gaffect(y,res);avma=av;
  520.     }
  521.   return res;
  522. }
  523.  
  524. GEN  izeta(x,prec)
  525.      GEN x;
  526.      long prec;
  527. {
  528.   long av=avma,tetpil,k,fl,n,s;
  529.   GEN y,p1,kk,pitemp,qn,z,q;
  530.  
  531.   if(typ(x)!=1) err(zetaer2);
  532.   if(gcmp1(x)) err(zetaer1);
  533.   s=signe(x);if(!s) return gneg(ghalf);
  534.   if(s<0)
  535.     {
  536.       k=itos(x);
  537.       if(k&1)
  538.     {
  539.       y=bernreal(1-k,prec);tetpil=avma;
  540.       return gerepile(av,tetpil,gdivgs(y,k-1));
  541.     }
  542.       else return gzero;
  543.     }
  544.   else
  545.     {
  546.       if(gcmpgs(x,((prec-2)<<5)+1)>0) {affsr(1,y=cgetr(prec));return y;}
  547.       else
  548.     {
  549.       k=itos(x);pitemp=mppi(prec);setexpo(pitemp,2);
  550.       if(k&1)
  551.         {
  552.           if((k&3)==3)
  553.         {
  554.           y=gzero;kk=stoi(k+1);
  555.           for(n=0;n<=((k+1)>>1);n+=2)
  556.             {
  557.               p1=gmul(binome(kk,n),gmul(bernreal(k+1-n,prec),bernreal(n,prec)));
  558.               if(n==((k+1)>>1)) p1=gmul2n(p1,-1);
  559.               y=((n>>1)&1)?gsub(y,p1):gadd(y,p1);
  560.             }
  561.           y=gmul(divrr(gpuigs(pitemp,k),mpfactr(k+1,prec)),y);
  562.           q=mpexp(pitemp,prec);z=divsr(1,addrs(q,-1));qn=q;fl=1;
  563.           for(n=2;fl;n++)
  564.             {
  565.               qn=mulrr(qn,q);p1=divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1)));
  566.               z=gadd(z,p1);fl=gexpo(p1)>= -1-((prec-2)<<5);
  567.             }
  568.           y=gadd(y,gmul2n(z,1));
  569.           tetpil=avma;return gerepile(av,tetpil,gneg(y));
  570.         }
  571.           else
  572.         {
  573.           y=gzero;kk=stoi(k+1);
  574.           for(n=0;n<=((k-1)>>1);n+=2)
  575.             {
  576.               p1=gmulsg(k+1-(n<<1),gmul(binome(kk,n),gmul(bernreal(k+1-n,prec),bernreal(n,prec))));
  577.               y=((n>>1)&1)?gsub(y,p1):gadd(y,p1);
  578.             }
  579.           y=divrs(gmul(divrr(gpuigs(pitemp,k),mpfactr(k+1,prec)),y),k-1);
  580.           q=mpexp(pitemp,prec);z=gzero;affsr(1,qn=cgetg(prec));fl=1;
  581.           for(n=1;fl;n++)
  582.             {
  583.               qn=mulrr(qn,q);p1=mulir(gpuigs(stoi(n),k),gsqr(addrs(qn,-1)));
  584.               p1=divrr(addrs(mulrr(qn,addsr(1,divrs(mulsr(n<<1,pitemp),k-1))),-1),p1);
  585.               z=gadd(z,p1);fl=gexpo(p1)>= -1-((prec-2)<<5);
  586.             }
  587.           z=gmul2n(z,1);
  588.           tetpil=avma;return gerepile(av,tetpil,gsub(y,z));
  589.         }
  590.         }
  591.       else
  592.         {
  593.           y=divrr(mulrr(gpuigs(pitemp,k),absr(bernreal(k,prec))),mpfactr(k,prec));
  594.           tetpil=avma;return gerepile(av,tetpil,gmul2n(y,-1));
  595.         }
  596.     }
  597.     }
  598. }
  599.       
  600.  
  601. GEN     gzeta(x,prec)
  602.     
  603.      GEN     x;
  604.      long    prec;
  605.     
  606. {
  607.   long    i,lx;
  608.   GEN     y;
  609.  
  610.     switch(typ(x))
  611.       {
  612.       case 1 : y=izeta(x,prec);break;
  613.       case 2 : 
  614.       case 6 : y=czeta(x,prec);break;
  615.       case 3 : 
  616.       case 7 : err(zetaer2);
  617.       case 11: err(impl,"zeta of power series");
  618.       case 17:
  619.       case 18:
  620.       case 19: lx=lg(x);y=cgetg(lx,typ(x));
  621.     for(i=1;i<lx;i++)
  622.       y[i]=(long)gzeta(x[i],prec);
  623.     break;
  624.       default: y=transc(gzeta,x,prec);
  625.       }
  626.   return y;
  627. }
  628.  
  629. void gzetaz(x,y)
  630.     
  631.      GEN     x,y;
  632.     
  633. {
  634.   long    av,prec;
  635.   GEN     p;
  636.  
  637.   prec=precision(y);
  638.   if(!prec) err(zetaer3);
  639.   av=avma;p=gzeta(x,prec);
  640.   gaffect(p,y);avma=av;
  641. }
  642.  
  643. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  644. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  645. /*~                                                                     ~*/
  646. /*~                     FONCTION POLYLOGARITHME                         ~*/
  647. /*~                                                                     ~*/
  648. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  649. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  650.  
  651. GEN cxpolylog(m,x,prec)
  652.      GEN x;
  653.      long m,prec; 
  654.  
  655. /* Ici, x peut etre complexe et le domaine de validite contient .005<|z|<230 */
  656.  
  657. {
  658.   long av=avma,tetpil,fl,i,n;
  659.   GEN p1,z,z2,h,q,s;
  660.  
  661.   if(gcmp1(x)) return izeta(stoi(m),prec);
  662.   z=glog(x,prec);h=gneg(glog(gneg(z),prec));for(i=1;i<=m-1;i++) h=gadd(h,gdivsg(1,stoi(i)));
  663.   q=gun;s=izeta(stoi(m),prec);
  664.   for(n=1;n<=m+1;n++) {q=gdivgs(gmul(q,z),n);s=gadd(s,gmul((n==(m-1))?h:izeta(stoi(m-n),prec),q));}
  665.   fl=1;z2=gmul(z,z);
  666.   for(n=m+3;fl;n+=2)
  667.     {
  668.       q=gdivgs(gmul(q,z2),(n-1)*n);p1=gmul(izeta(stoi(m-n),prec),q);
  669.       tetpil=avma;s=gadd(s,p1);fl=gexpo(p1)>= -((prec-2)<<5)-1;
  670.     }
  671.   return gerepile(av,tetpil,s);
  672. }
  673.  
  674. GEN polylog(m,x,prec)
  675.      GEN x;
  676.      long m,prec;
  677.  
  678. {
  679.   long av,tetpil,l,e,sx,i;
  680.   GEN p1,p2,unreel,n,y,logx,logx2;
  681.   
  682.   if(m<0) err(polyloger1);
  683.   if(!m) return gneg(ghalf);
  684.   if(gcmp0(x)) return gcopy(x);
  685.   av=avma;
  686.   if(m==1)
  687.     {
  688.       p1=glog(gsubsg(1,x),prec);tetpil=avma;return gerepile(av,tetpil,gneg(p1));
  689.     }
  690.   if(!(l=precision(x))) {l=prec;affsr(1,unreel=cgetr(l));x=gmul(unreel,x);}
  691.   e=gexpo(gnorm(x));if((!e)||(e== -1)) return cxpolylog(m,x,prec);
  692.   if(e>0) 
  693.     {
  694.       logx=glog(x,l);sx=gsigne(gimag(x));
  695.       if(!sx) {if(m%2) sx=gsigne(greal(gsubsg(1,x)));else sx= -gsigne(greal(x));}
  696.       x=ginv(x,l);
  697.     }
  698.   y=x;n=gun;p1=x;
  699.   do
  700.     {
  701.       n=addsi(1,n);p1=gmul(x,p1);p2=gdiv(p1,gpuigs(n,m));
  702.       tetpil=avma;y=gadd(y,p2);
  703.     }
  704.   while(gexpo(p2)>-32*(l-2));
  705.   if(e<=0) return gerepile(av,tetpil,y);
  706.   logx2=gmul(logx,logx);
  707.   p1=gneg(ghalf);
  708.   for(i=m-2;i>=0;i-=2)
  709.     {
  710.       p1=gadd(izeta(stoi(m-i),l),gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
  711.     }
  712.   if(m%2) p1=gmul(logx,p1);
  713.   p2=cgetg(3,6);p2[1]=zero;p2[2]=ldiv(gmulsg(sx,mppi(l)),mpfact(m-1));
  714.   p1=gadd(gmul2n(p1,1),gmul(p2,gpuigs(logx,m-1)));tetpil=avma;
  715.   y=(m%2)?gadd(p1,y):gsub(p1,y);
  716.   return gerepile(av,tetpil,y);
  717. }
  718.  
  719. GEN dilog(x,prec)
  720.      GEN x;
  721.      long prec;
  722. {
  723.   GEN p1,p2,p3,y;
  724.   long av,tetpil;
  725.  
  726.   av=avma;if(gcmpgs(gnorm(x,prec),1)<1)
  727.     {
  728.       tetpil=avma;return gerepile(av,tetpil,polylog(2,x,prec));
  729.     }
  730.   y=gneg(polylog(2,ginv(x),prec));p3=mppi(prec);p2=gdivgs(gmul(p3,p3),6);
  731.   p1=glog(gneg(x),prec);p1=gadd(p2,gmul2n(gmul(p1,p1),-1));
  732.   tetpil=avma;return gerepile(av,tetpil,gsub(y,p1));
  733. }
  734.  
  735. GEN polylogd(m,x,prec)
  736.      GEN x;
  737.      long m,prec;
  738.  
  739. {
  740.   long k,l,av,tetpil,fl,m2;
  741.   GEN p1,p2,p3,y,unreel;
  742.   
  743.   m2=m&1;av=avma;
  744.   if(gcmp0(x)) return gcopy(x);
  745.   if(gcmp1(x)&&(m>=2)) return m2?izeta(stoi(m),prec):gzero;
  746.   if(!(l=precision(x))) {l=prec;affsr(1,unreel=cgetr(l));x=gmul(unreel,x);}
  747.   p1=gabs(x,prec);fl=0;
  748.   if(gcmpgs(p1,1)>0) {x=ginv(x);p1=gabs(x,prec);fl=!m2;}
  749.   p1=gneg(glog(p1,prec));
  750.   y=m2?greal(polylog(m,x,prec)):gimag(polylog(m,x,prec));p2=gun;
  751.   for(k=1;k<=(m-1);k++)
  752.     {
  753.       p2=gdivgs(gmul(p2,p1),k);
  754.       p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
  755.       tetpil=avma;y=gadd(y,gmul(p2,p3));
  756.     }
  757.   if(m2)
  758.     {
  759.       p2=gdivgs(gmul(glog(gnorm(gsubsg(1,x)),prec),p1),2*m);tetpil=avma;y=gadd(y,p2);
  760.     }
  761.   if(fl) {tetpil=avma;return gerepile(av,tetpil,gneg(y));}
  762.   else {return gerepile(av,tetpil,y);}
  763. }
  764.  
  765. GEN polylogdold(m,x,prec)
  766.      GEN x;
  767.      long m,prec;
  768.  
  769. {
  770.   long k,l,av,tetpil,fl,m2;
  771.   GEN p1,p2,p3,y,unreel;
  772.   
  773.   m2=m&1;av=avma;
  774.   if(gcmp0(x)) return gcopy(x);
  775.   if(gcmp1(x)&&(m>=2)) return m2?izeta(stoi(m),prec):gzero;
  776.   if(!(l=precision(x))) {l=prec;affsr(1,unreel=cgetr(l));x=gmul(unreel,x);}
  777.   p1=gabs(x,prec);fl=0;
  778.   if(gcmpgs(p1,1)>0) {x=ginv(x);p1=gabs(x,prec);fl=!m2;}
  779.   p1=gneg(glog(p1,prec));
  780.   y=m2?greal(polylog(m,x,prec)):gimag(polylog(m,x,prec));p2=gun;
  781.   for(k=1;k<=(m-1);k++)
  782.     {
  783.       p2=gdivgs(gmul(p2,p1),k);
  784.       p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
  785.       tetpil=avma;y=gadd(y,gmul(p2,p3));
  786.     }
  787.   if(m2)
  788.     {
  789.       p2=gdivgs(gmul(p2,p1),-2*m);tetpil=avma;y=gadd(y,p2);
  790.     }
  791.   if(fl) {tetpil=avma;return gerepile(av,tetpil,gneg(y));}
  792.   else {return gerepile(av,tetpil,y);}
  793. }
  794.  
  795.  
  796. GEN polylogp(m,x,prec)
  797.      GEN x;
  798.      long m,prec;
  799.  
  800. {
  801.   long k,l,av,tetpil,fl,m2;
  802.   GEN p1,p2,p3,p4,p5,p51,y,unreel;
  803.   
  804.   m2=m&1;av=avma;
  805.   if(gcmp0(x)) return gcopy(x);
  806.   if(gcmp1(x)&&(m>=2)) return m2?izeta(stoi(m),prec):gzero;
  807.   if(!(l=precision(x))) {l=prec;affsr(1,unreel=cgetr(l));x=gmul(unreel,x);}
  808.   p1=gabs(x,prec);fl=0;
  809.   if(gcmpgs(p1,1)>0) {x=ginv(x);p1=gabs(x,prec);fl=!m2;}
  810.   p1=gmul2n(glog(p1,prec),1);mpbern(m>>1,prec);p51=cgetr(prec);
  811.   y=m2?greal(polylog(m,x,prec)):gimag(polylog(m,x,prec));p2=gun;
  812.   if(m==1)
  813.     {
  814.       p2=gmul2n(p1,-2);tetpil=avma;y=gadd(y,p2);
  815.     }
  816.   else
  817.     for(k=1;k<=(m-1);k++)
  818.       {
  819.     p2=gdivgs(gmul(p2,p1),k);
  820.     if((!(k&1))||(k==1))
  821.       {
  822.         if(k!=1) 
  823.           {
  824.         if(bernzone[2]>prec) {affrr(bern(k>>1),p51);p5=p51;}
  825.         else p5=(GEN)bern(k>>1);
  826.         p4=gmul(p2,p5);
  827.           }
  828.         else p4=gneg(gmul2n(p2,-1));
  829.         p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
  830.         tetpil=avma;y=gadd(y,gmul(p4,p3));
  831.       }
  832.       }
  833.   if(fl) {tetpil=avma;return gerepile(av,tetpil,gneg(y));}
  834.   else {return gerepile(av,tetpil,y);}
  835. }
  836.  
  837. GEN     gpolylog(m,x,prec)
  838.     
  839.      GEN     x;
  840.      long    m,prec;
  841.     
  842. {
  843.   long    i,lx,av=avma,tetpil,v,n;
  844.   GEN     y,p1,p2;
  845.  
  846.   if(m<=0)
  847.     {
  848.       p1=polx[0];p2=gsubsg(1,p1);
  849.       for(i=1;i<=(-m);i++) p1=gmul(polx[0],gadd(gmul(p2,deriv(p1,0)),gmulsg(i,p1)));
  850.       p1=gdiv(p1,gpuigs(p2,1-m));tetpil=avma;
  851.       y=gerepile(av,tetpil,gsubst(p1,0,x));
  852.     }
  853.   else
  854.     {
  855.       switch(typ(x))
  856.     {
  857.     case 1 : 
  858.     case 2 : 
  859.     case 4 :
  860.     case 5 :
  861.     case 6 :
  862.     case 8 : y=polylog(m,x,prec);break;
  863.     case 3 : 
  864.     case 7 : 
  865.     case 9 : p1=roots(x[1],prec);lx=lg(p1);p2=cgetg(lx,18);
  866.       for(i=1;i<lx;i++) p2[i]=lpoleval(x[2],p1[i]);
  867.       tetpil=avma;y=cgetg(lx,18);
  868.       for(i=1;i<lx;i++) y[i]=(long)polylog(m,p2[i],prec);
  869.       y=gerepile(av,tetpil,y);break;
  870.     case 10:
  871.     case 13:
  872.     case 14: p1=tayl(x,gvar(x),precdl);tetpil=avma;
  873.       y=gerepile(av,tetpil,gpolylog(m,p1,prec));
  874.       break;
  875.     case 11: if(m<0) err(polyloger1);
  876.       if(!m) return gneg(ghalf);
  877.       if(m==1)
  878.         {
  879.           p1=glog(gsubsg(1,x),prec);tetpil=avma;return gerepile(av,tetpil,gneg(p1));
  880.         }
  881.       if(valp(x)<=0) err(impl,"polylog around a!=0");
  882.       v=varn(x);n=(lg(x)-2)/valp(x);y=ggrando(polx[v],lg(x)-2);
  883.       for(i=n;i>=1;i--)
  884.         {
  885.           p1=gadd(gpuigs(stoi(i),-m),y);tetpil=avma;
  886.           y=gmul(x,p1);
  887.         }
  888.       y=gerepile(av,tetpil,y);break;
  889.     case 17:
  890.     case 18:
  891.     case 19: lx=lg(x);y=cgetg(lx,typ(x));
  892.       for(i=1;i<lx;i++)
  893.         y[i]=(long)gpolylog(m,x[i],prec);
  894.       break;
  895.     default: err(zetaer2);
  896.     }
  897.     }
  898.   return y;
  899. }
  900.  
  901. void gpolylogz(m,x,y)
  902.     
  903.      GEN     x,y;
  904.      long m;
  905. {
  906.   long    av,prec;
  907.   GEN     p;
  908.  
  909.   prec=precision(y);
  910.   if(!prec) err(zetaer3);
  911.   av=avma;p=gpolylog(m,x,prec);
  912.   gaffect(p,y);avma=av;
  913. }
  914.  
  915.  
  916. GEN qq(x,prec)
  917.     GEN x;
  918.     long prec;
  919. {
  920.   long av=avma,tetpil,l,tx=typ(x);
  921.   GEN p1,p2,q;
  922.   
  923.   if(tx==7) return gcopy(x);
  924.   if(tx<10)
  925.   {
  926.     if(!(l=precision(x))) l=prec;
  927.     p1=mppi(l);setexpo(p1,2);p2=cgetg(3,6);p2[1]=zero;p2[2]=(long)p1;
  928.     q=gmul(x,p2);tetpil=avma;
  929.     return gerepile(av,tetpil,gexp(q,prec));
  930.   }
  931.   else return tayl(x,gvar(x),precdl);
  932. }
  933.  
  934. GEN inteta(q)
  935.     GEN q;
  936. {
  937.   long av=avma,tetpil,tx=typ(q),l,n,f,v;
  938.   GEN p1,ps,qn,y0,y;
  939.  
  940.   y=gun;n=0;qn=gun;ps=gun;
  941.   if(tx==7)
  942.     {
  943.       do
  944.     {
  945.       n++;p1=gneg(gmul(ps,gmul(q,gmul(qn,qn))));y=gadd(y,p1);qn=gmul(qn,q);
  946.       ps=gmul(p1,qn);tetpil=avma;y0=y;y=gadd(y,ps);
  947.     }
  948.       while(!gegal(y0,y));
  949.     }
  950.   else
  951.     {
  952.       if(tx<10) l=precision(q);else v=gvar(q);
  953.       do
  954.     {
  955.       n++;p1=gneg(gmul(ps,gmul(q,gmul(qn,qn))));y=gadd(y,p1);qn=gmul(qn,q);
  956.       ps=gmul(p1,qn);tetpil=avma;y=gadd(y,ps);
  957.       f=(tx<10)?(gexpo(ps)-gexpo(y)>=-32*(l-2)):((!gcmp0(ps))&&(gval(ps,v)<precdl));
  958.     }
  959.       while(f);
  960.     }
  961.   return gerepile(av,tetpil,y);
  962. }
  963.  
  964.    
  965. GEN eta(x,prec)
  966.     GEN x;
  967.     long prec;
  968. {
  969.   long av=avma,tetpil;
  970.   GEN q;
  971.  
  972.   q=qq(x,prec);tetpil=avma;
  973.   return gerepile(av,tetpil,inteta(q));
  974. }
  975.  
  976. GEN jell(x,prec)
  977.      GEN x;
  978.      long prec;
  979. {
  980.   long av=avma,tetpil;
  981.   GEN p1,p2,q;
  982.   
  983.   q=qq(x,prec);p1=gdiv(inteta(gmul(q,q)),inteta(q));
  984.   p1=gmul2n(gmul(p1,p1),1);p1=gmul(q,gpuigs(p1,12));p2=gaddsg(768,gadd(gmul(p1,p1),gdivsg(4096,p1)));
  985.   p1=gmulsg(48,p1);tetpil=avma;
  986.   return gerepile(av,tetpil,gadd(p2,p1));
  987. }
  988.  
  989. GEN wf2(x,prec)
  990.      GEN x;
  991.      long prec;
  992. {
  993.   long av=avma,tetpil;
  994.   GEN p1,p2,q;
  995.   
  996.   q=qq(x,prec);p1=gmul(gdiv(inteta(gmul(q,q)),inteta(q)),gsqrt(deux,prec));
  997.   p2=cgetg(3,6);p2[1]=zero;p2[2]=ldivrs(mppi(prec),12);p2=gexp(gmul(x,p2),prec);tetpil=avma;
  998.   return gerepile(av,tetpil,gmul(p1,p2));
  999. }
  1000.  
  1001. GEN wf(x,prec)
  1002.      GEN x;
  1003.      long prec;
  1004. {
  1005.   long av=avma,tetpil;
  1006.   GEN p1,p2,q;
  1007.   
  1008.   q=qq(gmul2n(gaddgs(x,1),-1),prec);p1=gdiv(inteta(q),inteta(gmul(q,q)));
  1009.   p2=cgetg(3,6);p2[1]=zero;p2[2]=ldivrs(mppi(prec),-24);p2=gexp(gmul(p2,x),prec);tetpil=avma;
  1010.   return gerepile(av,tetpil,gmul(p1,p2));
  1011. }
  1012.  
  1013. GEN sagm(x,prec)
  1014.      GEN x;
  1015.      long prec;
  1016. {
  1017.   GEN z,p1,a,b,a1,b1;
  1018.   long av,tetpil,pp,ep;
  1019.  
  1020.   if(gcmp0(x)) return gcopy(x);
  1021.   switch(typ(x))
  1022.     {
  1023.     case 2:
  1024.     case 6: av=avma;if(pp=precision(x)) prec=pp;
  1025.       a1=x;b1=gun;
  1026.       do
  1027.     {
  1028.       a=a1;b=b1;
  1029.       a1=gmul2n(gadd(a,b),-1);
  1030.       b1=gsqrt(gmul(a,b),prec);
  1031.     }
  1032.       while(gexpo(gsub(b1,a1))>=gexpo(b1)-((prec-2)<<5)+5);
  1033.       tetpil=avma;z=gerepile(av,tetpil,gcopy(a1));break;
  1034.     case 7: av=avma;a1=x;b1=gun;pp=precp(x);
  1035.       do
  1036.     {
  1037.       a=a1;b=b1;
  1038.       a1=gmul2n(gadd(a,b),-1);b1=gsqrt(gmul(a,b));
  1039.       p1=gsub(b1,a1);ep=valp(p1)-valp(b1);
  1040.       if(ep<=0) {b1=gneg(b1);p1=gsub(b1,a1);ep=valp(p1)-valp(b1);}
  1041.     }
  1042.       while((ep<pp)&&(!gcmp0(p1)));
  1043.       tetpil=avma;z=gerepile(av,tetpil,gcopy(a1));break;
  1044.     case 11: av=avma;a1=x;b1=gun;pp=lg(x)-2;
  1045.       do
  1046.     {
  1047.       a=a1;b=b1;
  1048.       a1=gmul2n(gadd(a,b),-1);b1=gsqrt(gmul(a,b));
  1049.       p1=gsub(b1,a1);ep=valp(p1)-valp(b1);
  1050.       if(ep<=0) {b1=gneg(b1);p1=gsub(b1,a1);ep=valp(p1)-valp(b1);}
  1051.     }
  1052.       while(ep<pp);
  1053.       tetpil=avma;z=gerepile(av,tetpil,gcopy(a1));break;
  1054.     case 3: err(impl,"agm of mod");
  1055.     default: z=transc(sagm,x,prec);
  1056.     }
  1057.   return z;
  1058. }
  1059.  
  1060. GEN agm(x,y,prec)
  1061.      GEN x,y;
  1062.      long prec;
  1063. {
  1064.   GEN z;
  1065.   long av,tetpil;
  1066.  
  1067.   if(typ(y)>=17)
  1068.     {
  1069.       if(typ(x)>=17) err(agmer1);
  1070.       {z=x;x=y;y=z;}
  1071.     }
  1072.   if(gcmp0(y)) return gcopy(y);
  1073.   av=avma;z=sagm(gdiv(x,y),prec);tetpil=avma;
  1074.   return gerepile(av,tetpil,gmul(y,z));
  1075. }
  1076.   
  1077. GEN logagm(q)
  1078.      GEN q;
  1079. {
  1080.   long av=avma,prec=lg(q),tetpil,s,n,lim;
  1081.   GEN y,q4,q1,pitemp;
  1082.  
  1083.   if(typ(q)!=2) err(loger1);
  1084.   if(signe(q)<=0) err(loger2);
  1085.   lim= -((prec-2)<<4);n=0;
  1086.   if(expo(q)>=0) {q=ginv(q);s=0;} else s=1;
  1087.   while(expo(q)>=lim) {q1=q;q=mulrr(q,q);n=n+1;}
  1088.   q4=gmul2n(q,2);pitemp=mppi(prec);
  1089.   if(!n) y=divrr(pitemp,agm(addsr(1,q4),gmul2n(gsqrt(q,prec),2),prec));
  1090.   else y=divrr(pitemp,agm(addsr(1,q4),gmul2n(q1,2),prec));
  1091.   tetpil=avma;y=gmul2n(y,-n);if(s) setsigne(y,-1);
  1092.   return gerepile(av,tetpil,y);
  1093. }
  1094.  
  1095. GEN glogagm(x,prec)
  1096.   
  1097.    GEN   x;
  1098.    long  prec;
  1099.   
  1100. {
  1101.   long    av,tetpil,v;
  1102.   GEN y,p1,p2;
  1103.  
  1104.   switch(typ(x))
  1105.   {
  1106.   case 2 : if(signe(x)>=0) y=logagm(x);
  1107.   else
  1108.     {
  1109.     y=cgetg(3,6);y[2]=lmppi(lg(x));
  1110.     setsigne(x,1);y[1]=(long)logagm(x);
  1111.     setsigne(x,-1);
  1112.     }
  1113.     break;
  1114.   
  1115.   case 6 : y=cgetg(3,6);y[2]=larg(x,prec);
  1116.     av=avma;p1=glogagm(gnorm(x),prec);tetpil=avma;
  1117.     y[1]=lpile(av,tetpil,gmul2n(p1,-1));
  1118.     break;
  1119.   
  1120.   case 7 : y = palog(x); break;
  1121.   case 3 : err(loger3);
  1122.   
  1123.   case 11: av=avma;if(valp(x)) err(loger5);
  1124.     v=varn(x);p1=gdiv(deriv(x,v),x);
  1125.     if(gcmp1(x[2]))
  1126.     {
  1127.       tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));
  1128.     }
  1129.     else
  1130.     {
  1131.       p1=integ(p1,v);p2=glogagm(x[2],prec);
  1132.       tetpil=avma;y=gerepile(av,tetpil,gadd(p1,p2));
  1133.     }
  1134.     break;
  1135.   
  1136.   default: y=transc(glogagm,x,prec);
  1137.   }
  1138.   return y;
  1139. }
  1140.  
  1141. GEN theta(q,z,prec)
  1142.      GEN q,z;
  1143.      long prec;
  1144. {
  1145.   long av=avma,tetpil,l,n;
  1146.   GEN p1,ps,qn,qnold,y,zy,lq,ps2,unreel,k,zold;
  1147.  
  1148.   if(gexpo(q)>=0) err(thetaer1);
  1149.   if(l=precision(q)) prec=l;
  1150.   affsr(1,unreel=cgetr(prec));z=gmul(unreel,z);
  1151.   if(!l) q=gmul(unreel,q);
  1152.   if(gcmp0(zy=gimag(z))) k=gzero;
  1153.   else
  1154.     {
  1155.       lq=glog(q,prec);k=ground(gdiv(zy,greal(lq)));
  1156.       if(!gcmp0(k)) {zold=z;z=gadd(z,gdiv(gmul(lq,k),gi));}
  1157.     }
  1158.   y=gsin(z,prec);n=0;qn=gun;ps=gneg(ps2=gmul(q,q));
  1159.   do
  1160.     {
  1161.       n++;p1=gsin(gmulsg(n+n+1,z),prec);qnold=qn;qn=gmul(qn,ps);
  1162.       ps=gmul(ps,ps2);p1=gmul(p1,qn);
  1163.       y=gadd(y,p1);
  1164.     }
  1165.   while(gexpo(qnold)>= -((prec-2)<<5));
  1166.   if(!gcmp0(k))
  1167.     {
  1168.       y=gmul(y,gmul(gpui(q,gmul(k,k)),gexp(gmul2n(gmul(gmul(gi,zold),k),1),prec)));
  1169.       if(mpodd(k)) y=gneg(y);
  1170.     }
  1171.   p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1);tetpil=avma;
  1172.   return gerepile(av,tetpil,gmul(p1,y));
  1173. }
  1174.  
  1175. GEN thetanullk(q,k,prec)
  1176.      GEN q;
  1177.      long k,prec;
  1178. {
  1179.   long av=avma,tetpil,l,n;
  1180.   GEN p1,ps,qn,y,ps2,unreel;
  1181.  
  1182.   if(gexpo(q)>=0) err(thetaer1);
  1183.   if(!(k&1)) return gzero;
  1184.   n=0;qn=gun;ps=gneg(ps2=gmul(q,q));
  1185.   y=gun;if(!(l=precision(q))) 
  1186.     {
  1187.       l=prec;affsr(1,unreel=cgetr(prec));q=gmul(unreel,q);
  1188.     }
  1189.   do
  1190.     {
  1191.       n++;p1=gpuigs(stoi(n+n+1),k);qn=gmul(qn,ps);
  1192.       ps=gmul(ps,ps2);p1=gmul(p1,qn);
  1193.       y=gadd(y,p1);
  1194.     }
  1195.   while(gexpo(p1)>= -((l-2)<<5));
  1196.   p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1);tetpil=avma;
  1197.   if(k&2) {p1=gneg(p1);tetpil=avma;}
  1198.   return gerepile(av,tetpil,gmul(p1,y));
  1199. }
  1200.  
  1201. GEN jbesselh(n,z,prec)
  1202.      GEN n,z;
  1203.      long prec;
  1204. {
  1205.   long av,tetpil,k,l,i,lz;
  1206.   GEN y,p1,p2,zinv,p0,s,c;
  1207.  
  1208.   if(typ(n)!=1) err(jbesselher1);
  1209.   k=itos(n);if(k<0) err(impl,"ybesselh");
  1210.   
  1211.   switch(typ(z))
  1212.     {
  1213.     case 2 : 
  1214.     case 6 :
  1215.       if(gcmp0(z)) return gzero;
  1216.       av=avma;zinv=ginv(z);
  1217.       l=precision(z);if(l>prec) prec=l;
  1218.       gsincos(z,&s,&c,prec);
  1219.       p0=gmul(zinv,s);p1=gmul(zinv,gsub(p0,c));
  1220.       if(!k) p1=p0;
  1221.       else
  1222.     {
  1223.       for(i=2;i<=k;i++)
  1224.         {
  1225.           p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0);p0=p1;p1=p2;
  1226.         }
  1227.     }
  1228.       p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec);
  1229.       tetpil=avma;y=gerepile(av,tetpil,gmul(p2,p1));
  1230.       break;
  1231.  
  1232.     case 7 : err(impl,"p-adic jbessel function");
  1233.     case 3 : err(gamer3);
  1234.     case 11: err(impl,"jbessel of power series");
  1235.     case 17:
  1236.     case 18:
  1237.     case 19: lz=lg(z);y=cgetg(lz,typ(z));
  1238.       for(i=1;i<lz;i++)
  1239.         y[i]=(long)jbesselh(n,z[i],prec);
  1240.       break;
  1241.     case 1 :
  1242.     case 4 :
  1243.     case 5 : av=avma;p1=cgetr(prec);gaffect(z,p1);tetpil=avma;
  1244.       y=gerepile(av,tetpil,jbesselh(n,p1,prec));break;
  1245.     case 8 : av=avma;p1=cgetr(prec);affsr(1,p1);
  1246.       p1=gmul(z,p1);tetpil=avma;
  1247.       y=gerepile(av,tetpil,jbesselh(n,p1,prec));break;
  1248.     case 10:
  1249.     case 13:
  1250.     case 14: av=avma;p1=tayl(z,gvar(z),precdl);tetpil=avma;
  1251.       y=gerepile(av,tetpil,jbesselh(n,p1,prec));
  1252.       break;
  1253.     case 9 : av=avma;p1=roots(z[1],prec);lz=lg(p1);p2=cgetg(lz,18);
  1254.       for(i=1;i<lz;i++) p2[i]=lpoleval(z[2],p1[i]);
  1255.       tetpil=avma;y=cgetg(lz,18);
  1256.       for(i=1;i<lz;i++) y[i]=(long)jbesselh(n,p2[i],prec);
  1257.       y=gerepile(av,tetpil,y);break;
  1258.     case 15:
  1259.     case 16: err(transcer1);
  1260.     }
  1261.   return y;
  1262. }
  1263.  
  1264.  
  1265.  
  1266.