home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / arith1 < prev    next >
Encoding:
Text File  |  1991-11-28  |  39.7 KB  |  1,527 lines
  1. /*********************************************************************/
  2. /*********************************************************************/
  3. /**                                                                 **/
  4. /**                    FONCTIONS ARITHMETIQUES                      **/
  5. /**                                                                 **/
  6. /**                       (premiere partie)                         **/
  7. /**                                                                 **/
  8. /**                      copyright Babe Cool                        **/
  9. /**                                                                 **/
  10. /**                                                                 **/
  11. /*********************************************************************/
  12. /*********************************************************************/
  13.  
  14. #include  "genpari.h"
  15.  
  16. /*********************************************************************/
  17. /*       Decomposition binaire d'un I < 2^32 dans un vecteur         */
  18. /*********************************************************************/
  19.  
  20. GEN binaire(x)
  21.      GEN x;
  22.      
  23. {
  24.   unsigned m = 0x80000000, u,lx=lgef(x),ly=33;
  25.   GEN y;
  26.   
  27.   if (lx==2) {y=cgetg(2,17);y[1]=zero;}
  28.   else
  29.     {
  30.       u=x[--lx];
  31.       while( !(m&u) ) {m >>= 1 ;ly--;}
  32.       y=cgetg(ly,17);ly=1;
  33.       do {y[ly]=m&u ? un : zero;ly++;} while( m>>=1 ) ;
  34.     }
  35.   return y;
  36. }
  37.  
  38. /* Renvoie 0 ou 1 selon que le bit numero n de x est a 0 ou 1 */
  39.  
  40. int bittest(x,n)
  41.      GEN x;
  42.      long n;
  43.      
  44. {
  45.   long l,n2;
  46.   
  47.   if(typ(x)!=1) err(arither1);
  48.   if((!signe(x))||(n<0)) return 0;
  49.   l=lgef(x)-1-(n>>5);
  50.   if(l<=1) return 0;
  51.   n2=(1<<(n&31))&x[l];return n2 ? 1 : 0;
  52. }
  53.  
  54. /*********************************************************************/
  55. /**                                                                 **/
  56. /**            ORDRE DE x entier MODULO n dans (Z/nZ)*              **/
  57. /**                                                                 **/
  58. /*********************************************************************/
  59.  
  60. GEN order(x) 
  61.      GEN x;
  62.      
  63. {
  64.   long av=avma,av1,m,i,p,e,u;
  65.   GEN y,t,o,o1;
  66.   
  67.   if(typ(x)!=3) err(orderer);
  68.   if(!gcmp1(mppgcd(m=x[1],x[2]))) err(orderer);
  69.   t=decomp(o=phi(m));
  70.   for(i=lg(t[1])-1;i;i--)
  71.     {
  72.       p=coeff(t,i,1);e=itos(coeff(t,i,2));
  73.       do
  74.     { 
  75.       y=gpui(x,o1=divii(o,p));e--;
  76.       if(u=gcmp1(y[2])) o=o1;
  77.     }
  78.       while(e&&u);
  79.     }
  80.   av1=avma;
  81.   return gerepile(av,av1,gcopy(o));
  82. }
  83.  
  84. /******************************************************************/
  85. /**               GENERATEUR DE (Z/mZ)* lorsque m=p^e             */
  86. /******************************************************************/
  87.  
  88. GEN gener(m) 
  89.      GEN m;
  90.      
  91. {
  92.   long av=avma,av1,k,i,p,e,u;
  93.   GEN fi,x,xi,y,t,q;
  94.   
  95.   if(typ(m)!=1) err(arither1);
  96.   t=decomp(m);
  97.   if(lg(t[1])!=2) err(generer);
  98.   p=coeff(t,1,1);e=itos(coeff(t,1,2));
  99.   q=fi=subis(p,1);
  100.   if(e>1) fi=mulii(q,gpuigs(p,e-1));
  101.   t=decomp(q);k=lg(t[1])-1;
  102.   xi=stoi(1);
  103.   do
  104.     {
  105.       xi[2]++;
  106.       i=k;u=1;
  107.       if(gcmp1(mppgcd(xi,m)))
  108.     {
  109.       x=gmodulcp(xi,m);
  110.       if(e>1)
  111.         {
  112.           y=gpui(x,divii(fi,p));
  113.           if(gcmp1(y[2])) {i=u=0;}
  114.         }
  115.       while(i)
  116.         {
  117.           y=gpui(x,divii(fi,coeff(t,i,1)));
  118.           if (gcmp1(y[2])) {i=u=0;} else i--;
  119.         }
  120.     }
  121.       else u=0;
  122.     }
  123.   while(!u);
  124.   av1=avma;
  125.   return gerepile(av,av1,gcopy(x));
  126. }
  127.  
  128. /*********************************************************************/
  129. /*********************************************************************/
  130. /**                                                                 **/
  131. /**                     FONCTION RACINE                             **/
  132. /**                                                                 **/
  133. /*********************************************************************/
  134. /*********************************************************************/
  135.  
  136. GEN racine(a)
  137.      GEN a;
  138. {
  139.   GEN x,y,z;
  140.   long av,av2,k,count=0;
  141.   if (typ(a)!=1) err(arither1);
  142.   switch (signe(a))
  143.     {
  144.       case-1:
  145.       x=cgetg(3,6); x[1]=zero;
  146.       setsigne(a, 1); x[2]=(long) racine(a); setsigne(a, -1);
  147.       return x;
  148.     case 0 : return gzero;
  149.     case 1 :
  150.       k=sqrt((double)(unsigned long)mant(a,1));
  151.       x=shifts(k+1,(lgef(a)-3)*16);
  152.       av=avma;
  153.       y=cgeti(lgef(x));
  154.       av2=avma;
  155.       do
  156.     {
  157.       divisz(addii(x,divii(a,x)),2,y);
  158.       z=y;y=x;x=z;
  159.       avma=av2;
  160.       count++;
  161.     }
  162.       while (cmpii(x,y)<0);
  163.       if (odd(count)) x=y;else affii(y,x);
  164.       avma=av;
  165.       return x;
  166.     }
  167. }
  168.  
  169. /*********************************************************************/
  170. /*********************************************************************/
  171. /**                                                                 **/
  172. /**                     FONCTION CARRE PARFAIT                      **/
  173. /**                                                                 **/
  174. /*********************************************************************/
  175. /*********************************************************************/
  176.  
  177.  
  178. long carrecomplet(x,pt)
  179.      
  180.      GEN x,*pt;
  181. {
  182.   
  183.   static carresmod64[]={1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,
  184.               1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  185.               0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,
  186.               0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0};
  187.   static carresmod63[]={1,1,0,0,1,0,0,1,0,1,0,0,0,0,0,0,1,0,
  188.               1,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,
  189.               1,1,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,
  190.               0,0,0,0,1,0,0,0,0};
  191.   static carresmod65[]={1,1,0,0,1,0,0,0,0,1,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,
  192.               1,0,0,1,1,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,1,0,1,
  193.               0,0,0,1,1,0,0,0,0,1,0,0,1};
  194.   static carresmod11[]={1,1,0,1,1,1,0,0,0,1,0};
  195.   
  196.   GEN y;
  197.   long av,t,result,tetpil;
  198.   if (typ(x)!=1) err(arither1);
  199.   switch(signe(x))
  200.     {
  201.     case -1: return 0;
  202.     case 0: {*pt=gzero;return 1;}
  203.     case 1:  if (!carresmod64[63&(mant(x,lgef(x)-2))]) return 0;
  204.     }
  205.   av=avma;
  206.   t=mant(modis(x,63),1);if (!carresmod63[t]) {avma=av;return 0;}
  207.   t=mant(modis(x,65),1);if (!carresmod65[t]) {avma=av;return 0;}
  208.   t=mant(modis(x,11),1);if (!carresmod11[t]) {avma=av;return 0;}
  209.   y=racine(x);
  210.   result=cmpii(mulii(y,y),x);
  211.   if(result) {avma=av;return 0;} 
  212.   else {tetpil=avma;*pt=gerepile(av,tetpil,gcopy(y));return 1;}
  213. }
  214.  
  215. long carreparfait(x)
  216.      GEN x;
  217. {
  218.   GEN p1;
  219.   long av=avma,f;
  220.   
  221.   f=carrecomplet(x,&p1);
  222.   avma=av;return f;
  223. }
  224.  
  225. /*********************************************************************/
  226. /*********************************************************************/
  227. /**                                                                 **/
  228. /**                     SYMBOLE DE KRONECKER                        **/
  229. /**                                                                 **/
  230. /*********************************************************************/
  231. /*********************************************************************/
  232.  
  233. long kronecker(x,y)
  234.      GEN x,y;
  235. {
  236.   GEN x1,y1,z;
  237.   long av,r,s=1;
  238.   if ((typ(x)!=1) || (typ(y)!=1)) err(arither1);
  239.   av=avma;
  240.   switch (signe(y))
  241.     {
  242.       case-1: y1=negi(y);if (signe(x)<0) s= -1;break;
  243.     case 0: return (lgef(x)==3)&&(x[2]==1);
  244.     case 1: y1=y;
  245.     }
  246.   if (r=vali(y1))
  247.     if (mpodd(x))
  248.       {
  249.     if (odd(r)&&(abs((x[lgef(x)-1]&7)-4)==1)) s= -s;
  250.     y1=shifti(y1,-r);
  251.       }
  252.     else {avma=av;return 0;}
  253.   x1=modii(x,y1);
  254.   while (signe(x1))
  255.     {
  256.       if (r=vali(x1))
  257.     {
  258.       if (odd(r)&&(abs((y1[lgef(y1)-1]&7)-4)==1)) s= -s;
  259.       x1=shifti(x1,-r);
  260.     }
  261.       if ((y1[lgef(y1)-1]&2)&&(x1[lgef(x1)-1]&2)) s= -s;
  262.       z=resii(y1,x1);y1=x1;x1=z;
  263.     }
  264.   avma=av;
  265.   return cmpsi(1,y1) ? 0 : s;
  266. }
  267.  
  268. long  kro8(x,y)
  269.      
  270.      GEN  x,y;
  271.      
  272.      /*  a usage interne: aucune verification de types  */
  273.      
  274. {
  275.   GEN  p1,p2;
  276.   long k,av=avma;
  277.   
  278.   p1=(GEN)(x[1]);p2=(GEN)(p1[3]);
  279.   p1=gsub(gmul(p2,p2),gmul2n(p1[2],2));
  280.   k=kronecker(p1,y);
  281.   avma=av;return k;
  282. }
  283.  
  284. long krogs(x,y)
  285.      GEN x;
  286.      long y;
  287. {
  288.   long av,r,s=1,x1,z;
  289.   
  290.   if (typ(x)!=1) err(arither1);
  291.   av=avma;
  292.   if(y<=0)
  293.     {
  294.       if(y) {y= -y;if(signe(x)<0) s= -1;}
  295.       else  return (lgef(x)==3)&&(x[2]==1);
  296.     }
  297.   if (r=vals(y))
  298.     if (mpodd(x))
  299.       {
  300.     if (odd(r)&&(abs((x[lgef(x)-1]&7)-4)==1)) s= -s;
  301.     y>>=r;
  302.       }
  303.     else return 0;
  304.   x1=itos(modis(x,y));
  305.   while (x1)
  306.     {
  307.       if (r=vals(x1))
  308.     {
  309.       if (odd(r)&&(abs((y&7)-4)==1)) s= -s;
  310.       x1>>=r;
  311.     }
  312.       if ((y&2)&&(x1&2)) s= -s;
  313.       z=y%x1;y=x1;x1=z;
  314.     }
  315.   avma=av;
  316.   return (y==1) ? s : 0;
  317. }
  318.  
  319. long  krosg(s,x)
  320.      
  321.      GEN  x;
  322.      long s;
  323.      
  324. {
  325.   long av,y;
  326.   
  327.   av=avma;y=kronecker(stoi(s),x);
  328.   avma=av;return y;
  329. }
  330.  
  331. long kross(x,y)
  332.      long x,y;
  333. {
  334.   long r,s=1,x1,z;
  335.   
  336.   if(y<=0)
  337.     {
  338.       if(y) {y= -y;if(x<0) s= -1;}
  339.       else  return (abs(x)==1);
  340.     }
  341.   if (r=vals(y))
  342.     if (odd(x))
  343.       {
  344.     if (odd(r)&&(abs((x&7)-4)==1)) s= -s;
  345.     y>>=r;
  346.       }
  347.     else return 0;
  348.   x1=x%y;if(x1<0) x1+=y;
  349.   while (x1)
  350.     {
  351.       if (r=vals(x1))
  352.     {
  353.       if (odd(r)&&(abs((y&7)-4)==1)) s= -s;
  354.       x1>>=r;
  355.     }
  356.       if ((y&2)&&(x1&2)) s= -s;
  357.       z=y%x1;y=x1;x1=z;
  358.     }
  359.   return (y==1) ? s : 0;
  360. }
  361.  
  362. long  hil(x,y,p)
  363.      
  364.      GEN x,y,p;
  365.      
  366. #define eps(t) (((signe(t)*(t[lgef(t)-1]))&3)==3)
  367. #define ome(t) ((((t[lgef(t)-1])&7)==3)||(((t[lgef(t)-1])&7)==5))
  368.      
  369. {
  370.   long a,b,av,tx=typ(x),ty=typ(y),z;
  371.   GEN p1,p2,u,v;
  372.   
  373.   if(gcmp0(x)||gcmp0(y)) return 0;
  374.   if(tx>ty) {p1=x;x=y;y=p1;av=tx,tx=ty;ty=av;}
  375.   av=avma;
  376.   switch(tx)
  377.     {
  378.     case 1: switch(ty)
  379.       {
  380.       case 1: if(signe(p)<=0) {z=((signe(x)<0)&&(signe(y)<0))?-1:1; return z;}
  381.     a=pvaluation(x,p,&u);b=pvaluation(y,p,&v);
  382.     if(cmpsi(2,p))
  383.           {
  384.             z=((odd(a)&&odd(b)&&eps(p)))? -1:1;
  385.             if(odd(a)&&(kronecker(v,p)<0)) z= -z;
  386.             if(odd(b)&&(kronecker(u,p)<0)) z= -z;
  387.           }
  388.     else
  389.           {
  390.             z=(eps(u)&&eps(v))? -1:1;
  391.             if(odd(a)&&ome(v)) z= -z;
  392.             if(odd(b)&&ome(u)) z= -z;
  393.           }
  394.     avma=av;return z;
  395.       case 2: z=((signe(x)<0)&&(signe(y)<0))?-1:1; return z;
  396.       case 3: if(!cmpsi(2,y[1])) err(hiler1);
  397.     return hil(x,y[2],y[1]);
  398.       case 4:
  399.       case 5: p1=mulii(y[1],y[2]);z=hil(x,p1,p);
  400.     avma=av;return z;
  401.       case 7: 
  402.     if((!cmpsi(2,y[2]))&&precp(y)<=2) err(hiler1);
  403.     p1=(odd(valp(y)))?mulii(y[2],y[4]):(GEN)y[4];
  404.     z=hil(x,p1,y[2]);avma=av;return z;
  405.       default: err(hiler2);
  406.       }
  407.     case 2: if((ty<4)||(ty>5)) err(hiler2);
  408.       if(signe(x)>0) return 1;
  409.       return signe(y[1])*signe(y[2]);
  410.     case 3: if(!cmpsi(2,y[1])) err(hiler1);
  411.       switch(ty)
  412.     {
  413.         case 3: if(cmpii(x[1],y[1])) err(hiler2);
  414.           return hil(x[2],y[2],x[1]);
  415.         case 4:
  416.         case 5: return hil(x[2],y,x[1]);
  417.         case 7: if(cmpii(x[1],y[2])) err(hiler2);
  418.           return hil(x[2],y);
  419.         default: err(hiler2);
  420.     }
  421.     case 4:
  422.     case 5: p1=mulii(x[1],x[2]);switch(ty)
  423.       {
  424.       case 4:
  425.       case 5: p2=mulii(y[1],y[2]);z=hil(p1,p2,p);avma=av;return z;
  426.       case 7: z=hil(p1,y);avma=av;return z;
  427.       default: err(hiler2);
  428.       }
  429.     case 7: if((ty>7)||(cmpii(x[2],y[2]))) err(hiler2);
  430.       p1=(odd(valp(x)))?mulii(x[2],x[4]):(GEN)x[4];
  431.       p2=(odd(valp(y)))?mulii(y[2],y[4]):(GEN)y[4];z=hil(p1,p2,x[2]);
  432.       avma=av;return z;
  433.     default: err(hiler2);
  434.     }
  435. }
  436.  
  437. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  438. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  439. /*                                                                 */
  440. /*                      RACINE CARREE MODULO                       */
  441. /*                                                                 */
  442. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  443. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  444.  
  445. #define sqrmod(x,p) (modii(mulii(x,x),p))
  446.  
  447. long mpsqrtmod(a,p,pr)
  448.      GEN     a,p,*pr;
  449.      
  450. {
  451.   long  l,av0,av1,av3,f,i,k,e,r;
  452.   GEN   p1,p2,p3,m,m1,v,y,w,n;
  453.   
  454.   if ((typ(a)!=1) || (typ(p)!=1)) err(arither1);
  455.   av0=avma;
  456.   l=lg(p);v=cgeti(l);av1=avma;
  457.   w=cgeti(l);y=cgeti(l);
  458.   p1=addsi(-1,p);e=vali(p1);
  459.   p2=shifti(p1,-e);n=stoi(2);av3=avma;
  460.   if(e==1) affii(p1,y);
  461.   else do
  462.     {
  463.       m=puissmodulo(n,p2,p);m1=m;
  464.       for(i=1,f=0;(i<e)&&!f;i++)
  465.     f=gcmp1(m=sqrmod(m,p));
  466.       if(f) addsiz(1,n,n);else  affii(m1,y);
  467.       avma=av3;
  468.     }
  469.   while(f);
  470.   
  471.   /*      y contient un generateur de (Z/pZ)* eleve a la puis (p-1)/2^e   */
  472.   
  473.   p1=shifti(p2,-1);
  474.   p2=puissmodulo(a,p1,p);
  475.   if(!signe(p2)) {avma=av0;*pr=gzero;return 1;}
  476.   p3=mulii(p2,a);modiiz(p3,p,v);
  477.   p1=mulii(mulii(p2,p2),a);modiiz(p1,p,w);
  478.   avma=av3;
  479.   r=e;
  480.   while(!gcmp1(w))
  481.     {
  482.       for(k=1 ,p1=w;!gcmp1(p1=sqrmod(p1,p));k++);
  483.       avma=av3;
  484.       if(k==r) {avma=av0;return 0;}
  485.       for(i=1,p1=y;i<r-k;i++)  p1=sqrmod(p1,p);
  486.       modiiz(mulii(p1,v),p,v);
  487.       p1=sqrmod(p1,p);modiiz(mulii(p1,w),p, w);
  488.       affii(p1,y);avma=av3;r=k;
  489.     }
  490.   p1=subii(p,y);if(cmpii(y,p1)>0) affii(p1,y);
  491.   *pr=v;avma=av1;return 1;
  492. }
  493.  
  494.  
  495. /*********************************************************************/
  496. /*********************************************************************/
  497. /**                                                                 **/
  498. /**                     FONCTION PGCD                               **/
  499. /**                                                                 **/
  500. /*********************************************************************/
  501. /*********************************************************************/
  502.  
  503. GEN mppgcd1(a,b)
  504.      GEN a,b;
  505. {
  506.   GEN d,d1;
  507.   long av,count=0;
  508.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  509.   if (lgef(a)<lgef(b))
  510.     {
  511.       d=b;b=a;a=d;
  512.     }
  513.   d=gcopy(a);
  514.   av=avma;
  515.   d1=gcopy(b);
  516.   while (signe(d1))
  517.     {
  518.       resiiz(d,d1,d);
  519.       a=d;d=d1;d1=a;
  520.       count++;
  521.     }
  522.   if (odd(count))
  523.     {
  524.       affii(d,d1);d=d1;
  525.     }
  526.   avma=av;
  527.   if (signe(d)== -1) setsigne(d,1);
  528.   return d;
  529. }
  530.  
  531. GEN mppgcd2(a,b)
  532.      GEN a,b;
  533. {
  534.   GEN r;
  535.   long av=avma,tetpil;
  536.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  537.   if (lgef(a)<lgef(b))
  538.     {
  539.       r=b;b=a;a=r;
  540.     }
  541.   while(signe(b)) {r=resii(a,b);a=b;b=r;}
  542.   tetpil=avma;return gerepile(av,tetpil,gabs(a));
  543. }
  544.  
  545. GEN mppgcd(a,b)
  546.      GEN a,b;
  547. {
  548.   GEN t;
  549.   long av=avma,tetpil,st,k,va,vb;
  550.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  551.   if (lgef(a)<lgef(b))
  552.     {
  553.       t=b;b=a;a=t;
  554.     }
  555.   b=absi(b);tetpil=avma;a=absi(a);
  556.   if(signe(b)) {t=resii(a,b);a=b;b=t;} else return(gerepile(av,tetpil,a));
  557.   if(!signe(b)) {avma=tetpil;return a;}
  558.   va=vali(a);vb=vali(b);k=min(va,vb);a=shifti(a,-va);b=shifti(b,-vb);
  559.   t=subii(a,b);st=signe(t);if(st) setsigne(t,1);
  560.   while(st)
  561.     {
  562.       t=shifti(t,-vali(t));
  563.       if(st>0) a=t;else b=t;
  564.       t=subii(a,b);st=signe(t);if(st) setsigne(t,1);
  565.     }
  566.   tetpil=avma;return gerepile(av,tetpil,shifti(a,k));
  567. }
  568.  
  569. /*********************************************************************/
  570. /*********************************************************************/
  571. /**                                                                 **/
  572. /**                     FONCTION BEZOUT                             **/
  573. /**                                                                 **/
  574. /*********************************************************************/
  575. /*********************************************************************/
  576.  
  577. GEN bezout(a,b,u,v)
  578.      GEN a,b,*u,*v;
  579.      
  580. {
  581.   GEN p,u1,v1,d1,d,q,uu,vv,*bof;
  582.   long av,av2,count=0;
  583.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  584.   if (lgef(b)>lgef(a))
  585.     {
  586.       u1=b;b=a;a=u1;bof=u;u=v;v=bof;
  587.     }
  588.   d=gcopy(a);
  589.   if(!signe(b))
  590.     {
  591.       *u=gun;*v=gzero;
  592.     }
  593.   else
  594.     {
  595.       *u=uu=cgeti(lgef(b));
  596.       *v=vv=cgeti(lgef(a));affsi(0,vv);
  597.       av=avma;
  598.       v1=cgeti(lgef(a));affsi(1,v1);
  599.       d1=gcopy(b);
  600.       q=cgeti(lgef(a));
  601.       av2=avma;
  602.       
  603.       while (signe(d1))
  604.     {
  605.       dvmdiiz(d,d1,q,d);
  606.       subiiz(vv,mulii(q,v1),vv);
  607.       p=vv;vv=v1;v1=p;
  608.       p=d;d=d1;d1=p;
  609.       avma=av2;
  610.       count++;
  611.     }
  612.       if (odd(count))
  613.     {
  614.       affii(vv,v1);vv=v1;affii(d,d1);d=d1;
  615.     }
  616.       
  617.       diviiz(subii(d,mulii(b,vv)),a,uu);
  618.       avma=av;
  619.       if (signe(d)== -1)
  620.     {
  621.       setsigne(d,1);setsigne(uu,-signe(uu));setsigne(vv,-signe(vv));
  622.     }
  623.     }
  624.   return d;
  625. }
  626.  
  627. GEN bezout1(a,b,u,v)
  628.      GEN a,b,*u,*v;
  629. {
  630.   long av=avma,tetpil,sa,sb,va,vb,vt,f1,f2,st,i,dec,k;
  631.   GEN u1,v1,t1,v3,t3,q,r,d;
  632.  
  633.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  634.   if (lgef(b)>lgef(a)) {u1=b;b=a;a=u1;f1=1;} else f1=0;
  635.   sb=signe(b);b=absi(b);tetpil=avma;sa=signe(a);a=absi(a);
  636.   if(!sb)
  637.     {
  638.       a=gerepile(av,tetpil,a);
  639.       if(f1) {*u=gzero;*v=(sb>=0)?gun:negi(gun);} 
  640.       else {*u=(sa>=0)?gun:negi(gun);*v=gzero;}
  641.       return a;
  642.     }
  643.   q=dvmdii(a,b,&r);a=b;b=r;
  644.   if(!signe(b)) 
  645.     {
  646.       avma=tetpil;
  647.       if(f1) {*u=(sa>=0)?gun:negi(gun);*v=gzero;}
  648.       else {*u=gzero;*v=(sb>=0)?gun:negi(gun);}
  649.       return a;
  650.     }
  651.   va=vali(a);vb=vali(b);k=min(va,vb);if(k) {a=shifti(a,-k);b=shifti(b,-k);}
  652.   if(mpodd(b)) f2=0;
  653.   else {f2=1;u1=b;b=a;a=u1;}
  654.   u1=gun;d=a;v1=v3=b;
  655.   if(mpodd(a)) {t1=gzero;t3=b;st= -1;}
  656.   else {t1=shifti(addsi(1,b),-1);t3=shifti(a,-1);st=1;}
  657.   while(st)
  658.     {
  659.       vt=vali(t3);
  660.       if(vt)
  661.     {
  662.       t3=shifti(t3,-vt);
  663.       for(i=1;i<=vt;i++)
  664.         t1=mpodd(t1)?shifti(addii(t1,b),-1):shifti(t1,-1);
  665.     }
  666.       if(st>0) {u1=t1;d=t3;} else {v1=subii(b,t1);v3=t3;}
  667.       t1=subii(u1,v1);t3=subii(d,v3);if(signe(t1)<0) t1=addii(t1,b);
  668.       st=signe(t3);if(st) setsigne(t3,1);
  669.     }
  670.   v1=divii(subii(d,mulii(a,u1)),b);if(f2) {t1=u1;u1=v1;v1=t1;}
  671.   u1=subii(u1,mulii(v1,q));if(!f1){t1=u1;u1=v1;v1=t1;}
  672.   tetpil=avma;u1=(sa>=0)?gcopy(u1):negi(u1);v1=(sb>=0)?gcopy(v1):negi(v1);
  673.   d=shifti(d,k);dec=lpile(av,tetpil,0)>>2;u1+=dec;v1+=dec;d+=dec;
  674.   *u=u1;*v=v1;return d;
  675. }
  676.  
  677. /*********************************************************************/
  678. /*********************************************************************/
  679. /**                                                                 **/
  680. /**                     FONCTION CHINOISE                           **/
  681. /**                                                                 **/
  682. /*********************************************************************/
  683. /*********************************************************************/
  684.  
  685. GEN chinois(x,y)
  686.      GEN x,y;
  687. {
  688.   GEN z,p1,p2,p3,p4,d,*u,*v;
  689.   long av,av1,tetpil,dec;
  690.   
  691.   z=cgetg(3,typ(x));
  692.   av=avma;
  693.   d=gbezout(x[1],y[1],&u,&v);
  694.   if(gcmp(gmod(x[2],d),gmod(y[2],d))) err(chiner);
  695.   p1=gdiv(x[1],d);
  696.   p2=gdiv(y[1],d);
  697.   p3=gmul(y[2],gmul(u,p1));
  698.   p4=gmul(x[2],gmul(v,p2));
  699.   p3=gadd(p3,p4);
  700.   tetpil=avma;
  701.   p1=gmul(x[1],p2);p2=gmod(p3,p1);
  702.   av1=avma;dec=lpile(av,tetpil,0)>>2;
  703.   z[1]=adecaler(p1,tetpil,av1)?(long)(p1+dec):(long)p1;
  704.   z[2]=adecaler(p2,tetpil,av1)?(long)(p2+dec):(long)p2;
  705.   return z;
  706. }
  707.  
  708. /*********************************************************************/
  709. /*********************************************************************/
  710. /**                                                                 **/
  711. /**                     FONCTION INVERSE MODULO                     **/
  712. /**                                                                 **/
  713. /**    si a est inversible,on renvoie TRUE et l'inverse dans res   **/
  714. /**    dans le cas contraire,on renvoie FALSE et le pgcd dans res  **/
  715. /**                                                                 **/
  716. /*********************************************************************/
  717. /*********************************************************************/
  718.  
  719. long inversemodulo(a,b,res)
  720.      GEN a,b,*res;
  721.      
  722. {
  723.   GEN u,u1,d1,q;
  724.   long av,av2,count=0;
  725.   if ((typ(a)!=1) || (typ(b)!=1)) err(arither1);
  726.   if (!signe(b))  {*res=mpabs(a);return 0;}
  727.   *res=mpabs(b);
  728.   av=avma;
  729.   u=cgeti(lgef(b));affsi(0,u);
  730.   u1=cgeti(lgef(b));affsi(1,u1);
  731.   d1=cgeti(lgef(b));
  732.   modiiz(a,*res,d1);
  733.   q=cgeti(lgef(b));
  734.   av2=avma;
  735.   
  736.   while (signe(d1))
  737.     {
  738.       dvmdiiz(*res,d1,q,*res);
  739.       subiiz(u,mulii(q,u1),u);
  740.       a=u;u=u1;u1=a;
  741.       a= *res;*res=d1;d1=a;
  742.       avma=av2;
  743.       count++;
  744.     }
  745.   if (odd(count))
  746.     {
  747.       affii(*res,d1);*res=d1;affii(u,u1);u=u1;
  748.     }
  749.   
  750.   if (cmpis(*res,1)) {avma=av;return 0;}
  751.   modiiz(u,b,*res);avma=av;return 1;
  752. }
  753.  
  754. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  755. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  756. /*                                                                 */
  757. /*                      INVERSE MODULO                             */
  758. /*                                                                 */
  759. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  760. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  761.  
  762. GEN     mpinvmod(a,m)
  763.      GEN     a,m;
  764.      
  765. {
  766.   GEN   res;
  767.   
  768.   if (inversemodulo(a,m,&res)) return res;
  769.   err(invmoder,"",gmodulcp(res,m));
  770. }
  771.  
  772.  
  773. /*********************************************************************/
  774. /*********************************************************************/
  775. /**                                                                 **/
  776. /**                     PUISSANCE MODULO                            **/
  777. /**     le resultat occupe autant de place que le module            **/
  778. /**                                                                 **/
  779. /*********************************************************************/
  780. /*********************************************************************/
  781.  
  782. GEN puissmodulo(a,n,m)
  783.      GEN a,n,m;
  784. {
  785.   GEN res,aux;
  786.   long av,*p,man,k,nb;
  787.   
  788.   if ((typ(a)!=1) || (typ(n)!=1) || (typ(m)!=1)) err(arither1);
  789.   res=cgeti(lgef(m));
  790.   av=avma;
  791.   switch (signe(n))
  792.     {
  793.     case -1:
  794.       if (inversemodulo(a,m,&aux))
  795.     {
  796.       affii(puissmodulo(aux,mpabs(n),m),res);
  797.       avma=av;
  798.     }
  799.       else err(puissmoder);
  800.       break;
  801.     case 0: if(divise(a,m)) affsi(0,res); else affsi(1,res);break;
  802.     case 1: modiiz(a,m,res); p=n+2;
  803.       for (man= *p,k=31;man>0;man<<=1) k--;
  804.       man<<=1;/* le premier bit est implicite */
  805.       for (nb=lgef(n)-2;nb;man= *++p,k=32,nb--)
  806.         for(;k;man<<=1,k--)
  807.       {
  808.         modiiz(mulii(res,res),m,res);
  809.         if (man<0) modiiz(mulii(res,a),m,res);
  810.         avma=av;
  811.       }
  812.     }
  813.   return res;
  814. }
  815.  
  816. /*********************************************************************/
  817. /*********************************************************************/
  818. /**                                                                 **/
  819. /**                     PSEUDO PRIMALITE                            **/
  820. /**             n est-il pseudo-premier fort en base a ?            **/
  821. /**                                                                 **/
  822. /*********************************************************************/
  823. /*********************************************************************/
  824.  
  825. long pseudopremier(n,a)
  826.      GEN n,a;
  827. {
  828.   long r,av,av2,result;
  829.   GEN t,z,c,b;
  830.   if ((typ(a)!=1) || (typ(n)!=1)) err(arither1);
  831.   av=avma;
  832.   t=subis(mpabs(n),1);
  833.   if ((r=vali(t))>0)
  834.     {
  835.       b=puissmodulo(a,shifti(t,-r),n);
  836.       if (cmpsi(1,b))
  837.     {
  838.       c=cgeti(lgef(n));
  839.       av2=avma;
  840.       for(;r;r--)
  841.         {
  842.           modiiz(mulii(b,b),n,c);
  843.           z=b;b=c;c=z;
  844.           avma=av2;
  845.           if (!cmpsi(1,b)) break;
  846.         }
  847.       if (r) result=!cmpii(c,t);
  848.       else result=0;
  849.     }
  850.       else result=1;
  851.     }
  852.   else result=!cmpsi(1,t);/* t impair ou nul,tester n=2 ou-2 */
  853.   avma=av;
  854.   return result;
  855. }
  856.  
  857. /*********************************************************************/
  858. /*********************************************************************/
  859. /**                                                                 **/
  860. /**                     MILLER-RABIN                                **/
  861. /**                                                                 **/
  862. /**             tester k fois la primalite de n                     **/
  863. /**                                                                 **/
  864. /*********************************************************************/
  865. /*********************************************************************/
  866.  
  867. millerrabin(n,k)
  868.      
  869.      GEN n;
  870.      long k;
  871. {
  872.   long r,r1,av,av2,result;
  873.   GEN t,t1,z,c,b;
  874.   if ((typ(n)!=1)) err(arither1);
  875.   av=avma;
  876.   t=subis(mpabs(n),1);
  877.   if ((r1=vali(t))>0)
  878.     {
  879.       b=cgeti(lgef(n));
  880.       c=cgeti(lgef(n));
  881.       t1=shifti(t,-r1);
  882.       av2=avma;
  883.       result=1;
  884.       for(;k;k--)
  885.     {
  886.       do z=modii(stoi(rand()),n);
  887.       while(!signe(z));
  888.       affii(puissmodulo(z,t1,n),b);
  889.       if (cmpsi(1,b))
  890.         {
  891.           for(r=r1;r;r--)
  892.         {
  893.           modiiz(mulii(b,b),n,c);
  894.           z=b;b=c;c=z;
  895.           avma=av2;
  896.           if (!cmpsi(1,b)) break;
  897.         }
  898.           if (r) result=!cmpii(c,t);
  899.           else result=0;
  900.         }
  901.       else result=1;
  902.       if (!result) break;
  903.     }
  904.     }
  905.   else result=!cmpsi(1,t);/* t impair ou nul,tester n=2 ou-2 */
  906.   avma=av;
  907.   return result;
  908. }
  909.  
  910. GEN     bigprem(n)
  911.      GEN     n;
  912.      
  913. {
  914.   long    av1,av2,av3;
  915.   
  916.   av1=avma;
  917.   if (typ(n)!=1) err(arither1);
  918.   if(gcmp(n,deux)<=0) return gdeux;
  919.   if(!mpodd(n)) n=addsi(1,n);else n=gcopy(n);
  920.   av3=av2=avma;
  921.   while(!millerrabin(n,10)) {av2=avma;n=addsi(2,n);}
  922.   if (av2!=av3) return gerepile(av1,av2,n);
  923.   return n;
  924. }
  925.  
  926. long isprime(x)
  927.      GEN x;
  928. {
  929.   long av=avma;
  930.   
  931.   if (typ(x)!=1) err(arither1);
  932.   x=absi(x);if(gcmpgs(x,3)<=0) {avma=av;return !gcmp1(x);}
  933.   if(!mpodd(x)) {avma=av;return 0;}
  934.   avma=av;return millerrabin(x,10);
  935. }
  936.  
  937. long ispsp(x)
  938.      GEN x;
  939. {
  940.   long av=avma;
  941.   
  942.   if (typ(x)!=1) err(arither1);
  943.   x=absi(x);if(gcmpgs(x,3)<=0) {avma=av;return !gcmp1(x);}
  944.   if(!mpodd(x)) {avma=av;return 0;}
  945.   avma=av;return millerrabin(x,1);
  946. }
  947.  
  948. long issquarefree(x)
  949.      GEN x;
  950. {
  951.   long av=avma,lx,i;
  952.   GEN p1;
  953.   
  954.   if(gcmp0(x)) return 0;
  955.   if(typ(x)==1) 
  956.     {
  957.       p1=(GEN)(auxdecomp(x,1)[2]);
  958.       lx=lg(p1);
  959.       for(i=1;(i<lx)&&gcmp1(p1[i]);i++);
  960.       avma=av;return (i==lx);
  961.     }
  962.   else
  963.     {
  964.       if(typ(x)!=10) err(issquer1);
  965.       p1=ggcd(x,deriv(x,varn(x)));
  966.       avma=av;return (lgef(p1)==3);
  967.     }
  968. }
  969.  
  970. long isfundamental(x)
  971.      GEN x;
  972. {
  973.   long av,f,r;
  974.   GEN p1;
  975.   
  976.   if(gcmp0(x)) return 0;
  977.   r=x[lgef(x)-1]&3;
  978.   if(r==0) 
  979.     {
  980.       av=avma;p1=shifti(x,-2);r=p1[lgef(p1)-1]&3;
  981.       if(r==0) return 0;
  982.       if(signe(x)<0) r=4-r;
  983.       f=(r==1)?0:issquarefree(p1);avma=av;return f;
  984.     }
  985.   if(signe(x)<0) r=4-r;
  986.   return (r==1)?issquarefree(x):0;
  987. }
  988.  
  989.  
  990. /*********************************************************************/
  991. /*********************************************************************/
  992. /**                                                                 **/
  993. /**                     FONCTION FACTORIELLE                        **/
  994. /**                                                                 **/
  995. /*********************************************************************/
  996. /*********************************************************************/
  997.  
  998. GEN mpfact(n)
  999.      long n;
  1000.      
  1001. {
  1002.   long av = avma, tetpil, limite = (bot + avma) / 2, k;
  1003.   GEN f = gun;
  1004.   if (n < 2) if (n < 0) err(facter); else return f;
  1005.   for (k = 2; k < n; k++)
  1006.     if (avma < limite)
  1007.       {
  1008.     tetpil = avma;
  1009.     f = gerepile(av, tetpil, mulsi(k,f));
  1010.       }
  1011.     else f = mulsi(k,f);
  1012.   tetpil = avma;
  1013.   return gerepile(av, tetpil, mulsi(k,f));
  1014. }
  1015.  
  1016. GEN mpfactr(n,prec)
  1017.      long n,prec;
  1018.      
  1019. {
  1020.   long av = avma, tetpil, limite = (bot + avma) / 2, k;
  1021.   GEN f;
  1022.  
  1023.   affsr(1,f=cgetr(prec));
  1024.   if (n < 2) if (n < 0) err(facter); else return f;
  1025.   for (k = 2; k < n; k++)
  1026.     if (avma < limite)
  1027.       {
  1028.     tetpil = avma;
  1029.     f = gerepile(av, tetpil, mulsr(k,f));
  1030.       }
  1031.     else f = mulsr(k,f);
  1032.   tetpil = avma;
  1033.   return gerepile(av, tetpil, mulsr(k,f));
  1034. }
  1035.  
  1036. /*********************************************************************/
  1037. /*********************************************************************/
  1038. /**                                                                 **/
  1039. /**                     LUCAS ET FIBONACCI                          **/
  1040. /**                                                                 **/
  1041. /*********************************************************************/
  1042. /*********************************************************************/
  1043.  
  1044. void lucas(n,ln,ln1)
  1045.      long n;
  1046.      GEN *ln,*ln1;
  1047.      
  1048. {
  1049.   long taille,av;
  1050.   GEN z,t;
  1051.   if (n)
  1052.     {
  1053.       taille=C3*(1+abs(n))+3;
  1054.       *ln=cgeti(taille);
  1055.       *ln1=cgeti(taille);
  1056.       av=avma;
  1057.       lucas(n/2,&z,&t);
  1058.       switch(n % 4)
  1059.     {
  1060.     case -3:
  1061.       addsiz(2,mulii(z,z),*ln1);
  1062.       subiiz(addsi(1,mulii(z,t)),*ln1,*ln);break;
  1063.     case -2:
  1064.       addsiz(2,mulii(z,z),*ln);addsiz(1,mulii(z,t),*ln1);break;
  1065.     case -1:
  1066.       subisz(mulii(z,z),2,*ln1);
  1067.       subiiz(subis(mulii(z,t),1),*ln1,*ln);break;
  1068.     case  0: subisz(mulii(z,z),2,*ln);subisz(mulii(z,t),1,*ln1);break;
  1069.     case  1: subisz(mulii(z,t),1,*ln);addsiz(2,mulii(t,t),*ln1);break;
  1070.     case  2: addsiz(2,mulii(z,z),*ln);addsiz(1,mulii(z,t),*ln1);break;
  1071.     case  3: addsiz(1,mulii(z,t),*ln);subisz(mulii(t,t),2,*ln1);
  1072.     }
  1073.       avma=av;
  1074.     }
  1075.   else
  1076.     {
  1077.       *ln=stoi(2);
  1078.       *ln1=stoi(1);
  1079.     }
  1080. }
  1081.  
  1082. GEN fibo(n)
  1083.      long n;
  1084.      
  1085. {
  1086.   long taille,av;
  1087.   GEN ln,ln1,f;
  1088.   taille=C3*abs(n)+3;
  1089.   f=cgeti(taille);
  1090.   av=avma;
  1091.   lucas(n-1,&ln,&ln1);
  1092.   divisz(addii(mulsi(2,ln),ln1),5,f);
  1093.   avma=av;
  1094.   return f;
  1095. }
  1096.  
  1097. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1098. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1099. /*                                                                 */
  1100. /*                      FRACTIONS CONTINUES                        */
  1101. /*                                                                 */
  1102. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1103. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1104.  
  1105. GEN  gcf(x)
  1106.      GEN x;
  1107. {
  1108.   return sfcont(x,x,0);
  1109. }
  1110.  
  1111. GEN  gboundcf(x,k)
  1112.      GEN x;
  1113.      long k;
  1114. {
  1115.   return sfcont(x,x,k);
  1116. }
  1117.  
  1118. GEN     sfcont(x,x1,k)
  1119.      GEN     x,x1;
  1120.      long    k;
  1121. {
  1122.   long  av,lx=lg(x),tx=typ(x),e,i,j,l,l1,l2,lx1,tetpil,f,f1;
  1123.   GEN   y,p1,p2,p3,p4,yp;
  1124.   
  1125.   if(tx<10)
  1126.     {
  1127.       if((tx>=6)||(tx==3)) err(sfconter1);
  1128.       if(gcmp0(x)) {y=cgetg(2,17);y[1]=zero;}
  1129.       else switch(tx)
  1130.         {
  1131.         case 1 : y=cgetg(2,17);y[1]=lcopy(x);break;
  1132.         case 2 : l=avma;p1=cgetg(3,5);
  1133.           p2=gcopy(x);settyp(p2,1);setlgef(p2,lx);
  1134.           p1[1]=(long)p2;
  1135.           e=32*(lx-2)-1-expo(x);
  1136.           if(e<0) err(sfconter2);
  1137.           l1=e/32+3;p2=cgeti(l1);
  1138.           p2[1]=0x1000000+l1;
  1139.           p2[2]=(1<<(e&31));
  1140.           for(i=3;i<l1;p2[i]=0,i++);
  1141.           p1[2]=(long)p2;p3=cgetg(3,5);
  1142.           p3[2]=lcopy(p2);
  1143.           p3[1]=laddsi(signe(x),p1[1]);
  1144.           p1=sfcont(p1,p1,k);tetpil=avma;
  1145.           p3=sfcont(p3,p1,k);
  1146.           y=gerepile(l,tetpil,p3);
  1147.           break;
  1148.         case 4 :
  1149.         case 5 : l=avma;l1=46.093443*((lx1=lgef(x[2]))-2)+3;
  1150.       if((k>0)&&(l1>k+1)) l1=k+1;
  1151.           if(lgef(x[1])>=lx1) p1=gcopy(x[1]);
  1152.       else affii(x[1],p1=cgeti(lx1));
  1153.       p2=gcopy(x[2]);l2=avma;lx1=lg(x1);
  1154.           y=cgetg(l1,17);f1=1;f=(x!=x1);i=0;
  1155.           
  1156.           while(f1&&(!gcmp0(p2))&&(i<=l1-2))
  1157.             {
  1158.               i++;y[i]=ldvmdii(p1,p2,&p3);
  1159.               if(signe(p3)<0)
  1160.                 {
  1161.                   p4=addii(p3,p2);affii(p4,p1);
  1162.                   cgiv(p4);cgiv(p3);cgiv(y[1]);
  1163.                   y[1]=laddsi(-1,y[1]);
  1164.                 }
  1165.               else {affii(p3,p1);cgiv(p3);}
  1166.               p4=p1;p1=p2;p2=p4;
  1167.               if(f)
  1168.                 f1=(i<lx1)&&(!cmpii(y[i],x1[i]));
  1169.             }
  1170.           if(!f1)--i;
  1171.           setlg(y,i+1);l2=l2-((l1-i-1)<<2);
  1172.           y=gerepile(l,l2,y);
  1173.         }
  1174.     }
  1175.   else
  1176.     {
  1177.       switch(tx)
  1178.     {
  1179.     case 10: y=cgetg(2,17);y[1]=lcopy(x);break;
  1180.     case 11: av=avma;p1=gtrunc(x);tetpil=avma;
  1181.       y=gerepile(av,tetpil,sfcont(p1,p1,k));break;
  1182.     case 13:
  1183.     case 14:
  1184.       av=avma;l1=lgef(x[1]);if(lgef(x[2])>l1) l1=lgef(x[2]);
  1185.       if((k>0)&&(l1>k+1)) l1=k+1;
  1186.       yp=cgetg(l1,17);p1=(GEN)x[1];i=0;p2=(GEN)x[2];
  1187.       while((!gcmp0(p2))&&(i<=l1-2))
  1188.         {
  1189.           i++;yp[i]=ldivres(p1,p2,&p3);
  1190.           p1=p2;p2=p3;
  1191.         }
  1192.       tetpil=avma;y=cgetg(i+1,17);
  1193.       for(j=1;j<=i;j++) y[j]=lcopy(yp[j]);
  1194.       y=gerepile(av,tetpil,y);break;
  1195.     default: err(sfconter1);
  1196.     }
  1197.     }
  1198.   return y;
  1199. }
  1200.  
  1201. GEN     gcf2(b,x)
  1202.      GEN b,x;
  1203. {
  1204.   long lb,tb=typ(b),i,av,tetpil;
  1205.   GEN y;
  1206.   
  1207.   if(tb<17) err(sfconter1);
  1208.   lb=lg(b);if(lb==1) return cgetg(1,17);
  1209.   if(tb==19)
  1210.     {
  1211.       if(lg(b[1])==1) return sfcont(x,x,0);
  1212.       av=avma;y=cgetg(lb,17);
  1213.       for(i=1;i<lb;i++) y[i]=coeff(b,1,i);
  1214.       tetpil=avma;return gerepile(av,tetpil,sfcont2(y,x));
  1215.     }
  1216.   else return sfcont2(b,x);
  1217. }
  1218.  
  1219. GEN  sfcont2(b,x)
  1220.      GEN b,x;
  1221.      
  1222. {
  1223.   long  av=avma,tx=typ(x),l1=lg(b),i,j,tetpil,f;
  1224.   GEN   y,z,p1;
  1225.   
  1226.   l1=lg(b);y=cgetg(l1,17);
  1227.   if(l1==1) return y;
  1228.   if(tx<10)
  1229.     {
  1230.       if((tx>=6)||(tx==3)) err(sfconter1);
  1231.     }
  1232.   else if(tx==11) x=gtrunc(x);
  1233.   if(!gcmp1(b[1])) x=gmul(b[1],x);
  1234.   y[1]=lfloor(x);p1=gsub(x,y[1]);f=!gcmp0(p1);i=2;
  1235.   for(;(i<l1)&&f;i++)
  1236.     {
  1237.       x=gdiv(b[i],p1);y[i]=lfloor(x);p1=gsub(x,y[i]);f=!gcmp0(p1);
  1238.     }
  1239.   tetpil=avma;z=cgetg(i,17);for(j=1;j<i;j++) z[j]=lcopy(y[j]);
  1240.   return gerepile(av,tetpil,z);
  1241. }
  1242.  
  1243. GEN  pnqn(x)
  1244.      GEN x;
  1245. {
  1246.   long av=avma,tetpil,lx,ly,tx=typ(x),i;
  1247.   GEN y,p0,p1,q0,q1,a,b,p2,q2;
  1248.   
  1249.   if(tx<17) err(pnqner1);
  1250.   lx=lg(x);if(lx==1) return idmat(2);
  1251.   if(tx<19)
  1252.     {
  1253.       p0=q1=gun;q0=gzero;p1=(GEN)x[1];
  1254.       for(i=2;i<lx;i++)
  1255.     {
  1256.       a=(GEN)x[i];
  1257.       p2=gadd(gmul(a,p1),p0);p0=p1;p1=p2;
  1258.       q2=gadd(gmul(a,q1),q0);q0=q1;q1=q2;
  1259.     }
  1260.       tetpil=avma;y=cgetg(3,19);
  1261.       p2=cgetg(3,18);y[1]=(long)p2;p2[1]=lcopy(p1);p2[2]=lcopy(q1);
  1262.       p2=cgetg(3,18);y[2]=(long)p2;p2[1]=lcopy(p0);p2[2]=lcopy(q0);
  1263.       return gerepile(av,tetpil,y);
  1264.     }
  1265.   else
  1266.     {
  1267.       ly=lg(x[1]);if((ly==1)||(ly>3)) err(pnqner2);
  1268.       if(ly==2) 
  1269.     {
  1270.       p1=cgetg(lx,17);for(i=1;i<lx;i++) p1[i]=(long)(((GEN)x[i])[1]);
  1271.       tetpil=avma;return gerepile(av,tetpil,pnqn(p1));
  1272.     }
  1273.       else
  1274.     {
  1275.       p0=gun;q0=gzero;p1=(GEN)coeff(x,2,1);q1=(GEN)coeff(x,1,1);
  1276.       for(i=2;i<lx;i++)
  1277.         {
  1278.           a=(GEN)coeff(x,2,i);b=(GEN)coeff(x,1,i);
  1279.           p2=gadd(gmul(a,p1),gmul(b,p0));p0=p1;p1=p2;
  1280.           q2=gadd(gmul(a,q1),gmul(b,q0));q0=q1;q1=q2;
  1281.         }
  1282.       tetpil=avma;y=cgetg(3,19);
  1283.       p2=cgetg(3,18);y[1]=(long)p2;p2[1]=lcopy(p1);p2[2]=lcopy(q1);
  1284.       p2=cgetg(3,18);y[2]=(long)p2;p2[1]=lcopy(p0);p2[2]=lcopy(q0);
  1285.       return gerepile(av,tetpil,y);
  1286.     }
  1287.     }
  1288. }
  1289.  
  1290.  
  1291.  
  1292.  
  1293.  
  1294.  
  1295.  
  1296. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1297. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1298. /*~                                                                     ~*/
  1299. /*~                  UNITE FONDAMENTALE ET REGULATEUR                   ~*/
  1300. /*~                                                                     ~*/
  1301. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1302. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1303.  
  1304. GEN   fundunit(x)
  1305.      GEN   x;
  1306.      
  1307. {
  1308.   GEN  p1,q1,y,u,v,a,u1,v1,sqd,f,m,c;
  1309.   long av,tetpil,r,p4,lx,flp,flq,av2;
  1310.   
  1311.   if(typ(x)!=1) err(arither1);
  1312.   if (signe(x)<=0) err(funder1);
  1313.   r=(x[lg(x)-1])&3;
  1314.   if((r==2)||(r==3)) err(funder2);
  1315.   p4=lclone(quadpoly(x));
  1316.   av=avma;sqd=racine(x);lx=lgef(sqd)+1;
  1317.   affsi(2,v=cgeti(lx));affsi(r,u=cgeti(lx));
  1318.   affii(shifti(addsi(r,sqd),-1),a=cgeti(lx));
  1319.   m=cgetg(3,19);
  1320.   c=cgetg(3,18);m[1]=(long)c;
  1321.   c[1]=(long)a;c[2]=un;
  1322.   c=cgetg(3,18);m[2]=(long)c;
  1323.   c[1]=un;c[2]=zero;
  1324.   av2=avma;
  1325.   f=cgetg(3,19);
  1326.   c=cgetg(3,18);f[1]=(long)c;
  1327.   c[1]=lcopy(a);c[2]=un;
  1328.   c=cgetg(3,18);f[2]=(long)c;
  1329.   c[1]=un;c[2]=zero;
  1330.   do
  1331.     {
  1332.       u1=subii(mulii(a,v),u);flp=cmpii(u,u1);affii(u1,u);
  1333.       v1=divii(subii(x,mulii(u,u)),v);flq=cmpii(v,v1);affii(v1,v);
  1334.       diviiz(addii(sqd,u),v,a);
  1335.       tetpil=avma;if(flp&&flq) f=gerepile(av2,tetpil,gmul(f,m));
  1336.     }
  1337.   while(flp&&flq);
  1338.   c=(GEN)(f[2]);p1=(GEN)(c[1]);q1=(GEN)(c[2]);
  1339.   y=cgetg(4,8);y[1]=p4;
  1340.   if(r) { y[2]=lsubii(p1,q1);y[3]=lcopy(q1);}
  1341.   else {y[2]=lcopy(p1);y[3]=lcopy(q1);}
  1342.   if(!flq)
  1343.     {
  1344.       f=gmul(f,m);
  1345.       c=(GEN)(f[2]);p1=(GEN)(c[1]);q1=(GEN)(c[2]);
  1346.       v1=cgetg(4,8);v1[1]=p4;
  1347.       if(r) { v1[2]=lsubii(p1,q1);v1[3]=lcopy(q1);}
  1348.       else {v1[2]=lcopy(p1);v1[3]=lcopy(q1);}
  1349.       u1=gconj(y);tetpil=avma;y=gdiv(v1,u1);
  1350.     }
  1351.   else
  1352.     {
  1353.       u1=gconj(y);tetpil=avma;y=gdiv(y,u1);
  1354.     }
  1355.   if(signe(y[3])<0) {tetpil=avma;y=gneg(y);}
  1356.   return gerepile(av,tetpil,y);
  1357. }
  1358.  
  1359. GEN   regula(x,prec)
  1360.      GEN   x;
  1361.      long  prec;
  1362.      
  1363. {
  1364.   GEN  reg,reg1,rsqd;
  1365.   GEN  rexp,y,u,v,a,u1,v1,sqd;
  1366.   long av,tetpil,r,lx,flp,flq,av2;
  1367.   
  1368.   if(typ(x)!=1) err(arither1);
  1369.   if (signe(x)<=0) err(funder1);
  1370.   r=(x[lg(x)-1])&3;
  1371.   if((r==2)||(r==3)) err(funder2);
  1372.   av=avma;sqd=racine(x);if(gegal(mulii(sqd,sqd),x)) err(reguler1);
  1373.   rsqd=gsqrt(x,prec);affsi(0,rexp=cgeti(4));
  1374.   lx=lgef(sqd)+1;affsr(2,reg=cgetr(prec));
  1375.   affsi(2,v=cgeti(lx));affsi(r,u=cgeti(lx));
  1376.   affii(shifti(addsi(r,sqd),-1),a=cgeti(lx));
  1377.   av2=avma;
  1378.   do
  1379.     {
  1380.       u1=subii(mulii(a,v),u);flp=cmpii(u,u1);affii(u1,u);
  1381.       reg1=divri(addir(u,rsqd),v);
  1382.       v1=divii(subii(x,mulii(u,u)),v);flq=cmpii(v,v1);
  1383.       if(flp&&flq) {affii(v1,v);diviiz(addii(sqd,u),v,a);}
  1384.       tetpil=avma;
  1385.       if(flp&&flq)
  1386.         {
  1387.           reg=gerepile(av2,tetpil,mulrr(reg,reg1));
  1388.           r=expo(reg);addsiz(r,rexp,rexp);setexpo(reg,0);
  1389.         }
  1390.     }
  1391.   while(flp&&flq);
  1392.   reg=shiftr(mulrr(reg,reg),-1);
  1393.   if(!flq)
  1394.     {
  1395.       u1=subii(mulii(a,v),u);reg1=divri(addir(u,rsqd),v);
  1396.       reg=divri(mulrr(reg,reg1),v);
  1397.     }
  1398.   else reg=divri(reg,v);
  1399.   y=glog(reg,prec);u1=mpshift(mulir(rexp,glog(gdeux,prec)),1);
  1400.   tetpil=avma;return gerepile(av,tetpil,gadd(y,u1));
  1401. }
  1402.  
  1403.  
  1404. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1405. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1406. /*~                                                                     ~*/
  1407. /*~                         NOMBRE DE CLASSES                           ~*/
  1408. /*~                                                                     ~*/
  1409. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1410. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1411.  
  1412. GEN classno2(x)
  1413.      GEN x;
  1414.      
  1415. {
  1416.   long av=avma,tetpil,n,i,k,s=signe(x),fl2,ex;
  1417.   GEN p1,p2,p3,p4,p5,p6,p7,p9,pi4,d,reg,logd;
  1418.   
  1419.   if (typ(x)!=1) err(arither1);
  1420.   if (!s) err(arither2);
  1421.   p1=auxdecomp(absi(x),1);p2=(GEN)p1[2];p1=(GEN)p1[1];n=lg(p1);p3=gun;fl2=0;
  1422.   for(i=1;i<n;i++)
  1423.     {
  1424.       ex=itos(p2[i]);
  1425.       if(ex>=2)
  1426.     {
  1427.       p4=(GEN)p1[i];
  1428.       if(ex&1) p3=mulii(p3,p4);
  1429.       if(gegal(p4,gdeux)) fl2=1;
  1430.     }
  1431.       else p3=mulii(p3,p1[i]);
  1432.     }
  1433.   if((p3[lgef(p3)-1]&3)!=(2-s))
  1434.     {
  1435.       if(fl2) p3=shifti(p3,2);
  1436.       else err(classer2);
  1437.     }
  1438.   else fl2=0;
  1439.   p9=gun;x=(s<0) ? gneg(p3) : p3;
  1440.   for(i=1;i<n;i++)
  1441.     {
  1442.       ex=itos(p2[i]);p4=(GEN)p1[i];
  1443.       if(gegal(p4,deux)&&fl2) ex-=2;
  1444.       if(ex>=2)
  1445.     {
  1446.       p9=mulii(p9,subis(p4,kronecker(x,p4)));
  1447.       if(ex>=4) p9=mulii(p9,gpuigs(p4,(ex>>1)-1));
  1448.     }
  1449.     }
  1450.   pi4=mppi(4);
  1451.   if(s>0)
  1452.     {
  1453.       reg=regula(x,4);logd=glog(x,4);
  1454.       p1=gsqrt(gdiv(gmul(x,logd),gmul2n(pi4,1)),4);
  1455.       p2=gsubsg(1,gmul2n(gdiv(glog(reg,4),logd),1));
  1456.       p3=gsqrt(gdivsg(2,logd),4);
  1457.       if(gcmp(p2,p3)>=0) p1=gmul(p2,p1);
  1458.       p1=gtrunc(p1);
  1459.       if(lgef(p1)!=3) err(classer1);
  1460.       n=p1[2];if(n<0) err(classer1);
  1461.       p1=gsqrt(x,4);p4=divri(pi4,x);p3=gzero;
  1462.       p7=divsr(1,mpsqrt(pi4));
  1463.       for(i=1;i<=n;i++)
  1464.         {
  1465.           k=krogs(x,i);
  1466.           if(k)
  1467.             {
  1468.               p6=mulir(mulss(i,i),p4);p5=subsr(1,mulrr(p7,incgam3(ghalf,p6,4)));
  1469.               p5=addrr(divrs(mulrr(p1,p5),i),eint1(p6,4));
  1470.               if(k>0) p3=addrr(p3,p5);else p3=subrr(p3,p5);
  1471.             }
  1472.         }
  1473.       p3=shiftr(divrr(p3,reg),-1);
  1474.     }
  1475.   else
  1476.     {
  1477.       d=p3;
  1478.       if(gcmpgs(x,-11)>=0) {tetpil=avma;return gerepile(av,tetpil,gcopy(p9));}
  1479.       p1=gtrunc(gsqrt(gdiv(gmul(d,glog(d,4)),gmul2n(pi4,1)),4));
  1480.       if(lgef(p1)!=3) err(classer1);
  1481.       n=p1[2];if(n<0) err(classer1);
  1482.       p1=gdiv(gsqrt(d,4),pi4);p4=divri(pi4,d);p3=gzero;
  1483.       p7=divsr(1,mpsqrt(pi4));
  1484.       for(i=1;i<=n;i++)
  1485.         {
  1486.           k=krogs(x,i);
  1487.           if(k)
  1488.             {
  1489.               p6=mulir(mulss(i,i),p4);p5=subsr(1,mulrr(p7,incgam3(ghalf,p6,4)));
  1490.               p5=addrr(p5,divrr(divrs(p1,i),mpexp(mulir(mulss(i,i),p4))));
  1491.               if(k>0) p3=addrr(p3,p5);else p3=subrr(p3,p5);
  1492.             }
  1493.         }
  1494.     }
  1495.   p3=ground(p3);tetpil=avma;return gerepile(av,tetpil,gmul(p9,p3));
  1496. }
  1497.  
  1498. GEN classno3(x)
  1499.      GEN x;
  1500.      
  1501. {
  1502.   long d,a,b,h,b2,f,av,tetpil;
  1503.   GEN y;
  1504.   
  1505.   d= -itos(x);if((d>0)||((d&3)>1)) return gzero;
  1506.   if(!d) return gdivgs(un,-12);
  1507.   h=0;b=d&1;b2=(b-d)>>2;f=0;
  1508.   if(!b)
  1509.     {
  1510.       for(a=1;a*a<b2;a++) {if(!(b2%a)) h++;}
  1511.       f=(a*a==b2);b=2;b2=(4-d)>>2;
  1512.     }
  1513.   while(b2*3+d<0)
  1514.     {
  1515.       if(!(b2%b)) h++;
  1516.       for(a=b+1;a*a<b2;a++) {if(!(b2%a)) h+=2;}
  1517.       if(a*a==b2) h++;
  1518.       b+=2;b2=(b*b-d)>>2;
  1519.     }
  1520.   if(b2*3+d==0) 
  1521.     {
  1522.       av=avma;y=gdivgs(gun,3);tetpil=avma;
  1523.       return gerepile(av,tetpil,gaddsg(h,y));
  1524.     }
  1525.   if(f) return gaddsg(h,ghalf);else return stoi(h);
  1526. }
  1527.