home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / trans1.c < prev    next >
Text File  |  1992-05-20  |  33KB  |  1,419 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                +++++++++++++++++++++++++++++++                 **/
  5. /**                +                             +                 **/
  6. /**                +  FONCTIONS TRANSCENDANTES   +                 **/
  7. /**                +                             +                 **/
  8. /**                +     copyright Babe Cool     +                 **/
  9. /**                +                             +                 **/
  10. /**                +++++++++++++++++++++++++++++++                 **/
  11. /**                                                                **/
  12. /**                                                                **/
  13. /********************************************************************/
  14. /********************************************************************/
  15.  
  16. # include "genpari.h"
  17.  
  18. /********************************************************************/
  19. /********************************************************************/
  20. /**                                                                **/
  21. /**                        FONCTION PI                             **/
  22. /**                                                                **/
  23. /********************************************************************/
  24. /********************************************************************/
  25.  
  26. void constpi(prec)
  27.   
  28.    long  prec;
  29.   
  30. {
  31.   long    l,n,n1,av1,av2;
  32.   double  alpha;
  33.   GEN p1,p2,p3, tokill;
  34.  
  35. #define k1     545140134
  36. #define k2    13591409
  37. #define k3      640320
  38. #define alpha2 1.4722004  /*   3*log(k3/12)/(32*log(2))   */
  39.  
  40.   if ((tokill = gpi) && (lg(gpi)>=prec)) return;
  41.  
  42.   gpi=newbloc(prec); *gpi = 0x2010000+prec;
  43.   av1=avma;
  44.  
  45.   n=1+(prec-2)/alpha2;
  46.   n1=6*n-1;
  47.   p1=cgetr(prec);
  48.   p2=addsi(k2,mulss(n,k1));
  49.   affir(p2,p1);
  50.   /*initialisation longueur mantisse*/
  51.   if (prec>=4) l=4;else l=prec;alpha=l;
  52.   setlg(p1,l);
  53.  
  54.   while (n)
  55.   {
  56.     av2=avma;
  57.     if(n>1290)
  58.     {
  59.       if(n1>46340)
  60.       p3=divrs(divrs(mulsr(n1-4,mulsr(n1,mulsr(n1-2,p1))),n*n),n);
  61.       else
  62.       p3=divrs(divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n),n);
  63.     }
  64.     else p3=divrs(mulsr(n1-4,mulsr(n1*(n1-2),p1)),n*n*n);
  65.     p3=divrs(divrs(p3,100100025),327843840);
  66.     subisz(p2,k1,p2);subirz(p2,p3,p1);
  67.     avma=av2;
  68.     alpha+=alpha2;l=1+alpha;/*nouvelle longueur mantisse*/
  69.     if (l>prec) l=prec;
  70.     setlg(p1,l);
  71.     n--; n1-=6;
  72.   }
  73.  
  74.   p1=divsr(53360,p1);
  75.   mulrrz(p1,gsqrt(stoi(k3),prec),gpi);
  76.   avma=av1;
  77.   killbloc(tokill);
  78. }
  79.  
  80. GEN mppi(prec)
  81.   
  82.    long  prec;
  83.   
  84. {
  85.   GEN x;
  86.  
  87.   constpi(prec);
  88.   affrr(gpi,x=cgetr(prec));
  89.   return x;
  90. }
  91.  
  92. /********************************************************************/
  93. /********************************************************************/
  94. /**                                                                **/
  95. /**                      FONCTION EULER                            **/
  96. /**                                                                **/
  97. /********************************************************************/
  98. /********************************************************************/
  99.  
  100. void consteuler(prec)
  101.   
  102.    long  prec;
  103.   
  104. {
  105.   long    l,n,k,x,xx,av1,av2;
  106.   GEN u,v,a,b,tokill;
  107.  
  108.   if ((tokill = geuler) && (lg(geuler)>=prec)) return;
  109.   l=prec+1;
  110.   geuler=newbloc(prec); *geuler = 0x2010000 + prec;
  111.   x=1+8*(l-2)*LOG2;xx=x*x;
  112.   n=1+3.591*x;     /*    a1=3.591  :   a1*[ ln(a1)-1 ]=1 */
  113.   av1=avma;
  114.   affsr(x,a=cgetr(l));
  115.   u=mplog(a);setsigne(u,-1);affrr(u,a);
  116.   affsr(1,b=cgetr(l));
  117.   affsr(1,v=cgetr(l));
  118.   for (k=1;k<=n;k++)
  119.   {
  120.     av2=avma;
  121.     divrsz(mulsr(xx,b),k*k,b);
  122.     divrsz(addrr(divrs(mulsr(xx,a),k),b),k,a);
  123.     addrrz(u,a,u);addrrz(v,b,v);
  124.     avma=av2;
  125.   }
  126.   divrrz(u,v,geuler);
  127.   avma=av1;
  128.   killbloc(tokill);
  129. }
  130.  
  131. GEN mpeuler(prec)
  132.   
  133.    long  prec;
  134.   
  135. {
  136.   GEN x;
  137.  
  138.   consteuler(prec);
  139.   affrr(geuler,x=cgetr(prec));
  140.   return x;
  141. }
  142.  
  143. /********************************************************************/
  144. /********************************************************************/
  145. /**                                                                **/
  146. /**                     CONVERSION DE TYPES                        **/
  147. /**                 POUR FONCTIONS TRANSCENDANTES                  **/
  148. /**                                                                **/
  149. /********************************************************************/
  150. /********************************************************************/
  151.  
  152. GEN transc(f,x,prec)
  153.   
  154.    GEN    x,(*f)();
  155.    long   prec;
  156.   
  157. {
  158.   GEN p1,p2,y;
  159.   long  lx,i,av,tetpil;
  160.  
  161.   av=avma;
  162.   switch(typ(x))
  163.   {
  164.   case 1 :
  165.   case 4 :
  166.   case 5 : p1=cgetr(prec);gaffect(x,p1);tetpil=avma;
  167.     y=gerepile(av,tetpil,(*f)(p1,prec));break;
  168.   
  169.   case 6 :
  170.   case 8 : av=avma;p1=cgetr(prec);affsr(1,p1);
  171.     p1=gmul(x,p1);tetpil=avma;
  172.     y=gerepile(av,tetpil,(*f)(p1,prec));break;
  173.   
  174.   case 10:
  175.   case 13:
  176.   case 14: p1=tayl(x,gvar(x),precdl);tetpil=avma;
  177.     y=gerepile(av,tetpil,(*f)(p1,prec));
  178.     break;
  179.   
  180.   case 17:
  181.   case 18:
  182.   case 19: lx=lg(x);y=cgetg(lx,typ(x));
  183.     for(i=1;i<lx;i++)
  184.     y[i]=(long)(*f)(x[i],prec);break;
  185.   case 9 : av=avma;p1=roots(x[1],prec);lx=lg(p1);p2=cgetg(lx,18);
  186.     for(i=1;i<lx;i++) p2[i]=lpoleval(x[2],p1[i]);
  187.     tetpil=avma;y=cgetg(lx,18);
  188.     for(i=1;i<lx;i++) y[i]=(long)(*f)(p2[i],prec);
  189.     y=gerepile(av,tetpil,y);break;
  190.   case 15:
  191.   case 16: err(transcer1);
  192.   default: y=(*f)(x,prec);
  193.   }
  194.   return y;
  195. }
  196.  
  197.  
  198. /********************************************************************/
  199. /********************************************************************/
  200. /**                                                                **/
  201. /**                        RACINE CARREE                           **/
  202. /**                                                                **/
  203. /********************************************************************/
  204. /********************************************************************/
  205.  
  206. GEN mpsqrt(x)
  207.   
  208.    GEN   x;
  209.   
  210. {
  211.   long    l,l1,l0,l2,s,eps,n,i,alpha,ex,av;
  212.   double  beta;
  213.   GEN y,p1,p2,p3;
  214.  
  215.   if (typ(x)!=2) err(sqrter1);
  216.   s=signe(x);
  217.   if (s<0) err(sqrter2);
  218.   if (!s)
  219.   {
  220.     y=cgetr(3);y[1]=0x800000+(expo(x)>>1);
  221.     y[2]=0;
  222.   }
  223.   else
  224.   {
  225.     l=lg(x);y=cgetr(l);
  226.     av=avma;
  227.     p1=cgetr(l+1);
  228.     mpaff(x,p1);
  229.     ex=expo(x);eps=ex % 2;ex /=2;
  230.     if (eps<0)
  231.     {
  232.       eps+=2;ex-=1;
  233.     }
  234.     setexpo(p1,eps);setlg(p1,3);
  235.     l1=1;
  236.     beta=sqrt((eps+1)*(2+mant(p1,1)/C31));
  237.     n=2+log((double) (l-2))/LOG2;
  238.     p2=cgetr(l+1);
  239.     setexpo(p2,0);setlg(p2,3);setsigne(p2,1);
  240.     l2=3;
  241.     alpha=(beta-2)*C31;
  242.     setmant(p2,1,alpha);
  243.     p3=cgetr(l+1);l-=2;
  244.     if (!alpha)
  245.     {
  246.       setmant(p2,1,0x80000000);
  247.       setexpo(p2,1);
  248.     }
  249.   
  250.     for (i=1;i<=n;i++)
  251.     {
  252.       l0=l1+l1;
  253.       if (l0<=l)
  254.       {
  255.         l2=l2+l1;
  256.         setlg(p2,l2);
  257.         setlg(p3,l0+2);
  258.         setlg(p1,l0+2);
  259.         l1=l0;
  260.       }
  261.       else
  262.       {
  263.         l2=l2-l1+l+1;
  264.         setlg(p2,l2);
  265.         setlg(p3,l+3);setlg(p1,l+3);
  266.         l1=l+1;
  267.       }
  268.       divrrz(p1,p2,p3);
  269.       addrrz(p2,p3,p2);
  270.       setexpo(p2,expo(p2)-1);
  271.     }
  272.     affrr(p2,y);
  273.     setexpo(y,expo(y)+ex);
  274.     avma=av;
  275.   }
  276.   return y;
  277. }
  278.  
  279. GEN   pasqrt(x)
  280.   GEN x;
  281. {
  282.   long    av,av2,l3,e,lp,lpp,pp,fl,re;
  283.   GEN y,p1,p2;
  284.  
  285.   e=valp(x);
  286.   y=cgetg(5,7);y[2]=x[2];
  287.   if(gcmp0(x))
  288.   {
  289.     y[4]=zero;e=(e+1)>>1;
  290.     y[3]=lpuigs(x[2],e);
  291.     setvalp(y,e);setprecp(y,precp(x));
  292.   }
  293.   else
  294.   {
  295.     if(e&1) err(sqrter6);l3=lgef(x[3]);
  296.     e>>=1;fl=cmpis(x[2],2);pp=precp(x);
  297.     y[4]=lgeti(l3);av=avma;
  298.     if(fl)
  299.     {
  300.       if(!mpsqrtmod(x[4],x[2],&p1)) err(sqrter5);
  301.       lp=1;affii(p1,y[4]);avma=av;
  302.     }
  303.     else
  304.     {
  305.       affsi(1,y[4]);
  306.       re=(((GEN)x[4])[lgef(x[4])-1])&7;
  307.       if((re!=1)&&(pp>=2))
  308.       {if((pp!=2) || (re!=5)) err(sqrter5);}
  309.       lp=3;
  310.     }
  311.     if(pp<=lp)
  312.     {
  313.       setprecp(y,1);y[3]=x[2];
  314.     }
  315.     else
  316.     {
  317.       y[3]=lgeti(l3);av=avma;
  318.       p2=gcopy(x);setvalp(p2,0);av2=avma;
  319.       setvalp(y,0);if(!fl) affsi(8,y[3]);else affii(x[2],y[3]);
  320.       while(lp<pp)
  321.       {
  322.     if(fl)
  323.       {
  324.         lp<<=1;if(lp>=pp) {lp=pp;affii(x[3],y[3]);}
  325.         else muliiz(y[3],y[3],y[3]);
  326.       }
  327.     else
  328.       {
  329.         lpp=lp;lp=(lp<<1)-1;if(lp>=pp) {lp=pp;affii(x[3],y[3]);}
  330.         else mpshiftz(y[3],lpp-1,y[3]);
  331.       }
  332.         setprecp(y,lp);p1=gadd(y,gdiv(p2,y));
  333.         gdivz(p1,gdeux,y);
  334.     if(!fl) 
  335.       {
  336.         setprecp(y,lp-1);if(lp<pp) lp--;
  337.         mpshiftz(y[3],-1,y[3]);
  338.         p1=subii(y[4],y[3]);if(signe(p1)>=0) affii(p1,y[4]);
  339.       }
  340.     avma=av2;
  341.       }
  342.       avma=av;
  343.     }
  344.     setvalp(y,e);
  345.   }
  346.   return y;
  347. }
  348.  
  349. GEN gsqrt(x,prec)
  350.   
  351.    GEN   x;
  352.    long  prec;
  353.   
  354. {
  355.   long    av,tetpil,e;
  356.   GEN y,p1,p2;
  357.  
  358.   switch(typ(x))
  359.   {
  360.   case 2 : if (signe(x)>=0) y =mpsqrt(x);
  361.   else
  362.     {
  363.     y=cgetg(3,6);y[1]=zero;
  364.     av=avma;p1=mpneg(x);tetpil=avma;
  365.     y[2]=lpile(av,tetpil,mpsqrt(p1));
  366.     }
  367.     break;
  368.   
  369.   case 3 : y=cgetg(3,3);y[1]=copyifstack(x[1]);
  370.     if(!mpsqrtmod(x[2],y[1],y+2)) err(sqrter5);
  371.     break;
  372.   
  373.   case 6 : y=cgetg(3,6);av=avma;
  374.     if(gcmp0(x[2]))
  375.     {
  376.       if(signe(x[1])>=0)
  377.       {
  378.         y[1]=lsqrt(x[1],prec);
  379.         y[2]=zero;
  380.       }
  381.       else
  382.       {
  383.         y[1]=zero;p1=mpneg(x[1]);
  384.         tetpil=avma;
  385.         y[2]=lpile(av,tetpil,gsqrt(p1,prec));
  386.       }
  387.     }
  388.     else
  389.     {
  390.       p1=gmul(x[1],x[1]);
  391.       p2=gmul(x[2],x[2]);
  392.       p1=gsqrt(gadd(p1,p2),prec);
  393.       if (gcmp0(p1))
  394.       {
  395.         y[1]=lsqrt(p1);
  396.         y[2]=lcopy(y[1]);
  397.       }
  398.       else
  399.       {
  400.         if(signe(x[1])<0)
  401.         {
  402.           p2=gsub(p1,x[1]);p1=gmul2n(p2,-1);
  403.           tetpil=avma;y[2]=lpile(av,tetpil,gsqrt(p1,prec));
  404.           av=avma;p1=gmul2n(y[2],1);tetpil=avma;
  405.           y[1]=lpile(av,tetpil,gdiv(x[2],p1));
  406.         }
  407.         else
  408.         {
  409.           p1=gadd(p1,x[1]);p1=gmul2n(p1,-1);
  410.           tetpil=avma;y[1]=lpile(av,tetpil,gsqrt(p1,prec));
  411.           av=avma;p1=gmul2n(y[1],1);tetpil=avma;
  412.           y[2]=lpile(av,tetpil,gdiv(x[2],p1));
  413.         }
  414.       }
  415.     }
  416.     break;
  417.   
  418.   case 7 : y = pasqrt(x); break;
  419.   
  420.   case 11: e=valp(x);
  421.     if(gcmp0(x))
  422.     {y=cgetg(3,11);y[1]=0x8000+((e+1)>>1);}
  423.     else
  424.     {
  425.       av=avma;
  426.       if(e&1) err(sqrter6);
  427.       e>>=1;p2=gcopy(x);setvalp(p2,0);
  428.       tetpil=avma;
  429.       y=gerepile(av,tetpil,gpui(p2,ghalf,prec));
  430.       setvalp(y,e);
  431.     }
  432.     setvarn(y,varn(x));
  433.     break;
  434.   
  435.   default: y=transc(gsqrt,x,prec);
  436.   }
  437.   return y;
  438. }
  439.  
  440. void gsqrtz(x,y)
  441.   
  442.    GEN   x,y;
  443.   
  444. {
  445.   long    av,prec;
  446.   GEN p;
  447.  
  448.   prec=precision(y);
  449.   if(!prec) err(sqrter4);
  450.   av=avma;p=gsqrt(x,prec);
  451.   gaffect(p,y);avma=av;
  452. }
  453.  
  454. /********************************************************************/
  455. /********************************************************************/
  456. /**                                                                **/
  457. /**                    FONCTION EXPONENTIELLE-1                    **/
  458. /**                                                                **/
  459. /********************************************************************/
  460. /********************************************************************/
  461.  
  462. GEN mpexp1(x)
  463.    GEN   x;
  464.   
  465. {
  466.   long    l,l0,l1,l2,i,n,m,ex,s;
  467.   long    sgn,av;
  468.   double  alpha,beta,gama=2.0 /*optimise pour SUN3*/,aa,bb;
  469.   GEN y,p1,p2,p3,p4;
  470.  
  471.   if (typ(x)!=2) err(exp1er);
  472.   else
  473.   {
  474.     sgn=signe(x);
  475.     if (!sgn)
  476.     {
  477.       y=cgetr(3);y[1]=x[1];
  478.       y[2]=0;
  479.     }
  480.     else
  481.     {
  482.       l=lg(x);
  483.       y=cgetr(l);     /*reservation pour resultat*/
  484.       av=avma;
  485.       p1=cgetr(l+1);
  486.       mpabsz(x,p1);   /*p1 recoit |x|  */
  487.       ex=expo(x);
  488.       alpha= -1-log(2+mant(x,1)/C31)-ex*LOG2;
  489.       beta=5+32*(l-2)*LOG2;
  490.       aa=sqrt((beta)/(gama*LOG2));
  491.       bb=(alpha+0.5*log(beta*gama/LOG2))/LOG2;
  492.       if (aa>=bb)
  493.       {
  494.         n=1+sqrt(beta*gama/LOG2);
  495.         m=1+aa-bb;
  496.         setexpo(p1,ex-m);
  497.       }
  498.       else
  499.       {
  500.         n=1+beta/alpha;
  501.         m=0;
  502.       }
  503.       l2=l+1+m/32;
  504.       p2=cgetr(l2);
  505.       affsr(1,p2);setlg(p2,4);
  506.       s=0;
  507.       p3=cgetr(l2);setlg(p3,4);
  508.       p4=cgetr(l2);affrr(p1,p4);setlg(p4,4);
  509.       l1=2;l2-=2;
  510.     
  511.       for (i=n;i>=2;--i)
  512.       {
  513.         divrsz(p4,i,p3);
  514.         mulrrz(p3,p2,p2);
  515.         ex=expo(p3);s=s-ex;l0=s/32;l1+=l0;
  516.         if (l1>l2) l1=l2;
  517.         s %=32;
  518.         setlg(p2,l1+2);setlg(p3,l1+2);
  519.         setlg(p4,l1+2);
  520.         addsrz(1,p2,p2);
  521.       }
  522.       l2+=2;
  523.       setlg(p2,l2);setlg(p3,l2);setlg(p4,l2);
  524.       mulrrz(p4,p2,p2);
  525.     
  526.       for (i=1;i<=m;i++)
  527.       {
  528.         addsrz(2,p2,p3);
  529.         mulrrz(p3,p2,p2);
  530.       }
  531.     
  532.       if (sgn== -1)
  533.       {
  534.         addsrz(1,p2,p2);divsrz(1,p2,p2);
  535.         subrsz(p2,1,y);
  536.       }
  537.       else
  538.       affrr(p2,y);
  539.       avma=av;
  540.     }
  541.   }
  542.   return y;
  543. }
  544.  
  545. /********************************************************************/
  546. /********************************************************************/
  547. /**                                                                **/
  548. /**                   FONCTION EXPONENTIELLE                       **/
  549. /**                                                                **/
  550. /********************************************************************/
  551. /********************************************************************/
  552.  
  553. GEN mpexp(x)
  554.    GEN   x;
  555.   
  556. {
  557.   GEN y;
  558.   long av,tetpil;
  559.  
  560.   if(gcmp0(x)) return addsr(1,x);
  561.   if(signe(x)>=0) return addsr(1,mpexp1(x));
  562.   av=avma;y=addsr(1,mpexp1(negr(x)));tetpil=avma;
  563.   return gerepile(av,tetpil,divsr(1,y));
  564. }
  565.  
  566. GEN paexp(x)
  567.   GEN x;
  568.  
  569. {
  570.   GEN y, r, p1;
  571.   long k, av = avma, tetpil;
  572.   long e = valp(x), pp = precp(x), n = e + pp;
  573.   
  574.   if (!signe(x[4])) return gaddgs(x, 1);
  575.   if ((e <= 0) || (!cmpis(x[2], 2) && (e == 1))) err(paexper1);
  576.   av = avma;
  577.   if (cmpis(x[2], 2))
  578.   {
  579.     p1 = subis(x[2], 1);
  580.     k = itos(dvmdii(subis(mulis(p1, n), 1), subis(mulis(p1, e), 1), &r));
  581.     if (!signe(r)) k--;
  582.   }
  583.   else
  584.   {
  585.     k = (n - 1) / (e - 1); if (!((n - 1)%(e-1))) k--;
  586.   }
  587.   for(y = gun; k; k--)
  588.   {
  589.     tetpil = avma; y = gaddsg(1, gdivgs(gmul(y, x), k));
  590.   }
  591.   return gerepile(av, tetpil, y);
  592. }
  593.  
  594. GEN gexp(x,prec)
  595.   
  596.    GEN x;
  597.    long  prec;
  598.   
  599. {
  600.   long    av,tetpil,i,ly,j,ex;
  601.   GEN r,u,v;
  602.   GEN y,p1,p2;
  603.  
  604.   switch(typ(x))
  605.   {
  606.     case 2 : y=mpexp(x); break;
  607.     
  608.     case 6 : y=cgetg(3,6);av=avma;
  609.       r=gexp(x[1],prec);
  610.       gsincos(x[2],&u,&v,prec);
  611.       tetpil=avma;
  612.       y[1]=lmul(r,v);y[2]=lmul(r,u);
  613.       gerepile(av,tetpil,1);
  614.       break;
  615.     
  616.     case 7 : y = paexp(x); break;
  617.     
  618.     case 11: if(gcmp0(x)) y=gaddsg(1,x);
  619.       else
  620.       {
  621.         ex=valp(x);if(ex<0) err(exper3);
  622.         if(ex)
  623.         {
  624.           ly=lg(x)+ex;y=cgetg(ly,11); 
  625.           y[1]=0x1008000;y[2]=un;setvarn(y,varn(x));
  626.           for(i=3;i<ex+2;i++) y[i]=zero;
  627.           for(i=ex+2;i<ly;i++)
  628.           {
  629.             av=avma;p1=gzero;
  630.             for(j=ex;j<i-1;j++)
  631.               p1=gadd(p1,gmulgs(gmul(x[j-ex+2],y[i-j]),j));
  632.             tetpil=avma;y[i]=lpile(av,tetpil,gdivgs(p1,i-2));
  633.           }
  634.         }
  635.         else
  636.         {
  637.           av=avma;p1=gcopy(x);p1[2]=zero;
  638.           normalize(&p1);p2=gexp(p1,prec);
  639.           p1=gexp(x[2],prec);tetpil=avma;
  640.           y=gerepile(av,tetpil,gmul(p1,p2));
  641.         }
  642.       }
  643.     break;
  644.     
  645.     case 3 : err(exper1);
  646.     
  647.     default: y = transc(gexp,x,prec);
  648.   }
  649.   return y;
  650. }
  651.  
  652. void gexpz(x,y)
  653.   
  654.    GEN   x,y;
  655.   
  656. {
  657.   long    av,prec;
  658.   GEN p;
  659.  
  660.   prec=precision(y);
  661.   if(!prec) err(exper2);
  662.   av=avma;p=gexp(x,prec);
  663.   gaffect(p,y);avma=av;
  664. }
  665.  
  666. /********************************************************************/
  667. /********************************************************************/
  668. /**                                                                **/
  669. /**                      FONCTION LOGARITHME                       **/
  670. /**                                                                **/
  671. /********************************************************************/
  672. /********************************************************************/
  673.  
  674. GEN mplog(x)
  675.   
  676.    GEN   x;
  677.   
  678. {
  679.   long    l,l0,l1,l2 ,m,m1,n,i,ex,s;
  680.   long    sgn,avmacourant,av;
  681.   double  alpha,beta,aa,bb,p1copie;
  682.   GEN y,p1,p2,p3,p4,p5;
  683.  
  684.   if (typ(x)!=2) err(loger1);
  685.   if (signe(x)<=0) err(loger2);
  686.   l=lg(x);
  687.   y=cgetr(l);
  688.   avmacourant=avma;
  689.   p1=cgetr(l+1);
  690.   affrr(x,p1);
  691.   sgn=cmpsr (1,p1);
  692.   if (!sgn) affsr(0,y);
  693.   else
  694.   {
  695.     if (sgn>0) divsrz(1,p1,p1);/* x<1 changer x en 1/x*/
  696.     m1=0;
  697.     while (expo(p1)>=1)
  698.     {
  699.       av=avma;affrr(mpsqrt(p1),p1);
  700.       avma=av;m1+=1;
  701.     }
  702.     p1copie=2+mant(p1,1)/C31;
  703.     if (p1copie==1) p1copie=p1copie+0.00000001;
  704.     l-=2;
  705.     alpha= -log(p1copie-1);
  706.     beta=  16*l*LOG2;
  707.     aa=alpha/LOG2;
  708.     bb=4*sqrt(l/3.0);
  709.     if (aa<=bb)
  710.     {
  711.       n=1+4*sqrt(3.0*l);
  712.       m=1+bb-aa;
  713.     }
  714.     else
  715.     {
  716.       n=1+beta/alpha;
  717.       m=0;
  718.     }
  719.     l2=l+3+m/32;
  720.     p2=cgetr(l2);p3=cgetr(l2);
  721.     p4=cgetr(l2);p5=cgetr(l2);
  722.     affrr(p1,p4);
  723.   
  724.     for (i=1;i<=m;i++)
  725.     {
  726.       av=avma;affrr(mpsqrt(p4),p4);
  727.       avma=av;
  728.     }
  729.   
  730.     subrsz(p4,1,p2);
  731.     addsrz(1,p4,p3);
  732.     divrrz(p2,p3,p2);
  733.     mulrrz(p2,p2 ,p3);
  734.     l1=2;
  735.     setlg(p4,l1+2);setlg(p5,l1+2);
  736.     divssz(1,2*n+1,p4);
  737.     s=0;setlg(p3,4);
  738.     ex=expo(p3);
  739.     l2-=2;
  740.   
  741.     for (i=n;i>=1;--i)
  742.     {
  743.       mulrrz(p4,p3,p4);
  744.       divssz(1,2*i-1,p5);
  745.       s-=ex;l0=s/32;l1+=l0;
  746.       if (l1>l2) l1=l2;
  747.       s %=32;
  748.       setlg(p3,l1+2);
  749.       setlg(p4,l1+2);
  750.       setlg(p5,l1+2);
  751.       addrrz(p4,p5,p4);
  752.     }
  753.   
  754.     l2+=2;setlg(p4,l2);
  755.     mulrrz(p2,p4,y);
  756.     setexpo(y,expo(y)+m+1+m1);
  757.     if (sgn==1) mpnegz(y,y);
  758.   }
  759.   avma=avmacourant;
  760.   return y;
  761. }
  762.  
  763. GEN   teich(x)
  764.   GEN x;
  765.   
  766. {
  767.   GEN aux,y,z,p1;
  768.   long av,n,k;
  769.   
  770.   if(typ(x)!=7) err(teicher1);
  771.   if(!signe(x[4])) return gcopy(x);
  772.   y=cgetp(x); setvalp(y, 0);
  773.   if(!cmpis(x[2],2))
  774.   {
  775.     if(((GEN)x[4])[lgef(x[4])-1]&2) subisz(x[3],1,y[4]); 
  776.     else affsi(1,y[4]);
  777.   }
  778.   else
  779.   {
  780.     av=avma;p1=addsi(-1,x[2]);aux=divii(addsi(-1,x[3]),p1);
  781.     z=(GEN)x[4];n=precp(x);
  782.     for(k=1;k<n;k<<=1)
  783.       z=modii(mulii(z,addsi(1,mulii(aux,addsi(-1,puissmodulo(z,p1,x[3]))))),x[3]);
  784.     affii(z,y[4]);avma=av;
  785.   }
  786.   return y;
  787. }
  788.  
  789. GEN palogaux(x)
  790.   GEN x;
  791. {
  792.   long av,av1=avma,tetpil,k,e,pp,al;
  793.   GEN y,s,p1,y2;
  794.   
  795.   if(!cmpis(x[4],1)) 
  796.   {
  797.     y=gaddgs(x,-1);
  798.     if(!cmpis(x[2],2))
  799.     {
  800.       setvalp(y,valp(y)-1);mpshiftz(y[3],-1,y[3]);
  801.     }
  802.     return y;
  803.   }
  804.   y=gdiv(gaddgs(x,-1),gaddgs(x,1));e=valp(y);pp=precp(y);
  805.   if(cmpis(x[2],2))
  806.   {
  807.     av=avma;for(al=0,p1=stoi(e);cmpis(p1,pp+al+e)<0;al++) p1=mulii(p1,x[2]);
  808.     avma=av;
  809.   }
  810.   else al=1;
  811.   k=(pp+al+e-2)/e;if(!odd(k)) k--;
  812.   y2=gmul(y,y);s=gdivgs(gun,k);
  813.   while(k>=3)
  814.   {
  815.     k-=2;s=gadd(gmul(y2,s),gdivgs(gun,k));
  816.   }
  817.   tetpil=avma;return gerepile(av1,tetpil,gmul(s,y));
  818. }
  819.   
  820. GEN   palog(x)
  821.   GEN x;
  822. {
  823.   long av=avma,tetpil;
  824.   GEN p1,y;
  825.   
  826.   if (!signe(x[4])) err(loger6);
  827.   if(cmpis(x[2],2))
  828.   {
  829.     y=cgetp(x);setvalp(y,0);
  830.     p1=gsubgs(x[2],1);affii(puissmodulo(x[4],p1,x[3]),y[4]);
  831.     y=gmulgs(palogaux(y),2);tetpil=avma;
  832.     return gerepile(av,tetpil,gdiv(y,p1));
  833.   }
  834.   else
  835.   {
  836.     y=gsqr(x);setvalp(y,0);tetpil=avma;
  837.     return gerepile(av,tetpil,palogaux(y));
  838.   }
  839. }
  840.  
  841. GEN glog(x,prec)
  842.   
  843.    GEN   x;
  844.    long  prec;
  845.   
  846. {
  847.   long    av,tetpil,v;
  848.   GEN y,p1,p2;
  849.  
  850.   switch(typ(x))
  851.   {
  852.   case 2 : if(signe(x)>=0) y=mplog(x);
  853.   else
  854.     {
  855.     y=cgetg(3,6);y[2]=lmppi(lg(x));
  856.     setsigne(x,1);y[1]=lmplog(x);
  857.     setsigne(x,-1);
  858.     }
  859.     break;
  860.   
  861.   case 6 : y=cgetg(3,6);y[2]=larg(x,prec);
  862.     av=avma;p1=glog(gnorm(x),prec);tetpil=avma;
  863.     y[1]=lpile(av,tetpil,gmul2n(p1,-1));
  864.     break;
  865.   
  866.   case 7 : y = palog(x); break;
  867.   case 3 : err(loger3);
  868.   
  869.   case 11: av=avma;if(valp(x)) err(loger5);
  870.     v=varn(x);p1=gdiv(deriv(x,v),x);
  871.     if(gcmp1(x[2]))
  872.     {
  873.       tetpil=avma;y=gerepile(av,tetpil,integ(p1,v));
  874.     }
  875.     else
  876.     {
  877.       p1=integ(p1,v);p2=glog(x[2],prec);
  878.       tetpil=avma;y=gerepile(av,tetpil,gadd(p1,p2));
  879.     }
  880.     break;
  881.   
  882.   default: y=transc(glog,x,prec);
  883.   }
  884.   return y;
  885. }
  886.  
  887. void glogz(x,y)
  888.   
  889.    GEN   x,y;
  890.   
  891. {
  892.   long    av,prec;
  893.   GEN p;
  894.  
  895.   prec=precision(y);
  896.   if(!prec) err(loger4);
  897.   av=avma;p=glog(x,prec);
  898.   gaffect(p,y);avma=av;
  899. }
  900.  
  901. /********************************************************************/
  902. /********************************************************************/
  903. /**                                                                **/
  904. /**                      FONCTION SINCOS-1                         **/
  905. /**                                                                **/
  906. /********************************************************************/
  907. /********************************************************************/
  908.  
  909. GEN mpsc1(x,ptmod8)
  910.      long  *ptmod8;
  911.      GEN   x;
  912. {
  913.   long l,l0,l1,l2,l4,ee,i,n,m,s,t,avmacourant,av,mmax,mod4;
  914.   double alpha,beta,aa,bb,cc,dd;
  915.   GEN y,p1,p2,p3,p4,pitemp;
  916.   
  917.   mmax=23170;
  918.   if(typ(x)!=2) err(sc1er1);
  919.   else
  920.   {
  921.     if(!signe(x))
  922.     {
  923.       *ptmod8=0;y=cgetr(3);y[2]=0;y[1]=0x7fffff+(expo(x)<<1);
  924.     }
  925.     else
  926.     {
  927.       l=lg(x);y=cgetr(l);avmacourant=avma;pitemp=mppi(l+1);setexpo(pitemp,-1);
  928.       p1=cgetr(l+1);addrrz(x,pitemp,p1);setexpo(pitemp,0);
  929.       if(expo(p1)>=(32*(l-2)+3)) err(sc1er2);
  930.       av=avma;p3=divrr(p1,pitemp);p2=mpent(p3);
  931.       if(signe(p2)) subrrz(x,mulir(p2,pitemp),p1);else affrr(x,p1);
  932.       *ptmod8=(signe(p1)==-1)*4;mod4=0;
  933.       if(signe(p2))
  934.       {
  935.         mod4=mant(p2,lgef(p2)-2)&3;if((signe(p2)<0)&&(mod4)) mod4=4-mod4;
  936.         *ptmod8+=mod4;
  937.       }
  938.       avma=av;
  939.       if(gcmp0(p1)) alpha=1000000.0;
  940.       else {m=expo(p1);alpha=(m< -1023) ? -1-m*LOG2 : -1-log(fabs(rtodbl(p1)));}
  941.       beta=5+32*(l-2)*LOG2;aa=0.5/LOG2;bb=0.5*aa;
  942.       cc=aa+sqrt((beta+bb)/LOG2);
  943.       dd=((beta/cc)-alpha-log(cc))/LOG2;
  944.       if(dd>=0) {m=1+dd;n=(1+cc)/2.0;setexpo(p1,expo(p1)-m);}
  945.       else {m=0;n=(1+beta/alpha)/2.0;}
  946.       l2=l+1+m/32;
  947.       p2=cgetr(l2);p3=cgetr(l2);p4=cgetr(l2);
  948.       s=0;l1=2;
  949.       affrr(p1,p4);affsr(1,p2);mulrrz(p4,p4,p4);
  950.       setlg(p2,l1+2);setlg(p3,l1+2);setlg(p4,l1+2);
  951.       l4=l1+2;l2-=2;
  952.       if(n<=mmax)
  953.       {
  954.     setlg(p4,4);setlg(p3,4);divrsz(p4,(2*n+1)*(2*n+2),p3);ee= -expo(p3);
  955.     s=0;l1=1+(ee>>5);l4=l1+2;
  956.     setlg(p2,l4);setlg(p3,l4);setlg(p4,l4);
  957.         for(i=n;i>=2;--i)
  958.         {
  959.           divrsz(p4,2*i*(2*i-1),p3);mulrrz(p3,p2,p2);s-=expo (p3);
  960.       t=s&31;if(!t) l0=(s>>5);else l0=(s>>5)+1;l1+=l0;
  961.           if(l1>l2) {l0+=l2-l1;l1=l2;}l4+=l0;
  962.           setlg(p2,l4);setlg(p3,l4);setlg(p4,l4);
  963.           subsrz(1,p2,p2);
  964.         }
  965.       }
  966.       else
  967.       {
  968.     setlg(p4,4);setlg(p3,4);divrsz(p4,2*n+1,p3);divrsz(p3,2*n+2,p3);ee= -expo(p3);
  969.     s=0;l1=1+(ee>>5);l4=l1+2;
  970.     setlg(p2,l4);setlg(p3,l4);setlg(p4,l4);
  971.         for(i=n;i>=2;--i)
  972.         {
  973.           divrsz(p4,2*i,p3);divrsz(p3,2*i-1,p3);mulrrz(p3,p2,p2);s-=expo (p3);
  974.       t=s&31;if(!t) l0=(s>>5);else l0=(s>>5)+1;l1+=l0;
  975.           if(l1>l2) {l0+=l2-l1;l1=l2;}l4+=l0;
  976.           setlg(p2,l4);setlg(p3,l4);setlg(p4,l4);
  977.           subsrz(1,p2,p2);
  978.         }
  979.       }
  980.       setlg(p2,l4);setlg(p3,l4);setlg(p4,l4);
  981.       setexpo(p4,expo(p4)-1);setsigne(p4,-signe (p4));
  982.       mulrrz(p4,p2,p2);
  983.       for(i=1;i<=m;i++)
  984.       {
  985.         addsrz(2,p2,p3);mulrrz(p2,p3,p2);setexpo(p2,expo(p2)+1);
  986.       }
  987.       mpaff(p2,y);avma=avmacourant;
  988.     }
  989.   }
  990.   return y;
  991. }
  992.  
  993. /********************************************************************/
  994. /********************************************************************/
  995. /**                                                                **/
  996. /**                      FONCTION sqrt(-x*(x+2))                   **/
  997. /**                                                                **/
  998. /********************************************************************/
  999. /********************************************************************/
  1000.  
  1001. GEN mpaut(x)
  1002.   
  1003.    GEN   x;
  1004.   
  1005. {
  1006.   long    av,tetpil;
  1007.   GEN y,p1;
  1008.  
  1009.   av=avma;
  1010.   p1=addsr(2,x);
  1011.   p1=mulrr(x,p1);
  1012.   setsigne(p1,-signe(p1));tetpil=avma;
  1013.   y=gerepile(av,tetpil,mpsqrt(p1));
  1014.   return y;
  1015. }
  1016.  
  1017. /********************************************************************/
  1018. /********************************************************************/
  1019. /**                                                                **/
  1020. /**                       FONCTION COSINUS                         **/
  1021. /**                                                                **/
  1022. /********************************************************************/
  1023. /********************************************************************/
  1024.  
  1025. GEN mpcos(x)
  1026.   
  1027.    GEN   x;
  1028.   
  1029. {
  1030.   long    mod8,av,tetpil;
  1031.   GEN y,p1;
  1032.  
  1033.   if(typ(x)!=2) err(coser1);
  1034.   if(!signe(x))
  1035.   {
  1036.     y=addsr(1,x);
  1037.   }
  1038.   else
  1039.   {
  1040.     av=avma;
  1041.     p1=mpsc1(x,&mod8);tetpil=avma;
  1042.     switch(mod8)
  1043.     {
  1044.     case 0: y=addsr(1,p1);break;
  1045.     case 1: y=mpaut(p1);setsigne(y,-signe(y));break;
  1046.     case 2: y=subsr(-1,p1);break;
  1047.     case 3: y=mpaut(p1);break;
  1048.     case 4: y=addsr(1,p1);break;
  1049.     case 5: y=mpaut(p1);break;
  1050.     case 6: y=subsr(-1,p1);break;
  1051.     case 7: y=mpaut(p1);setsigne(y,-signe(y));break;
  1052.     default:;
  1053.     }
  1054.     y=gerepile(av,tetpil,y);
  1055.   }
  1056.   return y;
  1057. }
  1058.  
  1059. GEN gcos(x,prec)
  1060.   
  1061.    GEN   x;
  1062.    long  prec;
  1063.   
  1064. {
  1065.   long    av,tetpil;
  1066.   GEN r,u,v;
  1067.   GEN y,p1,p2;
  1068.  
  1069.   switch(typ(x))
  1070.   {
  1071.   case 2 : y=mpcos(x);break;
  1072.   case 6 : y=cgetg(3,6);av=avma;
  1073.     r=gexp(x[2],prec);p1=ginv(r);
  1074.     p2=gmul2n(mpadd(p1,r),-1);
  1075.     p1=mpsub(p2,r);
  1076.     gsincos(x[1],&u,&v,prec);
  1077.     tetpil=avma;
  1078.     y[1]=lmul(p2,v);y[2]=lmul(p1,u);
  1079.     gerepile(av,tetpil,1);
  1080.     break;
  1081.   
  1082.   case 3 :
  1083.   case 7 : err(coser2);
  1084.   
  1085.   case 11: if(gcmp0(x)) y=gaddsg(1,x);
  1086.   else
  1087.     {
  1088.     av=avma;if(valp(x)<0) err(coser4);
  1089.     gsincos(x,&u,&v,prec);tetpil=avma;
  1090.     y=gerepile(av,tetpil,gcopy(v));
  1091.     }
  1092.     break;
  1093.   
  1094.   default: y=transc(gcos,x,prec);
  1095.   }
  1096.   return y;
  1097. }
  1098.  
  1099. void gcosz(x,y)
  1100.   
  1101.    GEN   x,y;
  1102.   
  1103. {
  1104.   long    av,prec;
  1105.   GEN p;
  1106.  
  1107.   prec=precision(y);
  1108.   if(!prec) err(coser3);
  1109.   av=avma;p=gcos(x,prec);
  1110.   gaffect(p,y);avma=av;
  1111. }
  1112.  
  1113. /********************************************************************/
  1114. /********************************************************************/
  1115. /**                                                                **/
  1116. /**                       FONCTION SINUS                           **/
  1117. /**                                                                **/
  1118. /********************************************************************/
  1119. /********************************************************************/
  1120.  
  1121. GEN mpsin(x)
  1122.   
  1123.    GEN   x;
  1124.   
  1125. {
  1126.   long    mod8,av,tetpil;
  1127.   GEN y,p1;
  1128.  
  1129.   if(typ(x)!=2) err(siner1);
  1130.   if(!signe(x))
  1131.   {
  1132.     y=cgetr(3);y[1]=x[1];y[2]=0;
  1133.   }
  1134.   else
  1135.   {
  1136.     y=cgetr(lg(x));
  1137.     av=avma;
  1138.     p1=mpsc1(x,&mod8);tetpil=avma;
  1139.     switch(mod8)
  1140.     {
  1141.     case 0: y=mpaut(p1);break;
  1142.     case 1: y=addsr(1,p1);break;
  1143.     case 2: y=mpaut(p1);setsigne(y,-signe(y));break;
  1144.     case 3: y=subsr(-1,p1);break;
  1145.     case 4: y=mpaut(p1);setsigne(y,-signe(y));break;
  1146.     case 5: y=addsr(1,p1);break;
  1147.     case 6: y=mpaut(p1);break;
  1148.     case 7: y=subsr(-1,p1);break;
  1149.     default:;
  1150.     }
  1151.     y=gerepile(av,tetpil,y);
  1152.   }
  1153.   return y;
  1154. }
  1155.  
  1156. GEN gsin(x,prec)
  1157.   
  1158.    GEN   x;
  1159.    long  prec;
  1160.   
  1161. {
  1162.   long    av,tetpil;
  1163.   GEN r,u,v;
  1164.   GEN y,p1,p2;
  1165.  
  1166.   switch(typ(x))
  1167.   {
  1168.   case 2 : y=mpsin(x);break;
  1169.   case 6 : y=cgetg(3,6);av=avma;
  1170.     r=gexp(x[2],prec);p1=ginv(r);
  1171.     p2=gmul2n(mpadd(p1,r),-1);
  1172.     p1=mpsub(p2,p1);
  1173.     gsincos(x[1],&u,&v,prec);
  1174.     tetpil=avma;
  1175.     y[1]=lmul(p2,u);y[2]=lmul(p1,v);
  1176.     gerepile(av,tetpil,1);
  1177.     break;
  1178.   
  1179.   case 3 :
  1180.   case 7 : err(siner2);
  1181.   
  1182.   case 11: if(gcmp0(x)) y=gcopy(x);
  1183.   else
  1184.     {
  1185.     av=avma;if(valp(x)<0) err(siner4);
  1186.     gsincos(x,&u,&v,prec);tetpil=avma;
  1187.     y=gerepile(av,tetpil,gcopy(u));
  1188.     }
  1189.     break;
  1190.   
  1191.   default: y=transc(gsin,x,prec);
  1192.   }
  1193.   return y;
  1194. }
  1195.  
  1196. void gsinz(x,y)
  1197.   
  1198.    GEN   x,y;
  1199.   
  1200. {
  1201.   long    av,prec;
  1202.   GEN p;
  1203.  
  1204.   prec=precision(y);
  1205.   if(!prec) err(siner3);
  1206.   av=avma;p=gsin(x,prec);
  1207.   gaffect(p,y);avma=av;
  1208. }
  1209.  
  1210. /********************************************************************/
  1211. /********************************************************************/
  1212. /**                                                                **/
  1213. /**                    PROCEDURE SINUS,COSINUS                     **/
  1214. /**                                                                **/
  1215. /********************************************************************/
  1216. /********************************************************************/
  1217.  
  1218. void mpsincos(x,s,c)
  1219.   
  1220.    GEN   x,*s,*c;
  1221.   
  1222. {
  1223.   long    av,av1,tetpil,dec,mod8;
  1224.   GEN p1;
  1225.  
  1226.   if(typ(x)!=2) err(sicoer1);
  1227.   if(!signe(x))
  1228.   {
  1229.     p1=cgetr(3);*s=p1;p1[1]=x[1];
  1230.     p1[2]=0;*c=addsr(1,x);
  1231.   }
  1232.   else
  1233.   {
  1234.     av=avma;
  1235.     p1=mpsc1(x,&mod8);tetpil=avma;
  1236.   
  1237.     switch(mod8)
  1238.     {
  1239.       case 0: *c=addsr(1,p1);*s=mpaut(p1);break;
  1240.       case 1: *s=addsr(1,p1);*c=mpaut(p1);
  1241.            setsigne(*c,-signe(*c)); break;
  1242.       case 2: *c=subsr(-1,p1);*s=mpaut(p1);
  1243.            setsigne(*s,-signe(*s)); break;
  1244.       case 3: *s=subsr(-1,p1);*c=mpaut(p1); break;
  1245.       case 4: *c=addsr(1,p1);*s=mpaut(p1);
  1246.            setsigne(*s,-signe(*s)); break;
  1247.       case 5: *s=addsr(1,p1);*c=mpaut(p1); break;
  1248.       case 6: *c=subsr(-1,p1);*s=mpaut(p1); break;
  1249.       case 7: *s=subsr(-1,p1);*c=mpaut(p1);
  1250.            setsigne(*c,-signe(*c)); break;
  1251.     }
  1252.     av1=avma;dec=lpile(av,tetpil,0)>>2;
  1253.     if(adecaler(*s,tetpil,av1)) (*s)+=dec;
  1254.     if(adecaler(*c,tetpil,av1)) (*c)+=dec;
  1255.   }
  1256. }
  1257.  
  1258. void gsincos(x,s,c,prec)
  1259.   
  1260.    GEN   x,*s,*c;
  1261.    long  prec;
  1262.   
  1263. {
  1264.   long    av,av1,tetpil,dec,ii,i,j,ex,ex2,lx,ly;
  1265.   GEN r,u,v,u1,v1,p1,p2,p3,p4,ps,pc;
  1266.  
  1267.   switch(typ(x))
  1268.   {
  1269.     case 1 :
  1270.     case 4 :
  1271.     case 5 : av=avma;p1=cgetr(prec);gaffect(x,p1);
  1272.       tetpil=avma;mpsincos(p1,s,c);av1=avma;
  1273.       dec=lpile(av,tetpil,0)>>2;
  1274.       if(adecaler(*s,tetpil,av1)) *s+=dec;
  1275.       if(adecaler(*c,tetpil,av1)) *c+=dec;
  1276.       break;
  1277.     case 2 : mpsincos(x,s,c);break;
  1278.     case 6 : ps=cgetg(3,6);pc=cgetg(3,6);
  1279.       *s=ps;*c=pc;av=avma;
  1280.       r=gexp(x[2],prec);p1=ginv(r);
  1281.       p2=gmul2n(mpadd(p1,r),-1);
  1282.       p1=mpsub(p2,p1);r=mpsub(p2,r);
  1283.       gsincos(x[1],&u,&v,prec);
  1284.       tetpil=avma;
  1285.       ps[1]=lmul(p2,u);ps[2]=lmul(p1,v);
  1286.       pc[1]=lmul(p2,v);pc[2]=lmul(r,u);
  1287.       av1=avma;dec=lpile(av,tetpil,0);
  1288.       if(adecaler(ps[1],tetpil,av1)) ps[1]+=dec;
  1289.       if(adecaler(ps[2],tetpil,av1)) ps[2]+=dec;
  1290.       if(adecaler(pc[1],tetpil,av1)) pc[1]+=dec;
  1291.       if(adecaler(pc[2],tetpil,av1)) pc[2]+=dec;
  1292.       break;
  1293.     
  1294.     case 11: if(gcmp0(x)) {*c=gaddsg(1,x);*s=gcopy(x);}
  1295.     else
  1296.     {
  1297.       ex=valp(x);lx=lg(x);ex2=2*ex+2;
  1298.       if(ex<0) err(sicoer3);
  1299.       if(ex2>lx)
  1300.       {
  1301.         *s=gcopy(x);av=avma;p1=gdivgs(gmul(x,x),2);
  1302.         tetpil=avma;*c=gerepile(av,tetpil,gsubsg(1,p1));
  1303.       }
  1304.       if(!ex)
  1305.       {
  1306.         av=avma;p1=gcopy(x);p1[2]=zero;
  1307.         normalize(&p1);gsincos(p1,&u,&v,prec);
  1308.         gsincos(x[2],&u1,&v1,prec);
  1309.         p1=gmul(v1,v);p2=gmul(u1,u);p3=gmul(v1,u);
  1310.         p4=gmul(u1,v);tetpil=avma;
  1311.         *c=gsub(p1,p2);*s=gadd(p3,p4);
  1312.         av1=avma;dec=lpile(av,tetpil,0)>>2;
  1313.         if(adecaler(*c,tetpil,av1)) (*c)+=dec;
  1314.     if(adecaler(*s,tetpil,av1)) (*s)+=dec;
  1315.       }
  1316.       else
  1317.       {
  1318.         ly=lx+2*ex;
  1319.         pc=cgetg(ly,11);ps=cgetg(lx,11);
  1320.         *c=pc;*s=ps;pc[1]=0x1008000;setvarn(pc,varn(x));
  1321.         pc[2]=un;ps[1]=x[1];
  1322.         for(i=2;i<ex+2;i++) ps[i]=lcopy(x[i]);
  1323.         for(i=3;i<ex2;i++) pc[i]=zero;
  1324.         for(i=ex2;i<ly;i++)
  1325.         {
  1326.           ii=i-ex;av=avma;p1=gzero;
  1327.           for(j=ex;j<ii-1;j++)
  1328.             p1=gadd(p1,gmulgs(gmul(x[j-ex+2],ps[ii-j]),j));
  1329.           tetpil=avma;
  1330.           pc[i]=lpile(av,tetpil,gdivgs(p1,2-i));
  1331.           if(ii<lx)
  1332.           {
  1333.             av=avma;p1=gzero;
  1334.             for(j=ex;j<=i-ex2;j++)
  1335.               p1=gadd(p1,gmulgs(gmul(x[j-ex+2],pc[i-j]),j));
  1336.             p1=gdivgs(p1,i-2);tetpil=avma;
  1337.             ps[i-ex]=lpile(av,tetpil,gadd(p1,x[i-ex]));
  1338.           }
  1339.         }
  1340.       }
  1341.     }
  1342.     break;
  1343.     
  1344.     default: err(sicoer2);
  1345.   }
  1346.   return;
  1347. }
  1348.  
  1349. /********************************************************************/
  1350. /********************************************************************/
  1351. /**                                                                **/
  1352. /**                      FONCTION TANGENTE                         **/
  1353. /**                                                                **/
  1354. /********************************************************************/
  1355. /********************************************************************/
  1356.  
  1357. GEN mptan(x)
  1358.   
  1359.    GEN   x;
  1360.   
  1361. {
  1362.   long    av,tetpil;
  1363.   GEN s,c,y;
  1364.  
  1365.   av=avma;
  1366.   mpsincos(x,&s,&c);tetpil=avma;
  1367.   y=gerepile(av,tetpil,divrr(s,c));
  1368.   return y;
  1369. }
  1370.  
  1371. GEN gtan(x,prec)
  1372.   
  1373.    GEN   x;
  1374.    long  prec;
  1375.   
  1376. {
  1377.   long    av,tetpil;
  1378.   GEN s,c;
  1379.   GEN y;
  1380.  
  1381.   switch(typ(x))
  1382.   {
  1383.   case 2 : y=mptan(x);break;
  1384.   case 6 : av=avma;gsincos(x,&s,&c,prec);
  1385.     tetpil=avma;
  1386.     y=gerepile(av,tetpil,gdiv(s,c));
  1387.     break;
  1388.   
  1389.   case 3 :
  1390.   case 7 : err(taner1);
  1391.   
  1392.   case 11: if(gcmp0(x)) y=gcopy(x);
  1393.     else
  1394.     {
  1395.       av=avma;if(valp(x)<0) err(taner3);
  1396.       gsincos(x,&s,&c,prec);tetpil=avma;
  1397.       y=gerepile(av,tetpil,gdiv(s,c));
  1398.     }
  1399.     break;
  1400.   
  1401.   default: y=transc(gtan,x,prec);
  1402.   }
  1403.   return y;
  1404. }
  1405.  
  1406. void gtanz(x,y)
  1407.   
  1408.    GEN   x,y;
  1409.   
  1410. {
  1411.   long    av,prec;
  1412.   GEN p;
  1413.  
  1414.   prec=precision(y);
  1415.   if(!prec) err(taner2);
  1416.   av=avma;p=gtan(x,prec);
  1417.   gaffect(p,y);avma=av;
  1418. }
  1419.