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