home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / arith2.c < prev    next >
C/C++ Source or Header  |  1992-05-20  |  34KB  |  1,271 lines
  1. /*********************************************************************/
  2. /*********************************************************************/
  3. /**                                                                 **/
  4. /**                    FONCTIONS ARITHMETIQUES                      **/
  5. /**                                                                 **/
  6. /**                       (deuxieme partie)                         **/
  7. /**                                                                 **/
  8. /**                      copyright Babe Cool                        **/
  9. /**                                                                 **/
  10. /**                                                                 **/
  11. /*********************************************************************/
  12. /*********************************************************************/
  13.  
  14. #include  "genpari.h"
  15.  
  16. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  17. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  18. /*~                                                                     ~*/
  19. /*~                     NOMBRES PREMIERS                                ~*/
  20. /*~                                                                     ~*/
  21. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  22. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  23.  
  24.  
  25. GEN prime (n)
  26.      long n;
  27.      
  28. {
  29.   byteptr p=diffptr;
  30.   long c,prime=0;
  31.   while (n--) if (c= *p++) prime += c; else err(primer1);
  32.   return stoi(prime);
  33. }
  34.  
  35. GEN primes (n)
  36.      long n;
  37.      
  38. {
  39.   GEN y,z;
  40.   byteptr p=diffptr;
  41.   long c,prime=0;
  42.   z=y=cgetg(n+1,17);
  43.   while (n--)
  44.     {
  45.       if (c= *p++) prime+=c; else err(primer1);
  46.       *++z=(long)stoi(prime);
  47.     }
  48.   return y;
  49. }
  50.  
  51. byteptr initprimes(maxnum)
  52.      long maxnum;
  53. {
  54.   register long k,size=(maxnum+257)>>1;
  55.   byteptr p=(byteptr)calloc(size,1);
  56.   register byteptr q,r,s,fin=p+size;
  57.   for(r=q=p,k=1;r<fin;)
  58.     {
  59.       do {r+=k; k+=2; r+=k;} while (*++q);
  60.       for(s=r;s<fin;s+=k) *s=1;
  61.     }
  62.   r=p; *r++=2; *r++=1;
  63.   for(s=q=r-1;;)
  64.     {
  65.       while (*++q);
  66.       if (q>=fin) break;
  67.       *r++=(q-s)<<1;
  68.       s=q;
  69.     }
  70.   *r++=0;
  71.   return (byteptr)realloc(p,r-p);
  72. }
  73.  
  74. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  75. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  76. /*~                                                                     ~*/
  77. /*~                   FONCTIONS ARITHMETIQUES DE BASE                   ~*/
  78. /*~                                                                     ~*/
  79. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  80. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  81.  
  82. /* Attention, le parametre doit etre une variable. */
  83. #define pseudoprime(p) millerrabin(p,5*lgef(p))
  84.  
  85. GEN gmu(n)
  86.      GEN n;
  87. {
  88.   long l,tx,i;GEN y;
  89.  
  90.   if((tx=typ(n))>=17)
  91.     {
  92.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)gmu(n[i]);
  93.       return y;
  94.     }
  95.   if(tx!=1) err(arither1);
  96.   return stoi(mu(n));
  97. }
  98.  
  99. long mu(n)
  100.      GEN n;
  101.      
  102. {
  103.   byteptr d=diffptr;
  104.   unsigned char c;
  105.   GEN p, f;
  106.   long av=avma,s=1;
  107.   if (typ(n)!=1) err(arither1);
  108.   if (!signe(n)) err(arither2);
  109.   if ((n[2]==1)&&(lgef(n)==3)) return 1;
  110.   n=absi(n);
  111.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  112.   while (c= *d++)
  113.     {
  114.       p[2]+=c;
  115.       if (!mpdivis(n,p,n)) continue;
  116.       if (divise(n,p)) {avma=av;return 0;}
  117.       s= -s;
  118.       if ((n[2]==1)&&(lgef(n)==3)) break;
  119.     }
  120.   while ((n[2]!=1)||(lgef(n)!=3))
  121.     {
  122.       f=gcopy(n);
  123.       while (!((cmpii(mulii(p,p),f)>0)||pseudoprime(f))) f=ellfacteur(f);
  124.       diviiz(n,f,n);
  125.       if (divise(n,f)) {avma=av;return 0;}
  126.       s= -s;
  127.     }
  128.   avma=av;
  129.   return s;
  130. }
  131.  
  132. GEN phi(n)
  133.      GEN n;
  134.      
  135. {
  136.   byteptr d=diffptr;
  137.   unsigned char c;
  138.   GEN f,p,m;
  139.   long av1,av2,tx,i,l;
  140.   
  141.   if((tx=typ(n))>=17)
  142.     {
  143.       l=lg(n);p=cgetg(l,tx);for(i=1;i<l;i++) p[i]=(long)phi(n[i]);
  144.       return p;
  145.     }
  146.   if (tx!=1) err(arither1);
  147.   if (!signe(n)) err(arither2);
  148.   if ((n[2]==1)&&(lgef(n)==3)) return stoi(1);
  149.   m=cgeti(lgef(n));affsi(1,m);
  150.   av1=avma;
  151.   n=absi(n);
  152.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  153.   av2=avma;
  154.   while (c= *d++)
  155.     {
  156.       p[2]+=c;
  157.       if (!mpdivis(n,p,n)) continue;
  158.       muliiz(m,subis(p,1),m);
  159.       while (mpdivis(n,p,n)) muliiz(m,p,m);
  160.       if ((n[2]==1)&&(lgef(n)==3)) break;
  161.       avma=av2;
  162.     }
  163.   while ((n[2]!=1) || (lgef(n)!=3))
  164.     {
  165.       f=gcopy(n);
  166.       while (!((cmpii(mulii(p,p),f)>0) || pseudoprime(f))) f = ellfacteur(f);
  167.       mpdivis(n,f,n);
  168.       muliiz(m,subis(f,1),m);
  169.       while (mpdivis(n,f,n)) muliiz(m,f,m);
  170.       avma=av2;
  171.     }
  172.   avma=av1;
  173.   return m;
  174. }
  175.  
  176. GEN auxdecomp(n,all)
  177.      GEN n; long all;
  178.      
  179. {
  180.   byteptr d=diffptr;
  181.   unsigned char c;
  182.   GEN p,z,p1,p2,f;
  183.   long nb=0,j,k,lim,av1,av2,av3,av4,sn;
  184.   
  185.   if (typ(n)!=1) err(arither1);
  186.   if (!(sn=signe(n))) err(arither2);
  187.   if(sn<0) {stoi(-1);stoi(1);nb++;}
  188.   if ((n[2]!=1) || (lgef(n)!=3))
  189.     {
  190.       av1=avma;
  191.       n=absi(n);
  192.       p=cgeti(3);p[1]=0x1000003;p[2]=0;
  193.       av2=avma;lim=(all<=1)?VERYBIGINT:all;
  194.       while ((c= *d++)&&(p[2]<=lim))
  195.         {
  196.           p[2]+=c;
  197.           if (!mpdivis(n,p,n)) continue;
  198.           nb++;
  199.           for (k=1;mpdivis(n,p,n);k++);
  200.           gcopy(p);
  201.           stoi(k);
  202.           if ((n[2]==1)&&(lgef(n)==3)) break;
  203.         }
  204.       while ((n[2]!=1)||(lgef(n)!=3))
  205.         {
  206.           av3=avma;
  207.           f=n;
  208.       if(all==1)
  209.         while (!((cmpii(mulii(p,p),f)>0)||pseudoprime(f))) f=ellfacteur(f);
  210.           av4=avma;
  211.           f=gerepile(av3,av4,gcopy(f));
  212.           nb++;
  213.           for (k=0;mpdivis(n,f,n);k++);
  214.           stoi(k);
  215.         }
  216.       gerepile(av1,av2,0);
  217.     }
  218.   p=(GEN)avma;
  219.   z=cgetg(3,19);
  220.   p1=cgetg(nb+1,18);z[1]=(long)p1;
  221.   p2=cgetg(nb+1,18);z[2]=(long)p2;
  222.   for (j=nb;j;j--)
  223.     {
  224.       p2[j]=(long)p;
  225.       p+=lg(p);
  226.       p1[j]=(long)p;
  227.       p+=lg(p);
  228.     }
  229.   return z;
  230. }
  231.  
  232. GEN decomp(n)
  233.      GEN n;
  234. {
  235.   return auxdecomp(n,1);
  236. }
  237.  
  238. GEN smallfact(n)
  239.      GEN n;
  240. {
  241.   GEN p1,p2,p3,p4,p5,y;
  242.   long av,tetpil,tx,i,l;
  243.  
  244.   if((tx=typ(n))>=17)
  245.     {
  246.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)smallfact(n[i]);
  247.       return y;
  248.     }
  249.   switch(tx)
  250.     {
  251.     case 1:return auxdecomp(n,0);
  252.     case 5:av=avma;n=gred(n);
  253.     case 4:if(typ(n)==4) av=avma;p1=auxdecomp(n[1],0);
  254.       p2=auxdecomp(n[2],0);p4=concat(p1[1],p2[1]);
  255.       p5=concat(p1[2],gneg(p2[2]));p3=indexsort(p4);
  256.       tetpil=avma;y=cgetg(3,19);y[1]=(long)extract(p4,p3);
  257.       y[2]=(long)extract(p5,p3);return gerepile(av,tetpil,y);
  258.     default: err(arither1);
  259.     }
  260. }
  261.  
  262. GEN boundfact(n,lim)
  263.      GEN n;
  264.      long lim;
  265. {
  266.   GEN p1,p2,p3,p4,p5,y;
  267.   long av,tetpil,tx,i,l;
  268.  
  269.   if(lim<=1) lim=0;
  270.   if((tx=typ(n))>=17)
  271.     {
  272.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)boundfact(n[i],lim);
  273.       return y;
  274.     }
  275.   switch(tx)
  276.     {
  277.     case 1:return auxdecomp(n,lim);
  278.     case 5:av=avma;n=gred(n);
  279.     case 4:if(typ(n)==4) av=avma;p1=auxdecomp(n[1],lim);
  280.       p2=auxdecomp(n[2],lim);p4=concat(p1[1],p2[1]);
  281.       p5=concat(p1[2],gneg(p2[2]));p3=indexsort(p4);
  282.       tetpil=avma;y=cgetg(3,19);y[1]=(long)extract(p4,p3);
  283.       y[2]=(long)extract(p5,p3);return gerepile(av,tetpil,y);
  284.     default: err(arither1);
  285.     }
  286. }
  287.  
  288. GEN divisors(n)
  289.      GEN n;
  290. {
  291.   long tetpil,av=avma,i,j,l,start;
  292.   GEN p,t=decomp(n),p1=(GEN)t[1],p2=(GEN)t[2],t1,t2,t3,nbdiv=gun;
  293.   l=lg(p1);
  294.   start=1+((l>1)&&(signe(p1[1])<0));
  295.   for(i=start;i<l;i++)nbdiv=mulis(nbdiv,1+itos(p2[i]));
  296.   p=t=cgetg(itos(nbdiv)+1,17);
  297.   *++p=lstoi(1);
  298.   for(i=start;i<l;i++)
  299.     for(t1=t,j=itos(p2[i]);j;j--)
  300.       {
  301.     t2=p;
  302.     for(t3=t1;t3<t2;)
  303.       *++p=lmulii(*++t3,p1[i]);
  304.     t1=t2;
  305.       }
  306.   tetpil=avma;
  307.   return gerepile(av,tetpil,sort(t));
  308. }
  309.  
  310.  
  311. GEN gomega(n)
  312.      GEN n;
  313. {
  314.   long l,tx,i;GEN y;
  315.  
  316.   if((tx=typ(n))>=17)
  317.     {
  318.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)gomega(n[i]);
  319.       return y;
  320.     }
  321.   if(tx!=1) err(arither1);
  322.   return stoi(omega(n));
  323. }
  324.  
  325. long omega(n)
  326.      GEN n;
  327.      
  328. {
  329.   byteptr d=diffptr;
  330.   unsigned char c;
  331.   GEN p,f;
  332.   long nb=0,av=avma,av2;
  333.   if (typ(n)!=1) err(arither1);
  334.   if (!signe(n)) err(arither2);
  335.   if ((n[2]==1)&&(lgef(n)==3)) return 0;
  336.   n=absi(n);
  337.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  338.   av2=avma;
  339.   while (c= *d++)
  340.     {
  341.       p[2]+=c;
  342.       if (!mpdivis(n,p,n)) continue;
  343.       nb++;
  344.       while (mpdivis(n,p,n));
  345.       if ((n[2]==1)&&(lgef(n)==3)) break;
  346.     }
  347.   while ((n[2]!=1) || (lgef(n)!=3))
  348.     {
  349.       f=gcopy(n);
  350.       while (!((cmpii(mulii(p,p),f)>0)||pseudoprime(f))) f=ellfacteur(f);
  351.       nb++;
  352.       while (mpdivis(n,f,n));
  353.       avma=av2;
  354.     }
  355.   avma=av;
  356.   return nb;
  357. }
  358.  
  359.  
  360. GEN gbigomega(n)
  361.      GEN n;
  362. {
  363.   long l,tx,i;GEN y;
  364.  
  365.   if((tx=typ(n))>=17)
  366.     {
  367.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)gbigomega(n[i]);
  368.       return y;
  369.     }
  370.   if(tx!=1) err(arither1);
  371.   return stoi(bigomega(n));
  372. }
  373.  
  374. long bigomega(n)
  375.      GEN n;
  376.      
  377. {
  378.   byteptr d=diffptr;
  379.   unsigned char c;
  380.   GEN p,f;
  381.   long nb=0,av=avma,av2;
  382.   if (typ(n)!=1) err(arither1);
  383.   if (!signe(n)) err(arither2);
  384.   if ((n[2]==1)&&(lgef(n)==3)) return 0;
  385.   n=absi(n);
  386.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  387.   av2=avma;
  388.   while (c= *d++)
  389.     {
  390.       p[2]+=c;
  391.       if (!mpdivis(n,p,n)) continue;
  392.       do nb++;while (mpdivis(n,p,n));
  393.       if ((n[2]==1)&&(lgef(n)==3)) break;
  394.     }
  395.   while ((n[2]!=1) || (lgef(n)!=3))
  396.     {
  397.       f=gcopy(n);
  398.       while (!((cmpii(mulii(p,p),f)>0) || pseudoprime(f))) f = ellfacteur(f);
  399.       while (mpdivis(n,f,n)) nb++;
  400.       avma=av2;
  401.     }
  402.   avma=av;
  403.   return nb;
  404. }
  405.  
  406. GEN numbdiv(n)
  407.      GEN n;
  408.      
  409. {
  410.   byteptr d=diffptr;
  411.   unsigned char c;
  412.   GEN p,f,m,m1,y;
  413.   long l,av=avma,av2,av3,i,tx;
  414.  
  415.   if((tx=typ(n))>=17)
  416.     {
  417.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)numbdiv(n[i]);
  418.       return y;
  419.     }
  420.   if (tx!=1) err(arither1);
  421.   if (!signe(n)) err(arither2);
  422.   if ((n[2]==1)&&(lgef(n)==3)) return stoi(1);
  423.   n=absi(n);
  424.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  425.   av2=avma;
  426.   m=stoi(1);
  427.   while (c= *d++)
  428.     {
  429.       p[2]+=c;
  430.       if (!mpdivis(n,p,n)) continue;
  431.       for (l=2;mpdivis(n,p,n);l++);
  432.       av3=avma;
  433.       m1=mulsi(l,m);
  434.       if (lgef (m1) ==lgef(m)) affii(m1,m);
  435.       else m=gerepile(av2,av3,m1);
  436.       if ((n[2]==1)&&(lgef(n)==3)) break;
  437.     }
  438.   while ((n[2]!=1) || (lgef(n)!=3))
  439.     {
  440.       f=gcopy(n);
  441.       while (!((cmpii(mulii(p,p),f)>0) || pseudoprime(f))) f = ellfacteur(f);
  442.       for (l=2;mpdivis(n,f,n);l++);
  443.       av3=avma;
  444.       m=gerepile(av2,av3,mulsi(l,m));
  445.     }
  446.   return gerepile(av,av2,m);
  447. }
  448.  
  449. GEN sumdiv(n)
  450.      GEN n;
  451.      
  452. {
  453.   byteptr d=diffptr;
  454.   unsigned char c;
  455.   GEN p,f,m,m1,y;
  456.   long av1=avma,av2,av3,i,tx,l;
  457.  
  458.   if((tx=typ(n))>=17)
  459.     {
  460.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)sumdiv(n[i]);
  461.       return y;
  462.     }
  463.   if (tx!=1) err(arither1);
  464.   if (!signe(n)) err(arither2);
  465.   if ((n[2]==1)&&(lgef(n)==3)) return stoi(1);
  466.   n=absi(n);
  467.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  468.   av2=avma;
  469.   m=gun;
  470.   while (c= *d++)
  471.     {
  472.       p[2]+=c;
  473.       if (!mpdivis(n,p,n)) continue;
  474.       m1=addsi(1,p);
  475.       while (mpdivis(n,p,n)) m1=addsi(1,mulii(p,m1));
  476.       av3=avma;m=gerepile(av2,av3,mulii(m1,m));
  477.       if ((n[2]==1)&&(lgef(n)==3)) break;
  478.     }
  479.   while ((n[2]!=1)||(lgef(n)!=3))
  480.     {
  481.       f=gcopy(n);
  482.       while (!((cmpii(mulii(p,p),f)>0)||pseudoprime(f))) f=ellfacteur(f);
  483.       m1=gun;
  484.       while (mpdivis(n,f,n)) m1=addsi(1,mulii(f,m1));
  485.       av3=avma;m=gerepile(av2,av3,mulii(m1,m));
  486.     }
  487.   return gerepile(av1,av2,m);
  488. }
  489.  
  490. GEN sumdivk(k,n)
  491.      long k;
  492.      GEN n;
  493.      
  494. {
  495.   byteptr d=diffptr;
  496.   unsigned char c;
  497.   GEN p,p1,f,m,m1,pk,y;
  498.   long av1=avma,av2,av3,k1,tx,i,l;
  499.  
  500.   if((tx=typ(n))>=17)
  501.     {
  502.       l=lg(n);y=cgetg(l,tx);for(i=1;i<l;i++) y[i]=(long)sumdivk(k,n[i]);
  503.       return y;
  504.     }
  505.   if (tx!=1) err(arither1);
  506.   if (!signe(n)) err(arither2);
  507.   if ((n[2]==1)&&(lgef(n)==3)) return stoi(1);
  508.   n=absi(n);k1=k;if(k<0) {k= -k;p1=gpuigs(n,k1);}
  509.   p=cgeti(3);p[1]=0x1000003;p[2]=0;
  510.   av2=avma;
  511.   m=gun;
  512.   while (c= *d++)
  513.     {
  514.       p[2]+=c;
  515.       if (!mpdivis(n,p,n)) continue;
  516.       pk=gpuigs(p,k);
  517.       m1=addsi(1,pk);
  518.       while (mpdivis(n,p,n)) m1=addsi(1,mulii(pk,m1));
  519.       av3=avma;m=gerepile(av2,av3,mulii(m1,m));
  520.       if ((n[2]==1)&&(lgef(n)==3)) break;
  521.     }
  522.   while ((n[2]!=1)||(lgef(n)!=3))
  523.     {
  524.       f=gcopy(n);
  525.       while (!((cmpii(mulii(p,p),f)>0)||pseudoprime(f))) f=ellfacteur(f);
  526.       pk=gpuigs(f,k);
  527.       m1=gun;
  528.       while (mpdivis(n,f,n)) m1=addsi(1,mulii(pk,m1));
  529.       av3=avma;m=gerepile(av2,av3,mulii(m1,m));
  530.     }
  531.   av3=avma;return (k1>=0) ? gerepile(av1,av2,m) : gerepile(av1,av3,gmul(p1,m));
  532. }
  533.  
  534. /* decomposition de nombres par la methode des courbes elliptiques.
  535.    La seule fonction exportee est ellfacteur.
  536.    Cette fonction renvoie un facteur non trivial de n.
  537.    On suppose en entree que n n'est pas premier, et n'est divisible par
  538.    aucun petit nombre premier (pas de facteur<65536,en tout cas.)        */
  539.  
  540. static GEN x1,y11,x2,y2,aux,w,n,global;
  541. static long nombre,taillef;
  542.  
  543. static GEN cree(nb)
  544.      long nb;
  545. {
  546.   GEN z=cgetg(nb+1,17);
  547.   long i;
  548.   for(i=1;i<=nb;i++) z[i]=lgeti(taillef);
  549.   return z;
  550. }
  551.  
  552. static long pur(x,p)
  553.      long x,*p;
  554. {
  555.   byteptr d=diffptr;
  556.   long m=0;
  557.   do m+= *d++;while (x % m);
  558.   do x /=m;while (!(x % m));
  559.   *p=m;
  560.   return x==1;
  561. }
  562.  
  563. static GEN elladd() 
  564. {
  565.   GEN v1,glob,lambda;
  566.   long i,av=avma;
  567.   for(i=1;i<=nombre;i++)
  568.     {subiiz(x1[i],x2[i],aux[i]); if (signe(aux[i])<0) addiiz(aux[i],n,aux[i]);}
  569.   for(i=1;i<=nombre;i++)
  570.     {modiiz(mulii(aux[i],w[i]),n,w[i+1]);avma=av;}
  571.   if (!inversemodulo(w[nombre+1],n,&glob)) return glob;
  572.   affii(glob,global);
  573.   for(i=nombre;i;i--)
  574.     {
  575.       v1=modii(mulii(global,w[i]),n);
  576.       modiiz(mulii(global,aux[i]),n,global);
  577.       lambda=modii(mulii(subii(y11[i],y2[i]),v1),n);
  578.       modiiz(subii(mulii(lambda,lambda),addii(x2[i],x1[i])),n,x2[i]);
  579.       modiiz(subii(mulii(lambda,subii(x1[i],x2[i])),y11[i]),n,y2[i]);
  580.       avma=av;
  581.     }
  582.   return (GEN)0;
  583. }
  584.  
  585. static GEN elldouble()
  586. {
  587.   GEN v1,v2,glob,lambda;
  588.   long i,av=avma;
  589.   for(i=1;i<=nombre;i++) {modiiz(shifti(y2[i],1),n,aux[i]);avma=av;}
  590.   for(i=1;i<=nombre;i++) {modiiz(mulii(aux[i],w[i]),n,w[i+1]);avma=av;}
  591.   if (!inversemodulo(w[nombre+1],n,&glob)) return glob;
  592.   affii(glob,global);
  593.   for(i=nombre;i;i--)
  594.     {
  595.       v1=modii(mulii(global,w[i]),n);
  596.       modiiz(mulii(global,aux[i]),n,global);
  597.       lambda=modii(mulii(addsi(1,mulsi(3,mulii(x2[i],x2[i]))),v1),n);
  598.       v2=modii(subii(mulii(lambda,lambda),shifti(x2[i],1)),n);
  599.       modiiz(subii(mulii(lambda,subii(x2[i],v2)),y2[i]),n,y2[i]);
  600.       affii(v2,x2[i]);
  601.       avma=av;
  602.     }
  603.   return (GEN)0;
  604. }
  605.  
  606. static GEN ellmult(k) 
  607.      long k;
  608. {
  609.   GEN result = (GEN)0;
  610.   if (k>1)
  611.     {
  612.       if (result = ellmult(k/2)) return result;
  613.       if (result = elldouble()) return result;
  614.       if (k&1) result = elladd();
  615.     }
  616.   return result;
  617. }
  618.  
  619. GEN ellfacteur(n1)
  620.      GEN n1;
  621. {
  622.   static long i,k,m,av;
  623.   static GEN glob;
  624.   taillef=lgef(n1);
  625.   nombre=10*lgef(n1);
  626.   global=cgeti(taillef);
  627.   av=avma;
  628.   n=absi(n1);
  629.   x1=cree(nombre);
  630.   x2=cree(nombre);
  631.   y11=cree(nombre);
  632.   y2=cree(nombre);
  633.   aux=cree(nombre);
  634.   w=cree(nombre+1);
  635.   w[1]=un;
  636.   for(i=nombre;i;i--) {affsi(rand(),x2[i]);affsi(rand(),y2[i]);}
  637.   for (m=2;;m++)
  638.     if (pur(m,&k))
  639.       {
  640.         for(i=nombre;i;i--) {affii(x2[i],x1[i]);affii(y2[i],y11[i]);}
  641.         if(glob = ellmult(k))
  642.           {
  643.             if (cmpii(mulii(glob,glob),n)>0) diviiz(n,glob,global);
  644.             else affii(glob,global);
  645.             avma=av;
  646.             return global;
  647.           }
  648.       }
  649. }
  650.  
  651. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  652. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  653. /*~                                                                     ~*/
  654. /*~                   COMPOSITION DES FORMES QUADRATIQUES               ~*/
  655. /*~                                                                     ~*/
  656. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  657. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  658.  
  659. GEN  qfi(x, y, z)
  660.      GEN x, y, z;
  661.      
  662. {
  663.   GEN t;
  664.   
  665.   if((typ(x)!=1)||(typ(y)!=1)||(typ(z)!=1)) err(qfer1);
  666.   t=cgetg(4,16);
  667.   t[1] = lcopy(x); t[2] = lcopy(y); t[3] = lcopy(z);
  668.   return t;
  669. }
  670.  
  671.  
  672. GEN  qfr(x, y, z, d)
  673.      GEN x, y, z, d;
  674.      
  675. {
  676.   GEN t;
  677.   
  678.   if((typ(x)!=1)||(typ(y)!=1)||(typ(z)!=1)||(typ(d)!=2)) err(qfer1);
  679.   t=cgetg(5,15);
  680.   t[1]=lcopy(x);t[2]=lcopy(y);t[3]=lcopy(z);t[4]=lcopy(d);
  681.   return t;
  682. }
  683.  
  684. GEN compose(x,y)
  685.      GEN x,y;
  686.      
  687. {
  688.   long av,tetpil;
  689.   GEN s,n,d,d1,d2,x1,x2,y1,y2,v1,v2,v3,b3,c3,m,z,p1,r;
  690.   
  691.   if(cmpii(x[1],y[1])>0) {s=x;x=y;y=s;}
  692.   av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  693.   d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
  694.   v1=divii(x[1],d1);v2=divii(y[1],d1);
  695.   d2=(gcmp1(d1))?gun:mppgcd(d1,mppgcd(x[3],mppgcd(y[3],n)));
  696.   v3=mulii(d2,v1);
  697.   m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
  698.   r=modii(m,v3);b3=shifti((p1=mulii(v2,r)),1);
  699.   c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
  700.   z=cgetg(4,16);z[1]=lmulii(v3,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v3);
  701.   tetpil=avma;return gerepile(av,tetpil,redcomp(z));
  702. }
  703.  
  704. GEN compreal(x,y)
  705.      GEN x,y;
  706.      
  707. {
  708.   long av,tetpil;
  709.   GEN s,n,d,d1,d2,x1,x2,y1,y2,v1,v2,v3,b3,c3,m,z,p1,r;
  710.   
  711.   av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  712.   d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
  713.   v1=divii(x[1],d1);v2=divii(y[1],d1);
  714.   d2=(gcmp1(d1))?gun:mppgcd(d1,mppgcd(x[3],mppgcd(y[3],n)));
  715.   v3=mulii(d2,v1);
  716.   m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
  717.   r=modii(m,v3);b3=shifti((p1=mulii(v2,r)),1);
  718.   c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
  719.   z=cgetg(5,15);z[1]=lmulii(v3,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v3);
  720.   z[4]=laddrr(x[4],y[4]);tetpil=avma;return gerepile(av,tetpil,redreal(z));
  721. }
  722.  
  723. GEN comprealraw(x,y)
  724.      GEN x,y;
  725.      
  726. {
  727.   long av,tetpil;
  728.   GEN s,n,d,d1,d2,x1,x2,y1,y2,v1,v2,v3,b3,c3,m,z,p1,r;
  729.   
  730.   av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  731.   d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
  732.   v1=divii(x[1],d1);v2=divii(y[1],d1);
  733.   d2=(gcmp1(d1))?gun:mppgcd(d1,mppgcd(x[3],mppgcd(y[3],n)));
  734.   v3=mulii(d2,v1);
  735.   m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
  736.   r=modii(m,v3);b3=shifti((p1=mulii(v2,r)),1);
  737.   c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
  738.   z=cgetg(5,15);z[1]=lmulii(v3,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v3);
  739.   z[4]=laddrr(x[4],y[4]);tetpil=avma;return gerepile(av,tetpil,gcopy(z));
  740. }
  741.  
  742. GEN nucomp(x,y,l)
  743.      GEN x,y,l;
  744.      
  745.      /* l=floor((|d|/4)^(1/4)) */
  746.      
  747. {
  748.   long av=avma,tetpil,cz;
  749.   GEN s,n,d,d1,v,u1,u,v1,v2,z,p1,p2,p3;
  750.   GEN a,b,e,f,g,a1,a2,a3,b3,v3,q,t2,t3,c3,q1,q2,q3,q4;
  751.   
  752.   if((typ(x)!=16)||(typ(y)!=16)) err(nucomper1);
  753.   if(x==y) return nudupl(x,l);
  754.   if(cmpii(x[1],y[1])<0) {s=x;x=y;y=s;}
  755.   s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  756.   a1=(GEN)x[1];a2=(GEN)y[1];d=bezout(a2,a1,&u,&v);
  757.   if(gcmp1(d)) {a=negi(gmul(u,n));d1=d;}
  758.   else
  759.     if(divise(s,d)) 
  760.       {
  761.     a=negi(gmul(u,n));d1=d;a1=divii(a1,d1);a2=divii(a2,d1);s=divii(s,d1);
  762.       }
  763.     else
  764.       {
  765.     d1=bezout(s,d,&u1,&v1);
  766.     if(!gcmp1(d1)) 
  767.       {
  768.         a1=divii(a1,d1);a2=divii(a2,d1);s=divii(s,d1);d=divii(d,d1);
  769.       }
  770.     p1=resii(x[3],d);p2=resii(y[3],d);
  771.     p3=modii(negi(mulii(u1,addii(mulii(u,p1),mulii(v,p2)))),d);
  772.     a=subii(mulii(p3,divii(a1,d)),mulii(u,divii(n,d)));
  773.       }
  774.   a=modii(a,a1);p1=subii(a1,a);if(cmpii(a,p1)>0) a=negi(p1);
  775.   v=gzero;d=a1;v2=gun;v3=a;cz=0;
  776.   while(cmpii(v3,l)>=0)
  777.     {
  778.       q=dvmdii(d,v3,&t3);t2=subii(v,mulii(q,v2));
  779.       v=v2;d=v3;v2=t2;v3=t3;cz++;
  780.     }
  781.   if(!cz)
  782.     {
  783.       q1=mulii(a2,v3);q2=addii(q1,n);f=divii(q2,d);
  784.       g=divii(addii(mulii(v3,s),y[3]),d);
  785.       a3=mulii(d,a2);c3=addii(mulii(v3,f),mulii(g,d1));
  786.       b3=addii(shifti(q1,1),y[2]);
  787.     }
  788.   else
  789.     {
  790.       if(cz&1) {v3=negi(v3);v2=negi(v2);}
  791.       b=divii(addii(mulii(a2,d),mulii(n,v)),a1);q1=mulii(b,v3);
  792.       q2=addii(q1,n);f=divii(q2,d);
  793.       e=divii(addii(mulii(s,d),mulii(y[3],v)),a1);
  794.       q3=mulii(e,v2);q4=subii(q3,s);g=divii(q4,v);
  795.       if(!gcmp1(d1)) {v2=mulii(d1,v2);v=mulii(d1,v);}
  796.       a3=addii(mulii(d,b),mulii(e,v));c3=addii(mulii(v3,f),mulii(g,v2));
  797.       b3=addii(addii(q1,q2),mulii(d1,addii(q3,q4)));
  798.     }
  799.   z=cgetg(4,16);z[1]=(long)a3;z[2]=(long)b3;z[3]=(long)c3;  
  800.   tetpil=avma;return gerepile(av,tetpil,redcomp(z));
  801. }
  802.  
  803. GEN nudupl(x,l)
  804.      GEN x,l;
  805.      
  806.      /* l=floor((|d|/4)^(1/4)) */
  807.  
  808. {
  809.   long av=avma,tetpil,cz;
  810.   GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,q,t2,t3,e,g;
  811.  
  812.   d1=bezout(x[2],x[1],&u,&v);a=divii(x[1],d1);b=divii(x[2],d1);
  813.   c=modii(negi(mulii(u,x[3])),a);p1=subii(a,c);if(cmpii(c,p1)>0) c=negi(p1);
  814.   v=gzero;d=a;v2=gun;v3=c;cz=0;
  815.   while(cmpii(v3,l)>=0)
  816.     {
  817.       q=dvmdii(d,v3,&t3);t2=subii(v,mulii(q,v2));
  818.       v=v2;d=v3;v2=t2;v3=t3;cz++;
  819.     }
  820.   if(!cz)
  821.     {
  822.       g=divii(addii(mulii(b,v3),x[3]),d);
  823.       z=cgetg(4,16);z[1]=lmulii(d,d);
  824.       z[2]=laddii(x[2],shifti(mulii(d,v3),1));
  825.       z[3]=laddii(mulii(v3,v3),mulii(g,d1));
  826.     }
  827.   else
  828.     {
  829.       if(cz&1) {v=negi(v);d=negi(d);}
  830.       e=divii(addii(mulii(x[3],v),mulii(b,d)),a);
  831.       g=divii(subii(mulii(e,v2),b),v);b2=addii(mulii(e,v2),mulii(v,g));
  832.       if(!gcmp1(d1)) {b2=mulii(d1,b2);v=mulii(d1,v);v2=mulii(d1,v2);}
  833.       z=cgetg(4,16);z[1]=laddii(mulii(d,d),mulii(e,v));
  834.       z[2]=laddii(b2,shifti(mulii(d,v3),1));
  835.       z[3]=laddii(mulii(v3,v3),mulii(g,v2));
  836.     }
  837.   tetpil=avma;return gerepile(av,tetpil,redcomp(z));
  838. }
  839.  
  840. GEN sqcomp(x)
  841.      GEN x;
  842.      
  843. {
  844.   long av,tetpil;
  845.   GEN d1,d2,x2,y2,v1,v3,b3,c3,m,z,p1,r;
  846.   
  847.   av=avma;
  848.   d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
  849.   d2=(gcmp1(d1))?gun:mppgcd(d1,x[3]);v3=mulii(d2,v1);
  850.   m=mulii(x[3],x2);setsigne(m,-signe(m));
  851.   r=modii(m,v3);b3=shifti((p1=mulii(v1,r)),1);
  852.   c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
  853.   z=cgetg(4,16);z[1]=lmulii(v3,v1);
  854.   z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v3);
  855.   tetpil=avma;return gerepile(av,tetpil,redcomp(z));
  856. }
  857.  
  858. GEN sqcompreal(x)
  859.      GEN x;
  860.      
  861. {
  862.   long av,tetpil;
  863.   GEN d1,d2,x2,y2,v1,v3,b3,c3,m,z,p1,r;
  864.   
  865.   av=avma;
  866.   d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
  867.   d2=(gcmp1(d1))?gun:mppgcd(d1,x[3]);v3=mulii(d2,v1);
  868.   m=mulii(x[3],x2);setsigne(m,-signe(m));
  869.   r=modii(m,v3);b3=shifti((p1=mulii(v1,r)),1);
  870.   c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
  871.   z=cgetg(5,15);z[1]=lmulii(v3,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v3);
  872.   z[4]=lshiftr(x[4],1);tetpil=avma;return gerepile(av,tetpil,redreal(z));
  873. }
  874.  
  875. GEN sqcomprealraw(x)
  876.      GEN x;
  877.      
  878. {
  879.   long av,tetpil;
  880.   GEN d1,d2,x2,y2,v1,v3,b3,c3,m,z,p1,r;
  881.   
  882.   av=avma;
  883.   d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
  884.   d2=(gcmp1(d1))?gun:mppgcd(d1,x[3]);v3=mulii(d2,v1);
  885.   m=mulii(x[3],x2);setsigne(m,-signe(m));
  886.   r=modii(m,v3);b3=shifti((p1=mulii(v1,r)),1);
  887.   c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
  888.   z=cgetg(5,15);z[1]=lmulii(v3,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v3);
  889.   z[4]=lshiftr(x[4],1);tetpil=avma;return gerepile(av,tetpil,gcopy(z));
  890. }
  891.  
  892. GEN powrealraw(x,n,prec)
  893.      GEN x;
  894.      long n,prec;
  895. {
  896.   long av,tetpil;
  897.   GEN p1,p2,y,z;
  898.  
  899.   if(n<0) err(impl,"powrealraw for negative exponents");
  900.   if(n==1) return gcopy(x);
  901.   z=x;av=avma;y=cgetg(5,15);y[1]=un;
  902.   p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  903.   p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  904.   y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  905.   p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  906.   p1[2]=0;tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  907.   if(!n) return y;
  908.   else
  909.     {
  910.       for(;n>1;n>>=1)
  911.       {
  912.         if (n&1) y=comprealraw(y,z);
  913.         z=sqcomprealraw(z);
  914.       }
  915.       tetpil=avma;return gerepile(av,tetpil,comprealraw(y,z));
  916.     }
  917. }
  918.  
  919. GEN nupow(x,n)
  920.      GEN x,n;
  921. {
  922.   long av,tetpil,i,j;
  923.   unsigned long m;
  924.   GEN p1,p2,y,z,l;
  925.  
  926.   if(typ(n)!=1) err(nupower1);
  927.   if(gcmp1(n)) return gcopy(x);
  928.   z=x;av=avma;y=cgetg(4,16);y[1]=un;y[2]=mpodd(x[2]) ? un : zero;
  929.   p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);y[3]=lsubii(p1,p2);
  930.   if(!signe(n)) {tetpil=avma;return gerepile(av,tetpil,gcopy(y));}
  931.   else
  932.     {
  933.       l=racine(shifti(racine(y[3]),1));
  934.       for (i=lgef(n)-1;i>2;i--)
  935.       {
  936.         for (m=n[i],j=0;j<32;j++,m>>=1)
  937.         {
  938.           if (m&1) y=nucomp(y,z,l);
  939.           z=nudupl(z,l);
  940.         }
  941.       }
  942.       for (m=n[2];m>1;m>>=1)
  943.       {
  944.         if (m&1) y=nucomp(y,z,l);
  945.         z=nudupl(z,l);
  946.       }
  947.       tetpil=avma;y=nucomp(y,z,l);
  948.       if ((signe(n)<0)&&cmpii(y[1],y[2])&&cmpii(y[1],y[3]))
  949.     setsigne(y[2],-signe(y[2]));
  950.       return gerepile(av,tetpil,y);
  951.     }
  952. }
  953.  
  954. GEN redcomp(x)
  955.      GEN x;
  956.      
  957. {
  958.   long av,tetpil,s,fl,fg;
  959.   GEN b1,c1,p1,a,b,c,d,z;
  960.   
  961.   av=avma;a=(GEN)x[1];b=(GEN)x[2];c=(GEN)x[3];
  962.   fl=cmpii(a,c);s=signe(b);setsigne(b,abs(s));fg=cmpii(a,b);
  963.   while((fl>0)||(fg<0))
  964.     {
  965.       p1=shifti(c,1);d=dvmdii(b,p1,&b1);
  966.       setsigne(b,s);
  967.       if(s>=0)
  968.         {
  969.           if(cmpii(b1,c)<=0) {setsigne(d,-signe(d));setsigne(b1,-signe(b1));}
  970.           else {d=subsi(-1,d);b1=subii(p1,b1);}
  971.         }
  972.       else
  973.         {
  974.           if(cmpii(b1,c)>=0) {d=addsi(1,d);b1=subii(b1,p1);}
  975.         }
  976.       c1=addii(a,mulii(d,shifti(subii(b,b1),-1)));a=c;
  977.       b=b1;c=c1;
  978.       fl=cmpii(a,c);s=signe(b);setsigne(b,abs(s));fg=cmpii(a,b);
  979.     }
  980.   if(fl&&fg&&(s<0)) setsigne(b,s);tetpil=avma;
  981.   z=cgetg(4,16);z[1]=lcopy(a);z[2]=lcopy(b);z[3]=lcopy(c);
  982.   return gerepile(av,tetpil,z);
  983. }
  984.  
  985. GEN rhoreal(x)
  986.      GEN x;
  987. {
  988.   long av,tetpil,s,l;
  989.   GEN y,p1,nn,sqrtD,isqrtD,D;
  990.   
  991.   if(typ(x)!=15) err(rhoer1);
  992.   av=avma;y=cgetg(5,15);
  993.   l=max(lg(x[4]),((-expo(x[4]))>>5)+2);if(l<=2) l=3;
  994.   y[1]=lcopy(x[3]);sqrtD=gsqrt(D=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2)),l);
  995.   isqrtD=mptrunc(sqrtD);
  996.   s=signe(y[1]);setsigne(y[1],1);
  997.   if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
  998.   else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
  999.   p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
  1000.   setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
  1001.   y[4]=laddrr(shiftr(mplog(absr(divrr(addir(x[2],sqrtD),subir(x[2],sqrtD)))),-1),x[4]);
  1002.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1003. }
  1004.  
  1005. GEN redreal(x)
  1006.      GEN x;
  1007. {
  1008.   long fl,av=avma,tetpil,l,s;
  1009.   GEN y,p1,sqrtD,isqrtD,D,z,nn;
  1010.   
  1011.   if(typ(x)!=15) err(rhoer1);
  1012.   if(signe(x[4])) l=lg(x[4]);
  1013.   else {l=((-expo(x[4]))>>5)+2;if(l<=2) l=3;}
  1014.   sqrtD=gsqrt(D=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2)),l);
  1015.   isqrtD=mptrunc(sqrtD);
  1016.   y=cgetg(5,15);y[1]=x[1];y[2]=x[2];y[3]=x[3];y[4]=x[4];
  1017.   if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
  1018.   else
  1019.     {
  1020.       p1=subii(isqrtD,shifti(absi(x[1]),1));
  1021.       if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
  1022.     }
  1023.   while(fl)
  1024.     {
  1025.       z=cgetg(5,15);
  1026.       z[1]=y[3];s=signe(z[1]);setsigne(z[1],1);
  1027.       if(cmpii(isqrtD,z[1])>=0) nn=divii(addii(isqrtD,y[2]),p1=shifti(z[1],1));
  1028.       else nn=divii(addii(z[1],y[2]),p1=shifti(z[1],1));
  1029.       p1=mulii(nn,p1);z[2]=lsubii(p1,y[2]);
  1030.       setsigne(z[1],s);p1=shifti(subii(mulii(z[2],z[2]),D),-2);z[3]=ldivii(p1,z[1]);
  1031.       z[4]=laddrr(shiftr(mplog(absr(divrr(addir(y[2],sqrtD),subir(y[2],sqrtD)))),-1),y[4]);
  1032.       y=z;
  1033.       if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
  1034.       else
  1035.     {
  1036.       p1=subii(isqrtD,shifti(absi(y[1]),1));
  1037.       if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
  1038.     }
  1039.     }
  1040.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1041. }
  1042.  
  1043. GEN rhorealnod(x,isqrtD)
  1044.      GEN x,isqrtD;
  1045. {
  1046.   long av,tetpil,s;
  1047.   GEN y,p1,nn,D;
  1048.   
  1049.   if(typ(x)!=15) err(rhoer1);
  1050.   if(typ(isqrtD)!=1) err(arither1);
  1051.   av=avma;y=cgetg(5,15);
  1052.   y[1]=lcopy(x[3]);D=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1053.   s=signe(y[1]);setsigne(y[1],1);
  1054.   if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
  1055.   else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
  1056.   p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
  1057.   setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
  1058.   affsr(0,y[4]=lgetr(3));
  1059.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1060. }
  1061.  
  1062. GEN redrealnod(x,isqrtD)
  1063.      GEN x,isqrtD;
  1064. {
  1065.   long fl,av=avma,tetpil,s;
  1066.   GEN y,p1,D,z,nn;
  1067.   
  1068.   if(typ(x)!=15) err(rhoer1);
  1069.   if(typ(isqrtD)!=1) err(arither1);
  1070.   D=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1071.   y=cgetg(5,15);y[1]=x[1];y[2]=x[2];y[3]=x[3];
  1072.   if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
  1073.   else
  1074.     {
  1075.       p1=subii(isqrtD,shifti(absi(x[1]),1));
  1076.       if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
  1077.     }
  1078.   while(fl)
  1079.     {
  1080.       z=cgetg(5,15);
  1081.       z[1]=y[3];s=signe(z[1]);setsigne(z[1],1);
  1082.       if(cmpii(isqrtD,z[1])>=0) nn=divii(addii(isqrtD,y[2]),p1=shifti(z[1],1));
  1083.       else nn=divii(addii(z[1],y[2]),p1=shifti(z[1],1));
  1084.       p1=mulii(nn,p1);z[2]=lsubii(p1,y[2]);
  1085.       setsigne(z[1],s);p1=shifti(subii(mulii(z[2],z[2]),D),-2);z[3]=ldivii(p1,z[1]);
  1086.       y=z;
  1087.       if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
  1088.       else
  1089.     {
  1090.       p1=subii(isqrtD,shifti(absi(y[1]),1));
  1091.       if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
  1092.     }
  1093.     }
  1094.   affsr(0,y[4]=lgetr(3));
  1095.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1096. }
  1097.  
  1098. GEN primeform(x,p,prec)
  1099.      GEN x,p;
  1100.      long prec;
  1101. {
  1102.   long av,tetpil,s,t,sb,sx=signe(x);
  1103.   GEN y,b,c;
  1104.   
  1105.   if((typ(x)!=1)||(!sx)) err(arither1);
  1106.   if(gcmpgs(p,2))
  1107.     {
  1108.       y=(sx<0) ? cgetg(4,16) : cgetg(5,15);y[1]=lcopy(p);y[2]=lgeti(lgef(p));
  1109.       av=avma;
  1110.       if(!mpsqrtmod(x,p,&b)) err(sqrter5);
  1111.       s=b[lgef(b)-1]&1;t=x[lgef(x)-1]&1;
  1112.       if(odd(s+t)) subiiz(p,b,y[2]);else affii(b,y[2]);
  1113.       c=shifti(subii(mulii(y[2],y[2]),x),-2);tetpil=avma;
  1114.       y[3]=lpile(av,tetpil,divii(c,p));
  1115.     }
  1116.   else
  1117.     {
  1118.       s=x[lgef(x)-1]&7;if(signe(x)<0) s=8-s;
  1119.       switch(s)
  1120.         {
  1121.         case 0:
  1122.         case 8: sb=0;break;
  1123.         case 1: sb=1;break;
  1124.         case 4: sb=2;break;
  1125.         default: err(sqrter5);
  1126.         }
  1127.       y=(sx<0) ? cgetg(4,16) : cgetg(5,15);y[1]=lcopy(deux);y[2]=lstoi(sb);
  1128.       av=avma;c=gsubsg(sb*sb,x);tetpil=avma;
  1129.       y[3]=lpile(av,tetpil,shifti(c,-3));
  1130.     }
  1131.   if(sx>0) affsr(0,y[4]=lgetr(prec));
  1132.   return y;
  1133. }
  1134.  
  1135. GEN classno(x)
  1136.      GEN x;
  1137.      
  1138.      /* calcul de h(x) pour x<0 par la methode de Shanks */
  1139.      
  1140. {
  1141.   static long prem,ptforme;
  1142.   long av,av1,av2,av3,tetpil,k,k2,i,j,j1,lim,limite,succes,com,j2,s,tx;
  1143.   GEN tabla, tablb, count, index, hash;
  1144.   static long court[3]={0x1000003,0x1010003,0};
  1145.   GEN p1,p2,p3,hin,q,formes[100],l,h,hp,f,fh,fg,ftest,pm2;
  1146.   static byteptr p;
  1147.   
  1148.   if((tx=typ(x))>=17) 
  1149.     {
  1150.       k=lg(x);p1=cgetg(k,tx);for(i=1;i<k;i++) p1[i]=(long)classno(x[i]);
  1151.       return p1;
  1152.     }
  1153.   if(tx!=1) err(arither1);
  1154.   if (!(s=signe(x))) err(arither2);
  1155.   if(s>0) return classno2(x);
  1156.   k=x[lgef(x)-1]&3;
  1157.   if((k==1)||(k==2)) err(classer2);
  1158.   if(gcmpgs(x,-12)>=0) return gun;
  1159.   tabla = newbloc(10000);
  1160.   tablb = newbloc(10000);
  1161.   count = newbloc(256);
  1162.   index = newbloc(257);
  1163.   hash = newbloc(10000);
  1164.   prem=ptforme=0;p=diffptr;av=avma;limite=(av+bot)>>1;
  1165.   p1=divrr(gsqrt(absi(x),4),mppi(4));
  1166.   l=gtrunc(shiftr(gsqrt(gsqrt(absi(x),4),4),1));
  1167.   if(gcmpgs(l,1000)<0) affsi(1000,l);
  1168.   while((gcmpgs(l,prem)>=0)&&(*p))
  1169.     {
  1170.       prem+= *p++;
  1171.       k=krogs(x,prem); 
  1172.       if(k)
  1173.     {
  1174.       av2=avma;
  1175.       if(k>0)
  1176.         {
  1177.           divrsz(mulsr(prem,p1),prem-1,p1);avma=av2;
  1178.           if((ptforme<100)&&(prem>2))
  1179.         {
  1180.           court[2]=prem;
  1181.           formes[ptforme++]=redcomp(primeform(x,court,3));
  1182.         }
  1183.         }
  1184.       else {divrsz(mulsr(prem,p1),prem+1,p1);avma=av2;}
  1185.     }
  1186.     }
  1187.   hin=ground(p1);h=gcopy(hin);f=formes[0];fh=gpui(f,h);
  1188.   s=2*itos(gceil(gpui(p1,gdivgs(gun,4),4)));
  1189.   if(s>=10000) s=10000;
  1190.   p1=fh;av2=avma;
  1191.   for(i=0;i<=255;i++) count[i]=0;
  1192.   for(i=0;i<=s-1;i++)
  1193.     {
  1194.       p2=(GEN)p1[1];tabla[i]=p2[lgef(p2)-1];j=tabla[i]&255;count[j]++;
  1195.       p2=(GEN)p1[2];tablb[i]=p2[lgef(p2)-1];
  1196.       p1=compose(p1,f);
  1197.     }
  1198.   fg=gpuigs(f,s);ftest=gpuigs(p1,0);
  1199.   index[0]=0;for(i=0;i<=254;i++) index[i+1]=index[i]+count[i];
  1200.   for(i=0;i<=s-1;i++) hash[index[tabla[i]&255]++]=i;
  1201.   index[0]=0;for(i=0;i<=255;i++) index[i+1]=index[i]+count[i];
  1202.   succes=0;com=0;av3=avma;
  1203.   while(!succes)
  1204.     {
  1205.       p1=(GEN)ftest[1];k=p1[lgef(p1)-1];j=k&255;
  1206.       pm2=negi(ftest[2]);
  1207.       for(j1=index[j];(j1<index[j+1])&&(!succes);j1++)
  1208.     {
  1209.       j2=hash[j1];
  1210.       if(tabla[j2]==k)
  1211.         {
  1212.           p2=(GEN)ftest[2];k2=p2[lgef(p2)-1];
  1213.           if((tablb[j2]==k2)||(tablb[j2]== -k2))
  1214.         {
  1215.           p1=gmul(gpuigs(f,j2),fh);
  1216.           succes=(!cmpii(p1[1],ftest[1]))&&((!cmpii(p1[2],ftest[2]))||(!cmpii(p1[2],pm2)));
  1217.         }
  1218.         }
  1219.     }
  1220.       if(!succes)
  1221.     {
  1222.       com++;
  1223.       if(avma>=limite) ftest=gmul(ftest,fg);
  1224.       else {tetpil=avma;ftest=gerepile(av3,tetpil,gmul(ftest,fg));}
  1225.       if (gcmp1(ftest[1]))
  1226.         {
  1227.           err(impl, "classno with too small order");
  1228.         }
  1229.     }
  1230.     }
  1231.   h=addis(h,j2);p2=mulsi(s,stoi(com));tetpil=avma;
  1232.   h=(!cmpii(p1[2],ftest[2]))?subii(h,p2):addii(h,p2);
  1233.   p2=factor(h);
  1234.   p1=(GEN)p2[1];
  1235.   p2=(GEN)p2[2];
  1236.   for(i=1;i<lg(p1);i++)
  1237.     {
  1238.       p3=divii(h,p1[i]);fh=gpui(f,p3);lim=itos(p2[i]);
  1239.       for(j=1;(j<=lim)&&(!cmpii(fh[1],un));j++)
  1240.     {
  1241.       h=p3;
  1242.       if(j<lim) {p3=divii(h,p1[i]);fh=gpui(f,p3);}
  1243.     }
  1244.     }
  1245.   q=ground(gdiv(hin,h));tetpil=avma;
  1246.   hp=mulii(q,h);av1=avma;
  1247.   for(i=1;(i<=10)&&(i<ptforme);i++)
  1248.     {
  1249.       f=formes[i];com=0;
  1250.       fg=gpui(f,h);
  1251.       fh=gpui(fg,q);
  1252.       ftest=fg;
  1253.       if(cmpis(fh[1],1))
  1254.     {
  1255.       for(;;)
  1256.         {
  1257.           com++;
  1258.           if ((!cmpii(fh[1],ftest[1]))&&((!cmpii(fh[2],ftest[2]))||(!cmpii(fh[2],negi(ftest[2]))))) break;
  1259.           ftest=gmul(ftest,fg);
  1260.         }
  1261.       if(!cmpii(fh[2],ftest[2])) com= -com;
  1262.       p2=mulsi(com,h);q=addsi(com,q);tetpil=avma;
  1263.       hp=addii(hp,p2);av1=avma;
  1264.     }
  1265.     }
  1266.   avma=av1;
  1267.   killbloc(tabla); killbloc(tablb); killbloc(count);
  1268.   killbloc(index); killbloc(hash);
  1269.   return gerepile(av,tetpil,hp);
  1270. }
  1271.