home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / buch.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  46.3 KB  |  1,653 lines
  1.  /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*                                                                 */
  4. /*            ALGORITHMES SOUS-EXPONENTIELS DE CALCUL              */
  5. /*            DU GROUPE DE CLASSES ET DU REGULATEUR                */
  6. /*                    (McCURLEY, BUCHMANN)                         */
  7. /*                                                                 */
  8. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  9. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  10.  
  11. #include "genpari.h"
  12. #define C 15
  13. #define HASHT 1024
  14. #define EXP220 1048576
  15.  
  16. void hphaseimag1(),hphaseimag2(),hphasereal1(),hphasereal2();
  17. long *subfactorbaseimag(),*subfactorbasereal(),factorbase();
  18. long factorise(),*largeprime(),*largeprime2();
  19. GEN **powsubfactimag(),**powsubfactreal(),rangi(),clean(),hnfimag();
  20. GEN sqrealform3(),comprealform3(),rhorealform3(),redrealform3();
  21. GEN sqrealform5(),powrealform5(),comprealform5(),rhorealform5(),redrealform5();
  22. GEN redrealform(),gcdreal(),hnfreal(),lfunc();
  23.  
  24. GEN buchimag(D,gcbach,gCO)
  25.      GEN D,gcbach,gCO;
  26. {
  27.   long CO;
  28.   double cbach;
  29.   long limc,limc2,mglob,m,cp,nbram,auxrel,lo,lo1,ran,fl,fl2;
  30.   long *fpd,b1,b2,fpc,pp,ep,b,extrarel,c,ic,jc,limhash;
  31.   long kc2,nbcol1,nbcol2,nbcol3,nbrow1,nbrow2,nbrow3,*numbase,*base,*subbase;
  32.   long **mat,*vectprime;
  33.   long primfact[100],expoprimfact[100],primfact1[100],expoprimfact1[100];
  34.   long badprim[100],nbbadprim;
  35.   long av=avma,tetpil,kc,kcco,i,j,*pro,p2,*p1,*ex,q,s,nbtest,mm,av3;
  36.   double drc,lim,logd;
  37.   GEN dr,v,detmat,detmatsur2;
  38.   GEN matgen,**vp,form,form1,pc,mit,met,mot,p3,p4;
  39.   GEN matpro; 
  40.   long* hashtab[HASHT];
  41.  
  42.   if((typ(D)!=1)||(typ(gCO)!=1)) 
  43.     err(talker,"incorrect type in buch");
  44.   if(!signe(D)) err(talker,"zero discriminant in buch");
  45.   if(signe(D)>0) err(talker,"positive discriminant in buchimag");
  46.   s=D[lgef(D)-1]&3;
  47.   if((s==1)||(s==2)) err(talker,"discriminant not 0 or 1 mod 4 in buch");
  48.   if(signe(gCO)<=0) err(talker,"not enough relations in buch");
  49.   CO=itos(gCO);cbach=gtodouble(gcbach);
  50.   if((!cmpis(D,-3))||(!cmpis(D,-4)))
  51.     {
  52.       p3=cgetg(3,17);p3[1]=un;p3[2]=lgetg(1,17);return p3;
  53.     }
  54.   dr=cgetr(3);affir(D,dr);drc=rtodbl(dr);
  55.   lim=sqrt(fabs(drc)/3.0);
  56.   logd=log(fabs(drc));limc=cbach*logd*logd;cp=exp(sqrt(logd*log(logd)/8.0));
  57.   if(cp>limc) limc=cp;limc2=6*logd*logd;if(limc>limc2) limc2=limc;
  58.   subbase=subfactorbaseimag(D,&v,lim,&nbram);
  59.   mglob=m=subbase[0];
  60.   vp=powsubfactimag(v,m,C+7);
  61.   kc2=factorbase(D,limc2,limc,&numbase,&base,&kc,badprim,&nbbadprim);
  62.   pro=(long*)malloc(4*(kc2+1));vectprime=(long*)malloc(4*(kc2+1));
  63.   for(i=1;i<=kc2;i++)
  64.     {
  65.       p2=vectprime[i]=base[i];
  66.       pro[i]=(p2>0) ? p2 : -p2;
  67.     }
  68.   free(base);base=pro;
  69.   kcco=kc+CO;
  70.   mat=(long **)malloc(4*(kcco+1));
  71.   for(i=1;i<=kcco;i++)
  72.     {
  73.       p1=(long *)malloc(4*(kc+1));mat[i]=p1;
  74.       for(j=1;j<=kc;j++) p1[j]=0;
  75.     }
  76.   ex=(long*)malloc((m+1)*4);
  77.   q=31-(long)ceil(log((double)C)/log(2.0));
  78.   s=0;nbtest=0;
  79.   mm=m+nbram+CO;
  80.   limhash=(limc<32767)?limc*limc:(2<<30);
  81.   for(i=0;i<HASHT;i++) hashtab[i]=(long*)0;
  82.   auxrel=0;
  83.   while(s<mm)
  84.     {
  85.       for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  86.       av3=avma;
  87.       form=vp[1][ex[1]];
  88.       for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
  89.       fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  90.       nbtest++;
  91.       if(fpc>1)
  92.     {
  93.       fpd=largeprime(fpc,ex,0,0,mglob,hashtab);
  94.       if(fpd)
  95.         {
  96.           auxrel++;
  97.           s++;lo1=lo;
  98.           for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
  99.           form1=vp[1][fpd[2]];
  100.           for(i=2;i<=m;i++) 
  101.         form1=compose(form1,vp[i][fpd[i+1]]);
  102.           b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
  103.           factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  104.           if(b1==b2)
  105.         {
  106.           for(i=1;i<=m;i++)
  107.             mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
  108.           for(j=1;j<=lo;j++)
  109.             {
  110.               pp=primfact[j];ep=expoprimfact[j];
  111.               b1=itos(modis(form1[2],(pp<<1)));
  112.               if(b1>pp) ep= -ep;
  113.               mat[s][numbase[pp]]-=ep;
  114.             }
  115.           for(j=1;j<=lo1;j++)
  116.             {
  117.               pp=primfact1[j];ep=expoprimfact1[j];
  118.               b1=itos(modis(form[2],(pp<<1)));
  119.               if(b1>pp) ep= -ep;
  120.               mat[s][numbase[pp]]+=ep;
  121.             }
  122.         }
  123.           else
  124.         {
  125.           if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
  126.           else
  127.             {
  128.               for(i=1;i<=m;i++)
  129.             mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
  130.               for(j=1;j<=lo;j++)
  131.             {
  132.               pp=primfact[j];ep=expoprimfact[j];
  133.               b1=itos(modis(form1[2],(pp<<1)));
  134.               if(b1>pp) ep= -ep;
  135.               mat[s][numbase[pp]]+=ep;
  136.             }
  137.               for(j=1;j<=lo1;j++)
  138.             {
  139.               pp=primfact1[j];ep=expoprimfact1[j];
  140.               b1=itos(modis(form[2],(pp<<1)));
  141.               if(b1>pp) ep= -ep;
  142.               mat[s][numbase[pp]]+=ep;
  143.             }
  144.             }
  145.         }
  146.         }
  147.     }      
  148.       if(fpc==1)
  149.     {
  150.       s++;
  151.       for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
  152.       for(j=1;j<=lo;j++)
  153.         {
  154.           pp=primfact[j];ep=expoprimfact[j];
  155.           b=itos(modis(form[2],(pp<<1)));
  156.           if(b>pp) ep= -ep;
  157.           mat[s][numbase[pp]]+=ep;
  158.         }
  159.     }
  160.       avma=av3;
  161.     }
  162.   while(s<kcco)
  163.     {
  164.       for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  165.       av3=avma;
  166.       form=vp[1][ex[1]];
  167.       for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
  168.       pc=primeform(D,stoi(base[s+1-CO]));
  169.       form=compose(form,pc);
  170.       fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  171.       nbtest++;
  172.       if(fpc>1)
  173.     {
  174.       fpd=largeprime(fpc,ex,s+1-CO,0,mglob,hashtab);
  175.       if(fpd)
  176.         {
  177.           auxrel++;
  178.           s++;lo1=lo;
  179.           for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
  180.           form1=vp[1][fpd[2]];
  181.           for(i=2;i<=m;i++) 
  182.         form1=compose(form1,vp[i][fpd[i+1]]);
  183.           if(fpd[m+2])
  184.         {
  185.           pc=primeform(D,stoi(base[fpd[m+2]]));
  186.           form1=compose(form1,pc);
  187.         }
  188.           b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
  189.           factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  190.           if(b1==b2)
  191.         {
  192.           for(i=1;i<=m;i++)
  193.             mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
  194.           mat[s][s-CO]= -1;
  195.           if(fpd[m+2]) mat[s][fpd[m+2]]++;
  196.           for(j=1;j<=lo;j++)
  197.             {
  198.               pp=primfact[j];ep=expoprimfact[j];
  199.               b1=itos(modis(form1[2],(pp<<1)));
  200.               if(b1>pp) ep= -ep;
  201.               mat[s][numbase[pp]]-=ep;
  202.             }
  203.           for(j=1;j<=lo1;j++)
  204.             {
  205.               pp=primfact1[j];ep=expoprimfact1[j];
  206.               b1=itos(modis(form[2],(pp<<1)));
  207.               if(b1>pp) ep= -ep;
  208.               mat[s][numbase[pp]]+=ep;
  209.             }
  210.         }
  211.           else
  212.         {
  213.           if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
  214.           else
  215.            {
  216.               for(i=1;i<=m;i++)
  217.             mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
  218.               mat[s][s-CO]= -1;
  219.               if(fpd[m+2]) mat[s][fpd[m+2]]--;
  220.               for(j=1;j<=lo;j++)
  221.             {
  222.               pp=primfact[j];ep=expoprimfact[j];
  223.               b1=itos(modis(form1[2],(pp<<1)));
  224.               if(b1>pp) ep= -ep;
  225.               mat[s][numbase[pp]]+=ep;
  226.             }
  227.               for(j=1;j<=lo1;j++)
  228.             {
  229.               pp=primfact1[j];ep=expoprimfact1[j];
  230.               b1=itos(modis(form[2],(pp<<1)));
  231.               if(b1>pp) ep= -ep;
  232.               mat[s][numbase[pp]]+=ep;
  233.             }
  234.             }
  235.         }
  236.         }
  237.     }
  238.       if(fpc==1)
  239.     {
  240.       s++;
  241.       for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
  242.       mat[s][s-CO]= -1;
  243.       for(j=1;j<=lo;j++)
  244.         {
  245.           pp=primfact[j];ep=expoprimfact[j];
  246.           b=itos(modis(form[2],(pp<<1)));
  247.           if(b>pp) ep= -ep;
  248.           mat[s][numbase[pp]]+=ep;
  249.         }
  250.       if(!(mat[s][s-CO])) {for(i=1;i<=kc;i++) mat[s][i]=0;s--;}
  251.     }
  252.       avma=av3;
  253.     }
  254.   nbtest=auxrel=0;
  255.   while(s<kc2)
  256.     {
  257.       for (i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  258.       av3=avma;
  259.       form=vp[1][ex[1]];
  260.       for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
  261.       pp=base[s+1];
  262.       pc=primeform(D,stoi(pp));
  263.       form=compose(form,pc);
  264.       fpc=factorise(form,s,pp,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  265.       nbtest++;
  266.       if(fpc>1)
  267.     {
  268.       fpd=largeprime2(fpc,ex,s+1,mglob,hashtab);
  269.       if(fpd&&(fpc!=pp)) {auxrel++;s++;}
  270.     }
  271.       if(fpc==1) s++;
  272.       avma=av3;
  273.     }
  274.   hphaseimag1(mat,vectprime,kcco,kc,&nbcol1,&nbrow1);
  275.   hphaseimag2(mat,vectprime,nbcol1,nbrow1,&nbcol2,&nbrow2);
  276.   nbcol3=nbcol2;nbrow3=nbrow2;
  277.   matgen=cgetg(nbcol3+1,19);
  278.   for(i=1;i<=nbcol3;i++) matgen[i]=lgetg(nbrow3+1,18);
  279.   for(i=1;i<=nbcol3;i++)
  280.     {
  281.       for(j=1;j<=nbrow3;j++) ((GEN)(matgen[i]))[j]=lstoi(mat[i][j]);
  282.     }
  283.   detmat=absi(rangi(matgen,&ran));
  284.   detmatsur2=shifti(detmat,-1);
  285.   if(gcmp1(absi(detmat))) {mit=idmat(1);nbrow3=1;nbcol3=lg(matpro);}
  286.   else mit=hnfimag(matgen,detmat,detmatsur2);
  287.   extrarel=nbcol3-nbrow3;
  288.   p1=gun;
  289.   for(i=1;i<=nbrow3;i++)
  290.     {
  291.       p3=(GEN)coeff(mit,i,i);
  292.       if(signe(p3)) p1=mulii(p1,p3);
  293.     }
  294.   c=0;for(i=1;i<=nbrow3;i++) if(!gcmp1(coeff(mit,i,i))) c++;
  295.   if(!c) met=idmat(1);
  296.   else
  297.     {
  298.       met=cgetg(c+1,19);ic=0;
  299.       for(i=1;i<=nbrow3;i++)
  300.     {
  301.       if(!gcmp1(coeff(mit,i,i)))
  302.         {
  303.           ic++;p3=cgetg(c+1,18);met[ic]=(long)p3;jc=0;
  304.           for(j=1;j<=nbrow3;j++)
  305.         {
  306.           if(!gcmp1(coeff(mit,j,j)))
  307.             {
  308.               jc++;p3[jc]=coeff(mit,j,i);
  309.             }
  310.         }
  311.         }
  312.     }
  313.     }
  314.   met=smith(met);
  315.   lo=lg(met)-1;
  316.   c=0;for(i=1;i<=lo;i++) if(!gcmp1(met[i])) c++;
  317.   p3=gdiv(gmul(p1,mppi(4)),gmul(lfunc(D),gsqrt(absi(D),4)));
  318.   tetpil=avma;p4=cgetg(4,17);p4[1]=lcopy(p1);mot=cgetg(c+1,17);p4[2]=(long)mot;
  319.   for(i=1;i<=c;i++) mot[i]=lcopy(met[i]);p4[3]=lcopy(p3);
  320.   free(base);free(vectprime);free(ex);free(subbase);free(numbase);
  321.   for(i=1;i<=kcco;i++) free(mat[i]);free(mat);
  322.   for(i=1;i<=m;i++) free(vp[i]);free(vp);
  323.   return gerepile(av,tetpil,p4);
  324. }
  325.  
  326. long *subfactorbaseimag(d,w,ll,ptnbram)
  327.      long *ptnbram;
  328.      double ll;
  329.      GEN d,*w;
  330. {
  331.   byteptr delta=diffptr;
  332.   long av1,i,j,pp,tp,fl,*subbase,pro[100],toto=0;
  333.   double prod;
  334.   GEN p1;
  335.  
  336.   av1=avma;
  337.   i=1;pp=0;*ptnbram=0;fl=1;prod=1;
  338.   while(fl||(i==1))
  339.     {
  340.       pp+=*delta++;tp=krogs(d,pp);
  341.       if(tp==1) {pro[i]=pp;i++;prod*=pp;} else if(!tp) (*ptnbram)++;
  342.       if((prod>ll)&&(i>1)) fl=0;
  343.     }
  344.   avma=av1;i--;
  345.   *w=cgetg(i+1,18);
  346.   for(j=1;j<=i;j++) 
  347.     {
  348.       p1=primeform(d,stoi(pro[j]));(*w)[j]=(long)p1;
  349.     }
  350.   subbase=(long*)malloc(4*(i+1));subbase[0]=i;
  351.   for(j=1;j<=i;j++) subbase[j]=pro[j];
  352.   return subbase;
  353. }
  354.   
  355. GEN **powsubfactimag(w,n,a)
  356.      GEN w;
  357.      long n,a;
  358. {
  359.   long i,j,toto=0;
  360.   GEN **x;
  361.  
  362.   x=(GEN**)malloc(4*(n+1));
  363.   for(i=1;i<=n;i++) x[i]=(GEN*)malloc(4*(a+1));
  364.   for(i=1;i<=n;i++)
  365.     {
  366.       x[i][0]=gpuigs(w[1],0);
  367.       for(j=1;j<=a;j++) x[i][j]=compose(x[i][j-1],w[i]);
  368.     }
  369.   return x;
  370. }
  371.  
  372. long factorbase(d,n2,n,ptnum,ptbase,ptkc,badprim,nbbadprim)
  373.      GEN d;
  374.      long n2,n,**ptnum,**ptbase,*ptkc,*badprim,*nbbadprim;
  375. {
  376.   byteptr delta=diffptr;
  377.   long av2,i,pp,qq,fl,kr,r,*numbase,*base,kc;
  378.   GEN p1;
  379.  
  380.   numbase=(long*)malloc(4*(n2+1));*ptnum=numbase;
  381.   base=(long*)malloc(4*(n2+1));*ptbase=base;
  382.   *nbbadprim=0;
  383.   av2=avma;
  384.   i=0;pp=*delta++;qq=2;fl=1;
  385.   while(pp<=n2)
  386.     {
  387.       if((kr=krogs(d,pp))!=-1)
  388.     {
  389.       if(kr) {i++;numbase[pp]=i;base[i]=pp;}
  390.       else
  391.         {
  392.           p1=divis(d,pp);
  393.           if(signe(modis(p1,pp))) {i++;numbase[pp]=i;base[i]= -pp;}
  394.           else
  395.         {
  396.           if(pp==2)
  397.             {
  398.               r=p1[lgef(p1)-1]&7;if(signe(d)<0) r=8-r;
  399.               if(r>=4) {i++;numbase[pp]=i;base[i]= -pp;}
  400.               else badprim[++(*nbbadprim)]=pp;
  401.             }
  402.           else badprim[++(*nbbadprim)]=pp;
  403.         }
  404.         }
  405.     }
  406.       pp+=*delta++;
  407.       if((pp>n)&&fl) {kc=i;*ptkc=kc;fl=0;}
  408.     }
  409.   avma=av2;
  410.   return i;
  411. }
  412.  
  413. long factorise(f,n,limp,ptlo,primfact,expoprimfact,badprim,nbbadprim,base,limhash)
  414.      GEN f;
  415.      long n,limp,*ptlo,*primfact,*expoprimfact,*badprim,nbbadprim,*base,limhash;
  416.  
  417. {
  418.   long sr,i,p,k,fl=1,av=avma,q1,lo;
  419.   GEN x,q,r;
  420.  
  421.   lo=0;
  422.   x=(GEN)(f[1]);if(gcmp1(x)) {*ptlo=0;return 1;}
  423.   for(i=1;(i<=n)&&fl;i++)
  424.     {
  425.       p=base[i];q=dvmdis(x,p,&r);
  426.       if(sr=(!signe(r)))
  427.     {
  428.       primfact[++lo]=p;x=q;k=0;av=avma;
  429.       while(sr)
  430.         {
  431.           k++;q=dvmdis(x,p,&r);
  432.           if(sr=(!signe(r))) {x=q;av=avma;}
  433.         }
  434.       expoprimfact[lo]=k;
  435.     }
  436.       else avma=av;
  437.       fl=(cmpis(q,p)>0);
  438.     }
  439.   if(!fl)
  440.     {
  441.       if(gcmp1(x)) {avma=av;*ptlo=lo;return 1;}
  442.       else
  443.     {
  444.       if(cmpis(x,limp)<=0)
  445.         {
  446.           for(i=1;i<=nbbadprim;i++)
  447.         if(!signe(modis(x,badprim[i]))) {avma=av;*ptlo=lo;return 0;}
  448.           primfact[++lo]=itos(x);expoprimfact[lo]=1;
  449.           avma=av;*ptlo=lo;return 1;
  450.         }
  451.     }
  452.     }
  453.   *ptlo=lo;
  454.   if(cmpis(x,limhash)<=0)
  455.     {
  456.       q1=itos(x);
  457.       avma=av;return q1;
  458.     }
  459.   else {avma=av;return 0;}
  460.   avma=av;return 0;
  461. }
  462.  
  463. long *largeprime2(q1,ex,np,mglob,hashtab)
  464.      long q1,*ex,np,mglob,*hashtab[];
  465. {
  466.   long hashv,*pt,cpt,i,*p1;
  467.   hashv=((q1&2047)-1)>>1;
  468.   pt=hashtab[hashv];cpt=0;
  469.   while(pt&&(q1!=pt[1])) {pt=(long*)(*pt);cpt++;}
  470.   if(!pt)
  471.     {
  472.       if(!(p1=(long*)malloc((mglob+4)<<2))) err(talker,"debordement de pile de hash");
  473.       p1[1]=q1;for(i=2;i<mglob+2;i++) p1[i]=ex[i-1];p1[mglob+2]=np;p1[mglob+3]=0;
  474.       p1[0]=cpt ? (long)hashtab[hashv] : 0;
  475.       hashtab[hashv]=p1;return (long*)0;
  476.     }
  477.   else return (pt[mglob+2]==np) ? (long*)0 : pt;
  478. }
  479.     
  480. int compte(mat,row,longueur,firstnonzero)
  481.      long **mat,row,longueur,*firstnonzero;
  482. {
  483.   int n,j,p;
  484.   n=0;
  485.   for (j=1;j<=longueur;j++)
  486.     {
  487.       p=mat[j][row];
  488.       if(p) {if(abs(p)>=2) return 32767;else {n++;*firstnonzero=j;}}
  489.     }
  490.   return n;
  491. }
  492.  
  493. void hphaseimag1(mat,vectprime,kcco,kc,ptcol,ptrow)
  494.      long **mat,*vectprime,kcco,kc,*ptcol,*ptrow;
  495.  
  496. {
  497.   long numrow,i,i1,j,nb,n,extrarel;
  498.   extrarel=kcco-kc;
  499.   numrow=kc;
  500.   while(numrow>=1)
  501.     {
  502.       for(i=numrow;(i>=1)&&((nb=compte(mat,i,extrarel+numrow,&n))>1);i--);
  503.       if(!nb)
  504.     {
  505.       for(j=1;j<=extrarel+numrow;j++)
  506.         for(i1=i;i1<numrow;i1++) mat[j][i1]=mat[j][i1+1];
  507.       for(i1=i;i1<numrow;i1++) vectprime[i1]=vectprime[i1+1];
  508.       numrow--;extrarel++;
  509.     }
  510.       else if(nb==1)
  511.     {
  512.       for(j=1;j<=extrarel+numrow;j++) mat[j][i]=mat[j][numrow];
  513.       vectprime[i]=vectprime[numrow];
  514.       mat[n]=mat[extrarel+numrow];
  515.       numrow--;
  516.     }
  517.       else {*ptrow=numrow;*ptcol=numrow+extrarel;return;}
  518.     }
  519.   err(talker,"not enough relations or h=1 in hphaseimag1");
  520. }
  521.  
  522. void hphaseimag2(mat,vectprime,nbcol1,nbrow1,ptcol,ptrow)
  523.      long **mat,*vectprime,nbcol1,nbrow1,*ptcol,*ptrow;
  524.  
  525. {
  526.   long numrow,i,j,nb,n,sizemax,imin,nbmin,nmin,s,t,fl=1,extrarel;
  527.  
  528.   extrarel=nbcol1-nbrow1;
  529.   numrow=nbrow1;sizemax=0;
  530.   while(fl&&(numrow>=1)&&(sizemax<=1073741823))
  531.     {
  532.       imin=numrow;nbmin=compte(mat,numrow,extrarel+numrow,&n);nmin=n;
  533.       for(i=numrow-1;i>=1;i--)
  534.     {
  535.       nb=compte(mat,i,extrarel+numrow,&n);
  536.       if(nb<nbmin) {nbmin=nb;imin=i;nmin=n;}
  537.     }
  538.       if(nbmin<32767)
  539.     {
  540.       s=mat[nmin][imin];
  541.       for(j=1;j<nmin;j++)
  542.         {
  543.           sizemax=0;
  544.           if(t=mat[j][imin])
  545.         {
  546.           if(s==t)
  547.             {
  548.               for(i=1;i<=numrow;i++)
  549.             {
  550.               mat[j][i]-=mat[nmin][i];
  551.               sizemax=max(sizemax,abs(mat[j][i]));
  552.             }
  553.             }
  554.           else
  555.             {
  556.               for(i=1;i<=numrow;i++)
  557.             {
  558.               mat[j][i]+=mat[nmin][i];
  559.               sizemax=max(sizemax,abs(mat[j][i]));
  560.             }
  561.             }
  562.         }
  563.         }
  564.       for(j=1;j<=extrarel+numrow;j++) mat[j][imin]=mat[j][numrow];
  565.       vectprime[imin]=vectprime[numrow];
  566.       mat[nmin]=mat[extrarel+numrow];
  567.       numrow--;
  568.     }
  569.       else {*ptrow=numrow;*ptcol=numrow+extrarel;fl=0;}
  570.     }
  571.   if(!numrow) err(talker,"not enough relations or h=1 in hphaseimag2");
  572.   if(sizemax>1073741823) return;
  573.   return;
  574. }
  575.  
  576. GEN hnfimag(x,detmat,detmatsur2)
  577.      GEN x,detmat,detmatsur2;
  578. {
  579.   long li,co,av,tetpil,i,j,def,ldef,lim,dec;
  580.   GEN b,q,w,p1,y,d,u,v,p3,p4;
  581.  
  582.   lim=(avma+bot)>>1;
  583.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  584.   def=co;ldef=(li>co)?li-co+1:1;
  585.   for(i=li-1;i>=ldef;i--)
  586.     {
  587.       def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
  588.       while(j>1)
  589.     {
  590.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  591.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  592.       y[j]=lsub(gmul(p3=divii(coeff(y,i,j),d),y[j-1]),gmul(p4=divii(coeff(y,i,j-1),d),y[j]));
  593.       y[j]=(long)clean(y[j],detmat,detmatsur2);
  594.       y[j-1]=(long)p1;
  595.       y[j-1]=(long)clean(y[j-1],detmat,detmatsur2);
  596.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  597.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  598.     }
  599.       if(j==1)
  600.     {
  601.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  602.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  603.       y[1]=lsub(gmul(p3=divii(coeff(y,i,1),d),y[def]),gmul(p4=divii(coeff(y,i,def),d),y[1]));
  604.       y[1]=(long)clean(y[1],detmat,detmatsur2);
  605.       y[def]=(long)p1;
  606.       y[def]=(long)clean(y[def],detmat,detmatsur2);
  607.     }
  608.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  609.     }
  610.   tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  611.   b=gcopy(detmat);w=cgetg(li,19);def--;
  612.   for(i=li-1;i>=1;i--)
  613.     {
  614.       d=bezout(coeff(y,i,i+def),b,&u,&v);w[i]=lmod(gmul(u,y[i+def]),b);
  615.       if(!signe(coeff(w,i,i))) coeff(w,i,i)=(long)d;
  616.       if(i>1) b=divii(b,d);
  617.     }
  618.   for(i=li-2;i>=1;i--)
  619.     {
  620.       for(j=i+1;j<li;j++)
  621.     {
  622.       q=gdivent(coeff(w,i,j),coeff(w,i,i));w[j]=lsub(w[j],gmul(q,w[i]));
  623.     }
  624.       if(avma<lim) {tetpil=avma;w=gerepile(av,tetpil,gcopy(w));}
  625.     }
  626.   tetpil=avma;w=gerepile(av,tetpil,gcopy(w));
  627.   return w;
  628. }
  629.  
  630. GEN rangi(x,rr)
  631.      GEN  x;
  632.      long *rr;
  633. {
  634.   GEN pass,c,v,det1;
  635.   long i,j,k,rg,t,n,n1,m,m1,nbmot,av,av1,pivprec,piv,p1,p,vi,tetpil,cm=0;
  636.   
  637.   nbmot=20;n=(n1=lg(x))-1;m=(m1=lg(x[1]))-1;av=avma;
  638.   pivprec=lgeti(nbmot);affsi(0,det1=cgeti(nbmot));affsi(1,piv=lgeti(nbmot));pass=cgetg(m1,19);
  639.   for(j=1;j<=m;j++)
  640.     {
  641.       p=pass[j]=lgetg(m1,18);for(i=1;i<=m;i++) ((GEN)p)[i]=lgeti(nbmot);
  642.     }
  643.   c=cgeti(m1);for(k=1;k<=m;k++) c[k]=0;
  644.   v=cgetg(m1,18);av1=avma;k=1;rg=0;
  645.   while((k<=n)&&(rg<m))
  646.     {
  647.       for(t=0,i=1;i<=m;i++)
  648.     if(!c[i])
  649.       {
  650.         vi=lmulii(piv,coeff(x,i,k));
  651.         for(j=1;j<=m;j++)
  652.           if(c[j]) vi=laddii(vi,mulii(coeff(pass,i,j),coeff(x,j,k)));
  653.         v[i]=vi;
  654.         if(!t) if(signe(vi)) t=i;
  655.       }
  656.       if (t)
  657.     {
  658.       rg++;c[t]=k;
  659.       affii(piv,pivprec);affii(v[t],piv);
  660.       if(rg<m)
  661.         {
  662.           for(i=1;i<=m;i++)
  663.         if(!c[i])
  664.           for(j=1;j<=m;j++)
  665.             if(c[j]&&(j!=t))
  666.               {
  667.             p1=lsubii(mulii(piv,coeff(pass,i,j)),mulii(v[i],coeff(pass,t,j)));
  668.             if(rg>1) diviiz(p1,pivprec,coeff(pass,i,j));
  669.             else affii(p1,coeff(pass,i,j));
  670.               }
  671.           for(i=1;i<=m;i++) if(!c[i]) mpnegz(v[i],coeff(pass,i,t));
  672.         }
  673.       else c[t]=0;
  674.     }
  675.       if(rg==m) 
  676.     {
  677.       cm=1;affii(ggcd(piv,det1),det1);
  678.       affii(pivprec,piv);rg--;
  679.     }
  680.       avma=av1;k++;
  681.     }
  682.   *rr=rg+cm;tetpil=avma;return gerepile(av,tetpil,gcopy(det1));
  683. }
  684.  
  685. GEN clean(x,detmat,detmatsur2)
  686.      GEN x,detmat,detmatsur2;
  687. {
  688.   long lx=lg(x),i;
  689.   GEN y,p1;
  690.   y=cgetg(lx,18);
  691.   for(i=1;i<lx;i++)
  692.     {
  693.       p1=modii(x[i],detmat);
  694.       y[i]=(cmpii(p1,detmatsur2)>0) ? lsubii(p1,detmat) : (long)p1;
  695.     }
  696.   return y;
  697. }
  698.  
  699. GEN buchreal(D,gsens,gcbach,gCO)
  700.      GEN D,gsens,gcbach,gCO;
  701. {
  702.   long sens,CO;
  703.   double cbach;
  704.   long precreg,limc,limc2,mglob,m,cp,nbram,nbrederr,auxrel,lo,lo1,ran,fl,fl2;
  705.   long nrho,*fpd,b1,b2,fpc,pp,ep,b,extrarel,c,ic,jc,dec,limhash;
  706.   long kc2,nbcol1,nbcol2,nbcol3,nbrow1,nbrow2,nbrow3,*numbase,*base,*subbase;
  707.   long **mat,*vectprime;
  708.   long primfact[100],expoprimfact[100],primfact1[100],expoprimfact1[100];
  709.   long badprim[100],nbbadprim;
  710.   long av=avma,tetpil,kc,kcco,i,j,*pro,p2,*p1,*ex,q,s,nbtest,mm,av3;
  711.   double drc,lim,logd;
  712.   GEN dr,sqrtD,isqrtD,errglob,log2precis,v,detmat,detmatsur2;
  713.   GEN matgen,vecexpo,vecerr,ermax,**vp,form,form1,forminit,pc,mit,met,mot;
  714.   GEN wecexpo,wecerr,xecexpo,xecerr,matpro,p3,p4; 
  715.   long* hashtab[HASHT];
  716.  
  717.   if((typ(D)!=1)||(typ(gCO)!=1)||(typ(gsens)!=1)) 
  718.     err(talker,"incorrect type in buch");
  719.   if(!signe(D)) err(talker,"zero discriminant in buch");
  720.   if(signe(D)<0) err(talker,"negative discriminant in buchreal, try buchimag");
  721.   s=D[lgef(D)-1]&3;
  722.   if((s==2)||(s==3)) err(talker,"discriminant not 0 or 1 mod 4 in buch");
  723.   if(signe(gCO)<=0) err(talker,"not enough relations in buch");
  724.   CO=itos(gCO);cbach=gtodouble(gcbach);sens=abs(signe(gsens));
  725.   dr=cgetr(3);affir(D,dr);drc=rtodbl(dr);
  726.   lim=sqrt(fabs(drc)/3.0);
  727.   precreg=(gexpo(D)>>5)+5;
  728.   sqrtD=gsqrt(D,precreg);isqrtD=gfloor(sqrtD);
  729.   affsr(1,errglob=cgetr(4));setexpo(errglob,-((precreg-2)<<5));
  730.   log2precis=glog(gdeux,precreg);
  731.   logd=log(fabs(drc));limc=cbach*logd*logd;cp=exp(sqrt(logd*log(logd)/8.0));
  732.   if(cp>limc) limc=cp;limc2=6*logd*logd;if(limc>limc2) limc2=limc;
  733.   subbase=subfactorbasereal(D,&v,lim,precreg,&nbram,isqrtD,sqrtD,sens);
  734.   mglob=m=subbase[0];
  735.   vp=powsubfactreal(v,m,C+7,D,isqrtD,sqrtD,sens,precreg);
  736.   kc2=factorbase(D,limc2,limc,&numbase,&base,&kc,badprim,&nbbadprim);
  737.   pro=(long*)malloc(4*(kc2+1));vectprime=(long*)malloc(4*(kc2+1));
  738.   for(i=1;i<=kc2;i++)
  739.     {
  740.       p2=vectprime[i]=base[i];
  741.       pro[i]=(p2>0) ? p2 : -p2;
  742.     }
  743.   free(base);base=pro;
  744.   kcco=kc+CO;
  745.   mat=(long **)malloc(4*(kcco+1));
  746.   vecexpo=cgetg(kcco+1,17);vecerr=cgetg(kcco+1,17);
  747.   for(i=1;i<=kcco;i++) {vecexpo[i]=lgetr(precreg);vecerr[i]=lgetr(4);}
  748.   for(i=1;i<=kcco;i++)
  749.     {
  750.       p1=(long *)malloc(4*(kc+1));mat[i]=p1;
  751.       for(j=1;j<=kc;j++) p1[j]=0;
  752.     }
  753.   ex=(long*)malloc((m+1)*4);
  754.   q=31-(long)ceil(log((double)C)/log(2.0));
  755.   s=0;nbtest=0;
  756.   mm=m+nbram+CO;
  757.   limhash=(limc<32767)?limc*limc:(2<<30);
  758.   for(i=0;i<HASHT;i++) hashtab[i]=(long*)0;
  759.   auxrel=0;
  760.   while(s<mm)
  761.     {
  762.       for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  763.       av3=avma;nbrederr=0;
  764.       form=vp[1][ex[1]];
  765.       for (i=2;i<=m;i++) 
  766.     form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
  767.       fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  768.       nbtest++;
  769.       if(fpc>1)
  770.     {
  771.       fpd=largeprime(fpc,ex,0,0,mglob,hashtab);
  772.       if(fpd)
  773.         {
  774.           auxrel++;
  775.           s++;lo1=lo;
  776.           form=vp[1][ex[1]];
  777.           for (i=2;i<=m;i++) 
  778.         form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
  779.           for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
  780.           form1=vp[1][fpd[2]];
  781.           for(i=2;i<=m;i++) 
  782.         form1=comprealform5(form1,vp[i][fpd[i+1]],&nbrederr,D,isqrtD,sqrtD,sens);
  783.           b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
  784.           factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  785.           if(b1==b2)
  786.         {
  787.           for(i=1;i<=m;i++)
  788.             mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
  789.           for(j=1;j<=lo;j++)
  790.             {
  791.               pp=primfact[j];ep=expoprimfact[j];
  792.               b1=itos(modis(form1[2],(pp<<1)));
  793.               if(b1>pp) ep= -ep;
  794.               mat[s][numbase[pp]]-=ep;
  795.             }
  796.           for(j=1;j<=lo1;j++)
  797.             {
  798.               pp=primfact1[j];ep=expoprimfact1[j];
  799.               b1=itos(modis(form[2],(pp<<1)));
  800.               if(b1>pp) ep= -ep;
  801.               mat[s][numbase[pp]]+=ep;
  802.             }
  803.           affrr(shiftr(mpadd(mulir(mulsi(EXP220,subii(form[4],form1[4])),log2precis),
  804.                      mplog(absr(divrr(form[5],form1[5])))),-1),vecexpo[s]);
  805.           mulsrz(nbrederr,errglob,vecerr[s]);
  806.         }
  807.           else
  808.         {
  809.           if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
  810.           else
  811.             {
  812.               for(i=1;i<=m;i++)
  813.             mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
  814.               for(j=1;j<=lo;j++)
  815.             {
  816.               pp=primfact[j];ep=expoprimfact[j];
  817.               b1=itos(modis(form1[2],(pp<<1)));
  818.               if(b1>pp) ep= -ep;
  819.               mat[s][numbase[pp]]+=ep;
  820.             }
  821.               for(j=1;j<=lo1;j++)
  822.             {
  823.               pp=primfact1[j];ep=expoprimfact1[j];
  824.               b1=itos(modis(form[2],(pp<<1)));
  825.               if(b1>pp) ep= -ep;
  826.               mat[s][numbase[pp]]+=ep;
  827.             }
  828.               affrr(shiftr(mpadd(mulir(mulsi(EXP220,addii(form1[4],form[4])),log2precis),mplog(absr(mulrr(form1[5],form[5])))),-1),vecexpo[s]);
  829.               mulsrz(nbrederr,errglob,vecerr[s]);
  830.             }
  831.         }
  832.         }
  833.     }      
  834.       if(fpc==1)
  835.     {
  836.       s++;
  837.       form=vp[1][ex[1]];
  838.       for(i=2;i<=m;i++) 
  839.         form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
  840.       for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
  841.       for(j=1;j<=lo;j++)
  842.         {
  843.           pp=primfact[j];ep=expoprimfact[j];
  844.           b=itos(modis(form[2],(pp<<1)));
  845.           if(b>pp) ep= -ep;
  846.           mat[s][numbase[pp]]+=ep;
  847.         }
  848.       affrr(shiftr(mpadd(mulir(mulsi(EXP220,form[4]),log2precis),mplog(absr(form[5]))),-1),vecexpo[s]);
  849.       mulsrz(nbrederr,errglob,vecerr[s]);
  850.     }
  851.       avma=av3;
  852.     }
  853.   while(s<kcco)
  854.     {
  855.       for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  856.       av3=avma;
  857.       form=vp[1][ex[1]];
  858.       for (i=2;i<=m;i++) form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
  859.       pc=redrealform3(primeform(D,stoi(base[s+1-CO]),precreg),D,isqrtD,sens);
  860.       form=comprealform3(form,pc,D,isqrtD,sens);forminit=form;nrho=0;
  861.       fl2=0;
  862.       do
  863.     {
  864.       fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  865.       nbtest++;
  866.       if(fpc>1)
  867.         {
  868.           fpd=largeprime(fpc,ex,s+1-CO,nrho,mglob,hashtab);
  869.           if(fpd)
  870.         {
  871.           form=vp[1][ex[1]];
  872.           for (i=2;i<=m;i++) 
  873.             form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
  874.           pc=redrealform(primeform(D,stoi(base[s+1-CO]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
  875.           pc[4]=zero;affsr(1,pc[5]);
  876.           form=comprealform5(form,pc,&nbrederr,D,isqrtD,sqrtD,sens);
  877.           for(i=1;i<=nrho;i++) 
  878.             form=rhorealform5(form,&nbrederr,D,isqrtD,sqrtD);
  879.           auxrel++;
  880.           s++;lo1=lo;
  881.           for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
  882.           form1=vp[1][fpd[2]];
  883.           for(i=2;i<=m;i++) 
  884.             form1=comprealform5(form1,vp[i][fpd[i+1]],&nbrederr,D,isqrtD,sqrtD,sens);
  885.           if(fpd[m+2])
  886.             {
  887.               pc=redrealform(primeform(D,stoi(base[fpd[m+2]]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
  888.               pc[4]=zero;affsr(1,pc[5]);
  889.               form1=comprealform5(form1,pc,&nbrederr,D,isqrtD,sqrtD,sens);
  890.               for(i=1;i<=fpd[m+3];i++) 
  891.             form1=rhorealform5(form1,&nbrederr,D,isqrtD,sqrtD);
  892.               if((!sens)&&(signe(addii(form1[1],form1[3]))))
  893.             {setsigne(form1[1],1);setsigne(form1[3],-1);}
  894.             }
  895.           b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
  896.           factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  897.           if(b1==b2)
  898.             {
  899.               fl2=1;
  900.               for(i=1;i<=m;i++)
  901.             mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
  902.               mat[s][s-CO]= -1;
  903.               if(fpd[m+2]) mat[s][fpd[m+2]]++;
  904.               for(j=1;j<=lo;j++)
  905.             {
  906.               pp=primfact[j];ep=expoprimfact[j];
  907.               b1=itos(modis(form1[2],(pp<<1)));
  908.               if(b1>pp) ep= -ep;
  909.               mat[s][numbase[pp]]-=ep;
  910.             }
  911.               for(j=1;j<=lo1;j++)
  912.             {
  913.               pp=primfact1[j];ep=expoprimfact1[j];
  914.               b1=itos(modis(form[2],(pp<<1)));
  915.               if(b1>pp) ep= -ep;
  916.               mat[s][numbase[pp]]+=ep;
  917.             }
  918.               affrr(shiftr(mpadd(mulir(mulsi(EXP220,subii(form[4],form1[4])),log2precis),mplog(absr(divrr(form[5],form1[5])))),-1),vecexpo[s]);
  919.               mulsrz(nbrederr,errglob,vecerr[s]);
  920.             }
  921.           else
  922.             {
  923.               if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
  924.               else
  925.             {
  926.               fl2=1;
  927.               for(i=1;i<=m;i++)
  928.                 mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
  929.               mat[s][s-CO]= -1;
  930.               if(fpd[m+2]) mat[s][fpd[m+2]]--;
  931.               for(j=1;j<=lo;j++)
  932.                 {
  933.                   pp=primfact[j];ep=expoprimfact[j];
  934.                   b1=itos(modis(form1[2],(pp<<1)));
  935.                   if(b1>pp) ep= -ep;
  936.                   mat[s][numbase[pp]]+=ep;
  937.                 }
  938.               for(j=1;j<=lo1;j++)
  939.                 {
  940.                   pp=primfact1[j];ep=expoprimfact1[j];
  941.                   b1=itos(modis(form[2],(pp<<1)));
  942.                   if(b1>pp) ep= -ep;
  943.                   mat[s][numbase[pp]]+=ep;
  944.                 }
  945.               affrr(shiftr(mpadd(mulir(mulsi(EXP220,addii(form1[4],form[4])),log2precis),mplog(absr(mulrr(form1[5],form[5])))),-1),vecexpo[s]);
  946.               mulsrz(nbrederr,errglob,vecerr[s]);
  947.             }
  948.             }
  949.         }
  950.         }
  951.       if ((fpc!=1)&&(!fl2))
  952.         {
  953.           form=rhorealform3(form,D,isqrtD);nrho++;
  954.           if((sens)||(!signe(addii(form[1],form[3]))))
  955.         {form=rhorealform3(form,D,isqrtD);nrho++;}
  956.           else
  957.         {
  958.           setsigne(form[1],1);setsigne(form[3],-1);
  959.         }
  960.           fl=gegal(form[1],forminit[1]);
  961.           if(fl) fl=gegal(form[2],forminit[2]);
  962.           if(fl) fl=gegal(form[3],forminit[3]);
  963.         }
  964.     }
  965.       while((fpc!=1)&&(!fl)&&(!fl2));
  966.       if(fpc==1)
  967.     {
  968.       nbrederr=0;
  969.       form=vp[1][ex[1]];
  970.       for (i=2;i<=m;i++) 
  971.         form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
  972.       pc=redrealform(primeform(D,stoi(base[s+1-CO]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
  973.       pc[4]=zero;affsr(1,pc[5]);
  974.       form=comprealform5(form,pc,&nbrederr,D,isqrtD,sqrtD,sens);
  975.       for(i=1;i<=nrho;i++) 
  976.         form=rhorealform5(form,&nbrederr,D,isqrtD,sqrtD);
  977.       s++;
  978.       for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
  979.       mat[s][s-CO]= -1;
  980.       for(j=1;j<=lo;j++)
  981.         {
  982.           pp=primfact[j];ep=expoprimfact[j];
  983.           b=itos(modis(form[2],(pp<<1)));
  984.           if(b>pp) ep= -ep;
  985.           mat[s][numbase[pp]]+=ep;
  986.         }
  987.       affrr(shiftr(mpadd(mulir(mulsi(EXP220,form[4]),log2precis),mplog(absr(form[5]))),-1),vecexpo[s]);
  988.       mulsrz(nbrederr,errglob,vecerr[s]);
  989.       if(!(mat[s][s-CO])) {for(i=1;i<=kc;i++) mat[s][i]=0;s--;}
  990.     }
  991.       avma=av3;
  992.     }
  993.   nbtest=auxrel=0;
  994.   while(s<kc2)
  995.     {
  996.       for (i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
  997.       av3=avma;
  998.       form=vp[1][ex[1]];
  999.       for (i=2;i<=m;i++) form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
  1000.       pp=base[s+1];
  1001.       pc=redrealform3(primeform(D,stoi(pp),precreg),D,isqrtD,sens);
  1002.       form=comprealform3(form,pc,D,isqrtD,sens);forminit=form;fl2=0;
  1003.       do
  1004.     {
  1005.       fpc=factorise(form,s,pp,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
  1006.       nbtest++;
  1007.       if(fpc>1)
  1008.         {
  1009.           fpd=largeprime2(fpc,ex,s+1,mglob,hashtab);
  1010.           if(fpd&&(fpc!=pp)) {auxrel++;s++;fl2=1;}
  1011.         }
  1012.       if((fpc!=1)&&(!fl2))
  1013.         {
  1014.           form=rhorealform3(form,D,isqrtD);
  1015.           if((sens)||(!signe(addii(form[1],form[3])))) 
  1016.         form=rhorealform3(form,D,isqrtD);
  1017.           else
  1018.         {
  1019.           setsigne(form[1],1);setsigne(form[3],-1);
  1020.         }
  1021.           fl=gegal(form[1],forminit[1]);
  1022.           if(fl) fl=gegal(form[2],forminit[2]);
  1023.           if(fl) fl=gegal(form[3],forminit[3]);
  1024.         }
  1025.     }
  1026.       while((fpc!=1)&&(!fl)&&(!fl2));
  1027.       if(fpc==1) s++;
  1028.       avma=av3;
  1029.     }
  1030.   hphasereal1(mat,vectprime,vecexpo,vecerr,kcco,kc,&nbcol1,&nbrow1);
  1031.   hphasereal2(mat,vectprime,vecexpo,vecerr,nbcol1,nbrow1,&nbcol2,&nbrow2);
  1032.   nbcol3=nbcol2;nbrow3=nbrow2;
  1033.   matgen=cgetg(nbcol3+1,19);
  1034.   for(i=1;i<=nbcol3;i++) matgen[i]=lgetg(nbrow3+1,18);
  1035.   wecexpo=cgetg(nbcol3+1,17);wecerr=cgetg(nbcol3+1,17);
  1036.   for(i=1;i<=nbcol3;i++)
  1037.     {
  1038.       wecexpo[i]=vecexpo[i];wecerr[i]=vecerr[i];
  1039.       for(j=1;j<=nbrow3;j++) ((GEN)(matgen[i]))[j]=lstoi(mat[i][j]);
  1040.     }
  1041.   vecexpo=wecexpo;vecerr=wecerr;
  1042.   detmat=absi(rangi(matgen,&ran));
  1043.   detmatsur2=shifti(detmat,-1);
  1044.   matpro=kerint(matgen);
  1045.   xecexpo=gmul(vecexpo,matpro);
  1046.   xecerr=gmul(vecerr,gabs(matpro));
  1047.   if(gcmp1(absi(detmat))) {mit=idmat(1);nbrow3=1;nbcol3=lg(matpro);}
  1048.   else mit=hnfreal(matgen,detmat,detmatsur2,vecexpo,vecerr,&wecexpo,&wecerr);
  1049.   wecexpo=xecexpo;wecerr=xecerr;
  1050.   extrarel=nbcol3-nbrow3;
  1051.   for(i=2;i<=extrarel;i++)
  1052.     wecexpo[i]=(long)gcdreal(wecexpo[i-1],wecexpo[i],wecerr[i-1],wecerr[i],wecexpo+i-1);
  1053.   ermax=(GEN)wecexpo[1];for(i=2;i<extrarel;i++) ermax=gmax(ermax,wecexpo[i]);
  1054.   p1=gun;
  1055.   for(i=1;i<=nbrow3;i++)
  1056.     {
  1057.       p3=(GEN)coeff(mit,i,i);
  1058.       if(signe(p3)) p1=mulii(p1,p3);
  1059.     }
  1060.   c=0;for(i=1;i<=nbrow3;i++) if(!gcmp1(coeff(mit,i,i))) c++;
  1061.   if(!c) met=idmat(1);
  1062.   else
  1063.     {
  1064.       met=cgetg(c+1,19);ic=0;
  1065.       for(i=1;i<=nbrow3;i++)
  1066.     {
  1067.       if(!gcmp1(coeff(mit,i,i)))
  1068.         {
  1069.           ic++;p3=cgetg(c+1,18);met[ic]=(long)p3;jc=0;
  1070.           for(j=1;j<=nbrow3;j++)
  1071.         {
  1072.           if(!gcmp1(coeff(mit,j,j)))
  1073.             {
  1074.               jc++;p3[jc]=coeff(mit,j,i);
  1075.             }
  1076.         }
  1077.         }
  1078.     }
  1079.     }
  1080.   met=smith(met);
  1081.   lo=lg(met)-1;
  1082.   c=0;for(i=1;i<=lo;i++) if(!gcmp1(met[i])) c++;
  1083.   p2=(long)absr(wecexpo[extrarel]);
  1084.   p3=gdiv(gmul2n(gmul(p1,p2),1),gmul(lfunc(D),sqrtD));
  1085.   tetpil=avma;p4=cgetg(6,17);p4[1]=lcopy(p1);mot=cgetg(c+1,17);p4[2]=(long)mot;
  1086.   for(i=1;i<=c;i++) mot[i]=lcopy(met[i]);
  1087.   p4[3]=lcopy(p2);p4[4]=lcopy(ermax);p4[5]=lcopy(p3);
  1088.   free(base);free(vectprime);free(ex);free(subbase);free(numbase);
  1089.   for(i=1;i<=kcco;i++) free(mat[i]);free(mat);
  1090.   for(i=1;i<=m;i++) free(vp[i]);free(vp);
  1091.   return gerepile(av,tetpil,p4);
  1092. }
  1093.  
  1094. long *subfactorbasereal(d,w,ll,precreg,ptnbram,isqrtD,sqrtD,sens)
  1095.      long precreg,*ptnbram,sens;
  1096.      double ll;
  1097.      GEN d,*w;
  1098.      GEN isqrtD,sqrtD;
  1099. {
  1100.   byteptr delta=diffptr;
  1101.   long av1,i,j,pp,tp,fl,*subbase,pro[100],toto=0;
  1102.   double prod;
  1103.   GEN p1;
  1104.  
  1105.   av1=avma;
  1106.   i=1;pp=0;*ptnbram=0;fl=1;prod=1;
  1107.   while(fl||(i==1))
  1108.     {
  1109.       pp+=*delta++;tp=krogs(d,pp);
  1110.       if(tp==1) {pro[i]=pp;i++;prod*=pp;} else if(!tp) (*ptnbram)++;
  1111.       if((prod>ll)&&(i>1)) fl=0;
  1112.     }
  1113.   avma=av1;i--;
  1114.   *w=cgetg(i+1,18);
  1115.   for(j=1;j<=i;j++) 
  1116.     {
  1117.       p1=redrealform(primeform(d,stoi(pro[j]),precreg),&toto,d,isqrtD,sqrtD,sens,precreg);
  1118.       p1[4]=zero;affsr(1,p1[5]);
  1119.       (*w)[j]=(long)p1;
  1120.     }
  1121.   subbase=(long*)malloc(4*(i+1));subbase[0]=i;
  1122.   for(j=1;j<=i;j++) subbase[j]=pro[j];
  1123.   return subbase;
  1124. }
  1125.   
  1126. GEN **powsubfactreal(w,n,a,D,isqrtD,sqrtD,sens,precreg)
  1127.      GEN w,D,isqrtD,sqrtD;
  1128.      long n,a,sens,precreg;
  1129. {
  1130.   long i,j,toto=0;
  1131.   GEN **x;
  1132.  
  1133.   x=(GEN**)malloc(4*(n+1));
  1134.   for(i=1;i<=n;i++) x[i]=(GEN*)malloc(4*(a+1));
  1135.   for(i=1;i<=n;i++)
  1136.     {
  1137.       j=0;x[i][j]=powrealform5(w[1],0,D,isqrtD,sqrtD,sens,precreg);
  1138.       while(j<=a)
  1139.     {
  1140.       j++;
  1141.       x[i][j]=comprealform5(x[i][j-1],w[i],&toto,D,isqrtD,sqrtD,sens);
  1142.     }
  1143.     }
  1144.   fflush(stdout);
  1145.   return x;
  1146. }
  1147.  
  1148. long *largeprime(q1,ex,np,nrho,mglob,hashtab)
  1149.      long q1,*ex,np,nrho,mglob,*hashtab[];
  1150. {
  1151.   long hashv,*pt,cpt,i,*p1,fl;
  1152.   hashv=((q1&2047)-1)>>1;
  1153.   pt=hashtab[hashv];cpt=0;
  1154.   while(pt&&(q1!=pt[1])) {pt=(long*)(*pt);cpt++;}
  1155.   if(!pt)
  1156.     {
  1157.       if(!(p1=(long*)malloc((mglob+4)<<2))) err(talker,"debordement de pile de hash");
  1158.       p1[1]=q1;for(i=2;i<mglob+2;i++) p1[i]=ex[i-1];p1[mglob+2]=np;p1[mglob+3]=nrho;
  1159.       p1[0]=cpt ? (long)hashtab[hashv] : 0;
  1160.       hashtab[hashv]=p1;return (long*)0;
  1161.     }
  1162.   else
  1163.     {
  1164.       fl=1;i=2;while(fl&&(i<(mglob+2))) {fl=(pt[i]==ex[i-1]);i++;}
  1165.       if(fl) fl=(pt[i]==np);
  1166.       return fl ? (long*)0 : pt;
  1167.     }
  1168. }
  1169.  
  1170. void hphasereal1(mat,vectprime,vecexpo,vecerr,kcco,kc,ptcol,ptrow)
  1171.      long **mat,*vectprime,kcco,kc,*ptcol,*ptrow;
  1172.      GEN vecexpo,vecerr;
  1173.  
  1174. {
  1175.   long numrow,i,i1,j,nb,n,extrarel;
  1176.   extrarel=kcco-kc;
  1177.   numrow=kc;
  1178.   while(numrow>=1)
  1179.     {
  1180.       for(i=numrow;(i>=1)&&((nb=compte(mat,i,extrarel+numrow,&n))>1);i--);
  1181.       if(!nb)
  1182.     {
  1183.       for(j=1;j<=extrarel+numrow;j++)
  1184.         for(i1=i;i1<numrow;i1++) mat[j][i1]=mat[j][i1+1];
  1185.       for(i1=i;i1<numrow;i1++) vectprime[i1]=vectprime[i1+1];
  1186.       numrow--;extrarel++;
  1187.     }
  1188.       else if(nb==1)
  1189.     {
  1190.       for(j=1;j<=extrarel+numrow;j++) mat[j][i]=mat[j][numrow];
  1191.       vectprime[i]=vectprime[numrow];
  1192.       mat[n]=mat[extrarel+numrow];
  1193.       vecexpo[n]=vecexpo[extrarel+numrow];
  1194.       vecerr[n]=vecerr[extrarel+numrow];
  1195.       numrow--;
  1196.     }
  1197.       else {*ptrow=numrow;*ptcol=numrow+extrarel;return;}
  1198.     }
  1199.   err(talker,"not enough relations or h=1 in hphasereal1");
  1200. }
  1201.  
  1202. void hphasereal2(mat,vectprime,vecexpo,vecerr,nbcol1,nbrow1,ptcol,ptrow)
  1203.      long **mat,*vectprime,nbcol1,nbrow1,*ptcol,*ptrow;
  1204.      GEN vecexpo,vecerr;
  1205.  
  1206. {
  1207.   long numrow,i,j,nb,n,sizemax,imin,nbmin,nmin,s,t,fl=1,extrarel;
  1208.  
  1209.   extrarel=nbcol1-nbrow1;
  1210.   numrow=nbrow1;sizemax=0;
  1211.   while(fl&&(numrow>=1)&&(sizemax<=1073741823))
  1212.     {
  1213.       imin=numrow;nbmin=compte(mat,numrow,extrarel+numrow,&n);nmin=n;
  1214.       for(i=numrow-1;i>=1;i--)
  1215.     {
  1216.       nb=compte(mat,i,extrarel+numrow,&n);
  1217.       if(nb<nbmin) {nbmin=nb;imin=i;nmin=n;}
  1218.     }
  1219.       if(nbmin<32767)
  1220.     {
  1221.       s=mat[nmin][imin];
  1222.       for(j=1;j<nmin;j++)
  1223.         {
  1224.           sizemax=0;
  1225.           if(t=mat[j][imin])
  1226.         {
  1227.           if(s==t)
  1228.             {
  1229.               for(i=1;i<=numrow;i++)
  1230.             {
  1231.               mat[j][i]-=mat[nmin][i];
  1232.               sizemax=max(sizemax,abs(mat[j][i]));
  1233.             }
  1234.               subrrz(vecexpo[j],vecexpo[nmin],vecexpo[j]);
  1235.               addrrz(vecerr[j],vecerr[nmin],vecerr[j]);
  1236.             }
  1237.           else
  1238.             {
  1239.               for(i=1;i<=numrow;i++)
  1240.             {
  1241.               mat[j][i]+=mat[nmin][i];
  1242.               sizemax=max(sizemax,abs(mat[j][i]));
  1243.             }
  1244.               addrrz(vecexpo[j],vecexpo[nmin],vecexpo[j]);
  1245.               addrrz(vecerr[j],vecerr[nmin],vecerr[j]);
  1246.             }
  1247.         }
  1248.         }
  1249.       for(j=1;j<=extrarel+numrow;j++) mat[j][imin]=mat[j][numrow];
  1250.       vectprime[imin]=vectprime[numrow];
  1251.       mat[nmin]=mat[extrarel+numrow];
  1252.       vecexpo[nmin]=vecexpo[extrarel+numrow];
  1253.       vecerr[nmin]=vecerr[extrarel+numrow];
  1254.       numrow--;
  1255.     }
  1256.       else {*ptrow=numrow;*ptcol=numrow+extrarel;fl=0;}
  1257.     }
  1258.   if(!numrow) err(talker,"not enough relations or h=1 in hphasereal2");
  1259.   if(sizemax>1073741823) return;
  1260.   return;
  1261. }
  1262.  
  1263. GEN gcdreal(a,b,era,erb,erd)
  1264.      GEN a,b,era,erb,*erd;
  1265. {
  1266.   long av,tetpil,dec,e;
  1267.   GEN k1,r;
  1268.  
  1269.   if(expo(a)<0) {*erd=gcopy(a);return gcopy(b);}
  1270.   if(expo(b)<0) {*erd=gcopy(b);return gcopy(a);}
  1271.   av=avma;a=absr(a);b=absr(b);
  1272.   while((expo(b)>=0)&&(signe(b)))
  1273.     {
  1274.       k1=gcvtoi(divrr(a,b),&e);
  1275.       r=subrr(a,mulir(k1,b));a=b;b=r;
  1276.     }
  1277.   tetpil=avma;a=gcopy(a);b=gcopy(b);dec=lpile(av,tetpil,0)>>2;
  1278.   *erd=b+dec;return a+dec;
  1279. }
  1280.  
  1281. GEN hnfreal(x,detmat,detmatsur2,vlog,verr,wlog,werr)
  1282.      GEN x,detmat,detmatsur2,vlog,verr,*wlog,*werr;
  1283. {
  1284.   long li,co,av,tetpil,i,j,def,ldef,lim,dec;
  1285.   GEN b,q,w,p1,y,d,u,v,p3,p4;
  1286.  
  1287.   lim=(avma+bot)>>1;
  1288.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  1289.   def=co;ldef=(li>co)?li-co+1:1;
  1290.   for(i=li-1;i>=ldef;i--)
  1291.     {
  1292.       def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
  1293.       while(j>1)
  1294.     {
  1295.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  1296.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  1297.       y[j]=lsub(gmul(p3=divii(coeff(y,i,j),d),y[j-1]),gmul(p4=divii(coeff(y,i,j-1),d),y[j]));
  1298.       y[j]=(long)clean(y[j],detmat,detmatsur2);
  1299.       y[j-1]=(long)p1;
  1300.       y[j-1]=(long)clean(y[j-1],detmat,detmatsur2);
  1301.       p1=gadd(gmul(u,vlog[j]),gmul(v,vlog[j-1]));
  1302.       vlog[j]=lsub(gmul(p3,vlog[j-1]),gmul(p4,vlog[j]));
  1303.       vlog[j-1]=(long)p1;
  1304.       p1=gadd(gmul(absi(u),verr[j]),gmul(absi(v),verr[j-1]));
  1305.       verr[j]=ladd(gmul(gabs(p3),verr[j-1]),gmul(gabs(p4),verr[j]));
  1306.       verr[j-1]=(long)p1;
  1307.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  1308.       if(avma<lim)
  1309.         {
  1310.           tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
  1311.           dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
  1312.         }
  1313.     }
  1314.       if(j==1)
  1315.     {
  1316.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  1317.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  1318.       y[1]=lsub(gmul(p3=divii(coeff(y,i,1),d),y[def]),gmul(p4=divii(coeff(y,i,def),d),y[1]));
  1319.       y[1]=(long)clean(y[1],detmat,detmatsur2);
  1320.       y[def]=(long)p1;
  1321.       y[def]=(long)clean(y[def],detmat,detmatsur2);
  1322.       p1=gadd(gmul(u,vlog[1]),gmul(v,vlog[def]));
  1323.       vlog[1]=lsub(gmul(p3,vlog[def]),gmul(p4,vlog[1]));
  1324.       vlog[def]=(long)p1;
  1325.       p1=gadd(gmul(absi(u),verr[1]),gmul(absi(v),verr[def]));
  1326.       verr[1]=ladd(gmul(p3,verr[def]),gmul(p4,verr[1]));
  1327.       verr[def]=(long)p1;
  1328.     }
  1329.       if(avma<lim)
  1330.     {
  1331.       tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
  1332.       dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
  1333.     }
  1334.     }
  1335.   tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
  1336.   dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
  1337.   b=gcopy(detmat);w=cgetg(li,19);def--;
  1338.   for(i=li-1;i>=1;i--)
  1339.     {
  1340.       d=bezout(coeff(y,i,i+def),b,&u,&v);w[i]=lmod(gmul(u,y[i+def]),b);
  1341.       if(!signe(coeff(w,i,i))) coeff(w,i,i)=(long)d;
  1342.       vlog[i+def]=lmul(u,vlog[i+def]);verr[i+def]=labs(gmul(u,verr[i+def]));
  1343.       if(i>1) b=divii(b,d);
  1344.     }
  1345.   for(i=li-2;i>=1;i--)
  1346.     {
  1347.       for(j=i+1;j<li;j++)
  1348.     {
  1349.       q=gdivent(coeff(w,i,j),coeff(w,i,i));w[j]=lsub(w[j],gmul(q,w[i]));
  1350.       vlog[j+def]=lsub(vlog[j+def],gmul(q,vlog[i+def]));
  1351.       verr[j+def]=ladd(verr[j+def],gabs(gmul(q,verr[i+def])));
  1352.     }
  1353.       if(avma<lim)
  1354.     {
  1355.       tetpil=avma;w=gcopy(w);vlog=gcopy(vlog);verr=gcopy(verr);
  1356.       dec=lpile(av,tetpil,0)>>2;w+=dec;vlog+=dec;verr+=dec;
  1357.     }
  1358.     }
  1359.   tetpil=avma;w=gcopy(w);vlog=gcopy(vlog);verr=gcopy(verr);
  1360.   dec=lpile(av,tetpil,0)>>2;w+=dec;vlog+=dec;verr+=dec;
  1361.   *wlog=vlog;*werr=verr;return w;
  1362. }
  1363.  
  1364. GEN comprealform3(x,y,D,isqrtD,sens)
  1365.      GEN x,y,D,isqrtD;
  1366.      long sens;
  1367.      
  1368. {
  1369.   long av,tetpil;
  1370.   GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,b3,c3,m,z,p1,r;
  1371.   
  1372.   av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  1373.   d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
  1374.   v1=divii(x[1],d1);v2=divii(y[1],d1);
  1375.   m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
  1376.   r=modii(m,v1);b3=shifti((p1=mulii(v2,r)),1);
  1377.   c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
  1378.   z=cgetg(4,17);z[1]=lmulii(v1,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v1);
  1379.   tetpil=avma;return gerepile(av,tetpil,redrealform3(z,D,isqrtD,sens));
  1380. }
  1381.  
  1382. GEN comprealform5(x,y,pterr,D,isqrtD,sqrtD,sens)
  1383.      GEN x,y,D,isqrtD,sqrtD;
  1384.      long *pterr,sens;
  1385.      
  1386. {
  1387.   long av,tetpil,ss;
  1388.   GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,b3,c3,m,z,p1,r;
  1389.   
  1390.   av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
  1391.   d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
  1392.   v1=divii(x[1],d1);v2=divii(y[1],d1);
  1393.   m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
  1394.   r=modii(m,v1);b3=shifti((p1=mulii(v2,r)),1);
  1395.   c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
  1396.   z=cgetg(6,17);z[1]=lmulii(v1,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v1);
  1397.   z[5]=lmulrr(x[5],y[5]);(*pterr)++;
  1398.   if((ss=expo(z[5]))>=EXP220) {z[4]=laddii(addsi(1,x[4]),y[4]);setexpo(z[5],ss-EXP220);}
  1399.   else z[4]=laddii(x[4],y[4]);
  1400.   tetpil=avma;return gerepile(av,tetpil,redrealform5(z,pterr,D,isqrtD,sqrtD,sens));
  1401. }
  1402.  
  1403. GEN sqrealform3(x,D,isqrtD,sens)
  1404.      GEN x,D,isqrtD;
  1405.      long sens;
  1406.      
  1407. {
  1408.   long av,tetpil;
  1409.   GEN d1,x2,y2,v1,b3,c3,m,z,p1,r;
  1410.   
  1411.   av=avma;
  1412.   d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
  1413.   m=mulii(x[3],x2);setsigne(m,-signe(m));
  1414.   r=modii(m,v1);b3=shifti((p1=mulii(v1,r)),1);
  1415.   c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
  1416.   z=cgetg(4,17);z[1]=lmulii(v1,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v1);
  1417.   tetpil=avma;return gerepile(av,tetpil,redrealform3(z,D,isqrtD,sens));
  1418. }
  1419.  
  1420. GEN sqrealform5(x,pterr,D,isqrtD,sqrtD,sens)
  1421.      GEN x,D,isqrtD,sqrtD;
  1422.      long *pterr,sens;
  1423. {
  1424.   long av,tetpil,ss;
  1425.   GEN d1,x2,y2,v1,b3,c3,m,z,p1,r;
  1426.   
  1427.   av=avma;
  1428.   d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
  1429.   m=mulii(x[3],x2);setsigne(m,-signe(m));
  1430.   r=modii(m,v1);b3=shifti((p1=mulii(v1,r)),1);
  1431.   c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
  1432.   z=cgetg(6,17);z[1]=lmulii(v1,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v1);
  1433.   z[5]=lmulrr(x[5],x[5]);(*pterr)++;
  1434.   if((ss=expo(z[5]))>=EXP220) {z[4]=laddii(addsi(1,x[4]),x[4]);setexpo(z[5],ss-EXP220);}
  1435.   else z[4]=lshifti(x[4],1);
  1436.   tetpil=avma;return gerepile(av,tetpil,redrealform5(z,pterr,D,isqrtD,sqrtD,sens));
  1437. }
  1438.  
  1439. GEN rhorealform3(x,D,isqrtD)
  1440.      GEN x,D,isqrtD;
  1441. {
  1442.   long av,tetpil,s;
  1443.   GEN y,p1,nn;
  1444.   
  1445.   av=avma;y=cgetg(4,17);y[1]=lcopy(x[3]);
  1446.   s=signe(y[1]);setsigne(y[1],1);
  1447.   if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
  1448.   else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
  1449.   p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
  1450.   setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
  1451.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1452. }
  1453.  
  1454. GEN rhorealform5(x,pterr,D,isqrtD,sqrtD)
  1455.      GEN x,D,isqrtD,sqrtD;
  1456.      long *pterr;
  1457. {
  1458.   long av,tetpil,s,ss;
  1459.   GEN y,p1,nn;
  1460.   
  1461.   av=avma;y=cgetg(6,17);y[1]=lcopy(x[3]);
  1462.   s=signe(y[1]);setsigne(y[1],1);
  1463.   if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
  1464.   else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
  1465.   p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
  1466.   setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
  1467.   y[5]=lmulrr(divrr(addir(x[2],sqrtD),subir(x[2],sqrtD)),x[5]);(*pterr)++;
  1468.   if((ss=expo(y[5]))>=EXP220) {y[4]=laddsi(1,x[4]);y[5]=lshiftr(y[5],-EXP220);}
  1469.   else y[4]=x[4];
  1470.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1471. }
  1472.  
  1473. GEN redrealform3(x,D,isqrtD,sens)
  1474.      GEN x,D,isqrtD;
  1475.      long sens;
  1476. {
  1477.   long fl,av=avma,tetpil;
  1478.   GEN y,p1;
  1479.   
  1480.   y=cgetg(4,17);y[1]=x[1];y[2]=x[2];y[3]=x[3];
  1481.   if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
  1482.   else
  1483.     {
  1484.       p1=subii(isqrtD,shifti(absi(x[1]),1));
  1485.       if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
  1486.     }
  1487.   while(fl)
  1488.     {
  1489.       y=rhorealform3(y,D,isqrtD);
  1490.       if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
  1491.       else
  1492.     {
  1493.       p1=subii(isqrtD,shifti(absi(y[1]),1));
  1494.       if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
  1495.     }
  1496.     }
  1497.   if(signe(y[1])<0)
  1498.     {
  1499.       if(sens||(!signe(addii(y[1],y[3]))))
  1500.     {
  1501.       tetpil=avma;return gerepile(av,tetpil,rhorealform3(y,D,isqrtD));
  1502.     }
  1503.       else
  1504.     {
  1505.       tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  1506.       setsigne(y[1],1);setsigne(y[3],-1);return y;
  1507.     }
  1508.     }
  1509.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1510. }
  1511.  
  1512. GEN redrealform5(x,pterr,D,isqrtD,sqrtD,sens)
  1513.      GEN x,D,isqrtD,sqrtD;
  1514.      long *pterr,sens;
  1515. {
  1516.   long fl,av=avma,tetpil;
  1517.   GEN y,p1;
  1518.   
  1519.   y=cgetg(6,17);y[1]=x[1];y[2]=x[2];y[3]=x[3];y[4]=x[4];y[5]=x[5];
  1520.   if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
  1521.   else
  1522.     {
  1523.       p1=subii(isqrtD,shifti(absi(x[1]),1));
  1524.       if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
  1525.     }
  1526.   while(fl)
  1527.     {
  1528.       y=rhorealform5(y,pterr,D,isqrtD,sqrtD);
  1529.       if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
  1530.       else
  1531.     {
  1532.       p1=subii(isqrtD,shifti(absi(y[1]),1));
  1533.       if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
  1534.     }
  1535.     }
  1536.   if(signe(y[1])<0)
  1537.     {
  1538.       if(sens||(!signe(addii(y[1],y[3]))))
  1539.     {
  1540.       tetpil=avma;
  1541.       return gerepile(av,tetpil,rhorealform5(y,pterr,D,isqrtD,sqrtD));
  1542.     }
  1543.       else
  1544.     {
  1545.       tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  1546.       setsigne(y[1],1);setsigne(y[3],-1);return y;
  1547.     }
  1548.     }
  1549.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1550. }
  1551.  
  1552. GEN redrealform(x,pterr,D,isqrtD,sqrtD,sens,precreg)
  1553.      GEN x,D,isqrtD,sqrtD;
  1554.      long *pterr,sens,precreg;
  1555.  
  1556. {
  1557.   long fl,av=avma,tetpil;
  1558.   GEN y,p1;
  1559.   
  1560.   y=cgetg(6,17);y[1]=x[1];y[2]=x[2];y[3]=x[3];y[4]=zero;affsr(1,y[5]=lgetr(precreg));
  1561.   if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
  1562.   else
  1563.     {
  1564.       p1=subii(isqrtD,shifti(absi(x[1]),1));
  1565.       if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
  1566.     }
  1567.   while(fl)
  1568.     {
  1569.       y=rhorealform5(y,pterr,D,isqrtD,sqrtD);
  1570.       if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
  1571.       else
  1572.     {
  1573.       p1=subii(isqrtD,shifti(absi(y[1]),1));
  1574.       if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
  1575.     }
  1576.     }
  1577.   if(signe(y[1])<0)
  1578.     {
  1579.       if(sens||(!signe(addii(y[1],y[3]))))
  1580.     {
  1581.       tetpil=avma;
  1582.       return gerepile(av,tetpil,rhorealform5(y,pterr,D,isqrtD,sqrtD));
  1583.     }
  1584.       else
  1585.     {
  1586.       tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  1587.       setsigne(y[1],1);setsigne(y[3],-1);return y;
  1588.     }
  1589.     }
  1590.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1591. }
  1592.  
  1593. GEN powrealform5(x,n,D,isqrtD,sqrtD,sens,precreg)
  1594.      GEN x,D,isqrtD,sqrtD;
  1595.      long n,sens,precreg;
  1596. {
  1597.   GEN y,p1;
  1598.   long av,tetpil,fl,toto=0;
  1599.  
  1600.   if(!n)
  1601.     {
  1602.       y=cgetg(6,17);y[1]=un;
  1603.       if(mpodd(x[2])) y[2]=(mpodd(isqrtD)) ? lcopy(isqrtD) : laddsi(-1,isqrtD);
  1604.       else {y[2]=lcopy(isqrtD);(((GEN)y[2])[lgef(isqrtD)-1])&=0xfffffffe;}
  1605.       av=avma;p1=subii(mulii(y[2],y[2]),D);tetpil=avma;
  1606.       y[3]=lpile(av,tetpil,shifti(p1,-2));y[4]=zero;affsr(1,y[5]=lgetr(precreg));
  1607.       return y;
  1608.     }
  1609.   if(n<0)
  1610.     {
  1611.       p1=cgetg(6,17);p1[1]=x[1];p1[2]=lnegi(x[2]);p1[3]=x[3];
  1612.       p1[4]=lnegi(addsi(1,x[4]));p1[5]=lshiftr(divir(un,x[5]),EXP220);
  1613.       n= -n;x=p1;
  1614.     }
  1615.   if(n==1) 
  1616.     {
  1617.       tetpil=avma;
  1618.       return gerepile(av,tetpil,redrealform(x,&toto,D,isqrtD,sqrtD,sens,precreg));
  1619.     }
  1620.   for(fl=0;n>1;n>>=1)
  1621.     {
  1622.       if(n&1)
  1623.     {
  1624.       if(fl) y=comprealform5(y,x,&toto,D,isqrtD,sqrtD,sens);
  1625.       else {fl=1;y=x;}
  1626.     }
  1627.       x=sqrealform5(x,&toto,D,isqrtD,sqrtD,sens);
  1628.     }
  1629.   tetpil=avma;y=fl ? comprealform5(y,x,&toto,D,isqrtD,sqrtD,sens) : gcopy(x);
  1630.   return gerepile(av,tetpil,y);
  1631. }
  1632.  
  1633. #define LIMP 30000
  1634.  
  1635. GEN lfunc(D)
  1636.      GEN D;
  1637. {
  1638.   GEN y;
  1639.   long av=avma,tetpil,prime=0;
  1640.   byteptr p=diffptr;
  1641.   
  1642.   prime=*p++;affsr(1,y=cgetr(4));
  1643.   do
  1644.     {
  1645.       if(!*p) err(recprimer);
  1646.       y=mulsr(prime,divrs(y,prime-krogs(D,prime)));
  1647.       prime+=*p++;
  1648.     }
  1649.   while(prime<=LIMP);
  1650.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1651. }
  1652.  
  1653.