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