home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / trans2.c < prev    next >
Text File  |  1992-05-20  |  45KB  |  1,839 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                +++++++++++++++++++++++++++++++                 **/
  5. /**                +                             +                 **/
  6. /**                +  FONCTIONS TRANSCENDANTES   +                 **/
  7. /**                +     (deuxieme partie)       +                 **/
  8. /**                +                             +                 **/
  9. /**                +     copyright Babe Cool     +                 **/
  10. /**                +                             +                 **/
  11. /**                +++++++++++++++++++++++++++++++                 **/
  12. /**                                                                **/
  13. /**                                                                **/
  14. /********************************************************************/
  15. /********************************************************************/
  16.  
  17. # include "genpari.h"
  18.  
  19. /********************************************************************/
  20. /********************************************************************/
  21. /**                                                                **/
  22. /**                       FONCTION ARCTG                           **/
  23. /**                                                                **/
  24. /********************************************************************/
  25. /********************************************************************/
  26.  
  27. GEN     mpatan(x)
  28.     
  29.      GEN     x;
  30.     
  31. {
  32.   long    l,l1,l2,n,m,u,i,avmacourant,av,lp;
  33.   long    e,sgn,s;
  34.   double  alpha,beta,gama=1.0,delta,fi;
  35.   GEN     y,p1,p2,p3,p4 ,p5;
  36.  
  37.   if (typ(x)!=2) err(ataner1);
  38.   sgn=signe(x);
  39.   if (!sgn)
  40.     {
  41.       y=cgetr(3);y[1]=x[1];
  42.       y[2]=0;
  43.     }
  44.   else
  45.     {
  46.       l=lp=lg(x);if(expo(x)>0) lp+=(expo(x)>>5);
  47.       y=cgetr(lp);avmacourant=avma;
  48.       p1=cgetr(l+1);affrr(x,p1);
  49.       if (sgn== -1) setsigne(p1,1);
  50.       u=cmprs(p1,1);
  51.       if (!u)
  52.         {
  53.           y=mppi(l+1);setexpo(y,-1);
  54.         }
  55.       else
  56.         {
  57.           if (u==1) divsrz(1,p1,p1);
  58.           if(expo(p1)<-100)
  59.             {
  60.               alpha=log(PI)-expo(p1)*LOG2;
  61.             }
  62.           else
  63.             {
  64.               alpha=rtodbl(p1);
  65.               alpha=log(PI/atan(alpha));
  66.             }
  67.           beta =16*LOG2*(l-2);
  68.           delta=LOG2+beta-alpha/2;
  69.           fi=alpha-2*LOG2;
  70.           if (delta<=0)
  71.             {
  72.               n=1;m=0;
  73.             }
  74.           else
  75.             {
  76.               if (delta>=gama*fi*fi/LOG2)
  77.                 {
  78.                   n=1+sqrt(gama*delta/LOG2);
  79.                   m=1+sqrt(delta/(gama*LOG2))-fi/LOG2;
  80.                 }
  81.               else
  82.                 {
  83.                   n=1+beta/fi;m=0;
  84.                 }
  85.             }
  86.           l2=l+1+m/32;
  87.           p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);
  88.           p5=cgetr(l2);
  89.           affrr(p1,p4);
  90.         
  91.           for (i=1;i<=m;i++)
  92.             {
  93.               mulrrz(p4,p4,p5);
  94.               addsrz(1,p5,p5);
  95.               av=avma;affrr(mpsqrt(p5),p5);avma=av;
  96.               addsrz(1,p5,p5);
  97.               divrrz(p4,p5 ,p4);
  98.             }
  99.           affrr(p4,p2);
  100.           mulrrz(p4,p4 ,p3);
  101.           l1=2;l2-=2;
  102.           setlg(p4,4);setlg(p5,4);
  103.           divssz(1,2*n+1,p4);
  104.           s=0;
  105.           setlg(p3,4);
  106.           e=expo(p3);
  107.         
  108.           for (i=n;i>=1;i--)
  109.             {
  110.               mulrrz(p4,p3,p4);
  111.               divssz(1,2*i-1,p5);
  112.               s-=e;l1+=(s/32);
  113.               if (l1>l2) l1=l2;
  114.               s %=32;
  115.               setlg(p3,l1+2);
  116.               setlg(p4,l1+2);
  117.               setlg(p5,l1+2);
  118.               subrrz(p5,p4,p4);
  119.             }
  120.         
  121.           setlg(p4,l2+2);
  122.           setlg(p5,l2+2);
  123.           mulrrz(p2,p4,p4);
  124.           setexpo(p4,expo(p4)+m);
  125.           if (u==1)
  126.             {
  127.               p5=mppi(lp+1);setexpo(p5,0);
  128.               subrrz(p5,p4,y);
  129.             }
  130.           else affrr(p4,y);
  131.           avma=avmacourant;
  132.         }
  133.       if (sgn== -1) setsigne(y,-signe(y));
  134.     }
  135.   return y;
  136. }
  137.  
  138. GEN     gatan(x,prec)
  139.     
  140.      GEN     x;
  141.      long    prec;
  142.     
  143. {
  144.   long    av,tetpil,l,v;
  145.   GEN     y,p1;
  146.  
  147.   switch(typ(x))
  148.     {
  149.     case 2 : y=mpatan(x);break;
  150.     case 6 : av=avma;p1=cgetg(3,6);
  151.       p1[1]=lneg(x[2]);
  152.       p1[2]=x[1];tetpil=avma;
  153.       y=gerepile(av,tetpil,gath(p1,prec));
  154.       l=y[1];y[1]=y[2];
  155.       y[2]=l;gnegz(l,l);
  156.       break;
  157.     
  158.     case 3 :
  159.     case 7 : err(ataner2);
  160.     
  161.     case 11: av=avma;if(valp(x)<0) err(ataner4);
  162.       v=varn(x);p1=gdiv(deriv(x,v),gaddsg(1,gmul(x,x)));
  163.       if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  164.       else
  165.         {
  166.           y=integ(p1,v);p1=gatan(x[2],prec);tetpil=avma;
  167.           y=gerepile(av,tetpil,gadd(p1,y));
  168.         }
  169.       break;
  170.     
  171.     default: y=transc(gatan,x,prec);
  172.     }
  173.   return y;
  174. }
  175.  
  176. void gatanz(x,y)
  177.     
  178.      GEN     x,y;
  179.     
  180. {
  181.   long    av,prec;
  182.   GEN     p;
  183.  
  184.   prec=precision(y);
  185.   if(!prec) err(ataner3);
  186.   av=avma;p=gatan(x,prec);
  187.   gaffect(p,y);avma=av;
  188. }
  189.  
  190. /********************************************************************/
  191. /********************************************************************/
  192. /**                                                                **/
  193. /**                      FONCTION ARCSINUS                         **/
  194. /**                                                                **/
  195. /********************************************************************/
  196. /********************************************************************/
  197.  
  198. GEN     mpasin(x)
  199.     
  200.      GEN     x;
  201.     
  202. {
  203.   long    l,u,v,sgn,av;
  204.   GEN     y,p1,p2;
  205.  
  206.   if (typ(x)!=2) err(asiner1);
  207.   else
  208.     {
  209.       sgn=signe(x);
  210.       if (!sgn)
  211.         {
  212.           y=cgetr(3);y[1]=x[1];
  213.           y[2]=0;
  214.         }
  215.       else
  216.         {
  217.           u=cmprs(x,1);v=cmpsr(-1,x);
  218.           if (u==1 || v==1) err(asiner2);
  219.           if (!u || !v)
  220.             {
  221.               y=mppi(lg(x));setexpo(y,0);
  222.               if (sgn== -1) setsigne(y,-1);
  223.             }
  224.           else
  225.             {
  226.               l=lg(x);y=cgetr(l+1);
  227.               av=avma;
  228.               p1=cgetr(l+1);
  229.               mulrrz(x,x,p1);
  230.               subsrz(1,p1,p1);
  231.               p2=mpsqrt(p1);
  232.               divrrz(x,p2,p1);
  233.               affrr(mpatan(p1),y);
  234.               if (sgn== -1) setsigne(y,-1);
  235.               avma=av;
  236.             }
  237.         }
  238.     }
  239.   return y;
  240. }
  241.  
  242. GEN     gasin(x,prec)
  243.     
  244.      GEN     x;
  245.      long    prec;
  246.     
  247. {
  248.   long    av,tetpil,l,v;
  249.   GEN     y,p1,z;
  250.  
  251.   switch(typ(x))
  252.     {
  253.     case 2 : av=avma;
  254.       if(gcmpgs(z=mpabs(x),1)<=0)
  255.         {avma=av;y=mpasin(x);}
  256.       else
  257.         {
  258.           tetpil=avma;y=cgetg(3,6);p1=mpach(z);
  259.           setsigne(p1,-signe(p1));y[2]=(long)p1;
  260.           y[1]=lmppi(lg(x));
  261.           setexpo(y[1],0);
  262.           if (signe(x)<0) gnegz(y,y);
  263.           y=gerepile(av,tetpil,y);
  264.         }
  265.       break;
  266.     
  267.     case 6 : av=avma;p1=cgetg(3,6);
  268.       p1[1]=lneg(x[2]);
  269.       p1[2]=x[1];tetpil=avma;
  270.       y=gerepile(av,tetpil,gash(p1,prec));
  271.       l=y[1];y[1]=y[2];
  272.       y[2]=l;gnegz(l,l);
  273.       break;
  274.     
  275.     case 3 :
  276.     case 7 : err(asiner3);
  277.     
  278.     case 11: if(gcmp0(x)) y=gcopy(x);
  279.     else
  280.       {
  281.         av=avma;if(valp(x)<0) err(asiner5);
  282.         v=varn(x);p1=gdiv(deriv(x,v),gsqrt(gsubsg(1,gmul(x,x)),prec));
  283.         if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  284.         else
  285.           {
  286.             y=integ(p1,v);p1=gasin(x[2],prec);tetpil=avma;
  287.             y=gerepile(av,tetpil,gadd(p1,y));
  288.           }
  289.       }
  290.       break;
  291.     
  292.     default: y=transc(gasin,x,prec);
  293.     }
  294.   return y;
  295. }
  296.  
  297. void gasinz(x,y)
  298.     
  299.      GEN     x,y;
  300.     
  301. {
  302.   long    av,prec;
  303.   GEN     p;
  304.  
  305.   prec=precision(y);
  306.   if(!prec) err(asiner4);
  307.   av=avma;p=gasin(x,prec);
  308.   gaffect(p,y);avma=av;
  309. }
  310.  
  311. /********************************************************************/
  312. /********************************************************************/
  313. /**                                                                **/
  314. /**                      FONCTION ARCCOSINUS                       **/
  315. /**                                                                **/
  316. /********************************************************************/
  317. /********************************************************************/
  318.  
  319. GEN     mpacos(x)
  320.     
  321.      GEN     x;
  322.     
  323. {
  324.   long    l,u,v,e,av;
  325.   GEN     y,p1,p2,pitemp;
  326.  
  327.   if (typ(x)!=2) err(acoser1);
  328.   else
  329.     {
  330.       u=cmprs(x,1);v=cmpsr(-1,x);
  331.       if (u==1 || v==1) err(acoser2);
  332.       if (!signe(x))
  333.         {
  334.           e=expo(x)>>5;if (e>=0) e= -1;
  335.           y=mppi(2-e);setexpo(y,0);
  336.         }
  337.       else
  338.         {
  339.           l=lg(x);
  340.           if (!u)
  341.             {
  342.               y=cgetr(3);y[2]=0;
  343.               y[1]=0x800000-((l-2)<<4);
  344.             }
  345.           else
  346.             {
  347.               if (!v) y=mppi(lg(x));
  348.               else
  349.                 {
  350.                   e=expo(x);if (e<0) e-=1;
  351.                   y=cgetr(l);
  352.                   av=avma;
  353.                   p1=cgetr(l+1);
  354.                   if (e<= -2)
  355.                     {
  356.                       mulrrz(x,x,p1);
  357.                       subsrz(1,p1,p1);
  358.                       affrr(mpsqrt(p1),p1);
  359.                       divrrz(x,p1,p1);
  360.                       affrr(mpatan(p1),y);
  361.                       pitemp=mppi(l);
  362.                       setexpo(pitemp,0);
  363.                       subrrz(pitemp,y,y);
  364.                       avma=av;
  365.                     }
  366.                   else
  367.                     {
  368.                       p2=cgetr(l+1);
  369.                       if (signe(x)>0)
  370.                         {
  371.                           addsrz(1,x,p1);
  372.                           subsrz(2,p1,p2);
  373.                           mulrrz(p1,p2,p1);
  374.                           affrr(mpsqrt(p1),p1);
  375.                           divrrz(p1,x,p1);
  376.                           affrr(mpatan(p1),y);
  377.                         }
  378.                       else
  379.                         {
  380.                           subsrz(1,x,p1);
  381.                           subsrz(2,p1,p2);
  382.                           mulrrz(p1,p2,p1);
  383.                           affrr(mpsqrt(p1),p1);
  384.                           divrrz(p1,x,p1);
  385.                           affrr(mpatan(p1),y);
  386.                           pitemp=mppi(l);
  387.                           addrrz(pitemp,y,y);
  388.                         }
  389.                     }
  390.                   avma=av;
  391.                 }
  392.             }
  393.         }
  394.     }
  395.   return y;
  396. }
  397.  
  398. GEN     gacos(x,prec)
  399.     
  400.      GEN     x;
  401.      long    prec;
  402.     
  403. {
  404.   long    av,tetpil,l,v;
  405.   GEN     y,p1;
  406.  
  407.   switch(typ(x))
  408.     {
  409.     case 2 : av=avma;
  410.       if(gcmpgs(p1=mpabs(x),1)<=0)
  411.         {avma=av;y=mpacos(x);}
  412.       else
  413.         {
  414.           tetpil=avma;y=cgetg(3,6);y[2]=lmpach(p1);
  415.           if(signe(x)>=0)
  416.             {
  417.               y[1]=zero;
  418.               setsigne(y[2],-signe(y[2]));
  419.             }
  420.           else y[1]=lmppi(lg(x));
  421.           y=gerepile(av,tetpil,y);
  422.         }
  423.       break;
  424.     
  425.     case 6 : y=gach(x,prec);
  426.       l=y[1];y[1]=y[2];
  427.       y[2]=l;gnegz(l,l);
  428.       break;
  429.     
  430.     case 3 :
  431.     case 7 : err(acoser3);
  432.     
  433.     case 11: av=avma;if(valp(x)<0) err(acoser5);
  434.       v=varn(x);
  435.       p1=integ(gdiv(deriv(x,v),gsqrt(gsubsg(1,gmul(x,x)),prec)),v);
  436.       if(gcmp1(x[2])&&(!valp(x)))
  437.         {tetpil=avma;y=gerepile(av,tetpil,gneg(p1));}
  438.       else
  439.         {
  440.           if(valp(x))
  441.             {y=mppi(prec);setexpo(y,0);}
  442.           else y=gacos(x[2],prec);
  443.           tetpil=avma;y=gerepile(av,tetpil,gsub(y,p1));
  444.         }
  445.       break;
  446.     
  447.     default: y=transc(gacos,x,prec);
  448.     }
  449.   return y;
  450. }
  451.  
  452. void gacosz(x,y)
  453.     
  454.      GEN     x,y;
  455.     
  456. {
  457.   long    av,prec;
  458.   GEN     p;
  459.  
  460.   prec=precision(y);
  461.   if(!prec) err(acoser4);
  462.   av=avma;p=gacos(x,prec);
  463.   gaffect(p,y);avma=av;
  464. }
  465.  
  466. /********************************************************************/
  467. /********************************************************************/
  468. /**                                                                **/
  469. /**                      FONCTION ARGUMENT                         **/
  470. /**                                                                **/
  471. /********************************************************************/
  472. /********************************************************************/
  473.  
  474. GEN     mparg(x,y)
  475.     
  476.      GEN     x,y;
  477.     
  478. {
  479.  
  480.   long    l,l1,sgnx,sgny,av,tetpil;
  481.   GEN     theta,pitemp;
  482.  
  483.   if (typ(x)!=2 || typ(y)!=2) err(arger1);
  484.   sgnx=signe(x);sgny=signe(y);
  485.   if (!sgny)
  486.     {
  487.       if (!sgnx) err(arger2);
  488.       l=lg(x);
  489.       if (sgnx>0)
  490.         {
  491.           theta=cgetr(3);theta[1]=y[1]-expo(x);
  492.           theta[2]=0;
  493.         }
  494.       else theta=mppi(l);
  495.     }
  496.   else
  497.     {
  498.       l=lg(y);l1=lg(x);if(l1>l) l=l1;
  499.       if(!sgnx)
  500.         {
  501.           theta=mppi(l);setexpo(theta,0);
  502.           if(sgny<0) setsigne(theta,-1);
  503.         }
  504.       else
  505.         {
  506.           if((expo(x)-expo(y))>-2)
  507.             {
  508.               av=avma;theta=divrr(y,x);tetpil=avma;
  509.               theta=mpatan(theta);
  510.               if(sgnx>0) theta=gerepile(av,tetpil,theta);
  511.               else
  512.                 {
  513.                   pitemp=mppi(l);tetpil=avma;
  514.                   if(sgny>0) theta=gerepile(av,tetpil,gadd(pitemp,theta));
  515.                   else theta=gerepile(av,tetpil,gsub(theta,pitemp));
  516.                 }
  517.             }
  518.           else
  519.             {
  520.               av=avma;
  521.               pitemp=mppi(l);
  522.               theta=mpatan(divrr(x,y));
  523.               tetpil=avma;setexpo(pitemp,0);
  524.               if(sgny>0) theta=gerepile(av,tetpil,gsub(pitemp,theta));
  525.               else
  526.                 {
  527.                   theta=gerepile(av,tetpil,gadd(pitemp,theta));
  528.                   setsigne(theta,-signe(theta));
  529.                 }
  530.             }
  531.         }
  532.     }
  533.   return theta;
  534. }
  535.  
  536. GEN     sarg(x,y,prec)
  537.      GEN        x,y;
  538.      long prec;
  539.  
  540. {
  541.  
  542.   long tetpil,av=avma;
  543.   GEN p1;
  544.  
  545.   switch(typ(x))
  546.     {
  547.     case 1 :
  548.     case 4 :
  549.     case 5 : p1=cgetr(prec);gaffect(x,p1);
  550.       x=p1;break;
  551.     default:;
  552.     }
  553.   switch(typ(y))
  554.     {
  555.     case 1 :
  556.     case 4 :
  557.     case 5 : p1=cgetr(prec);gaffect(y,p1);
  558.       y=p1;break;
  559.     default:;
  560.     }
  561.   if (av==(tetpil=avma)) return mparg(x,y);
  562.   else return gerepile(av,tetpil,mparg(x,y));
  563. }
  564.  
  565. GEN     garg(x,prec)
  566.     
  567.      GEN   x;
  568.      long  prec;
  569. {
  570.   GEN  p1,y;
  571.   long av,tx=typ(x),tetpil;
  572.  
  573.   if(gcmp0(x)) err(arger2);
  574.   switch(tx)
  575.     {
  576.     case 2 : prec=lg(x);
  577.     case 1 :
  578.     case 4 :
  579.     case 5 :
  580.       if(signe(x)>0)
  581.         {
  582.           y=cgetr(3);y[1]=0x800000-((prec-2)<<5);
  583.           y[2]=zero;
  584.         }
  585.       else y=mppi(prec);
  586.       break;
  587.     case 8 : av=avma;gaffsg(1,p1=cgetr(prec));p1=gmul(p1,x);
  588.       tetpil=avma;y=gerepile(av,tetpil,garg(p1,prec));
  589.       break;
  590.     case 6 : y=sarg(x[1],x[2],prec);
  591.       break;
  592.     default: if(tx>=17) y=transc(garg,x,prec);else err(arger1);
  593.     }
  594.   return y;
  595. }
  596.  
  597. /********************************************************************/
  598. /********************************************************************/
  599. /**                                                                **/
  600. /**                 FONCTION COSINUS HYPERBOLIQUE                  **/
  601. /**                                                                **/
  602. /********************************************************************/
  603. /********************************************************************/
  604.  
  605. GEN     mpch(x)
  606.     
  607.      GEN     x;
  608.     
  609. {
  610.   long    l,av;
  611.   GEN     y,p1,p2;
  612.  
  613.   if (typ(x)!=2) err(cher3);
  614.   if(gcmp0(x)) y=gaddsg(1,x);
  615.   else
  616.     {
  617.       l=lg(x);y=cgetr(l);
  618.       av=avma;
  619.       p1=cgetr(l+1);p2=cgetr(l+1);
  620.       affrr(mpexp1(x),p1);
  621.       addsrz(1,p1,p2);
  622.       divsrz(1,p2,p2);
  623.       addrrz(p1,p2,p2);
  624.       addsrz(1,p2,y);
  625.       setexpo(y,expo(y)-1);
  626.       avma=av;
  627.     }
  628.   return y;
  629. }
  630.  
  631. GEN     gch(x,prec)
  632.     
  633.      GEN     x;
  634.      long    prec;
  635.     
  636. {
  637.   long    av,tetpil;
  638.   GEN     y,p1;
  639.  
  640.   switch(typ(x))
  641.     {
  642.     case 2 : y=mpch(x);break;
  643.     
  644.     case 6 : av=avma;p1=gexp(x,prec);
  645.       p1=gadd(p1,ginv(p1));tetpil=avma;
  646.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  647.       break;
  648.     
  649.     case 3 :
  650.     case 7 : err(cher1);
  651.     
  652.     case 11: av=avma;p1=gexp(x,prec);p1=gadd(p1,gdivsg(1,p1));
  653.       tetpil=avma;y=gerepile(av,tetpil,gmul2n(p1,-1));
  654.       break;
  655.     
  656.     default: y=transc(gch,x,prec);
  657.     }
  658.   return y;
  659. }
  660.  
  661. void gchz(x,y)
  662.     
  663.      GEN     x,y;
  664.     
  665. {
  666.   long    av,prec;
  667.   GEN     p;
  668.  
  669.   prec=precision(y);
  670.   if(!prec) err(cher2);
  671.   av=avma;p=gch(x,prec);
  672.   gaffect(p,y);avma=av;
  673. }
  674.  
  675. /********************************************************************/
  676. /********************************************************************/
  677. /**                                                                **/
  678. /**                 FONCTION SINUS HYPERBOLIQUE                    **/
  679. /**                                                                **/
  680. /********************************************************************/
  681. /********************************************************************/
  682.  
  683. GEN     mpsh(x)
  684.     
  685.      GEN     x;
  686.     
  687. {
  688.   long    l,av;
  689.   GEN     y,p1,p2;
  690.  
  691.   if (typ(x)!=2) err(sher3);
  692.   if(!signe(x))
  693.     {
  694.       y=cgetr(3);y[1]=x[1];y[2]=0;
  695.     }
  696.   else
  697.     {
  698.       l=lg(x);y=cgetr(l);
  699.       av=avma;
  700.       p1=cgetr(l+1);p2=cgetr(l+1);
  701.       affrr(mpexp1(x),p1);
  702.       addsrz(1,p1,p2);
  703.       divrrz(p1,p2,p2);
  704.       addrrz(p1,p2,y);
  705.       setexpo(y,expo(y)-1);
  706.       avma=av;
  707.     }
  708.   return y;
  709. }
  710.  
  711. GEN     gsh(x,prec)
  712.     
  713.      GEN     x;
  714.      long    prec;
  715.     
  716. {
  717.   long    av,tetpil;
  718.   GEN     y,p1;
  719.  
  720.   switch(typ(x))
  721.     {
  722.     case 2 : y=mpsh(x);break;
  723.     
  724.     case 6 : av=avma;p1=gexp(x,prec);
  725.       p1=gsub(p1,ginv(p1));tetpil=avma;
  726.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  727.       break;
  728.     
  729.     case 3 :
  730.     case 7 : err(sher1);
  731.     
  732.     case 11: if(gcmp0(x)) y=gcopy(x);
  733.     else
  734.       {
  735.         av=avma;p1=gexp(x,prec);p1=gsub(p1,gdivsg(1,p1));
  736.         tetpil=avma;y=gerepile(av,tetpil,gmul2n(p1,-1));
  737.       }
  738.       break;
  739.     
  740.     default: y=transc(gsh,x,prec);
  741.     }
  742.   return y;
  743. }
  744.  
  745. void gshz(x,y)
  746.     
  747.      GEN     x,y;
  748.     
  749. {
  750.   long    av,prec;
  751.   GEN     p;
  752.  
  753.   prec=precision(y);
  754.   if(!prec) err(sher2);
  755.   av=avma;p=gsh(x,prec);
  756.   gaffect(p,y);avma=av;
  757. }
  758.  
  759. /********************************************************************/
  760. /********************************************************************/
  761. /**                                                                **/
  762. /**                 FONCTION TANGENTE HYPERBOLIQUE                 **/
  763. /**                                                                **/
  764. /********************************************************************/
  765. /********************************************************************/
  766.  
  767. GEN     mpth(x)
  768.     
  769.      GEN     x;
  770.     
  771. {
  772.   long    l,av;
  773.   GEN     y,p1,p2;
  774.  
  775.   if (typ(x)!=2) err(ther1);
  776.   if (!signe(x))
  777.     {
  778.       y=cgetr(3);y[1]=x[1];y[2]=0;
  779.     }
  780.   else
  781.     {
  782.       l=lg(x);y=cgetr(l);
  783.       av=avma;
  784.       p1=cgetr(l+1);p2=cgetr(l+1);
  785.       affrr(x,p1);
  786.       setexpo(p1,expo(p1)+1);
  787.       affrr(mpexp1(p1),p1);
  788.       addsrz(2,p1,p2);
  789.       divrrz(p1,p2,y);
  790.       avma=av;
  791.     }
  792.   return y;
  793. }
  794.  
  795. GEN     gth(x,prec)
  796.     
  797.      GEN     x;
  798.      long    prec;
  799.     
  800. {
  801.   long    av,tetpil;
  802.   GEN     y,p1,p2,p3;
  803.  
  804.   switch(typ(x))
  805.     {
  806.     case 2 : y=mpth(x);break;
  807.     
  808.     case 6 : av=avma;p1=gexp(gmul2n(x,1),prec);
  809.       p1=gdivsg(-2,gaddgs(p1,1));tetpil=avma;
  810.       y=gerepile(av,tetpil,gaddsg(1,p1));
  811.       break;
  812.     
  813.     case 3 :
  814.     case 7 : err(ther2);
  815.     
  816.     case 11: if(gcmp0(x)) y=gcopy(x);
  817.     else
  818.       {
  819.         av=avma;p1=gexp(gmul2n(x ,1),prec);
  820.         p2=gsubgs(p1,1);p3=gaddgs(p1,1);
  821.         tetpil=avma;y=gerepile(av,tetpil,gdiv(p2,p3));
  822.       }
  823.       break;
  824.     
  825.     default: y=transc(gth,x,prec);
  826.     }
  827.   return y;
  828. }
  829.  
  830. void gthz(x,y)
  831.     
  832.      GEN     x,y;
  833.     
  834. {
  835.   long    av,prec;
  836.   GEN     p;
  837.  
  838.   prec=precision(y);
  839.   if(!prec) err(ther3);
  840.   av=avma;p=gth(x,prec);
  841.   gaffect(p,y);avma=av;
  842. }
  843.  
  844. /********************************************************************/
  845. /********************************************************************/
  846. /**                                                                **/
  847. /**             FONCTION ARGUMENT SINUS HYPERBOLIQUE               **/
  848. /**                                                                **/
  849. /********************************************************************/
  850. /********************************************************************/
  851.  
  852. GEN     mpash(x)
  853.     
  854.      GEN     x;
  855.     
  856. {
  857.   long    l,av;
  858.   GEN     y,p1;
  859.  
  860.   if (typ(x)!=2) err(asher1);
  861.   if (!signe(x))
  862.     {
  863.       y=cgetr(3);y[1]=x[1];y[2]=0;
  864.     }
  865.   else
  866.     {
  867.       l=lg(x);y=cgetr(l);
  868.       av=avma;
  869.       p1=cgetr(l+1);
  870.       affrr(x,p1);
  871.       mulrrz(p1,p1,p1);
  872.       addsrz(1,p1,p1);
  873.       affrr(mpsqrt(p1),p1);
  874.       addrrz(x,p1,p1);
  875.       affrr(mplog(p1),y);
  876.       avma=av;
  877.     }
  878.   return y;
  879. }
  880.  
  881. GEN     gash(x,prec)
  882.     
  883.      GEN     x;
  884.      long    prec;
  885.     
  886. {
  887.   long    av,tetpil,v;
  888.   GEN     y,p1;
  889.  
  890.   switch(typ(x))
  891.     {
  892.     case 2 : y=mpash(x);break;
  893.     
  894.     case 6 : av=avma;p1=gaddsg(1,gmul(x,x));
  895.       p1=gadd(x,gsqrt(p1,prec));
  896.       p1=glog(gmul(p1,p1),prec);tetpil=avma;
  897.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  898.       break;
  899.     
  900.     case 3 :
  901.     case 7 : err(asher2);
  902.     
  903.     case 11: if(gcmp0(x)) y=gcopy(x);
  904.     else
  905.       {
  906.         av=avma;if(valp(x)<0) err(asher4);
  907.         v=varn(x);p1=gdiv(deriv(x,v),gsqrt(gaddsg(1,gmul(x,x))));
  908.         if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  909.         else
  910.           {
  911.             y=integ(p1,v);p1=gash(x[2],prec);tetpil=avma;
  912.             y=gerepile(av,tetpil,gadd(p1,y));
  913.           }
  914.       }
  915.       break;
  916.     
  917.     default: y=transc(gash,x,prec);
  918.     }
  919.   return y;
  920. }
  921.  
  922. void gashz(x,y)
  923.     
  924.      GEN     x,y;
  925.     
  926. {
  927.   long    av,prec;
  928.   GEN     p;
  929.  
  930.   prec=precision(y);
  931.   if(!prec) err(asher3);
  932.   av=avma;p=gash(x,prec);
  933.   gaffect(p,y);avma=av;
  934. }
  935.  
  936. /********************************************************************/
  937. /********************************************************************/
  938. /**                                                                **/
  939. /**            FONCTION ARGUMENT COSINUS HYPERBOLIQUE              **/
  940. /**                                                                **/
  941. /********************************************************************/
  942. /********************************************************************/
  943.  
  944. GEN     mpach(x)
  945.     
  946.      GEN     x;
  947.     
  948. {
  949.   long    l,av;
  950.   GEN     y,p1;
  951.  
  952.   if ((typ(x)!=2) || (gcmpgs(x,1)<0)) err(acher1);
  953.   l=lg(x);y=cgetr(l);
  954.   av=avma;
  955.   p1=cgetr(l+1);
  956.   affrr(x,p1);
  957.   mulrrz(p1,p1,p1);
  958.   subrsz(p1,1,p1);
  959.   affrr(mpsqrt(p1),p1);
  960.   addrrz(x,p1,p1);
  961.   affrr(mplog(p1),y);
  962.   avma=av;
  963.   return y;
  964. }
  965.  
  966. GEN     gach(x,prec)
  967.     
  968.      GEN     x;
  969.      long    prec;
  970.     
  971. {
  972.   long    av,tetpil,v;
  973.   GEN     y,p1;
  974.  
  975.   switch(typ(x))
  976.     {
  977.     case 2 : if(gcmpgs(x,1)>=0) y=mpach(x);
  978.     else
  979.       {
  980.         y=cgetg(3,6);
  981.         if(gcmpgs(x,-1)>=0)
  982.           {
  983.             y[2]=lmpacos(x);y[1]=zero;
  984.           }
  985.         else
  986.           {
  987.             av=avma;p1=mpach(gneg(x));tetpil=avma;
  988.             y[1]=lpile(av,tetpil,gneg(p1));
  989.             y[2]=lmppi(lg(x));
  990.           }
  991.       }
  992.       break;
  993.     
  994.     case 6 : av=avma;p1=gaddsg(-1,gmul(x,x));
  995.       p1=gadd(x,gsqrt(p1,prec));tetpil=avma;
  996.       y=glog(p1,prec);
  997.       if(signe(y[2])<0)
  998.         {
  999.           tetpil=avma;y=gneg(y);
  1000.         }
  1001.       y=gerepile(av,tetpil,y);
  1002.       break;
  1003.     
  1004.     case 3 :
  1005.     case 7 : err(acher2);
  1006.     
  1007.     case 11: av=avma;if(valp(x)<0) err(acher4);
  1008.       v=varn(x);p1=gdiv(deriv(x,v),gsqrt(gsubgs(gmul(x,x),1),prec));
  1009.       if(gcmp1(x[2])&&(!valp(x)))
  1010.         {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  1011.       else
  1012.         {
  1013.           y=integ(p1,v);
  1014.           if(valp(x))
  1015.             {
  1016.               p1=cgetg(3,6);p1[1]=zero;p1[2]=lmppi(prec);
  1017.               setexpo(p1[2],0);
  1018.             }
  1019.           else p1=gach(x[2],prec);
  1020.           tetpil=avma;y=gerepile(av,tetpil,gadd(p1,y));
  1021.         }
  1022.       break;
  1023.     
  1024.     default: y=transc(gach,x,prec);
  1025.     }
  1026.   return y;
  1027. }
  1028.  
  1029. void gachz(x,y)
  1030.     
  1031.      GEN     x,y;
  1032.     
  1033. {
  1034.   long    av,prec;
  1035.   GEN     p;
  1036.  
  1037.   prec=precision(y);
  1038.   if(!prec) err(acher3);
  1039.   av=avma;p=gach(x,prec);
  1040.   gaffect(p,y);avma=av;
  1041. }
  1042.  
  1043. /********************************************************************/
  1044. /********************************************************************/
  1045. /**                                                                **/
  1046. /**           FONCTION ARGUMENT TANGENTE HYPERBOLIQUE              **/
  1047. /**                                                                **/
  1048. /********************************************************************/
  1049. /********************************************************************/
  1050.  
  1051. GEN     mpath(x)
  1052.     
  1053.      GEN     x;
  1054.     
  1055. {
  1056.   long    av,tetpil;
  1057.   GEN     y,p1;
  1058.  
  1059.   if (typ(x)!=2) err(ather1);
  1060.   if (!signe(x))
  1061.     {
  1062.       y=cgetr(3);y[1]=x[1];y[2]=0;
  1063.     }
  1064.   else
  1065.     {
  1066.       av=avma;
  1067.       p1=addrs(divsr(2,subsr(1,x)),-1);
  1068.       tetpil=avma;y=gerepile(av,tetpil,mplog(p1));
  1069.       setexpo(y,expo(y)-1);
  1070.     }
  1071.   return y;
  1072. }
  1073.  
  1074. GEN     gath(x,prec)
  1075.     
  1076.      GEN     x;
  1077.      long    prec;
  1078.     
  1079. {
  1080.   long    av,tetpil,v;
  1081.   GEN     y,p1;
  1082.  
  1083.   switch(typ(x))
  1084.     {
  1085.     case 2 : if(expo(x)<0) y=mpath(x);
  1086.     else
  1087.       {
  1088.         av=avma;p1=addrs(divsr(2,addsr(-1,x)),1);
  1089.         tetpil=avma;y=cgetg(3,6);
  1090.         p1=mplog(p1);setexpo(p1,expo(p1)-1);
  1091.         y[1]=(long)p1;
  1092.         y[2]=lmppi(lg(x));setexpo(y[2],0);
  1093.         y=gerepile(av,tetpil,y);
  1094.       }
  1095.       break;
  1096.     
  1097.     case 6 : av=avma;
  1098.       p1=gaddgs(gdivsg(2,gsubsg(1,x)),-1);
  1099.       p1=glog(p1,prec);tetpil=avma;
  1100.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  1101.       break;
  1102.     
  1103.     case 3 :
  1104.     case 7 : err(ather2);
  1105.     
  1106.     case 11: av=avma;if(valp(x)<0) err(ather4);
  1107.       v=varn(x);p1=gdiv(deriv(x,v),gsubsg(1,gmul(x,x)));
  1108.       if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  1109.       else
  1110.         {
  1111.           y=integ(p1,v);p1=gath(x[2],prec);tetpil=avma;
  1112.           y=gerepile(av,tetpil,gadd(p1,y));
  1113.         }
  1114.       break;
  1115.     
  1116.     default: y=transc(gath,x,prec);
  1117.     }
  1118.   return y;
  1119. }
  1120.  
  1121. void gathz(x,y)
  1122.     
  1123.      GEN     x,y;
  1124.     
  1125. {
  1126.   long    av,prec;
  1127.   GEN     p;
  1128.  
  1129.   prec=precision(y);
  1130.   if(!prec) err(ather3);
  1131.   av=avma;p=gath(x,prec);
  1132.   gaffect(p,y);avma=av;
  1133. }
  1134.  
  1135. /********************************************************************/
  1136. /********************************************************************/
  1137. /**                                                                **/
  1138. /**             FONCTION TABLEAU DES NOMBRES DE BERNOULLI          **/
  1139. /**                                                                **/
  1140. /********************************************************************/
  1141. /********************************************************************/
  1142.  
  1143. void mpbern(nomb,prec)
  1144.     
  1145.      long    nomb,prec; /* pour calculer B_0,B_2,...,B_2*nomb */
  1146.     
  1147. {
  1148.   long    n,m,i,j,d1,d2,av;
  1149.   GEN     p1,p2;
  1150.  
  1151.   if(nomb<0) nomb=0;
  1152.   if (bernzone)
  1153.     {
  1154.       if ((bernzone[1]>=nomb)&&(bernzone[2]>=prec)) return;
  1155.       killbloc(bernzone);
  1156.     }
  1157.   bernzone=newbloc(3+prec*(nomb+1));
  1158.   bernzone[1]=nomb;
  1159.   bernzone[2]=prec;
  1160.   av=avma;
  1161.   p1=cgetr(prec+1);p2=cgetr(prec+1);
  1162.   *((GEN)bern(0))=0x2000000+prec;
  1163.   affsr(1,bern(0));
  1164.  
  1165.   for (i=1;i<=nomb;i++)
  1166.     {
  1167.       *((GEN)bern(i))=0x2000000+prec;
  1168.       affsr(0,p1);affsr(4,p2);
  1169.       n=8;m=5;d1=i-1;d2=2*i-3;
  1170.       for (j=i-1;j>0;--j)
  1171.         {
  1172.           addrrz(bern(j),p1,p1);
  1173.           mulsrz(n*m,p1,p1);
  1174.           divrsz(p1,d1*d2,p1);
  1175.           mulsrz(4,p2,p2);
  1176.           n+=4;m+=2;d1--;d2-=2;
  1177.         }
  1178.       addsrz(1,p1,p1);
  1179.       divrsz(p1,2*i+1,p1);
  1180.       subsrz(1,p1,p1);
  1181.       divrrz(p1,p2,bern(i));
  1182.     }
  1183.   avma=av;
  1184. }
  1185.  
  1186. GEN bernreal(n,prec)
  1187.      long n,prec;
  1188.  
  1189. {
  1190.   long n1; GEN p1;
  1191.   
  1192.   if(n==1) {affsr(-1,p1=cgetr(prec));setexpo(p1,-1);return p1;}
  1193.   if((n<0)||(n&1)) return gzero;
  1194.   n1=n>>1;mpbern(n1+1,prec);
  1195.   p1=cgetr(prec);affrr(bern(n1),p1);
  1196.   return p1;
  1197. }
  1198.  
  1199. GEN bernvec(nomb)
  1200.     
  1201.      long    nomb; /* pour calculer le vecteur B_0,B_2,...,B_2*nomb */
  1202.     
  1203. {
  1204.   long    n,m,i,j,d1,d2,av,tetpil;
  1205.   GEN     p1,p2,y;
  1206.  
  1207.   y=cgetg(nomb+2,17);y[1]=un;
  1208.   for (i=1;i<=nomb;i++)
  1209.     {
  1210.       av=avma;
  1211.       p1=gzero;p2=stoi(4);
  1212.       n=8;m=5;d1=i-1;d2=2*i-3;
  1213.       for (j=i-1;j>0;--j)
  1214.         {
  1215.           p1=gdivgs(gmulsg(n*m,gadd(p1,y[j+1])),d1*d2);
  1216.           p2=shifti(p2,2);
  1217.           n+=4;m+=2;d1--;d2-=2;
  1218.         }
  1219.       p1=gsubsg(1,gdivgs(gaddsg(1,p1),2*i+1));
  1220.       tetpil=avma;y[i+1]=lpile(av,tetpil,gdiv(p1,p2));
  1221.     }
  1222.   return y;
  1223. }
  1224.  
  1225. /********************************************************************/
  1226. /********************************************************************/
  1227. /**                                                                **/
  1228. /**                      FONCTION GAMMA                            **/
  1229. /**                                                                **/
  1230. /********************************************************************/
  1231. /********************************************************************/
  1232.  
  1233. GEN     mpgamma(x)
  1234.     
  1235.      GEN     x;
  1236.     
  1237. {
  1238.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av,av1;
  1239.   double  alpha,beta,dk;
  1240.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p9,p91,pitemp;
  1241.  
  1242.   if (typ(x)!=2) err(gamer1);
  1243.   l=lg(x);y=cgetr(l);s1=signe(x);if (!s1) err(gamer2);
  1244.   av=avma;p1=cgetr(l+1);u=expo(x);
  1245.   if ((u<-1) || (s1<0))
  1246.     {
  1247.       av1=avma;p2=gfrac(x);
  1248.       if(gcmp0(p2)) err(gamer2);
  1249.       avma=av1;subsrz(1,x,p1);
  1250.     }
  1251.   else affrr(x,p1);
  1252.   alpha=rtodbl(p1);beta=(16*LOG2*(l-2)/PI)-alpha;
  1253.   if (beta>=0) n=1+K2*beta;else n=0;
  1254.   if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1255.   else
  1256.     {
  1257.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1258.       if(beta>1.) beta+=(log(beta)/LOG2);
  1259.       p=16*(l-2)/beta+1;l2=l+1;
  1260.     }
  1261.   mpbern(p,l2);p91=cgetr(l2);
  1262.   p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);p5=cgetr(l2);p6=cgetr(l2);
  1263.   if(n) addsrz(n,p1,p2);else p2=p1;affrr(mplog(p2),p3);
  1264.   affsr(1,p4);setexpo(p4,expo(p4)-1);
  1265.   subrrz(p2,p4,p4);
  1266.   mulrrz(p4,p3,p4);subrrz(p4,p2,p4);
  1267.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p3);setexpo(p3,-1);
  1268.   addrrz(p4,p3,p4);mulrrz(p2,p2,p3);divsrz(1,p3,p3);e=expo(p3);
  1269.   setlg(p3,4);setlg(p5,4);setlg(p6,4);
  1270.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1271.   divrsz(p9,2*p*(2*p-1),p5);
  1272.   s=0;l1=2;l2-=2;
  1273.   for (k=p;k>1;k--)
  1274.     {
  1275.       mulrrz(p3,p5,p5);
  1276.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1277.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1278.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1279.       s &=31;
  1280.       setlg(p3,l1+2);setlg(p5,l1+2);setlg(p6,l1+2);
  1281.       addrrz(p6,p5,p5);
  1282.     }
  1283.   setlg(p5,l2+2);
  1284.   divrrz(p5,p2,p5);addrrz(p4,p5,p4);affrr(mpexp(p4),p4);
  1285.   for (i=1;i<=n;i++)
  1286.     {
  1287.       subrsz(p2,1,p2);divrrz(p4,p2,p4);
  1288.     }
  1289.   if ((u<-1)||(s1<0))
  1290.     {
  1291.       pitemp=mppi(l+1);mulrrz(pitemp,x,p1);affrr(mpsin(p1),p1);
  1292.       mulrrz(p1,p4,p4);divrrz(pitemp,p4,y);
  1293.     }
  1294.   else affrr(p4,y);
  1295.   avma=av;return y;
  1296. }
  1297.  
  1298. GEN     cxgamma(x,prec)
  1299.      GEN     x;
  1300.      long prec;
  1301.  
  1302. {
  1303.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av;
  1304.   double  alpha,beta,dk;
  1305.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p9,p91,pitemp;
  1306.  
  1307.   if (typ(x)!=6) err(gamer1);
  1308.   l=16383;if(typ(x[1])==2) l=precision(x[1]);
  1309.   if(typ(x[2])==2) {l1=precision(x[2]);if(l1<l) l=l1;}
  1310.   if(l==16383) l=prec;
  1311.   y=cgetg(3,6);y[1]=lgetr(l);y[2]=lgetr(l);
  1312.   s1=gsigne(x[1]);av=avma;
  1313.   p1=cgetg(3,6);p1[1]=lgetr(l+1);p1[2]=lgetr(l+1);
  1314.   if(s1||(typ(x[1])==2)) u=gexpo(x[1]);else u= -2;
  1315.   if ((s1<=0)||(u<-1)) gsubsgz(1,x,p1);
  1316.   else gaffect(x,p1);
  1317.   alpha=rtodbl(gabs(p1,4));beta=(16*LOG2*(l-2)/PI)-alpha;
  1318.   if (beta>=0) n=1+K2*beta;else n=0;
  1319.   if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1320.   else
  1321.     {
  1322.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1323.       if(beta>1.) beta+=(log(beta)/LOG2);
  1324.       p=16*(l-2)/beta+1;l2=l+1;
  1325.     }
  1326.   mpbern(p,l2);p91=cgetr(l2);
  1327.   p2=cgetg(3,6);p2[1]=lgetr(l2);p2[2]=lgetr(l2);
  1328.   p3=cgetg(3,6);p3[1]=lgetr(l2);p3[2]=lgetr(l2);
  1329.   p4=cgetg(3,6);p4[1]=lgetr(l2);p4[2]=lgetr(l2);
  1330.   p5=cgetg(3,6);p5[1]=lgetr(l2);p5[2]=lgetr(l2);
  1331.   p6=cgetr(l2);
  1332.   if(n) {addsrz(n,p1[1],p2[1]);affrr(p1[2],p2[2]);} else p2=p1;
  1333.   gaffect(glog(p2,l2),p3);
  1334.   affsr(1,p4[1]);setexpo(p4[1],-1);
  1335.   subrrz(p2[1],p4[1],p4[1]);p4[2]=lcopy(p2[2]);
  1336.   gmulz(p4,p3,p4);gsubz(p4,p2,p4);
  1337.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p6);setexpo(p6,-1);
  1338.   addrrz(p4[1],p6,p4[1]);gmulz(p2,p2,p3);gdivsgz(1,p3,p3);e=gexpo(p3);
  1339.   setlg(p3[1],4);setlg(p5[1],4);setlg(p6,4);setlg(p3[2],4);setlg(p5[2],4);
  1340.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1341.   gdivgsz(p9,2*p*(2*p-1),p5);
  1342.   s=0;l1=2;l2-=2;
  1343.   for (k=p;k>1;k--)
  1344.     {
  1345.       gmulz(p3,p5,p5);
  1346.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1347.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1348.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1349.       s &=31;
  1350.       setlg(p3[1],l1+2);setlg(p5[1],l1+2);setlg(p6,l1+2);
  1351.       setlg(p3[2],l1+2);setlg(p5[2],l1+2);
  1352.       addrrz(p6,p5[1],p5[1]);
  1353.     }
  1354.   setlg(p5[1],l2+2);setlg(p5[2],l2+2);
  1355.   gdivz(p5,p2,p5);gaddz(p4,p5,p4);gaffect(gexp(p4),p4);
  1356.   for (i=1;i<=n;i++)
  1357.     {
  1358.       subrsz(p2[1],1,p2[1]);gdivz(p4,p2,p4);
  1359.     }
  1360.   if ((s1<=0)||(u<-1))
  1361.     {
  1362.       pitemp=mppi(l+1);gmulz(pitemp,x,p1);gaffect(gsin(p1),p1);
  1363.       gmulz(p1,p4,p4);gdivz(pitemp,p4,y);
  1364.     }
  1365.   else gaffect(p4,y);
  1366.   avma=av;return y;
  1367. }
  1368.     
  1369. GEN     ggamma(x,prec)
  1370.     
  1371.      GEN     x;
  1372.      long    prec;
  1373.     
  1374. {
  1375.   long    i,lx;
  1376.   GEN     y;
  1377.  
  1378.     switch(typ(x))
  1379.       {
  1380.       case 1: if(signe(x)<=0) err(gamer2);
  1381.     y=transc(ggamma,x,prec);break;
  1382.       case 2 : y=mpgamma(x);break;
  1383.       case 6 : y=(gcmp0(x[2])) ? ggamma(x[1],prec) : cxgamma(x,prec);break;
  1384.       case 7 : err(impl,"p-adic gamma function");
  1385.       case 3 : err(gamer3);
  1386.       case 11: y=gexp(glngamma(x,prec));break;
  1387.       case 17:
  1388.       case 18:
  1389.       case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1390.     for(i=1;i<lx;i++)
  1391.       y[i]=lgamma(x[i],prec);
  1392.     break;
  1393.       default: y=transc(ggamma,x,prec);
  1394.       }
  1395.   return y;
  1396. }
  1397.  
  1398. void ggammaz(x,y)
  1399.     
  1400.      GEN     x,y;
  1401.     
  1402. {
  1403.   long    av,prec;
  1404.   GEN     p;
  1405.  
  1406.   prec=precision(y);
  1407.   if(!prec) err(gamer4);
  1408.   av=avma;p=ggamma(x,prec);
  1409.   gaffect(p,y);avma=av;
  1410. }
  1411.  
  1412. GEN     mplngamma(x)
  1413.     
  1414.      GEN     x;
  1415.     
  1416. {
  1417.   long    l,l1,l2,u,i,k,e,f,s,s1,n,p,av,av1;
  1418.   double  alpha,beta,dk;
  1419.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p8,p9,p91,pitemp;
  1420.  
  1421.   if (typ(x)!=2) err(gamer1);
  1422.   l=lg(x);y=cgetr(l);s1=signe(x);if (!s1) err(gamer2);
  1423.   av=avma;p1=cgetr(l+1);u=expo(x);
  1424.   if ((u<-1) || (s1<0))
  1425.     {
  1426.       av1=avma;p2=gfrac(x);
  1427.       if(gcmp0(p2)) err(gamer2);
  1428.       avma=av1;subsrz(1,x,p1);
  1429.     }
  1430.   else affrr(x,p1);
  1431.   if(expo(p1)>1000) 
  1432.     {
  1433.       n=0;beta=log(K4/(l-2))/LOG2+expo(p1);beta+=(log(beta)/LOG2);
  1434.       p=16*(l-2)/beta+1;l2=l+1;
  1435.     }
  1436.   else
  1437.     {
  1438.       alpha=rtodbl(p1);beta=(16*LOG2*(l-2)/PI)-alpha;
  1439.       if (beta>=0) n=1+K2*beta;else n=0;
  1440.       if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1441.       else
  1442.     {
  1443.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1444.       if(beta>1.) beta+=(log(beta)/LOG2);
  1445.       p=16*(l-2)/beta+1;l2=l+1;
  1446.     }
  1447.     }
  1448.   mpbern(p,l2);p91=cgetr(l2);
  1449.   p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);p5=cgetr(l2);p6=cgetr(l2);
  1450.   p8=cgetr(l2);if(n) addsrz(n,p1,p2);else p2=p1;affrr(mplog(p2),p3);
  1451.   affsr(1,p4);affsr(1,p8);setexpo(p4,expo(p4)-1);
  1452.   subrrz(p2,p4,p4);mulrrz(p4,p3,p4);subrrz(p4,p2,p4);
  1453.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p3);setexpo(p3,-1);
  1454.   addrrz(p4,p3,p4);mulrrz(p2,p2,p3);divsrz(1,p3,p3);e=expo(p3);
  1455.   setlg(p3,4);setlg(p5,4);setlg(p6,4);
  1456.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1457.   divrsz(p9,2*p*(2*p-1),p5);
  1458.   s=0;l1=2;l2-=2;
  1459.   for (k=p;k>1;k--)
  1460.     {
  1461.       mulrrz(p3,p5,p5);
  1462.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1463.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1464.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1465.       s &=31;
  1466.       setlg(p3,l1+2);setlg(p5,l1+2);setlg(p6,l1+2);
  1467.       addrrz(p6,p5,p5);
  1468.     }
  1469.   setlg(p5,l2+2);
  1470.   divrrz(p5,p2,p5);addrrz(p4,p5,p4);
  1471.   for (i=1;i<=n;i++)
  1472.     {
  1473.       subrsz(p2,1,p2);mulrrz(p8,p2,p8);
  1474.     }
  1475.   f=signe(p8);subrrz(p4,mplog((f>0)?p8:negr(p8)),p4);
  1476.   if ((u<-1)||(s1<0))
  1477.     {
  1478.       pitemp=mppi(l+1);mulrrz(pitemp,x,p1);divrrz(pitemp,mpsin(p1),p1);
  1479.       f*=signe(p1);subrrz(mplog(absr(p1)),p4,y);
  1480.     }
  1481.   else affrr(p4,y);
  1482.   avma=av;if(f<0) {p2=cgetg(3,6);p2[1]=(long)y;p2[2]=(long)mppi(l);return p2;}
  1483.   else return y;
  1484. }
  1485.  
  1486. GEN     cxlngamma(x,prec)
  1487.      GEN     x;
  1488.      long prec;
  1489.  
  1490. {
  1491.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av,flag;
  1492.   double  alpha,beta,dk;
  1493.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p8,p9,p91,pitemp;
  1494.  
  1495.   if (typ(x)!=6) err(gamer1);
  1496.   l=16383;if(typ(x[1])==2) l=precision(x[1]);
  1497.   if(typ(x[2])==2) {l1=precision(x[2]);if(l1<l) l=l1;}
  1498.   if(l==16383) l=prec;
  1499.   y=cgetg(3,6);y[1]=lgetr(l);y[2]=lgetr(l);
  1500.   s1=gsigne(x[1]);av=avma;
  1501.   p1=cgetg(3,6);p1[1]=lgetr(l+1);p1[2]=lgetr(l+1);
  1502.   if(s1||(typ(x[1])==2)) u=gexpo(x[1]);else u= -2;
  1503.   if (((s1<=0)||(u<-1))&&(!gcmp0(x[2])&&(gexpo(x[2])<=16)))
  1504.     {gsubsgz(1,x,p1);flag=1;}
  1505.   else {gaffect(x,p1);flag=0;}
  1506.   p2=gabs(p1,4);
  1507.   if(expo(p2)>1000) 
  1508.     {
  1509.       n=0;beta=log(K4/(l-2))/LOG2+expo(p1);beta+=(log(beta)/LOG2);
  1510.       p=16*(l-2)/beta+1;l2=l+1;
  1511.     }
  1512.   else
  1513.     {
  1514.       alpha=rtodbl(p2);beta=(16*LOG2*(l-2)/PI)-alpha;
  1515.       if (beta>=0) n=1+K2*beta;else n=0;
  1516.       if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1517.       else
  1518.     {
  1519.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1520.       if(beta>1.) beta+=(log(beta)/LOG2);
  1521.       p=16*(l-2)/beta+1;l2=l+1;
  1522.     }
  1523.     }
  1524.   mpbern(p,l2);p91=cgetr(l2);
  1525.   p2=cgetg(3,6);p2[1]=lgetr(l2);p2[2]=lgetr(l2);
  1526.   p3=cgetg(3,6);p3[1]=lgetr(l2);p3[2]=lgetr(l2);
  1527.   p4=cgetg(3,6);p4[1]=lgetr(l2);p4[2]=lgetr(l2);
  1528.   p5=cgetg(3,6);p5[1]=lgetr(l2);p5[2]=lgetr(l2);
  1529.   p8=cgetg(3,6);p8[1]=lgetr(l2);p8[2]=lgetr(l2);
  1530.   gaffsg(1,p8);p6=cgetr(l2);
  1531.   if(n) {addsrz(n,p1[1],p2[1]);affrr(p1[2],p2[2]);} else p2=p1;
  1532.   gaffect(glog(p2,l2),p3);
  1533.   affsr(1,p4[1]);setexpo(p4[1],-1);
  1534.   subrrz(p2[1],p4[1],p4[1]);p4[2]=lcopy(p2[2]);
  1535.   gmulz(p4,p3,p4);gsubz(p4,p2,p4);
  1536.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p6);setexpo(p6,-1);
  1537.   addrrz(p4[1],p6,p4[1]);gmulz(p2,p2,p3);gdivsgz(1,p3,p3);e=gexpo(p3);
  1538.   setlg(p3[1],4);setlg(p5[1],4);setlg(p6,4);setlg(p3[2],4);setlg(p5[2],4);
  1539.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1540.   gdivgsz(p9,2*p*(2*p-1),p5);
  1541.   s=0;l1=2;l2-=2;
  1542.   for (k=p;k>1;k--)
  1543.     {
  1544.       gmulz(p3,p5,p5);
  1545.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1546.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1547.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1548.       s &=31;
  1549.       setlg(p3[1],l1+2);setlg(p5[1],l1+2);setlg(p6,l1+2);
  1550.       setlg(p3[2],l1+2);setlg(p5[2],l1+2);
  1551.       addrrz(p6,p5[1],p5[1]);
  1552.     }
  1553.   setlg(p5[1],l2+2);setlg(p5[2],l2+2);
  1554.   gdivz(p5,p2,p5);gaddz(p4,p5,p4);
  1555.   for (i=1;i<=n;i++)
  1556.     {
  1557.       subrsz(p2[1],1,p2[1]);gmulz(p8,p2,p8);
  1558.     }
  1559.   gsubz(p4,glog(p8),p4);
  1560.   if (flag)
  1561.     {
  1562.       pitemp=mppi(l+1);gmulz(pitemp,x,p1);gdivz(pitemp,gsin(p1),p1);
  1563.       gsubz(glog(p1),p4,y);
  1564.     }
  1565.   else gaffect(p4,y);
  1566.   avma=av;return y;
  1567. }
  1568.     
  1569. GEN     glngamma(x,prec)
  1570.     
  1571.      GEN     x;
  1572.      long    prec;
  1573.     
  1574. {
  1575.   long    i,lx,av,tetpil,v,n;
  1576.   GEN     y,p1;
  1577.  
  1578.     switch(typ(x))
  1579.       {
  1580.       case 1: if(signe(x)<=0) err(gamer2);
  1581.     y=transc(glngamma,x,prec);break;
  1582.       case 2 : y=mplngamma(x);break;
  1583.       case 6 : y=(gcmp0(x[2])) ? glngamma(x[1],prec) : cxlngamma(x,prec);break;
  1584.       case 7 : err(impl,"p-adic lngamma function");
  1585.       case 3 : err(gamer3);
  1586.       case 11: av=avma;if(valp(x)) err(loger5);
  1587.     v=varn(x);if(!gcmp1(x[2])) err(impl,"lngamma around a!=1");
  1588.     p1=gsubsg(1,x);n=(lg(x)-3)/valp(p1);
  1589.     y=ggrando(polx[v],lg(x)-2);
  1590.     for(i=n;i>=2;i--)
  1591.       {
  1592.         y=gmul(p1,gadd(gdivgs(izeta(stoi(i),prec),i),y));
  1593.       }
  1594.     y=gadd(mpeuler(prec),y);tetpil=avma;
  1595.     y=gerepile(av,tetpil,gmul(p1,y));break;
  1596.       case 17:
  1597.       case 18:
  1598.       case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1599.     for(i=1;i<lx;i++)
  1600.       y[i]=(long)glngamma(x[i],prec);
  1601.     break;
  1602.       default: y=transc(glngamma,x,prec);
  1603.       }
  1604.   return y;
  1605. }
  1606.  
  1607. void glngammaz(x,y)
  1608.     
  1609.      GEN     x,y;
  1610.     
  1611. {
  1612.   long    av,prec;
  1613.   GEN     p;
  1614.  
  1615.   prec=precision(y);
  1616.   if(!prec) err(gamer4);
  1617.   av=avma;p=glngamma(x,prec);
  1618.   gaffect(p,y);avma=av;
  1619. }
  1620.  
  1621. /********************************************************************/
  1622. /********************************************************************/
  1623. /**                                                                **/
  1624. /**             FONCTION GAMMA DES DEMI-ENTIERS                    **/
  1625. /**                                                                **/
  1626. /********************************************************************/
  1627. /********************************************************************/
  1628.  
  1629. GEN     mpgamd(x,prec)
  1630.     
  1631.      long    x,prec;
  1632.     
  1633. {
  1634.   long    i,j,a,l,av;
  1635.   GEN     y,p1;
  1636.  
  1637.   a=abs(x);
  1638.   l=prec+1+a/32;
  1639.   if (l>32767) err(gamder1);
  1640.   y=cgetr(prec);
  1641.   av=avma;
  1642.   p1=mpsqrt(mppi(l));
  1643.   if (x>=0)
  1644.     {
  1645.       j= -1;
  1646.       for (i=0;i<x;i++)
  1647.         {
  1648.           j+=2;mulsrz(j,p1,p1);
  1649.           setexpo(p1,expo(p1)-1);
  1650.         }
  1651.     }
  1652.   else
  1653.     {
  1654.       j=1;
  1655.       for (i=0;i<a;i++)
  1656.         {
  1657.           j-=2;divrsz(p1,j,p1);
  1658.           setexpo(p1,expo(p1)+1);
  1659.         }
  1660.     }
  1661.   affrr(p1,y);
  1662.   avma=av;
  1663.   return y;
  1664. }
  1665.  
  1666. void mpgamdz(s,y)
  1667.     
  1668.      long    s;
  1669.      GEN     y;
  1670.     
  1671. {
  1672.   long    l,av;
  1673.   GEN     p1;
  1674.  
  1675.   av=avma;
  1676.   l=lg(y);
  1677.   p1=mpgamd(s,l);
  1678.   affrr(p1,y);
  1679.   avma=av;
  1680. }
  1681.  
  1682. GEN     ggamd(x,prec)
  1683.     
  1684.      GEN     x;
  1685.      long    prec;
  1686.     
  1687. {
  1688.   long    av,tetpil,i,lx;
  1689.   GEN     y,p1;
  1690.  
  1691.   switch(typ(x))
  1692.     {
  1693.     case 1 : y=mpgamd(itos(x),prec);break;
  1694.     case 2 :
  1695.     case 4 :
  1696.     case 5 :
  1697.     case 6 :
  1698.     case 8 : av=avma;p1=gadd(x,ghalf);tetpil=avma;
  1699.       y=gerepile(av,tetpil,ggamma(p1,prec));
  1700.       break;
  1701.     case 17:
  1702.     case 18:
  1703.     case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1704.       for(i=1;i<lx;i++)
  1705.         y[i]=lgamd(x[i],prec);
  1706.       break;
  1707.     case 3 :
  1708.     case 7 : err(gamder2);
  1709.     case 11: err(impl,"gamd of a power series");
  1710.     default: y=transc(ggamd,x,prec);
  1711.     }
  1712.   return y;
  1713. }
  1714.  
  1715. void ggamdz(x,y)
  1716.     
  1717.      GEN     x,y;
  1718.     
  1719. {
  1720.   long    av,prec;
  1721.   GEN     p;
  1722.  
  1723.   prec=precision(y);
  1724.   if(!prec) err(gamder3);
  1725.   av=avma;p=ggamd(x,prec);
  1726.   gaffect(p,y);avma=av;
  1727. }
  1728.  
  1729.  
  1730. /********************************************************************/
  1731. /********************************************************************/
  1732. /**                                                                **/
  1733. /**                      FONCTION PSI                              **/
  1734. /**                                                                **/
  1735. /********************************************************************/
  1736. /********************************************************************/
  1737.  
  1738.  
  1739. GEN     mppsi(z)        /*   version p=2 */
  1740.     
  1741.      GEN        z;
  1742.     
  1743. {
  1744.   long  l,n,k,x,xx,av1,av2,tetpil;
  1745.   GEN   zk,u,v,a,b;
  1746.  
  1747.   av1=avma;l=lg(z);
  1748.   x=1+16*(l-2)*LOG2+1.58*rtodbl(z);xx=x*x;
  1749.   n=1+3.591*x;
  1750.   affsr(x,a=cgetr(l));
  1751.   a=mplog(a);
  1752.   gaffect(a,u=cgetr(l));
  1753.   gaffsg(1,b=cgetr(l));
  1754.   gaffsg(1,v=cgetr(l));
  1755.   for (k=1;k<=n;k++)
  1756.     {
  1757.       av2=avma;
  1758.       zk=(k>1) ? gaddsg(k-1,z) : z;
  1759.       gdivz(gmulsg(xx,b),gmul(zk,zk),b);
  1760.       gdivz(gsub(gdiv(gmulsg(xx,a),zk),b),zk,a);
  1761.       gaddz(u,a,u);gaddz(v,b,v);
  1762.       avma=av2;
  1763.     }
  1764.   tetpil=avma;return gerepile(av1,tetpil,gdiv(u,v));
  1765. }
  1766.  
  1767. GEN     cxpsi(z,prec)        /*   version p=2 */
  1768.      GEN        z;
  1769.      long prec;
  1770.  
  1771. {
  1772.   long  l,l1,n,k,x,xx,av1,av2,tetpil;
  1773.   GEN   zk,u,v,a,b;
  1774.  
  1775.   l=16383;if(typ(z[1])==2) l=precision(z[1]);
  1776.   if(typ(z[2])==2) {l1=precision(z[2]);if(l1<l) l=l1;}
  1777.   if(l==16383) l=prec;
  1778.   av1=avma;x=1+16*(l-2)*LOG2+1.58*rtodbl(gabs(z,4));xx=x*x;
  1779.   n=1+3.591*x;
  1780.   a=cgetg(3,6);a[1]=lgetr(l);a[2]=lgetr(l);gaffsg(x,a);
  1781.   b=cgetg(3,6);b[1]=lgetr(l);b[2]=lgetr(l);gaffsg(1,b);
  1782.   u=cgetg(3,6);u[1]=lgetr(l);u[2]=lgetr(l);
  1783.   v=cgetg(3,6);v[1]=lgetr(l);v[2]=lgetr(l);gaffsg(1,v);
  1784.   a=glog(a);gaffect(a,u);
  1785.   for (k=1;k<=n;k++)
  1786.     {
  1787.       av2=avma;
  1788.       zk=(k>1) ? gaddsg(k-1,z) : z;
  1789.       gdivz(gmulsg(xx,b),gmul(zk,zk),b);
  1790.       gdivz(gsub(gdiv(gmulsg(xx,a),zk),b),zk,a);
  1791.       gaddz(u,a,u);gaddz(v,b,v);
  1792.       avma=av2;
  1793.     }
  1794.   tetpil=avma;return gerepile(av1,tetpil,gdiv(u,v));
  1795. }
  1796.  
  1797.  
  1798.  
  1799. GEN     gpsi(x,prec)
  1800.     
  1801.      GEN     x;
  1802.      long    prec;
  1803.     
  1804. {
  1805.   long    i,lx;
  1806.   GEN     y;
  1807.  
  1808.   switch(typ(x))
  1809.     {
  1810.     case 2 : y=mppsi(x);break;
  1811.     case 6 : y=cxpsi(x,prec);break;
  1812.     case 3 :
  1813.     case 7 : err(psier1);
  1814.     case 11: err(impl,"psi of power series");
  1815.     case 17:
  1816.     case 18:
  1817.     case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1818.       for(i=1;i<lx;i++)
  1819.         y[i]=lpsi(x[i],prec);
  1820.       break;
  1821.     default: y=transc(gpsi,x,prec);
  1822.     }
  1823.   return y;
  1824. }
  1825.  
  1826. void gpsiz(x,y)
  1827.     
  1828.      GEN     x,y;
  1829.     
  1830. {
  1831.   long    av,prec;
  1832.   GEN     p;
  1833.  
  1834.   prec=precision(y);
  1835.   if(!prec) err(psier2);
  1836.   av=avma;p=gpsi(x,prec);
  1837.   gaffect(p,y);avma=av;
  1838. }
  1839.