home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / trans2 < prev    next >
Encoding:
Text File  |  1991-11-28  |  43.4 KB  |  1,838 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. 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. 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. 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.     
  538.      GEN        x,y;
  539. {
  540.  
  541.   long tetpil,av=avma;
  542.   GEN p1;
  543.  
  544.   switch(typ(x))
  545.     {
  546.     case 1 :
  547.     case 4 :
  548.     case 5 : p1=cgetr(prec);gaffect(x,p1);
  549.       x=p1;break;
  550.     default:;
  551.     }
  552.   switch(typ(y))
  553.     {
  554.     case 1 :
  555.     case 4 :
  556.     case 5 : p1=cgetr(prec);gaffect(y,p1);
  557.       y=p1;break;
  558.     default:;
  559.     }
  560.   if (av==(tetpil=avma)) return mparg(x,y);
  561.   else return gerepile(av,tetpil,mparg(x,y));
  562. }
  563.  
  564. GEN     garg(x,prec)
  565.     
  566.      GEN   x;
  567.      long  prec;
  568. {
  569.   GEN  p1,y;
  570.   long av,tetpil;
  571.  
  572.   if(gcmp0(x)) err(arger2);
  573.   switch(typ(x))
  574.     {
  575.     case 2 : prec=lg(x);
  576.     case 1 :
  577.     case 4 :
  578.     case 5 :
  579.       if(signe(x)>0)
  580.         {
  581.           y=cgetr(3);y[1]=0x800000-((prec-2)<<5);
  582.           y[2]=zero;
  583.         }
  584.       else y=mppi(prec);
  585.       break;
  586.     case 8 : av=avma;gaffsg(1,p1=cgetr(prec));p1=gmul(p1,x);
  587.       tetpil=avma;y=gerepile(av,tetpil,garg(p1,prec));
  588.       break;
  589.     case 6 : y=sarg(x[1],x[2],prec);
  590.       break;
  591.     default: err(arger1);
  592.     }
  593.   return y;
  594. }
  595.  
  596. /********************************************************************/
  597. /********************************************************************/
  598. /**                                                                **/
  599. /**                 FONCTION COSINUS HYPERBOLIQUE                  **/
  600. /**                                                                **/
  601. /********************************************************************/
  602. /********************************************************************/
  603.  
  604. GEN     mpch(x)
  605.     
  606.      GEN     x;
  607.     
  608. {
  609.   long    l,av;
  610.   GEN     y,p1,p2;
  611.  
  612.   if (typ(x)!=2) err(cher3);
  613.   if(gcmp0(x)) y=gaddsg(1,x);
  614.   else
  615.     {
  616.       l=lg(x);y=cgetr(l);
  617.       av=avma;
  618.       p1=cgetr(l+1);p2=cgetr(l+1);
  619.       affrr(mpexp1(x),p1);
  620.       addsrz(1,p1,p2);
  621.       divsrz(1,p2,p2);
  622.       addrrz(p1,p2,p2);
  623.       addsrz(1,p2,y);
  624.       setexpo(y,expo(y)-1);
  625.       avma=av;
  626.     }
  627.   return y;
  628. }
  629.  
  630. GEN     gch(x,prec)
  631.     
  632.      GEN     x;
  633.      long    prec;
  634.     
  635. {
  636.   long    av,tetpil;
  637.   GEN     y,p1;
  638.  
  639.   switch(typ(x))
  640.     {
  641.     case 2 : y=mpch(x);break;
  642.     
  643.     case 6 : av=avma;p1=gexp(x,prec);
  644.       p1=gadd(p1,ginv(p1));tetpil=avma;
  645.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  646.       break;
  647.     
  648.     case 3 :
  649.     case 7 : err(cher1);
  650.     
  651.     case 11: av=avma;p1=gexp(x,prec);p1=gadd(p1,gdivsg(1,p1));
  652.       tetpil=avma;y=gerepile(av,tetpil,gmul2n(p1,-1));
  653.       break;
  654.     
  655.     default: y=transc(gch,x,prec);
  656.     }
  657.   return y;
  658. }
  659.  
  660. gchz(x,y)
  661.     
  662.      GEN     x,y;
  663.     
  664. {
  665.   long    av,prec;
  666.   GEN     p;
  667.  
  668.   prec=precision(y);
  669.   if(!prec) err(cher2);
  670.   av=avma;p=gch(x,prec);
  671.   gaffect(p,y);avma=av;
  672. }
  673.  
  674. /********************************************************************/
  675. /********************************************************************/
  676. /**                                                                **/
  677. /**                 FONCTION SINUS HYPERBOLIQUE                    **/
  678. /**                                                                **/
  679. /********************************************************************/
  680. /********************************************************************/
  681.  
  682. GEN     mpsh(x)
  683.     
  684.      GEN     x;
  685.     
  686. {
  687.   long    l,av;
  688.   GEN     y,p1,p2;
  689.  
  690.   if (typ(x)!=2) err(sher3);
  691.   if(!signe(x))
  692.     {
  693.       y=cgetr(3);y[1]=x[1];y[2]=0;
  694.     }
  695.   else
  696.     {
  697.       l=lg(x);y=cgetr(l);
  698.       av=avma;
  699.       p1=cgetr(l+1);p2=cgetr(l+1);
  700.       affrr(mpexp1(x),p1);
  701.       addsrz(1,p1,p2);
  702.       divrrz(p1,p2,p2);
  703.       addrrz(p1,p2,y);
  704.       setexpo(y,expo(y)-1);
  705.       avma=av;
  706.     }
  707.   return y;
  708. }
  709.  
  710. GEN     gsh(x,prec)
  711.     
  712.      GEN     x;
  713.      long    prec;
  714.     
  715. {
  716.   long    av,tetpil;
  717.   GEN     y,p1;
  718.  
  719.   switch(typ(x))
  720.     {
  721.     case 2 : y=mpsh(x);break;
  722.     
  723.     case 6 : av=avma;p1=gexp(x,prec);
  724.       p1=gsub(p1,ginv(p1));tetpil=avma;
  725.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  726.       break;
  727.     
  728.     case 3 :
  729.     case 7 : err(sher1);
  730.     
  731.     case 11: if(gcmp0(x)) y=gcopy(x);
  732.     else
  733.       {
  734.         av=avma;p1=gexp(x,prec);p1=gsub(p1,gdivsg(1,p1));
  735.         tetpil=avma;y=gerepile(av,tetpil,gmul2n(p1,-1));
  736.       }
  737.       break;
  738.     
  739.     default: y=transc(gsh,x,prec);
  740.     }
  741.   return y;
  742. }
  743.  
  744. gshz(x,y)
  745.     
  746.      GEN     x,y;
  747.     
  748. {
  749.   long    av,prec;
  750.   GEN     p;
  751.  
  752.   prec=precision(y);
  753.   if(!prec) err(sher2);
  754.   av=avma;p=gsh(x,prec);
  755.   gaffect(p,y);avma=av;
  756. }
  757.  
  758. /********************************************************************/
  759. /********************************************************************/
  760. /**                                                                **/
  761. /**                 FONCTION TANGENTE HYPERBOLIQUE                 **/
  762. /**                                                                **/
  763. /********************************************************************/
  764. /********************************************************************/
  765.  
  766. GEN     mpth(x)
  767.     
  768.      GEN     x;
  769.     
  770. {
  771.   long    l,av;
  772.   GEN     y,p1,p2;
  773.  
  774.   if (typ(x)!=2) err(ther1);
  775.   if (!signe(x))
  776.     {
  777.       y=cgetr(3);y[1]=x[1];y[2]=0;
  778.     }
  779.   else
  780.     {
  781.       l=lg(x);y=cgetr(l);
  782.       av=avma;
  783.       p1=cgetr(l+1);p2=cgetr(l+1);
  784.       affrr(x,p1);
  785.       setexpo(p1,expo(p1)+1);
  786.       affrr(mpexp1(p1),p1);
  787.       addsrz(2,p1,p2);
  788.       divrrz(p1,p2,y);
  789.       avma=av;
  790.     }
  791.   return y;
  792. }
  793.  
  794. GEN     gth(x,prec)
  795.     
  796.      GEN     x;
  797.      long    prec;
  798.     
  799. {
  800.   long    av,tetpil;
  801.   GEN     y,p1,p2,p3;
  802.  
  803.   switch(typ(x))
  804.     {
  805.     case 2 : y=mpth(x);break;
  806.     
  807.     case 6 : av=avma;p1=gexp(gmul2n(x,1),prec);
  808.       p1=gdivsg(-2,gaddgs(p1,1));tetpil=avma;
  809.       y=gerepile(av,tetpil,gaddsg(1,p1));
  810.       break;
  811.     
  812.     case 3 :
  813.     case 7 : err(ther2);
  814.     
  815.     case 11: if(gcmp0(x)) y=gcopy(x);
  816.     else
  817.       {
  818.         av=avma;p1=gexp(gmul2n(x ,1),prec);
  819.         p2=gsubgs(p1,1);p3=gaddgs(p1,1);
  820.         tetpil=avma;y=gerepile(av,tetpil,gdiv(p2,p3));
  821.       }
  822.       break;
  823.     
  824.     default: y=transc(gth,x,prec);
  825.     }
  826.   return y;
  827. }
  828.  
  829. gthz(x,y)
  830.     
  831.      GEN     x,y;
  832.     
  833. {
  834.   long    av,prec;
  835.   GEN     p;
  836.  
  837.   prec=precision(y);
  838.   if(!prec) err(ther3);
  839.   av=avma;p=gth(x,prec);
  840.   gaffect(p,y);avma=av;
  841. }
  842.  
  843. /********************************************************************/
  844. /********************************************************************/
  845. /**                                                                **/
  846. /**             FONCTION ARGUMENT SINUS HYPERBOLIQUE               **/
  847. /**                                                                **/
  848. /********************************************************************/
  849. /********************************************************************/
  850.  
  851. GEN     mpash(x)
  852.     
  853.      GEN     x;
  854.     
  855. {
  856.   long    l,av;
  857.   GEN     y,p1;
  858.  
  859.   if (typ(x)!=2) err(asher1);
  860.   if (!signe(x))
  861.     {
  862.       y=cgetr(3);y[1]=x[1];y[2]=0;
  863.     }
  864.   else
  865.     {
  866.       l=lg(x);y=cgetr(l);
  867.       av=avma;
  868.       p1=cgetr(l+1);
  869.       affrr(x,p1);
  870.       mulrrz(p1,p1,p1);
  871.       addsrz(1,p1,p1);
  872.       affrr(mpsqrt(p1),p1);
  873.       addrrz(x,p1,p1);
  874.       affrr(mplog(p1),y);
  875.       avma=av;
  876.     }
  877.   return y;
  878. }
  879.  
  880. GEN     gash(x,prec)
  881.     
  882.      GEN     x;
  883.      long    prec;
  884.     
  885. {
  886.   long    av,tetpil,v;
  887.   GEN     y,p1;
  888.  
  889.   switch(typ(x))
  890.     {
  891.     case 2 : y=mpash(x);break;
  892.     
  893.     case 6 : av=avma;p1=gaddsg(1,gmul(x,x));
  894.       p1=gadd(x,gsqrt(p1,prec));
  895.       p1=glog(gmul(p1,p1),prec);tetpil=avma;
  896.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  897.       break;
  898.     
  899.     case 3 :
  900.     case 7 : err(asher2);
  901.     
  902.     case 11: if(gcmp0(x)) y=gcopy(x);
  903.     else
  904.       {
  905.         av=avma;if(valp(x)<0) err(asher4);
  906.         v=varn(x);p1=gdiv(deriv(x,v),gsqrt(gaddsg(1,gmul(x,x))));
  907.         if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  908.         else
  909.           {
  910.             y=integ(p1,v);p1=gash(x[2],prec);tetpil=avma;
  911.             y=gerepile(av,tetpil,gadd(p1,y));
  912.           }
  913.       }
  914.       break;
  915.     
  916.     default: y=transc(gash,x,prec);
  917.     }
  918.   return y;
  919. }
  920.  
  921. gashz(x,y)
  922.     
  923.      GEN     x,y;
  924.     
  925. {
  926.   long    av,prec;
  927.   GEN     p;
  928.  
  929.   prec=precision(y);
  930.   if(!prec) err(asher3);
  931.   av=avma;p=gash(x,prec);
  932.   gaffect(p,y);avma=av;
  933. }
  934.  
  935. /********************************************************************/
  936. /********************************************************************/
  937. /**                                                                **/
  938. /**            FONCTION ARGUMENT COSINUS HYPERBOLIQUE              **/
  939. /**                                                                **/
  940. /********************************************************************/
  941. /********************************************************************/
  942.  
  943. GEN     mpach(x)
  944.     
  945.      GEN     x;
  946.     
  947. {
  948.   long    l,av;
  949.   GEN     y,p1;
  950.  
  951.   if ((typ(x)!=2) || (gcmpgs(x,1)<0)) err(acher1);
  952.   l=lg(x);y=cgetr(l);
  953.   av=avma;
  954.   p1=cgetr(l+1);
  955.   affrr(x,p1);
  956.   mulrrz(p1,p1,p1);
  957.   subrsz(p1,1,p1);
  958.   affrr(mpsqrt(p1),p1);
  959.   addrrz(x,p1,p1);
  960.   affrr(mplog(p1),y);
  961.   avma=av;
  962.   return y;
  963. }
  964.  
  965. GEN     gach(x,prec)
  966.     
  967.      GEN     x;
  968.      long    prec;
  969.     
  970. {
  971.   long    av,tetpil,v;
  972.   GEN     y,p1;
  973.  
  974.   switch(typ(x))
  975.     {
  976.     case 2 : if(gcmpgs(x,1)>=0) y=mpach(x);
  977.     else
  978.       {
  979.         y=cgetg(3,6);
  980.         if(gcmpgs(x,-1)>=0)
  981.           {
  982.             y[2]=lmpacos(x);y[1]=zero;
  983.           }
  984.         else
  985.           {
  986.             av=avma;p1=mpach(gneg(x));tetpil=avma;
  987.             y[1]=lpile(av,tetpil,gneg(p1));
  988.             y[2]=lmppi(lg(x));
  989.           }
  990.       }
  991.       break;
  992.     
  993.     case 6 : av=avma;p1=gaddsg(-1,gmul(x,x));
  994.       p1=gadd(x,gsqrt(p1,prec));tetpil=avma;
  995.       y=glog(p1,prec);
  996.       if(signe(y[2])<0)
  997.         {
  998.           tetpil=avma;y=gneg(y);
  999.         }
  1000.       y=gerepile(av,tetpil,y);
  1001.       break;
  1002.     
  1003.     case 3 :
  1004.     case 7 : err(acher2);
  1005.     
  1006.     case 11: av=avma;if(valp(x)<0) err(acher4);
  1007.       v=varn(x);p1=gdiv(deriv(x,v),gsqrt(gsubgs(gmul(x,x),1),prec));
  1008.       if(gcmp1(x[2])&&(!valp(x)))
  1009.         {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  1010.       else
  1011.         {
  1012.           y=integ(p1,v);
  1013.           if(valp(x))
  1014.             {
  1015.               p1=cgetg(3,6);p1[1]=zero;p1[2]=lmppi(prec);
  1016.               setexpo(p1[2],0);
  1017.             }
  1018.           else p1=gach(x[2],prec);
  1019.           tetpil=avma;y=gerepile(av,tetpil,gadd(p1,y));
  1020.         }
  1021.       break;
  1022.     
  1023.     default: y=transc(gach,x,prec);
  1024.     }
  1025.   return y;
  1026. }
  1027.  
  1028. gachz(x,y)
  1029.     
  1030.      GEN     x,y;
  1031.     
  1032. {
  1033.   long    av,prec;
  1034.   GEN     p;
  1035.  
  1036.   prec=precision(y);
  1037.   if(!prec) err(acher3);
  1038.   av=avma;p=gach(x,prec);
  1039.   gaffect(p,y);avma=av;
  1040. }
  1041.  
  1042. /********************************************************************/
  1043. /********************************************************************/
  1044. /**                                                                **/
  1045. /**           FONCTION ARGUMENT TANGENTE HYPERBOLIQUE              **/
  1046. /**                                                                **/
  1047. /********************************************************************/
  1048. /********************************************************************/
  1049.  
  1050. GEN     mpath(x)
  1051.     
  1052.      GEN     x;
  1053.     
  1054. {
  1055.   long    av,tetpil;
  1056.   GEN     y,p1;
  1057.  
  1058.   if (typ(x)!=2) err(ather1);
  1059.   if (!signe(x))
  1060.     {
  1061.       y=cgetr(3);y[1]=x[1];y[2]=0;
  1062.     }
  1063.   else
  1064.     {
  1065.       av=avma;
  1066.       p1=addrs(divsr(2,subsr(1,x)),-1);
  1067.       tetpil=avma;y=gerepile(av,tetpil,mplog(p1));
  1068.       setexpo(y,expo(y)-1);
  1069.     }
  1070.   return y;
  1071. }
  1072.  
  1073. GEN     gath(x,prec)
  1074.     
  1075.      GEN     x;
  1076.      long    prec;
  1077.     
  1078. {
  1079.   long    av,tetpil,v;
  1080.   GEN     y,p1;
  1081.  
  1082.   switch(typ(x))
  1083.     {
  1084.     case 2 : if(expo(x)<0) y=mpath(x);
  1085.     else
  1086.       {
  1087.         av=avma;p1=addrs(divsr(2,addsr(-1,x)),1);
  1088.         tetpil=avma;y=cgetg(3,6);
  1089.         p1=mplog(p1);setexpo(p1,expo(p1)-1);
  1090.         y[1]=(long)p1;
  1091.         y[2]=lmppi(lg(x));setexpo(y[2],0);
  1092.         y=gerepile(av,tetpil,y);
  1093.       }
  1094.       break;
  1095.     
  1096.     case 6 : av=avma;
  1097.       p1=gaddgs(gdivsg(2,gsubsg(1,x)),-1);
  1098.       p1=glog(p1,prec);tetpil=avma;
  1099.       y=gerepile(av,tetpil,gmul2n(p1,-1));
  1100.       break;
  1101.     
  1102.     case 3 :
  1103.     case 7 : err(ather2);
  1104.     
  1105.     case 11: av=avma;if(valp(x)<0) err(ather4);
  1106.       v=varn(x);p1=gdiv(deriv(x,v),gsubsg(1,gmul(x,x)));
  1107.       if(valp(x)) {tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));}
  1108.       else
  1109.         {
  1110.           y=integ(p1,v);p1=gath(x[2],prec);tetpil=avma;
  1111.           y=gerepile(av,tetpil,gadd(p1,y));
  1112.         }
  1113.       break;
  1114.     
  1115.     default: y=transc(gath,x,prec);
  1116.     }
  1117.   return y;
  1118. }
  1119.  
  1120. gathz(x,y)
  1121.     
  1122.      GEN     x,y;
  1123.     
  1124. {
  1125.   long    av,prec;
  1126.   GEN     p;
  1127.  
  1128.   prec=precision(y);
  1129.   if(!prec) err(ather3);
  1130.   av=avma;p=gath(x,prec);
  1131.   gaffect(p,y);avma=av;
  1132. }
  1133.  
  1134. /********************************************************************/
  1135. /********************************************************************/
  1136. /**                                                                **/
  1137. /**             FONCTION TABLEAU DES NOMBRES DE BERNOULLI          **/
  1138. /**                                                                **/
  1139. /********************************************************************/
  1140. /********************************************************************/
  1141.  
  1142. void mpbern(nomb,prec)
  1143.     
  1144.      long    nomb,prec; /* pour calculer B_0,B_2,...,B_2*nomb */
  1145.     
  1146. {
  1147.   long    n,m,i,j,d1,d2,av;
  1148.   GEN     p1,p2;
  1149.  
  1150.   if(nomb<0) nomb=0;
  1151.   if (bernzone)
  1152.     {
  1153.       if ((bernzone[1]>=nomb)&&(bernzone[2]>=prec)) return;
  1154.       killbloc(bernzone);
  1155.     }
  1156.   bernzone=newbloc(3+prec*(nomb+1));
  1157.   bernzone[1]=nomb;
  1158.   bernzone[2]=prec;
  1159.   av=avma;
  1160.   p1=cgetr(prec+1);p2=cgetr(prec+1);
  1161.   *((GEN)bern(0))=0x2000000+prec;
  1162.   affsr(1,bern(0));
  1163.  
  1164.   for (i=1;i<=nomb;i++)
  1165.     {
  1166.       *((GEN)bern(i))=0x2000000+prec;
  1167.       affsr(0,p1);affsr(4,p2);
  1168.       n=8;m=5;d1=i-1;d2=2*i-3;
  1169.       for (j=i-1;j>0;--j)
  1170.         {
  1171.           addrrz(bern(j),p1,p1);
  1172.           mulsrz(n*m,p1,p1);
  1173.           divrsz(p1,d1*d2,p1);
  1174.           mulsrz(4,p2,p2);
  1175.           n+=4;m+=2;d1--;d2-=2;
  1176.         }
  1177.       addsrz(1,p1,p1);
  1178.       divrsz(p1,2*i+1,p1);
  1179.       subsrz(1,p1,p1);
  1180.       divrrz(p1,p2,bern(i));
  1181.     }
  1182.   avma=av;
  1183. }
  1184.  
  1185. GEN bernreal(n,prec)
  1186.      long n,prec;
  1187.  
  1188. {
  1189.   long n1; GEN p1;
  1190.   
  1191.   if(n==1) {affsr(-1,p1=cgetr(prec));setexpo(p1,-1);return p1;}
  1192.   if((n<0)||(n&1)) return gzero;
  1193.   n1=n>>1;mpbern(n1+1,prec);
  1194.   p1=cgetr(prec);affrr(bern(n1),p1);
  1195.   return p1;
  1196. }
  1197.  
  1198. GEN bernvec(nomb)
  1199.     
  1200.      long    nomb; /* pour calculer le vecteur B_0,B_2,...,B_2*nomb */
  1201.     
  1202. {
  1203.   long    n,m,i,j,d1,d2,av,tetpil;
  1204.   GEN     p1,p2,y;
  1205.  
  1206.   y=cgetg(nomb+2,17);y[1]=un;
  1207.   for (i=1;i<=nomb;i++)
  1208.     {
  1209.       av=avma;
  1210.       p1=gzero;p2=stoi(4);
  1211.       n=8;m=5;d1=i-1;d2=2*i-3;
  1212.       for (j=i-1;j>0;--j)
  1213.         {
  1214.           p1=gdivgs(gmulsg(n*m,gadd(p1,y[j+1])),d1*d2);
  1215.           p2=shifti(p2,2);
  1216.           n+=4;m+=2;d1--;d2-=2;
  1217.         }
  1218.       p1=gsubsg(1,gdivgs(gaddsg(1,p1),2*i+1));
  1219.       tetpil=avma;y[i+1]=lpile(av,tetpil,gdiv(p1,p2));
  1220.     }
  1221.   return y;
  1222. }
  1223.  
  1224. /********************************************************************/
  1225. /********************************************************************/
  1226. /**                                                                **/
  1227. /**                      FONCTION GAMMA                            **/
  1228. /**                                                                **/
  1229. /********************************************************************/
  1230. /********************************************************************/
  1231.  
  1232. GEN     mpgamma(x)
  1233.     
  1234.      GEN     x;
  1235.     
  1236. {
  1237.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av,av1;
  1238.   double  alpha,beta,dk;
  1239.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p9,p91,pitemp;
  1240.  
  1241.   if (typ(x)!=2) err(gamer1);
  1242.   l=lg(x);y=cgetr(l);s1=signe(x);if (!s1) err(gamer2);
  1243.   av=avma;p1=cgetr(l+1);u=expo(x);
  1244.   if ((u<-1) || (s1<0))
  1245.     {
  1246.       av1=avma;p2=gfrac(x);
  1247.       if(gcmp0(p2)) err(gamer2);
  1248.       avma=av1;subsrz(1,x,p1);
  1249.     }
  1250.   else affrr(x,p1);
  1251.   alpha=rtodbl(p1);beta=(16*LOG2*(l-2)/PI)-alpha;
  1252.   if (beta>=0) n=1+K2*beta;else n=0;
  1253.   if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1254.   else
  1255.     {
  1256.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1257.       if(beta>1.) beta+=(log(beta)/LOG2);
  1258.       p=16*(l-2)/beta+1;l2=l+1;
  1259.     }
  1260.   mpbern(p,l2);p91=cgetr(l2);
  1261.   p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);p5=cgetr(l2);p6=cgetr(l2);
  1262.   if(n) addsrz(n,p1,p2);else p2=p1;affrr(mplog(p2),p3);
  1263.   affsr(1,p4);setexpo(p4,expo(p4)-1);
  1264.   subrrz(p2,p4,p4);
  1265.   mulrrz(p4,p3,p4);subrrz(p4,p2,p4);
  1266.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p3);setexpo(p3,-1);
  1267.   addrrz(p4,p3,p4);mulrrz(p2,p2,p3);divsrz(1,p3,p3);e=expo(p3);
  1268.   setlg(p3,4);setlg(p5,4);setlg(p6,4);
  1269.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1270.   divrsz(p9,2*p*(2*p-1),p5);
  1271.   s=0;l1=2;l2-=2;
  1272.   for (k=p;k>1;k--)
  1273.     {
  1274.       mulrrz(p3,p5,p5);
  1275.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1276.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1277.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1278.       s &=31;
  1279.       setlg(p3,l1+2);setlg(p5,l1+2);setlg(p6,l1+2);
  1280.       addrrz(p6,p5,p5);
  1281.     }
  1282.   setlg(p5,l2+2);
  1283.   divrrz(p5,p2,p5);addrrz(p4,p5,p4);affrr(mpexp(p4),p4);
  1284.   for (i=1;i<=n;i++)
  1285.     {
  1286.       subrsz(p2,1,p2);divrrz(p4,p2,p4);
  1287.     }
  1288.   if ((u<-1)||(s1<0))
  1289.     {
  1290.       pitemp=mppi(l+1);mulrrz(pitemp,x,p1);affrr(mpsin(p1),p1);
  1291.       mulrrz(p1,p4,p4);divrrz(pitemp,p4,y);
  1292.     }
  1293.   else affrr(p4,y);
  1294.   avma=av;return y;
  1295. }
  1296.  
  1297. GEN     cxgamma(x,prec)
  1298.     
  1299.      GEN     x;
  1300.     
  1301. {
  1302.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av;
  1303.   double  alpha,beta,dk;
  1304.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p9,p91,pitemp;
  1305.  
  1306.   if (typ(x)!=6) err(gamer1);
  1307.   l=16383;if(typ(x[1])==2) l=precision(x[1]);
  1308.   if(typ(x[2])==2) {l1=precision(x[2]);if(l1<l) l=l1;}
  1309.   if(l==16383) l=prec;
  1310.   y=cgetg(3,6);y[1]=lgetr(l);y[2]=lgetr(l);
  1311.   s1=gsigne(x[1]);av=avma;
  1312.   p1=cgetg(3,6);p1[1]=lgetr(l+1);p1[2]=lgetr(l+1);
  1313.   if(s1||(typ(x[1])==2)) u=gexpo(x[1]);else u= -2;
  1314.   if ((s1<=0)||(u<-1)) gsubsgz(1,x,p1);
  1315.   else gaffect(x,p1);
  1316.   alpha=rtodbl(gabs(p1,4));beta=(16*LOG2*(l-2)/PI)-alpha;
  1317.   if (beta>=0) n=1+K2*beta;else n=0;
  1318.   if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1319.   else
  1320.     {
  1321.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1322.       if(beta>1.) beta+=(log(beta)/LOG2);
  1323.       p=16*(l-2)/beta+1;l2=l+1;
  1324.     }
  1325.   mpbern(p,l2);p91=cgetr(l2);
  1326.   p2=cgetg(3,6);p2[1]=lgetr(l2);p2[2]=lgetr(l2);
  1327.   p3=cgetg(3,6);p3[1]=lgetr(l2);p3[2]=lgetr(l2);
  1328.   p4=cgetg(3,6);p4[1]=lgetr(l2);p4[2]=lgetr(l2);
  1329.   p5=cgetg(3,6);p5[1]=lgetr(l2);p5[2]=lgetr(l2);
  1330.   p6=cgetr(l2);
  1331.   if(n) {addsrz(n,p1[1],p2[1]);affrr(p1[2],p2[2]);} else p2=p1;
  1332.   gaffect(glog(p2,l2),p3);
  1333.   affsr(1,p4[1]);setexpo(p4[1],-1);
  1334.   subrrz(p2[1],p4[1],p4[1]);p4[2]=lcopy(p2[2]);
  1335.   gmulz(p4,p3,p4);gsubz(p4,p2,p4);
  1336.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p6);setexpo(p6,-1);
  1337.   addrrz(p4[1],p6,p4[1]);gmulz(p2,p2,p3);gdivsgz(1,p3,p3);e=gexpo(p3);
  1338.   setlg(p3[1],4);setlg(p5[1],4);setlg(p6,4);setlg(p3[2],4);setlg(p5[2],4);
  1339.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1340.   gdivgsz(p9,2*p*(2*p-1),p5);
  1341.   s=0;l1=2;l2-=2;
  1342.   for (k=p;k>1;k--)
  1343.     {
  1344.       gmulz(p3,p5,p5);
  1345.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1346.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1347.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1348.       s &=31;
  1349.       setlg(p3[1],l1+2);setlg(p5[1],l1+2);setlg(p6,l1+2);
  1350.       setlg(p3[2],l1+2);setlg(p5[2],l1+2);
  1351.       addrrz(p6,p5[1],p5[1]);
  1352.     }
  1353.   setlg(p5[1],l2+2);setlg(p5[2],l2+2);
  1354.   gdivz(p5,p2,p5);gaddz(p4,p5,p4);gaffect(gexp(p4),p4);
  1355.   for (i=1;i<=n;i++)
  1356.     {
  1357.       subrsz(p2[1],1,p2[1]);gdivz(p4,p2,p4);
  1358.     }
  1359.   if ((s1<=0)||(u<-1))
  1360.     {
  1361.       pitemp=mppi(l+1);gmulz(pitemp,x,p1);gaffect(gsin(p1),p1);
  1362.       gmulz(p1,p4,p4);gdivz(pitemp,p4,y);
  1363.     }
  1364.   else gaffect(p4,y);
  1365.   avma=av;return y;
  1366. }
  1367.     
  1368. GEN     ggamma(x,prec)
  1369.     
  1370.      GEN     x;
  1371.      long    prec;
  1372.     
  1373. {
  1374.   long    i,lx;
  1375.   GEN     y;
  1376.  
  1377.     switch(typ(x))
  1378.       {
  1379.       case 1: if(signe(x)<=0) err(gamer2);
  1380.     y=transc(ggamma,x,prec);break;
  1381.       case 2 : y=mpgamma(x);break;
  1382.       case 6 : y=(gcmp0(x[2])) ? ggamma(x[1],prec) : cxgamma(x,prec);break;
  1383.       case 7 : err(impl,"p-adic gamma function");
  1384.       case 3 : err(gamer3);
  1385.       case 11: y=gexp(glngamma(x,prec));break;
  1386.       case 17:
  1387.       case 18:
  1388.       case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1389.     for(i=1;i<lx;i++)
  1390.       y[i]=lgamma(x[i],prec);
  1391.     break;
  1392.       default: y=transc(ggamma,x,prec);
  1393.       }
  1394.   return y;
  1395. }
  1396.  
  1397. ggammaz(x,y)
  1398.     
  1399.      GEN     x,y;
  1400.     
  1401. {
  1402.   long    av,prec;
  1403.   GEN     p;
  1404.  
  1405.   prec=precision(y);
  1406.   if(!prec) err(gamer4);
  1407.   av=avma;p=ggamma(x,prec);
  1408.   gaffect(p,y);avma=av;
  1409. }
  1410.  
  1411. GEN     mplngamma(x)
  1412.     
  1413.      GEN     x;
  1414.     
  1415. {
  1416.   long    l,l1,l2,u,i,k,e,f,s,s1,n,p,av,av1;
  1417.   double  alpha,beta,dk;
  1418.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p8,p9,p91,pitemp;
  1419.  
  1420.   if (typ(x)!=2) err(gamer1);
  1421.   l=lg(x);y=cgetr(l);s1=signe(x);if (!s1) err(gamer2);
  1422.   av=avma;p1=cgetr(l+1);u=expo(x);
  1423.   if ((u<-1) || (s1<0))
  1424.     {
  1425.       av1=avma;p2=gfrac(x);
  1426.       if(gcmp0(p2)) err(gamer2);
  1427.       avma=av1;subsrz(1,x,p1);
  1428.     }
  1429.   else affrr(x,p1);
  1430.   if(expo(p1)>1000) 
  1431.     {
  1432.       n=0;beta=log(K4/(l-2))/LOG2+expo(p1);beta+=(log(beta)/LOG2);
  1433.       p=16*(l-2)/beta+1;l2=l+1;
  1434.     }
  1435.   else
  1436.     {
  1437.       alpha=rtodbl(p1);beta=(16*LOG2*(l-2)/PI)-alpha;
  1438.       if (beta>=0) n=1+K2*beta;else n=0;
  1439.       if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1440.       else
  1441.     {
  1442.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1443.       if(beta>1.) beta+=(log(beta)/LOG2);
  1444.       p=16*(l-2)/beta+1;l2=l+1;
  1445.     }
  1446.     }
  1447.   mpbern(p,l2);p91=cgetr(l2);
  1448.   p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);p5=cgetr(l2);p6=cgetr(l2);
  1449.   p8=cgetr(l2);if(n) addsrz(n,p1,p2);else p2=p1;affrr(mplog(p2),p3);
  1450.   affsr(1,p4);affsr(1,p8);setexpo(p4,expo(p4)-1);
  1451.   subrrz(p2,p4,p4);mulrrz(p4,p3,p4);subrrz(p4,p2,p4);
  1452.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p3);setexpo(p3,-1);
  1453.   addrrz(p4,p3,p4);mulrrz(p2,p2,p3);divsrz(1,p3,p3);e=expo(p3);
  1454.   setlg(p3,4);setlg(p5,4);setlg(p6,4);
  1455.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1456.   divrsz(p9,2*p*(2*p-1),p5);
  1457.   s=0;l1=2;l2-=2;
  1458.   for (k=p;k>1;k--)
  1459.     {
  1460.       mulrrz(p3,p5,p5);
  1461.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1462.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1463.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1464.       s &=31;
  1465.       setlg(p3,l1+2);setlg(p5,l1+2);setlg(p6,l1+2);
  1466.       addrrz(p6,p5,p5);
  1467.     }
  1468.   setlg(p5,l2+2);
  1469.   divrrz(p5,p2,p5);addrrz(p4,p5,p4);
  1470.   for (i=1;i<=n;i++)
  1471.     {
  1472.       subrsz(p2,1,p2);mulrrz(p8,p2,p8);
  1473.     }
  1474.   f=signe(p8);subrrz(p4,mplog((f>0)?p8:negr(p8)),p4);
  1475.   if ((u<-1)||(s1<0))
  1476.     {
  1477.       pitemp=mppi(l+1);mulrrz(pitemp,x,p1);divrrz(pitemp,mpsin(p1),p1);
  1478.       f*=signe(p1);subrrz(mplog(absr(p1)),p4,y);
  1479.     }
  1480.   else affrr(p4,y);
  1481.   avma=av;if(f<0) {p2=cgetg(3,6);p2[1]=(long)y;p2[2]=(long)mppi(l);return p2;}
  1482.   else return y;
  1483. }
  1484.  
  1485. GEN     cxlngamma(x,prec)
  1486.     
  1487.      GEN     x;
  1488.     
  1489. {
  1490.   long    l,l1,l2,u,i,k,e,s,s1,n,p,av,flag;
  1491.   double  alpha,beta,dk;
  1492.   GEN     y,p1,p2,p3,p4,p5,p6,p7,p8,p9,p91,pitemp;
  1493.  
  1494.   if (typ(x)!=6) err(gamer1);
  1495.   l=16383;if(typ(x[1])==2) l=precision(x[1]);
  1496.   if(typ(x[2])==2) {l1=precision(x[2]);if(l1<l) l=l1;}
  1497.   if(l==16383) l=prec;
  1498.   y=cgetg(3,6);y[1]=lgetr(l);y[2]=lgetr(l);
  1499.   s1=gsigne(x[1]);av=avma;
  1500.   p1=cgetg(3,6);p1[1]=lgetr(l+1);p1[2]=lgetr(l+1);
  1501.   if(s1||(typ(x[1])==2)) u=gexpo(x[1]);else u= -2;
  1502.   if (((s1<=0)||(u<-1))&&(!gcmp0(x[2])&&(gexpo(x[2])<=16)))
  1503.     {gsubsgz(1,x,p1);flag=1;}
  1504.   else {gaffect(x,p1);flag=0;}
  1505.   p2=gabs(p1,4);
  1506.   if(expo(p2)>1000) 
  1507.     {
  1508.       n=0;beta=log(K4/(l-2))/LOG2+expo(p1);beta+=(log(beta)/LOG2);
  1509.       p=16*(l-2)/beta+1;l2=l+1;
  1510.     }
  1511.   else
  1512.     {
  1513.       alpha=rtodbl(p2);beta=(16*LOG2*(l-2)/PI)-alpha;
  1514.       if (beta>=0) n=1+K2*beta;else n=0;
  1515.       if(n) {p=1+PI*(alpha+n);l2=l+1+(n>>5);}
  1516.       else
  1517.     {
  1518.       dk=K4*alpha/(l-2);beta=log(dk)/LOG2;
  1519.       if(beta>1.) beta+=(log(beta)/LOG2);
  1520.       p=16*(l-2)/beta+1;l2=l+1;
  1521.     }
  1522.     }
  1523.   mpbern(p,l2);p91=cgetr(l2);
  1524.   p2=cgetg(3,6);p2[1]=lgetr(l2);p2[2]=lgetr(l2);
  1525.   p3=cgetg(3,6);p3[1]=lgetr(l2);p3[2]=lgetr(l2);
  1526.   p4=cgetg(3,6);p4[1]=lgetr(l2);p4[2]=lgetr(l2);
  1527.   p5=cgetg(3,6);p5[1]=lgetr(l2);p5[2]=lgetr(l2);
  1528.   p8=cgetg(3,6);p8[1]=lgetr(l2);p8[2]=lgetr(l2);
  1529.   gaffsg(1,p8);p6=cgetr(l2);
  1530.   if(n) {addsrz(n,p1[1],p2[1]);affrr(p1[2],p2[2]);} else p2=p1;
  1531.   gaffect(glog(p2,l2),p3);
  1532.   affsr(1,p4[1]);setexpo(p4[1],-1);
  1533.   subrrz(p2[1],p4[1],p4[1]);p4[2]=lcopy(p2[2]);
  1534.   gmulz(p4,p3,p4);gsubz(p4,p2,p4);
  1535.   p7=mppi(l2);setexpo(p7,2);affrr(mplog(p7),p6);setexpo(p6,-1);
  1536.   addrrz(p4[1],p6,p4[1]);gmulz(p2,p2,p3);gdivsgz(1,p3,p3);e=gexpo(p3);
  1537.   setlg(p3[1],4);setlg(p5[1],4);setlg(p6,4);setlg(p3[2],4);setlg(p5[2],4);
  1538.   if(bernzone[2]>l2) {affrr(bern(p),p91);p9=p91;} else p9=(GEN)bern(p);
  1539.   gdivgsz(p9,2*p*(2*p-1),p5);
  1540.   s=0;l1=2;l2-=2;
  1541.   for (k=p;k>1;k--)
  1542.     {
  1543.       gmulz(p3,p5,p5);
  1544.       if(bernzone[2]>l2) {affrr(bern(k-1),p91);p9=p91;} else p9=(GEN)bern(k-1);
  1545.       divrsz(p9,(2*k-2)*(2*k-3),p6);
  1546.       s-=e;l1+=(s>>5);if (l1>l2) l1=l2;
  1547.       s &=31;
  1548.       setlg(p3[1],l1+2);setlg(p5[1],l1+2);setlg(p6,l1+2);
  1549.       setlg(p3[2],l1+2);setlg(p5[2],l1+2);
  1550.       addrrz(p6,p5[1],p5[1]);
  1551.     }
  1552.   setlg(p5[1],l2+2);setlg(p5[2],l2+2);
  1553.   gdivz(p5,p2,p5);gaddz(p4,p5,p4);
  1554.   for (i=1;i<=n;i++)
  1555.     {
  1556.       subrsz(p2[1],1,p2[1]);gmulz(p8,p2,p8);
  1557.     }
  1558.   gsubz(p4,glog(p8),p4);
  1559.   if (flag)
  1560.     {
  1561.       pitemp=mppi(l+1);gmulz(pitemp,x,p1);gdivz(pitemp,gsin(p1),p1);
  1562.       gsubz(glog(p1),p4,y);
  1563.     }
  1564.   else gaffect(p4,y);
  1565.   avma=av;return y;
  1566. }
  1567.     
  1568. GEN     glngamma(x,prec)
  1569.     
  1570.      GEN     x;
  1571.      long    prec;
  1572.     
  1573. {
  1574.   long    i,lx,av,tetpil,v,n;
  1575.   GEN     y,p1;
  1576.  
  1577.     switch(typ(x))
  1578.       {
  1579.       case 1: if(signe(x)<=0) err(gamer2);
  1580.     y=transc(glngamma,x,prec);break;
  1581.       case 2 : y=mplngamma(x);break;
  1582.       case 6 : y=(gcmp0(x[2])) ? glngamma(x[1],prec) : cxlngamma(x,prec);break;
  1583.       case 7 : err(impl,"p-adic lngamma function");
  1584.       case 3 : err(gamer3);
  1585.       case 11: av=avma;if(valp(x)) err(loger5);
  1586.     v=varn(x);if(!gcmp1(x[2])) err(impl,"lngamma around a!=1");
  1587.     p1=gsubsg(1,x);n=(lg(x)-3)/valp(p1);
  1588.     y=ggrando(polx[v],lg(x)-2);
  1589.     for(i=n;i>=2;i--)
  1590.       {
  1591.         y=gmul(p1,gadd(gdivgs(izeta(stoi(i),prec),i),y));
  1592.       }
  1593.     y=gadd(mpeuler(prec),y);tetpil=avma;
  1594.     y=gerepile(av,tetpil,gmul(p1,y));break;
  1595.       case 17:
  1596.       case 18:
  1597.       case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1598.     for(i=1;i<lx;i++)
  1599.       y[i]=(long)glngamma(x[i],prec);
  1600.     break;
  1601.       default: y=transc(glngamma,x,prec);
  1602.       }
  1603.   return y;
  1604. }
  1605.  
  1606. glngammaz(x,y)
  1607.     
  1608.      GEN     x,y;
  1609.     
  1610. {
  1611.   long    av,prec;
  1612.   GEN     p;
  1613.  
  1614.   prec=precision(y);
  1615.   if(!prec) err(gamer4);
  1616.   av=avma;p=glngamma(x,prec);
  1617.   gaffect(p,y);avma=av;
  1618. }
  1619.  
  1620. /********************************************************************/
  1621. /********************************************************************/
  1622. /**                                                                **/
  1623. /**             FONCTION GAMMA DES DEMI-ENTIERS                    **/
  1624. /**                                                                **/
  1625. /********************************************************************/
  1626. /********************************************************************/
  1627.  
  1628. GEN     mpgamd(x,prec)
  1629.     
  1630.      long    x,prec;
  1631.     
  1632. {
  1633.   long    i,j,a,l,av;
  1634.   GEN     y,p1;
  1635.  
  1636.   a=abs(x);
  1637.   l=prec+1+a/32;
  1638.   if (l>32767) err(gamder1);
  1639.   y=cgetr(prec);
  1640.   av=avma;
  1641.   p1=mpsqrt(mppi(l));
  1642.   if (x>=0)
  1643.     {
  1644.       j= -1;
  1645.       for (i=0;i<x;i++)
  1646.         {
  1647.           j+=2;mulsrz(j,p1,p1);
  1648.           setexpo(p1,expo(p1)-1);
  1649.         }
  1650.     }
  1651.   else
  1652.     {
  1653.       j=1;
  1654.       for (i=0;i<a;i++)
  1655.         {
  1656.           j-=2;divrsz(p1,j,p1);
  1657.           setexpo(p1,expo(p1)+1);
  1658.         }
  1659.     }
  1660.   affrr(p1,y);
  1661.   avma=av;
  1662.   return y;
  1663. }
  1664.  
  1665. mpgamdz(s,y)
  1666.     
  1667.      long    s;
  1668.      GEN     y;
  1669.     
  1670. {
  1671.   long    l,av;
  1672.   GEN     p1;
  1673.  
  1674.   av=avma;
  1675.   l=lg(y);
  1676.   p1=mpgamd(s,l);
  1677.   affrr(p1,y);
  1678.   avma=av;
  1679. }
  1680.  
  1681. GEN     ggamd(x,prec)
  1682.     
  1683.      GEN     x;
  1684.      long    prec;
  1685.     
  1686. {
  1687.   long    av,tetpil,i,lx;
  1688.   GEN     y,p1;
  1689.  
  1690.   switch(typ(x))
  1691.     {
  1692.     case 1 : y=mpgamd(itos(x),prec);break;
  1693.     case 2 :
  1694.     case 4 :
  1695.     case 5 :
  1696.     case 6 :
  1697.     case 8 : av=avma;p1=gadd(x,ghalf);tetpil=avma;
  1698.       y=gerepile(av,tetpil,ggamma(p1,prec));
  1699.       break;
  1700.     case 17:
  1701.     case 18:
  1702.     case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1703.       for(i=1;i<lx;i++)
  1704.         y[i]=lgamd(x[i],prec);
  1705.       break;
  1706.     case 3 :
  1707.     case 7 : err(gamder2);
  1708.     case 11: err(impl,"gamd of a power series");
  1709.     default: y=transc(ggamd,x,prec);
  1710.     }
  1711.   return y;
  1712. }
  1713.  
  1714. ggamdz(x,y)
  1715.     
  1716.      GEN     x,y;
  1717.     
  1718. {
  1719.   long    av,prec;
  1720.   GEN     p;
  1721.  
  1722.   prec=precision(y);
  1723.   if(!prec) err(gamder3);
  1724.   av=avma;p=ggamd(x,prec);
  1725.   gaffect(p,y);avma=av;
  1726. }
  1727.  
  1728.  
  1729. /********************************************************************/
  1730. /********************************************************************/
  1731. /**                                                                **/
  1732. /**                      FONCTION PSI                              **/
  1733. /**                                                                **/
  1734. /********************************************************************/
  1735. /********************************************************************/
  1736.  
  1737.  
  1738. GEN     mppsi(z)        /*   version p=2 */
  1739.     
  1740.      GEN        z;
  1741.     
  1742. {
  1743.   long  l,n,k,x,xx,av1,av2,tetpil;
  1744.   GEN   zk,u,v,a,b;
  1745.  
  1746.   av1=avma;l=lg(z);
  1747.   x=1+16*(l-2)*LOG2+1.58*rtodbl(z);xx=x*x;
  1748.   n=1+3.591*x;
  1749.   affsr(x,a=cgetr(l));
  1750.   a=mplog(a);
  1751.   gaffect(a,u=cgetr(l));
  1752.   gaffsg(1,b=cgetr(l));
  1753.   gaffsg(1,v=cgetr(l));
  1754.   for (k=1;k<=n;k++)
  1755.     {
  1756.       av2=avma;
  1757.       zk=(k>1) ? gaddsg(k-1,z) : z;
  1758.       gdivz(gmulsg(xx,b),gmul(zk,zk),b);
  1759.       gdivz(gsub(gdiv(gmulsg(xx,a),zk),b),zk,a);
  1760.       gaddz(u,a,u);gaddz(v,b,v);
  1761.       avma=av2;
  1762.     }
  1763.   tetpil=avma;return gerepile(av1,tetpil,gdiv(u,v));
  1764. }
  1765.  
  1766. GEN     cxpsi(z,prec)        /*   version p=2 */
  1767.     
  1768.      GEN        z;
  1769.     
  1770. {
  1771.   long  l,l1,n,k,x,xx,av1,av2,tetpil;
  1772.   GEN   zk,u,v,a,b;
  1773.  
  1774.   l=16383;if(typ(z[1])==2) l=precision(z[1]);
  1775.   if(typ(z[2])==2) {l1=precision(z[2]);if(l1<l) l=l1;}
  1776.   if(l==16383) l=prec;
  1777.   av1=avma;x=1+16*(l-2)*LOG2+1.58*rtodbl(gabs(z,4));xx=x*x;
  1778.   n=1+3.591*x;
  1779.   a=cgetg(3,6);a[1]=lgetr(l);a[2]=lgetr(l);gaffsg(x,a);
  1780.   b=cgetg(3,6);b[1]=lgetr(l);b[2]=lgetr(l);gaffsg(1,b);
  1781.   u=cgetg(3,6);u[1]=lgetr(l);u[2]=lgetr(l);
  1782.   v=cgetg(3,6);v[1]=lgetr(l);v[2]=lgetr(l);gaffsg(1,v);
  1783.   a=glog(a);gaffect(a,u);
  1784.   for (k=1;k<=n;k++)
  1785.     {
  1786.       av2=avma;
  1787.       zk=(k>1) ? gaddsg(k-1,z) : z;
  1788.       gdivz(gmulsg(xx,b),gmul(zk,zk),b);
  1789.       gdivz(gsub(gdiv(gmulsg(xx,a),zk),b),zk,a);
  1790.       gaddz(u,a,u);gaddz(v,b,v);
  1791.       avma=av2;
  1792.     }
  1793.   tetpil=avma;return gerepile(av1,tetpil,gdiv(u,v));
  1794. }
  1795.  
  1796.  
  1797.  
  1798. GEN     gpsi(x,prec)
  1799.     
  1800.      GEN     x;
  1801.      long    prec;
  1802.     
  1803. {
  1804.   long    i,lx;
  1805.   GEN     y;
  1806.  
  1807.   switch(typ(x))
  1808.     {
  1809.     case 2 : y=mppsi(x);break;
  1810.     case 6 : y=cxpsi(x,prec);break;
  1811.     case 3 :
  1812.     case 7 : err(psier1);
  1813.     case 11: err(impl,"psi of power series");
  1814.     case 17:
  1815.     case 18:
  1816.     case 19: lx=lg(x);y=cgetg(lx,typ(x));
  1817.       for(i=1;i<lx;i++)
  1818.         y[i]=lpsi(x[i],prec);
  1819.       break;
  1820.     default: y=transc(gpsi,x,prec);
  1821.     }
  1822.   return y;
  1823. }
  1824.  
  1825. gpsiz(x,y)
  1826.     
  1827.      GEN     x,y;
  1828.     
  1829. {
  1830.   long    av,prec;
  1831.   GEN     p;
  1832.  
  1833.   prec=precision(y);
  1834.   if(!prec) err(psier2);
  1835.   av=avma;p=gpsi(x,prec);
  1836.   gaffect(p,y);avma=av;
  1837. }
  1838.