home *** CD-ROM | disk | FTP | other *** search
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* ALGORITHMES SOUS-EXPONENTIELS DE CALCUL */
- /* DU GROUPE DE CLASSES ET DU REGULATEUR */
- /* (McCURLEY, BUCHMANN) */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- #include "genpari.h"
- #define C 15
- #define HASHT 1024
- #define EXP220 1048576
-
- void hphaseimag1(),hphaseimag2(),hphasereal1(),hphasereal2();
- long *subfactorbaseimag(),*subfactorbasereal(),factorbase();
- long factorise(),*largeprime(),*largeprime2();
- GEN **powsubfactimag(),**powsubfactreal(),rangi(),clean(),hnfimag();
- GEN sqrealform3(),comprealform3(),rhorealform3(),redrealform3();
- GEN sqrealform5(),powrealform5(),comprealform5(),rhorealform5(),redrealform5();
- GEN redrealform(),gcdreal(),hnfreal(),lfunc();
-
- GEN buchimag(D,gcbach,gCO)
- GEN D,gcbach,gCO;
- {
- long CO;
- double cbach;
- long limc,limc2,mglob,m,cp,nbram,auxrel,lo,lo1,ran,fl,fl2;
- long *fpd,b1,b2,fpc,pp,ep,b,extrarel,c,ic,jc,limhash;
- long kc2,nbcol1,nbcol2,nbcol3,nbrow1,nbrow2,nbrow3,*numbase,*base,*subbase;
- long **mat,*vectprime;
- long primfact[100],expoprimfact[100],primfact1[100],expoprimfact1[100];
- long badprim[100],nbbadprim;
- long av=avma,tetpil,kc,kcco,i,j,*pro,p2,*p1,*ex,q,s,nbtest,mm,av3;
- double drc,lim,logd;
- GEN dr,v,detmat,detmatsur2;
- GEN matgen,**vp,form,form1,pc,mit,met,mot,p3,p4;
- GEN matpro;
- long* hashtab[HASHT];
-
- if((typ(D)!=1)||(typ(gCO)!=1))
- err(talker,"incorrect type in buch");
- if(!signe(D)) err(talker,"zero discriminant in buch");
- if(signe(D)>0) err(talker,"positive discriminant in buchimag");
- s=D[lgef(D)-1]&3;
- if((s==1)||(s==2)) err(talker,"discriminant not 0 or 1 mod 4 in buch");
- if(signe(gCO)<=0) err(talker,"not enough relations in buch");
- CO=itos(gCO);cbach=gtodouble(gcbach);
- if((!cmpis(D,-3))||(!cmpis(D,-4)))
- {
- p3=cgetg(3,17);p3[1]=un;p3[2]=lgetg(1,17);return p3;
- }
- dr=cgetr(3);affir(D,dr);drc=rtodbl(dr);
- lim=sqrt(fabs(drc)/3.0);
- logd=log(fabs(drc));limc=cbach*logd*logd;cp=exp(sqrt(logd*log(logd)/8.0));
- if(cp>limc) limc=cp;limc2=6*logd*logd;if(limc>limc2) limc2=limc;
- subbase=subfactorbaseimag(D,&v,lim,&nbram);
- mglob=m=subbase[0];
- vp=powsubfactimag(v,m,C+7);
- kc2=factorbase(D,limc2,limc,&numbase,&base,&kc,badprim,&nbbadprim);
- pro=(long*)malloc(4*(kc2+1));vectprime=(long*)malloc(4*(kc2+1));
- for(i=1;i<=kc2;i++)
- {
- p2=vectprime[i]=base[i];
- pro[i]=(p2>0) ? p2 : -p2;
- }
- free(base);base=pro;
- kcco=kc+CO;
- mat=(long **)malloc(4*(kcco+1));
- for(i=1;i<=kcco;i++)
- {
- p1=(long *)malloc(4*(kc+1));mat[i]=p1;
- for(j=1;j<=kc;j++) p1[j]=0;
- }
- ex=(long*)malloc((m+1)*4);
- q=31-(long)ceil(log((double)C)/log(2.0));
- s=0;nbtest=0;
- mm=m+nbram+CO;
- limhash=(limc<32767)?limc*limc:(2<<30);
- for(i=0;i<HASHT;i++) hashtab[i]=(long*)0;
- auxrel=0;
- while(s<mm)
- {
- for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
- fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime(fpc,ex,0,0,mglob,hashtab);
- if(fpd)
- {
- auxrel++;
- s++;lo1=lo;
- for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
- form1=vp[1][fpd[2]];
- for(i=2;i<=m;i++)
- form1=compose(form1,vp[i][fpd[i+1]]);
- b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
- factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- if(b1==b2)
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]-=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- }
- else
- {
- if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
- else
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- }
- }
- }
- }
- if(fpc==1)
- {
- s++;
- for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b=itos(modis(form[2],(pp<<1)));
- if(b>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- }
- avma=av3;
- }
- while(s<kcco)
- {
- for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
- pc=primeform(D,stoi(base[s+1-CO]));
- form=compose(form,pc);
- fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime(fpc,ex,s+1-CO,0,mglob,hashtab);
- if(fpd)
- {
- auxrel++;
- s++;lo1=lo;
- for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
- form1=vp[1][fpd[2]];
- for(i=2;i<=m;i++)
- form1=compose(form1,vp[i][fpd[i+1]]);
- if(fpd[m+2])
- {
- pc=primeform(D,stoi(base[fpd[m+2]]));
- form1=compose(form1,pc);
- }
- b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
- factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- if(b1==b2)
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
- mat[s][s-CO]= -1;
- if(fpd[m+2]) mat[s][fpd[m+2]]++;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]-=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- }
- else
- {
- if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
- else
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
- mat[s][s-CO]= -1;
- if(fpd[m+2]) mat[s][fpd[m+2]]--;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- }
- }
- }
- }
- if(fpc==1)
- {
- s++;
- for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
- mat[s][s-CO]= -1;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b=itos(modis(form[2],(pp<<1)));
- if(b>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- if(!(mat[s][s-CO])) {for(i=1;i<=kc;i++) mat[s][i]=0;s--;}
- }
- avma=av3;
- }
- nbtest=auxrel=0;
- while(s<kc2)
- {
- for (i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++) form=compose(form,vp[i][ex[i]]);
- pp=base[s+1];
- pc=primeform(D,stoi(pp));
- form=compose(form,pc);
- fpc=factorise(form,s,pp,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime2(fpc,ex,s+1,mglob,hashtab);
- if(fpd&&(fpc!=pp)) {auxrel++;s++;}
- }
- if(fpc==1) s++;
- avma=av3;
- }
- hphaseimag1(mat,vectprime,kcco,kc,&nbcol1,&nbrow1);
- hphaseimag2(mat,vectprime,nbcol1,nbrow1,&nbcol2,&nbrow2);
- nbcol3=nbcol2;nbrow3=nbrow2;
- matgen=cgetg(nbcol3+1,19);
- for(i=1;i<=nbcol3;i++) matgen[i]=lgetg(nbrow3+1,18);
- for(i=1;i<=nbcol3;i++)
- {
- for(j=1;j<=nbrow3;j++) ((GEN)(matgen[i]))[j]=lstoi(mat[i][j]);
- }
- detmat=absi(rangi(matgen,&ran));
- detmatsur2=shifti(detmat,-1);
- if(gcmp1(absi(detmat))) {mit=idmat(1);nbrow3=1;nbcol3=lg(matpro);}
- else mit=hnfimag(matgen,detmat,detmatsur2);
- extrarel=nbcol3-nbrow3;
- p1=gun;
- for(i=1;i<=nbrow3;i++)
- {
- p3=(GEN)coeff(mit,i,i);
- if(signe(p3)) p1=mulii(p1,p3);
- }
- c=0;for(i=1;i<=nbrow3;i++) if(!gcmp1(coeff(mit,i,i))) c++;
- if(!c) met=idmat(1);
- else
- {
- met=cgetg(c+1,19);ic=0;
- for(i=1;i<=nbrow3;i++)
- {
- if(!gcmp1(coeff(mit,i,i)))
- {
- ic++;p3=cgetg(c+1,18);met[ic]=(long)p3;jc=0;
- for(j=1;j<=nbrow3;j++)
- {
- if(!gcmp1(coeff(mit,j,j)))
- {
- jc++;p3[jc]=coeff(mit,j,i);
- }
- }
- }
- }
- }
- met=smith(met);
- lo=lg(met)-1;
- c=0;for(i=1;i<=lo;i++) if(!gcmp1(met[i])) c++;
- p3=gdiv(gmul(p1,mppi(4)),gmul(lfunc(D),gsqrt(absi(D),4)));
- tetpil=avma;p4=cgetg(4,17);p4[1]=lcopy(p1);mot=cgetg(c+1,17);p4[2]=(long)mot;
- for(i=1;i<=c;i++) mot[i]=lcopy(met[i]);p4[3]=lcopy(p3);
- free(base);free(vectprime);free(ex);free(subbase);free(numbase);
- for(i=1;i<=kcco;i++) free(mat[i]);free(mat);
- for(i=1;i<=m;i++) free(vp[i]);free(vp);
- return gerepile(av,tetpil,p4);
- }
-
- long *subfactorbaseimag(d,w,ll,ptnbram)
- long *ptnbram;
- double ll;
- GEN d,*w;
- {
- byteptr delta=diffptr;
- long av1,i,j,pp,tp,fl,*subbase,pro[100],toto=0;
- double prod;
- GEN p1;
-
- av1=avma;
- i=1;pp=0;*ptnbram=0;fl=1;prod=1;
- while(fl||(i==1))
- {
- pp+=*delta++;tp=krogs(d,pp);
- if(tp==1) {pro[i]=pp;i++;prod*=pp;} else if(!tp) (*ptnbram)++;
- if((prod>ll)&&(i>1)) fl=0;
- }
- avma=av1;i--;
- *w=cgetg(i+1,18);
- for(j=1;j<=i;j++)
- {
- p1=primeform(d,stoi(pro[j]));(*w)[j]=(long)p1;
- }
- subbase=(long*)malloc(4*(i+1));subbase[0]=i;
- for(j=1;j<=i;j++) subbase[j]=pro[j];
- return subbase;
- }
-
- GEN **powsubfactimag(w,n,a)
- GEN w;
- long n,a;
- {
- long i,j,toto=0;
- GEN **x;
-
- x=(GEN**)malloc(4*(n+1));
- for(i=1;i<=n;i++) x[i]=(GEN*)malloc(4*(a+1));
- for(i=1;i<=n;i++)
- {
- x[i][0]=gpuigs(w[1],0);
- for(j=1;j<=a;j++) x[i][j]=compose(x[i][j-1],w[i]);
- }
- return x;
- }
-
- long factorbase(d,n2,n,ptnum,ptbase,ptkc,badprim,nbbadprim)
- GEN d;
- long n2,n,**ptnum,**ptbase,*ptkc,*badprim,*nbbadprim;
- {
- byteptr delta=diffptr;
- long av2,i,pp,qq,fl,kr,r,*numbase,*base,kc;
- GEN p1;
-
- numbase=(long*)malloc(4*(n2+1));*ptnum=numbase;
- base=(long*)malloc(4*(n2+1));*ptbase=base;
- *nbbadprim=0;
- av2=avma;
- i=0;pp=*delta++;qq=2;fl=1;
- while(pp<=n2)
- {
- if((kr=krogs(d,pp))!=-1)
- {
- if(kr) {i++;numbase[pp]=i;base[i]=pp;}
- else
- {
- p1=divis(d,pp);
- if(signe(modis(p1,pp))) {i++;numbase[pp]=i;base[i]= -pp;}
- else
- {
- if(pp==2)
- {
- r=p1[lgef(p1)-1]&7;if(signe(d)<0) r=8-r;
- if(r>=4) {i++;numbase[pp]=i;base[i]= -pp;}
- else badprim[++(*nbbadprim)]=pp;
- }
- else badprim[++(*nbbadprim)]=pp;
- }
- }
- }
- pp+=*delta++;
- if((pp>n)&&fl) {kc=i;*ptkc=kc;fl=0;}
- }
- avma=av2;
- return i;
- }
-
- long factorise(f,n,limp,ptlo,primfact,expoprimfact,badprim,nbbadprim,base,limhash)
- GEN f;
- long n,limp,*ptlo,*primfact,*expoprimfact,*badprim,nbbadprim,*base,limhash;
-
- {
- long sr,i,p,k,fl=1,av=avma,q1,lo;
- GEN x,q,r;
-
- lo=0;
- x=(GEN)(f[1]);if(gcmp1(x)) {*ptlo=0;return 1;}
- for(i=1;(i<=n)&&fl;i++)
- {
- p=base[i];q=dvmdis(x,p,&r);
- if(sr=(!signe(r)))
- {
- primfact[++lo]=p;x=q;k=0;av=avma;
- while(sr)
- {
- k++;q=dvmdis(x,p,&r);
- if(sr=(!signe(r))) {x=q;av=avma;}
- }
- expoprimfact[lo]=k;
- }
- else avma=av;
- fl=(cmpis(q,p)>0);
- }
- if(!fl)
- {
- if(gcmp1(x)) {avma=av;*ptlo=lo;return 1;}
- else
- {
- if(cmpis(x,limp)<=0)
- {
- for(i=1;i<=nbbadprim;i++)
- if(!signe(modis(x,badprim[i]))) {avma=av;*ptlo=lo;return 0;}
- primfact[++lo]=itos(x);expoprimfact[lo]=1;
- avma=av;*ptlo=lo;return 1;
- }
- }
- }
- *ptlo=lo;
- if(cmpis(x,limhash)<=0)
- {
- q1=itos(x);
- avma=av;return q1;
- }
- else {avma=av;return 0;}
- avma=av;return 0;
- }
-
- long *largeprime2(q1,ex,np,mglob,hashtab)
- long q1,*ex,np,mglob,*hashtab[];
- {
- long hashv,*pt,cpt,i,*p1;
- hashv=((q1&2047)-1)>>1;
- pt=hashtab[hashv];cpt=0;
- while(pt&&(q1!=pt[1])) {pt=(long*)(*pt);cpt++;}
- if(!pt)
- {
- if(!(p1=(long*)malloc((mglob+4)<<2))) err(talker,"debordement de pile de hash");
- p1[1]=q1;for(i=2;i<mglob+2;i++) p1[i]=ex[i-1];p1[mglob+2]=np;p1[mglob+3]=0;
- p1[0]=cpt ? (long)hashtab[hashv] : 0;
- hashtab[hashv]=p1;return (long*)0;
- }
- else return (pt[mglob+2]==np) ? (long*)0 : pt;
- }
-
- int compte(mat,row,longueur,firstnonzero)
- long **mat,row,longueur,*firstnonzero;
- {
- int n,j,p;
- n=0;
- for (j=1;j<=longueur;j++)
- {
- p=mat[j][row];
- if(p) {if(abs(p)>=2) return 32767;else {n++;*firstnonzero=j;}}
- }
- return n;
- }
-
- void hphaseimag1(mat,vectprime,kcco,kc,ptcol,ptrow)
- long **mat,*vectprime,kcco,kc,*ptcol,*ptrow;
-
- {
- long numrow,i,i1,j,nb,n,extrarel;
- extrarel=kcco-kc;
- numrow=kc;
- while(numrow>=1)
- {
- for(i=numrow;(i>=1)&&((nb=compte(mat,i,extrarel+numrow,&n))>1);i--);
- if(!nb)
- {
- for(j=1;j<=extrarel+numrow;j++)
- for(i1=i;i1<numrow;i1++) mat[j][i1]=mat[j][i1+1];
- for(i1=i;i1<numrow;i1++) vectprime[i1]=vectprime[i1+1];
- numrow--;extrarel++;
- }
- else if(nb==1)
- {
- for(j=1;j<=extrarel+numrow;j++) mat[j][i]=mat[j][numrow];
- vectprime[i]=vectprime[numrow];
- mat[n]=mat[extrarel+numrow];
- numrow--;
- }
- else {*ptrow=numrow;*ptcol=numrow+extrarel;return;}
- }
- err(talker,"not enough relations or h=1 in hphaseimag1");
- }
-
- void hphaseimag2(mat,vectprime,nbcol1,nbrow1,ptcol,ptrow)
- long **mat,*vectprime,nbcol1,nbrow1,*ptcol,*ptrow;
-
- {
- long numrow,i,j,nb,n,sizemax,imin,nbmin,nmin,s,t,fl=1,extrarel;
-
- extrarel=nbcol1-nbrow1;
- numrow=nbrow1;sizemax=0;
- while(fl&&(numrow>=1)&&(sizemax<=1073741823))
- {
- imin=numrow;nbmin=compte(mat,numrow,extrarel+numrow,&n);nmin=n;
- for(i=numrow-1;i>=1;i--)
- {
- nb=compte(mat,i,extrarel+numrow,&n);
- if(nb<nbmin) {nbmin=nb;imin=i;nmin=n;}
- }
- if(nbmin<32767)
- {
- s=mat[nmin][imin];
- for(j=1;j<nmin;j++)
- {
- sizemax=0;
- if(t=mat[j][imin])
- {
- if(s==t)
- {
- for(i=1;i<=numrow;i++)
- {
- mat[j][i]-=mat[nmin][i];
- sizemax=max(sizemax,abs(mat[j][i]));
- }
- }
- else
- {
- for(i=1;i<=numrow;i++)
- {
- mat[j][i]+=mat[nmin][i];
- sizemax=max(sizemax,abs(mat[j][i]));
- }
- }
- }
- }
- for(j=1;j<=extrarel+numrow;j++) mat[j][imin]=mat[j][numrow];
- vectprime[imin]=vectprime[numrow];
- mat[nmin]=mat[extrarel+numrow];
- numrow--;
- }
- else {*ptrow=numrow;*ptcol=numrow+extrarel;fl=0;}
- }
- if(!numrow) err(talker,"not enough relations or h=1 in hphaseimag2");
- if(sizemax>1073741823) return;
- return;
- }
-
- GEN hnfimag(x,detmat,detmatsur2)
- GEN x,detmat,detmatsur2;
- {
- long li,co,av,tetpil,i,j,def,ldef,lim,dec;
- GEN b,q,w,p1,y,d,u,v,p3,p4;
-
- lim=(avma+bot)>>1;
- av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
- def=co;ldef=(li>co)?li-co+1:1;
- for(i=li-1;i>=ldef;i--)
- {
- def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
- while(j>1)
- {
- d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
- p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
- 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]));
- y[j]=(long)clean(y[j],detmat,detmatsur2);
- y[j-1]=(long)p1;
- y[j-1]=(long)clean(y[j-1],detmat,detmatsur2);
- j--;while(j&&(!signe(coeff(y,i,j)))) j--;
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- if(j==1)
- {
- d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
- p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
- y[1]=lsub(gmul(p3=divii(coeff(y,i,1),d),y[def]),gmul(p4=divii(coeff(y,i,def),d),y[1]));
- y[1]=(long)clean(y[1],detmat,detmatsur2);
- y[def]=(long)p1;
- y[def]=(long)clean(y[def],detmat,detmatsur2);
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
- b=gcopy(detmat);w=cgetg(li,19);def--;
- for(i=li-1;i>=1;i--)
- {
- d=bezout(coeff(y,i,i+def),b,&u,&v);w[i]=lmod(gmul(u,y[i+def]),b);
- if(!signe(coeff(w,i,i))) coeff(w,i,i)=(long)d;
- if(i>1) b=divii(b,d);
- }
- for(i=li-2;i>=1;i--)
- {
- for(j=i+1;j<li;j++)
- {
- q=gdivent(coeff(w,i,j),coeff(w,i,i));w[j]=lsub(w[j],gmul(q,w[i]));
- }
- if(avma<lim) {tetpil=avma;w=gerepile(av,tetpil,gcopy(w));}
- }
- tetpil=avma;w=gerepile(av,tetpil,gcopy(w));
- return w;
- }
-
- GEN rangi(x,rr)
- GEN x;
- long *rr;
- {
- GEN pass,c,v,det1;
- long i,j,k,rg,t,n,n1,m,m1,nbmot,av,av1,pivprec,piv,p1,p,vi,tetpil,cm=0;
-
- nbmot=20;n=(n1=lg(x))-1;m=(m1=lg(x[1]))-1;av=avma;
- pivprec=lgeti(nbmot);affsi(0,det1=cgeti(nbmot));affsi(1,piv=lgeti(nbmot));pass=cgetg(m1,19);
- for(j=1;j<=m;j++)
- {
- p=pass[j]=lgetg(m1,18);for(i=1;i<=m;i++) ((GEN)p)[i]=lgeti(nbmot);
- }
- c=cgeti(m1);for(k=1;k<=m;k++) c[k]=0;
- v=cgetg(m1,18);av1=avma;k=1;rg=0;
- while((k<=n)&&(rg<m))
- {
- for(t=0,i=1;i<=m;i++)
- if(!c[i])
- {
- vi=lmulii(piv,coeff(x,i,k));
- for(j=1;j<=m;j++)
- if(c[j]) vi=laddii(vi,mulii(coeff(pass,i,j),coeff(x,j,k)));
- v[i]=vi;
- if(!t) if(signe(vi)) t=i;
- }
- if (t)
- {
- rg++;c[t]=k;
- affii(piv,pivprec);affii(v[t],piv);
- if(rg<m)
- {
- for(i=1;i<=m;i++)
- if(!c[i])
- for(j=1;j<=m;j++)
- if(c[j]&&(j!=t))
- {
- p1=lsubii(mulii(piv,coeff(pass,i,j)),mulii(v[i],coeff(pass,t,j)));
- if(rg>1) diviiz(p1,pivprec,coeff(pass,i,j));
- else affii(p1,coeff(pass,i,j));
- }
- for(i=1;i<=m;i++) if(!c[i]) mpnegz(v[i],coeff(pass,i,t));
- }
- else c[t]=0;
- }
- if(rg==m)
- {
- cm=1;affii(ggcd(piv,det1),det1);
- affii(pivprec,piv);rg--;
- }
- avma=av1;k++;
- }
- *rr=rg+cm;tetpil=avma;return gerepile(av,tetpil,gcopy(det1));
- }
-
- GEN clean(x,detmat,detmatsur2)
- GEN x,detmat,detmatsur2;
- {
- long lx=lg(x),i;
- GEN y,p1;
- y=cgetg(lx,18);
- for(i=1;i<lx;i++)
- {
- p1=modii(x[i],detmat);
- y[i]=(cmpii(p1,detmatsur2)>0) ? lsubii(p1,detmat) : (long)p1;
- }
- return y;
- }
-
- GEN buchreal(D,gsens,gcbach,gCO)
- GEN D,gsens,gcbach,gCO;
- {
- long sens,CO;
- double cbach;
- long precreg,limc,limc2,mglob,m,cp,nbram,nbrederr,auxrel,lo,lo1,ran,fl,fl2;
- long nrho,*fpd,b1,b2,fpc,pp,ep,b,extrarel,c,ic,jc,dec,limhash;
- long kc2,nbcol1,nbcol2,nbcol3,nbrow1,nbrow2,nbrow3,*numbase,*base,*subbase;
- long **mat,*vectprime;
- long primfact[100],expoprimfact[100],primfact1[100],expoprimfact1[100];
- long badprim[100],nbbadprim;
- long av=avma,tetpil,kc,kcco,i,j,*pro,p2,*p1,*ex,q,s,nbtest,mm,av3;
- double drc,lim,logd;
- GEN dr,sqrtD,isqrtD,errglob,log2precis,v,detmat,detmatsur2;
- GEN matgen,vecexpo,vecerr,ermax,**vp,form,form1,forminit,pc,mit,met,mot;
- GEN wecexpo,wecerr,xecexpo,xecerr,matpro,p3,p4;
- long* hashtab[HASHT];
-
- if((typ(D)!=1)||(typ(gCO)!=1)||(typ(gsens)!=1))
- err(talker,"incorrect type in buch");
- if(!signe(D)) err(talker,"zero discriminant in buch");
- if(signe(D)<0) err(talker,"negative discriminant in buchreal, try buchimag");
- s=D[lgef(D)-1]&3;
- if((s==2)||(s==3)) err(talker,"discriminant not 0 or 1 mod 4 in buch");
- if(signe(gCO)<=0) err(talker,"not enough relations in buch");
- CO=itos(gCO);cbach=gtodouble(gcbach);sens=abs(signe(gsens));
- dr=cgetr(3);affir(D,dr);drc=rtodbl(dr);
- lim=sqrt(fabs(drc)/3.0);
- precreg=(gexpo(D)>>5)+5;
- sqrtD=gsqrt(D,precreg);isqrtD=gfloor(sqrtD);
- affsr(1,errglob=cgetr(4));setexpo(errglob,-((precreg-2)<<5));
- log2precis=glog(gdeux,precreg);
- logd=log(fabs(drc));limc=cbach*logd*logd;cp=exp(sqrt(logd*log(logd)/8.0));
- if(cp>limc) limc=cp;limc2=6*logd*logd;if(limc>limc2) limc2=limc;
- subbase=subfactorbasereal(D,&v,lim,precreg,&nbram,isqrtD,sqrtD,sens);
- mglob=m=subbase[0];
- vp=powsubfactreal(v,m,C+7,D,isqrtD,sqrtD,sens,precreg);
- kc2=factorbase(D,limc2,limc,&numbase,&base,&kc,badprim,&nbbadprim);
- pro=(long*)malloc(4*(kc2+1));vectprime=(long*)malloc(4*(kc2+1));
- for(i=1;i<=kc2;i++)
- {
- p2=vectprime[i]=base[i];
- pro[i]=(p2>0) ? p2 : -p2;
- }
- free(base);base=pro;
- kcco=kc+CO;
- mat=(long **)malloc(4*(kcco+1));
- vecexpo=cgetg(kcco+1,17);vecerr=cgetg(kcco+1,17);
- for(i=1;i<=kcco;i++) {vecexpo[i]=lgetr(precreg);vecerr[i]=lgetr(4);}
- for(i=1;i<=kcco;i++)
- {
- p1=(long *)malloc(4*(kc+1));mat[i]=p1;
- for(j=1;j<=kc;j++) p1[j]=0;
- }
- ex=(long*)malloc((m+1)*4);
- q=31-(long)ceil(log((double)C)/log(2.0));
- s=0;nbtest=0;
- mm=m+nbram+CO;
- limhash=(limc<32767)?limc*limc:(2<<30);
- for(i=0;i<HASHT;i++) hashtab[i]=(long*)0;
- auxrel=0;
- while(s<mm)
- {
- for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;nbrederr=0;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++)
- form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
- fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime(fpc,ex,0,0,mglob,hashtab);
- if(fpd)
- {
- auxrel++;
- s++;lo1=lo;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++)
- form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
- for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
- form1=vp[1][fpd[2]];
- for(i=2;i<=m;i++)
- form1=comprealform5(form1,vp[i][fpd[i+1]],&nbrederr,D,isqrtD,sqrtD,sens);
- b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
- factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- if(b1==b2)
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]-=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,subii(form[4],form1[4])),log2precis),
- mplog(absr(divrr(form[5],form1[5])))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- }
- else
- {
- if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
- else
- {
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,addii(form1[4],form[4])),log2precis),mplog(absr(mulrr(form1[5],form[5])))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- }
- }
- }
- }
- if(fpc==1)
- {
- s++;
- form=vp[1][ex[1]];
- for(i=2;i<=m;i++)
- form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
- for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b=itos(modis(form[2],(pp<<1)));
- if(b>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,form[4]),log2precis),mplog(absr(form[5]))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- }
- avma=av3;
- }
- while(s<kcco)
- {
- for(i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++) form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
- pc=redrealform3(primeform(D,stoi(base[s+1-CO]),precreg),D,isqrtD,sens);
- form=comprealform3(form,pc,D,isqrtD,sens);forminit=form;nrho=0;
- fl2=0;
- do
- {
- fpc=factorise(form,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime(fpc,ex,s+1-CO,nrho,mglob,hashtab);
- if(fpd)
- {
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++)
- form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
- pc=redrealform(primeform(D,stoi(base[s+1-CO]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
- pc[4]=zero;affsr(1,pc[5]);
- form=comprealform5(form,pc,&nbrederr,D,isqrtD,sqrtD,sens);
- for(i=1;i<=nrho;i++)
- form=rhorealform5(form,&nbrederr,D,isqrtD,sqrtD);
- auxrel++;
- s++;lo1=lo;
- for(j=1;j<=lo1;j++) {primfact1[j]=primfact[j];expoprimfact1[j]=expoprimfact[j];}
- form1=vp[1][fpd[2]];
- for(i=2;i<=m;i++)
- form1=comprealform5(form1,vp[i][fpd[i+1]],&nbrederr,D,isqrtD,sqrtD,sens);
- if(fpd[m+2])
- {
- pc=redrealform(primeform(D,stoi(base[fpd[m+2]]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
- pc[4]=zero;affsr(1,pc[5]);
- form1=comprealform5(form1,pc,&nbrederr,D,isqrtD,sqrtD,sens);
- for(i=1;i<=fpd[m+3];i++)
- form1=rhorealform5(form1,&nbrederr,D,isqrtD,sqrtD);
- if((!sens)&&(signe(addii(form1[1],form1[3]))))
- {setsigne(form1[1],1);setsigne(form1[3],-1);}
- }
- b1=itos(modis(form1[2],(fpc<<1)));b2=itos(modis(form[2],(fpc<<1)));
- factorise(form1,kc,limc,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- if(b1==b2)
- {
- fl2=1;
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]=fpd[i+1]-ex[i];
- mat[s][s-CO]= -1;
- if(fpd[m+2]) mat[s][fpd[m+2]]++;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]-=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,subii(form[4],form1[4])),log2precis),mplog(absr(divrr(form[5],form1[5])))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- }
- else
- {
- if((b1+b2)!=(fpc<<1)) {s--;auxrel--;}
- else
- {
- fl2=1;
- for(i=1;i<=m;i++)
- mat[s][numbase[subbase[i]]]= -fpd[i+1]-ex[i];
- mat[s][s-CO]= -1;
- if(fpd[m+2]) mat[s][fpd[m+2]]--;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b1=itos(modis(form1[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- for(j=1;j<=lo1;j++)
- {
- pp=primfact1[j];ep=expoprimfact1[j];
- b1=itos(modis(form[2],(pp<<1)));
- if(b1>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,addii(form1[4],form[4])),log2precis),mplog(absr(mulrr(form1[5],form[5])))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- }
- }
- }
- }
- if ((fpc!=1)&&(!fl2))
- {
- form=rhorealform3(form,D,isqrtD);nrho++;
- if((sens)||(!signe(addii(form[1],form[3]))))
- {form=rhorealform3(form,D,isqrtD);nrho++;}
- else
- {
- setsigne(form[1],1);setsigne(form[3],-1);
- }
- fl=gegal(form[1],forminit[1]);
- if(fl) fl=gegal(form[2],forminit[2]);
- if(fl) fl=gegal(form[3],forminit[3]);
- }
- }
- while((fpc!=1)&&(!fl)&&(!fl2));
- if(fpc==1)
- {
- nbrederr=0;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++)
- form=comprealform5(form,vp[i][ex[i]],&nbrederr,D,isqrtD,sqrtD,sens);
- pc=redrealform(primeform(D,stoi(base[s+1-CO]),precreg),&nbrederr,D,isqrtD,sqrtD,sens,precreg);
- pc[4]=zero;affsr(1,pc[5]);
- form=comprealform5(form,pc,&nbrederr,D,isqrtD,sqrtD,sens);
- for(i=1;i<=nrho;i++)
- form=rhorealform5(form,&nbrederr,D,isqrtD,sqrtD);
- s++;
- for(i=1;i<=m;i++) mat[s][numbase[subbase[i]]]= -ex[i];
- mat[s][s-CO]= -1;
- for(j=1;j<=lo;j++)
- {
- pp=primfact[j];ep=expoprimfact[j];
- b=itos(modis(form[2],(pp<<1)));
- if(b>pp) ep= -ep;
- mat[s][numbase[pp]]+=ep;
- }
- affrr(shiftr(mpadd(mulir(mulsi(EXP220,form[4]),log2precis),mplog(absr(form[5]))),-1),vecexpo[s]);
- mulsrz(nbrederr,errglob,vecerr[s]);
- if(!(mat[s][s-CO])) {for(i=1;i<=kc;i++) mat[s][i]=0;s--;}
- }
- avma=av3;
- }
- nbtest=auxrel=0;
- while(s<kc2)
- {
- for (i=1;i<=m;i++) ex[i]=(rand()>>q)+7;
- av3=avma;
- form=vp[1][ex[1]];
- for (i=2;i<=m;i++) form=comprealform3(form,vp[i][ex[i]],D,isqrtD,sens);
- pp=base[s+1];
- pc=redrealform3(primeform(D,stoi(pp),precreg),D,isqrtD,sens);
- form=comprealform3(form,pc,D,isqrtD,sens);forminit=form;fl2=0;
- do
- {
- fpc=factorise(form,s,pp,&lo,primfact,expoprimfact,badprim,nbbadprim,base,limhash);
- nbtest++;
- if(fpc>1)
- {
- fpd=largeprime2(fpc,ex,s+1,mglob,hashtab);
- if(fpd&&(fpc!=pp)) {auxrel++;s++;fl2=1;}
- }
- if((fpc!=1)&&(!fl2))
- {
- form=rhorealform3(form,D,isqrtD);
- if((sens)||(!signe(addii(form[1],form[3]))))
- form=rhorealform3(form,D,isqrtD);
- else
- {
- setsigne(form[1],1);setsigne(form[3],-1);
- }
- fl=gegal(form[1],forminit[1]);
- if(fl) fl=gegal(form[2],forminit[2]);
- if(fl) fl=gegal(form[3],forminit[3]);
- }
- }
- while((fpc!=1)&&(!fl)&&(!fl2));
- if(fpc==1) s++;
- avma=av3;
- }
- hphasereal1(mat,vectprime,vecexpo,vecerr,kcco,kc,&nbcol1,&nbrow1);
- hphasereal2(mat,vectprime,vecexpo,vecerr,nbcol1,nbrow1,&nbcol2,&nbrow2);
- nbcol3=nbcol2;nbrow3=nbrow2;
- matgen=cgetg(nbcol3+1,19);
- for(i=1;i<=nbcol3;i++) matgen[i]=lgetg(nbrow3+1,18);
- wecexpo=cgetg(nbcol3+1,17);wecerr=cgetg(nbcol3+1,17);
- for(i=1;i<=nbcol3;i++)
- {
- wecexpo[i]=vecexpo[i];wecerr[i]=vecerr[i];
- for(j=1;j<=nbrow3;j++) ((GEN)(matgen[i]))[j]=lstoi(mat[i][j]);
- }
- vecexpo=wecexpo;vecerr=wecerr;
- detmat=absi(rangi(matgen,&ran));
- detmatsur2=shifti(detmat,-1);
- matpro=kerint(matgen);
- xecexpo=gmul(vecexpo,matpro);
- xecerr=gmul(vecerr,gabs(matpro));
- if(gcmp1(absi(detmat))) {mit=idmat(1);nbrow3=1;nbcol3=lg(matpro);}
- else mit=hnfreal(matgen,detmat,detmatsur2,vecexpo,vecerr,&wecexpo,&wecerr);
- wecexpo=xecexpo;wecerr=xecerr;
- extrarel=nbcol3-nbrow3;
- for(i=2;i<=extrarel;i++)
- wecexpo[i]=(long)gcdreal(wecexpo[i-1],wecexpo[i],wecerr[i-1],wecerr[i],wecexpo+i-1);
- ermax=(GEN)wecexpo[1];for(i=2;i<extrarel;i++) ermax=gmax(ermax,wecexpo[i]);
- p1=gun;
- for(i=1;i<=nbrow3;i++)
- {
- p3=(GEN)coeff(mit,i,i);
- if(signe(p3)) p1=mulii(p1,p3);
- }
- c=0;for(i=1;i<=nbrow3;i++) if(!gcmp1(coeff(mit,i,i))) c++;
- if(!c) met=idmat(1);
- else
- {
- met=cgetg(c+1,19);ic=0;
- for(i=1;i<=nbrow3;i++)
- {
- if(!gcmp1(coeff(mit,i,i)))
- {
- ic++;p3=cgetg(c+1,18);met[ic]=(long)p3;jc=0;
- for(j=1;j<=nbrow3;j++)
- {
- if(!gcmp1(coeff(mit,j,j)))
- {
- jc++;p3[jc]=coeff(mit,j,i);
- }
- }
- }
- }
- }
- met=smith(met);
- lo=lg(met)-1;
- c=0;for(i=1;i<=lo;i++) if(!gcmp1(met[i])) c++;
- p2=(long)absr(wecexpo[extrarel]);
- p3=gdiv(gmul2n(gmul(p1,p2),1),gmul(lfunc(D),sqrtD));
- tetpil=avma;p4=cgetg(6,17);p4[1]=lcopy(p1);mot=cgetg(c+1,17);p4[2]=(long)mot;
- for(i=1;i<=c;i++) mot[i]=lcopy(met[i]);
- p4[3]=lcopy(p2);p4[4]=lcopy(ermax);p4[5]=lcopy(p3);
- free(base);free(vectprime);free(ex);free(subbase);free(numbase);
- for(i=1;i<=kcco;i++) free(mat[i]);free(mat);
- for(i=1;i<=m;i++) free(vp[i]);free(vp);
- return gerepile(av,tetpil,p4);
- }
-
- long *subfactorbasereal(d,w,ll,precreg,ptnbram,isqrtD,sqrtD,sens)
- long precreg,*ptnbram,sens;
- double ll;
- GEN d,*w;
- GEN isqrtD,sqrtD;
- {
- byteptr delta=diffptr;
- long av1,i,j,pp,tp,fl,*subbase,pro[100],toto=0;
- double prod;
- GEN p1;
-
- av1=avma;
- i=1;pp=0;*ptnbram=0;fl=1;prod=1;
- while(fl||(i==1))
- {
- pp+=*delta++;tp=krogs(d,pp);
- if(tp==1) {pro[i]=pp;i++;prod*=pp;} else if(!tp) (*ptnbram)++;
- if((prod>ll)&&(i>1)) fl=0;
- }
- avma=av1;i--;
- *w=cgetg(i+1,18);
- for(j=1;j<=i;j++)
- {
- p1=redrealform(primeform(d,stoi(pro[j]),precreg),&toto,d,isqrtD,sqrtD,sens,precreg);
- p1[4]=zero;affsr(1,p1[5]);
- (*w)[j]=(long)p1;
- }
- subbase=(long*)malloc(4*(i+1));subbase[0]=i;
- for(j=1;j<=i;j++) subbase[j]=pro[j];
- return subbase;
- }
-
- GEN **powsubfactreal(w,n,a,D,isqrtD,sqrtD,sens,precreg)
- GEN w,D,isqrtD,sqrtD;
- long n,a,sens,precreg;
- {
- long i,j,toto=0;
- GEN **x;
-
- x=(GEN**)malloc(4*(n+1));
- for(i=1;i<=n;i++) x[i]=(GEN*)malloc(4*(a+1));
- for(i=1;i<=n;i++)
- {
- j=0;x[i][j]=powrealform5(w[1],0,D,isqrtD,sqrtD,sens,precreg);
- while(j<=a)
- {
- j++;
- x[i][j]=comprealform5(x[i][j-1],w[i],&toto,D,isqrtD,sqrtD,sens);
- }
- }
- fflush(stdout);
- return x;
- }
-
- long *largeprime(q1,ex,np,nrho,mglob,hashtab)
- long q1,*ex,np,nrho,mglob,*hashtab[];
- {
- long hashv,*pt,cpt,i,*p1,fl;
- hashv=((q1&2047)-1)>>1;
- pt=hashtab[hashv];cpt=0;
- while(pt&&(q1!=pt[1])) {pt=(long*)(*pt);cpt++;}
- if(!pt)
- {
- if(!(p1=(long*)malloc((mglob+4)<<2))) err(talker,"debordement de pile de hash");
- p1[1]=q1;for(i=2;i<mglob+2;i++) p1[i]=ex[i-1];p1[mglob+2]=np;p1[mglob+3]=nrho;
- p1[0]=cpt ? (long)hashtab[hashv] : 0;
- hashtab[hashv]=p1;return (long*)0;
- }
- else
- {
- fl=1;i=2;while(fl&&(i<(mglob+2))) {fl=(pt[i]==ex[i-1]);i++;}
- if(fl) fl=(pt[i]==np);
- return fl ? (long*)0 : pt;
- }
- }
-
- void hphasereal1(mat,vectprime,vecexpo,vecerr,kcco,kc,ptcol,ptrow)
- long **mat,*vectprime,kcco,kc,*ptcol,*ptrow;
- GEN vecexpo,vecerr;
-
- {
- long numrow,i,i1,j,nb,n,extrarel;
- extrarel=kcco-kc;
- numrow=kc;
- while(numrow>=1)
- {
- for(i=numrow;(i>=1)&&((nb=compte(mat,i,extrarel+numrow,&n))>1);i--);
- if(!nb)
- {
- for(j=1;j<=extrarel+numrow;j++)
- for(i1=i;i1<numrow;i1++) mat[j][i1]=mat[j][i1+1];
- for(i1=i;i1<numrow;i1++) vectprime[i1]=vectprime[i1+1];
- numrow--;extrarel++;
- }
- else if(nb==1)
- {
- for(j=1;j<=extrarel+numrow;j++) mat[j][i]=mat[j][numrow];
- vectprime[i]=vectprime[numrow];
- mat[n]=mat[extrarel+numrow];
- vecexpo[n]=vecexpo[extrarel+numrow];
- vecerr[n]=vecerr[extrarel+numrow];
- numrow--;
- }
- else {*ptrow=numrow;*ptcol=numrow+extrarel;return;}
- }
- err(talker,"not enough relations or h=1 in hphasereal1");
- }
-
- void hphasereal2(mat,vectprime,vecexpo,vecerr,nbcol1,nbrow1,ptcol,ptrow)
- long **mat,*vectprime,nbcol1,nbrow1,*ptcol,*ptrow;
- GEN vecexpo,vecerr;
-
- {
- long numrow,i,j,nb,n,sizemax,imin,nbmin,nmin,s,t,fl=1,extrarel;
-
- extrarel=nbcol1-nbrow1;
- numrow=nbrow1;sizemax=0;
- while(fl&&(numrow>=1)&&(sizemax<=1073741823))
- {
- imin=numrow;nbmin=compte(mat,numrow,extrarel+numrow,&n);nmin=n;
- for(i=numrow-1;i>=1;i--)
- {
- nb=compte(mat,i,extrarel+numrow,&n);
- if(nb<nbmin) {nbmin=nb;imin=i;nmin=n;}
- }
- if(nbmin<32767)
- {
- s=mat[nmin][imin];
- for(j=1;j<nmin;j++)
- {
- sizemax=0;
- if(t=mat[j][imin])
- {
- if(s==t)
- {
- for(i=1;i<=numrow;i++)
- {
- mat[j][i]-=mat[nmin][i];
- sizemax=max(sizemax,abs(mat[j][i]));
- }
- subrrz(vecexpo[j],vecexpo[nmin],vecexpo[j]);
- addrrz(vecerr[j],vecerr[nmin],vecerr[j]);
- }
- else
- {
- for(i=1;i<=numrow;i++)
- {
- mat[j][i]+=mat[nmin][i];
- sizemax=max(sizemax,abs(mat[j][i]));
- }
- addrrz(vecexpo[j],vecexpo[nmin],vecexpo[j]);
- addrrz(vecerr[j],vecerr[nmin],vecerr[j]);
- }
- }
- }
- for(j=1;j<=extrarel+numrow;j++) mat[j][imin]=mat[j][numrow];
- vectprime[imin]=vectprime[numrow];
- mat[nmin]=mat[extrarel+numrow];
- vecexpo[nmin]=vecexpo[extrarel+numrow];
- vecerr[nmin]=vecerr[extrarel+numrow];
- numrow--;
- }
- else {*ptrow=numrow;*ptcol=numrow+extrarel;fl=0;}
- }
- if(!numrow) err(talker,"not enough relations or h=1 in hphasereal2");
- if(sizemax>1073741823) return;
- return;
- }
-
- GEN gcdreal(a,b,era,erb,erd)
- GEN a,b,era,erb,*erd;
- {
- long av,tetpil,dec,e;
- GEN k1,r;
-
- if(expo(a)<0) {*erd=gcopy(a);return gcopy(b);}
- if(expo(b)<0) {*erd=gcopy(b);return gcopy(a);}
- av=avma;a=absr(a);b=absr(b);
- while((expo(b)>=0)&&(signe(b)))
- {
- k1=gcvtoi(divrr(a,b),&e);
- r=subrr(a,mulir(k1,b));a=b;b=r;
- }
- tetpil=avma;a=gcopy(a);b=gcopy(b);dec=lpile(av,tetpil,0)>>2;
- *erd=b+dec;return a+dec;
- }
-
- GEN hnfreal(x,detmat,detmatsur2,vlog,verr,wlog,werr)
- GEN x,detmat,detmatsur2,vlog,verr,*wlog,*werr;
- {
- long li,co,av,tetpil,i,j,def,ldef,lim,dec;
- GEN b,q,w,p1,y,d,u,v,p3,p4;
-
- lim=(avma+bot)>>1;
- av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
- def=co;ldef=(li>co)?li-co+1:1;
- for(i=li-1;i>=ldef;i--)
- {
- def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
- while(j>1)
- {
- d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
- p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
- 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]));
- y[j]=(long)clean(y[j],detmat,detmatsur2);
- y[j-1]=(long)p1;
- y[j-1]=(long)clean(y[j-1],detmat,detmatsur2);
- p1=gadd(gmul(u,vlog[j]),gmul(v,vlog[j-1]));
- vlog[j]=lsub(gmul(p3,vlog[j-1]),gmul(p4,vlog[j]));
- vlog[j-1]=(long)p1;
- p1=gadd(gmul(absi(u),verr[j]),gmul(absi(v),verr[j-1]));
- verr[j]=ladd(gmul(gabs(p3),verr[j-1]),gmul(gabs(p4),verr[j]));
- verr[j-1]=(long)p1;
- j--;while(j&&(!signe(coeff(y,i,j)))) j--;
- if(avma<lim)
- {
- tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
- dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
- }
- }
- if(j==1)
- {
- d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
- p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
- y[1]=lsub(gmul(p3=divii(coeff(y,i,1),d),y[def]),gmul(p4=divii(coeff(y,i,def),d),y[1]));
- y[1]=(long)clean(y[1],detmat,detmatsur2);
- y[def]=(long)p1;
- y[def]=(long)clean(y[def],detmat,detmatsur2);
- p1=gadd(gmul(u,vlog[1]),gmul(v,vlog[def]));
- vlog[1]=lsub(gmul(p3,vlog[def]),gmul(p4,vlog[1]));
- vlog[def]=(long)p1;
- p1=gadd(gmul(absi(u),verr[1]),gmul(absi(v),verr[def]));
- verr[1]=ladd(gmul(p3,verr[def]),gmul(p4,verr[1]));
- verr[def]=(long)p1;
- }
- if(avma<lim)
- {
- tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
- dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
- }
- }
- tetpil=avma;y=gcopy(y);vlog=gcopy(vlog);verr=gcopy(verr);
- dec=lpile(av,tetpil,0)>>2;y+=dec;vlog+=dec;verr+=dec;
- b=gcopy(detmat);w=cgetg(li,19);def--;
- for(i=li-1;i>=1;i--)
- {
- d=bezout(coeff(y,i,i+def),b,&u,&v);w[i]=lmod(gmul(u,y[i+def]),b);
- if(!signe(coeff(w,i,i))) coeff(w,i,i)=(long)d;
- vlog[i+def]=lmul(u,vlog[i+def]);verr[i+def]=labs(gmul(u,verr[i+def]));
- if(i>1) b=divii(b,d);
- }
- for(i=li-2;i>=1;i--)
- {
- for(j=i+1;j<li;j++)
- {
- q=gdivent(coeff(w,i,j),coeff(w,i,i));w[j]=lsub(w[j],gmul(q,w[i]));
- vlog[j+def]=lsub(vlog[j+def],gmul(q,vlog[i+def]));
- verr[j+def]=ladd(verr[j+def],gabs(gmul(q,verr[i+def])));
- }
- if(avma<lim)
- {
- tetpil=avma;w=gcopy(w);vlog=gcopy(vlog);verr=gcopy(verr);
- dec=lpile(av,tetpil,0)>>2;w+=dec;vlog+=dec;verr+=dec;
- }
- }
- tetpil=avma;w=gcopy(w);vlog=gcopy(vlog);verr=gcopy(verr);
- dec=lpile(av,tetpil,0)>>2;w+=dec;vlog+=dec;verr+=dec;
- *wlog=vlog;*werr=verr;return w;
- }
-
- GEN comprealform3(x,y,D,isqrtD,sens)
- GEN x,y,D,isqrtD;
- long sens;
-
- {
- long av,tetpil;
- GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,b3,c3,m,z,p1,r;
-
- av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
- d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
- v1=divii(x[1],d1);v2=divii(y[1],d1);
- m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
- r=modii(m,v1);b3=shifti((p1=mulii(v2,r)),1);
- c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
- z=cgetg(4,17);z[1]=lmulii(v1,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v1);
- tetpil=avma;return gerepile(av,tetpil,redrealform3(z,D,isqrtD,sens));
- }
-
- GEN comprealform5(x,y,pterr,D,isqrtD,sqrtD,sens)
- GEN x,y,D,isqrtD,sqrtD;
- long *pterr,sens;
-
- {
- long av,tetpil,ss;
- GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,b3,c3,m,z,p1,r;
-
- av=avma;s=shifti(addii(x[2],y[2]),-1);n=subii(y[2],s);
- d=bezout(y[1],x[1],&y1,&x1);d1=bezout(s,d,&x2,&y2);
- v1=divii(x[1],d1);v2=divii(y[1],d1);
- m=addii(mulii(mulii(y1,y2),n),mulii(y[3],x2));setsigne(m,-signe(m));
- r=modii(m,v1);b3=shifti((p1=mulii(v2,r)),1);
- c3=addii(mulii(y[3],d1),mulii(r,addii(y[2],p1)));
- z=cgetg(6,17);z[1]=lmulii(v1,v2);z[2]=laddii(y[2],b3);z[3]=ldivii(c3,v1);
- z[5]=lmulrr(x[5],y[5]);(*pterr)++;
- if((ss=expo(z[5]))>=EXP220) {z[4]=laddii(addsi(1,x[4]),y[4]);setexpo(z[5],ss-EXP220);}
- else z[4]=laddii(x[4],y[4]);
- tetpil=avma;return gerepile(av,tetpil,redrealform5(z,pterr,D,isqrtD,sqrtD,sens));
- }
-
- GEN sqrealform3(x,D,isqrtD,sens)
- GEN x,D,isqrtD;
- long sens;
-
- {
- long av,tetpil;
- GEN d1,x2,y2,v1,b3,c3,m,z,p1,r;
-
- av=avma;
- d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
- m=mulii(x[3],x2);setsigne(m,-signe(m));
- r=modii(m,v1);b3=shifti((p1=mulii(v1,r)),1);
- c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
- z=cgetg(4,17);z[1]=lmulii(v1,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v1);
- tetpil=avma;return gerepile(av,tetpil,redrealform3(z,D,isqrtD,sens));
- }
-
- GEN sqrealform5(x,pterr,D,isqrtD,sqrtD,sens)
- GEN x,D,isqrtD,sqrtD;
- long *pterr,sens;
- {
- long av,tetpil,ss;
- GEN d1,x2,y2,v1,b3,c3,m,z,p1,r;
-
- av=avma;
- d1=bezout(x[2],x[1],&x2,&y2);v1=divii(x[1],d1);
- m=mulii(x[3],x2);setsigne(m,-signe(m));
- r=modii(m,v1);b3=shifti((p1=mulii(v1,r)),1);
- c3=addii(mulii(x[3],d1),mulii(r,addii(x[2],p1)));
- z=cgetg(6,17);z[1]=lmulii(v1,v1);z[2]=laddii(x[2],b3);z[3]=ldivii(c3,v1);
- z[5]=lmulrr(x[5],x[5]);(*pterr)++;
- if((ss=expo(z[5]))>=EXP220) {z[4]=laddii(addsi(1,x[4]),x[4]);setexpo(z[5],ss-EXP220);}
- else z[4]=lshifti(x[4],1);
- tetpil=avma;return gerepile(av,tetpil,redrealform5(z,pterr,D,isqrtD,sqrtD,sens));
- }
-
- GEN rhorealform3(x,D,isqrtD)
- GEN x,D,isqrtD;
- {
- long av,tetpil,s;
- GEN y,p1,nn;
-
- av=avma;y=cgetg(4,17);y[1]=lcopy(x[3]);
- s=signe(y[1]);setsigne(y[1],1);
- if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
- else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
- p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
- setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN rhorealform5(x,pterr,D,isqrtD,sqrtD)
- GEN x,D,isqrtD,sqrtD;
- long *pterr;
- {
- long av,tetpil,s,ss;
- GEN y,p1,nn;
-
- av=avma;y=cgetg(6,17);y[1]=lcopy(x[3]);
- s=signe(y[1]);setsigne(y[1],1);
- if(cmpii(isqrtD,y[1])>=0) nn=divii(addii(isqrtD,x[2]),p1=shifti(y[1],1));
- else nn=divii(addii(y[1],x[2]),p1=shifti(y[1],1));
- p1=mulii(nn,p1);y[2]=lsubii(p1,x[2]);
- setsigne(y[1],s);p1=shifti(subii(mulii(y[2],y[2]),D),-2);y[3]=ldivii(p1,y[1]);
- y[5]=lmulrr(divrr(addir(x[2],sqrtD),subir(x[2],sqrtD)),x[5]);(*pterr)++;
- if((ss=expo(y[5]))>=EXP220) {y[4]=laddsi(1,x[4]);y[5]=lshiftr(y[5],-EXP220);}
- else y[4]=x[4];
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN redrealform3(x,D,isqrtD,sens)
- GEN x,D,isqrtD;
- long sens;
- {
- long fl,av=avma,tetpil;
- GEN y,p1;
-
- y=cgetg(4,17);y[1]=x[1];y[2]=x[2];y[3]=x[3];
- if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(x[1]),1));
- if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
- }
- while(fl)
- {
- y=rhorealform3(y,D,isqrtD);
- if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(y[1]),1));
- if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
- }
- }
- if(signe(y[1])<0)
- {
- if(sens||(!signe(addii(y[1],y[3]))))
- {
- tetpil=avma;return gerepile(av,tetpil,rhorealform3(y,D,isqrtD));
- }
- else
- {
- tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
- setsigne(y[1],1);setsigne(y[3],-1);return y;
- }
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN redrealform5(x,pterr,D,isqrtD,sqrtD,sens)
- GEN x,D,isqrtD,sqrtD;
- long *pterr,sens;
- {
- long fl,av=avma,tetpil;
- GEN y,p1;
-
- y=cgetg(6,17);y[1]=x[1];y[2]=x[2];y[3]=x[3];y[4]=x[4];y[5]=x[5];
- if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(x[1]),1));
- if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
- }
- while(fl)
- {
- y=rhorealform5(y,pterr,D,isqrtD,sqrtD);
- if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(y[1]),1));
- if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
- }
- }
- if(signe(y[1])<0)
- {
- if(sens||(!signe(addii(y[1],y[3]))))
- {
- tetpil=avma;
- return gerepile(av,tetpil,rhorealform5(y,pterr,D,isqrtD,sqrtD));
- }
- else
- {
- tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
- setsigne(y[1],1);setsigne(y[3],-1);return y;
- }
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN redrealform(x,pterr,D,isqrtD,sqrtD,sens,precreg)
- GEN x,D,isqrtD,sqrtD;
- long *pterr,sens,precreg;
-
- {
- long fl,av=avma,tetpil;
- GEN y,p1;
-
- 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));
- if((signe(x[2])<=0)||(cmpii(x[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(x[1]),1));
- if(signe(p1)<0) fl=(cmpii(x[2],absi(p1))<0);else fl=(cmpii(x[2],p1)<=0);
- }
- while(fl)
- {
- y=rhorealform5(y,pterr,D,isqrtD,sqrtD);
- if((signe(y[2])<=0)||(cmpii(y[2],isqrtD)>0)) fl=1;
- else
- {
- p1=subii(isqrtD,shifti(absi(y[1]),1));
- if(signe(p1)<0) fl=(cmpii(y[2],absi(p1))<0);else fl=(cmpii(y[2],p1)<=0);
- }
- }
- if(signe(y[1])<0)
- {
- if(sens||(!signe(addii(y[1],y[3]))))
- {
- tetpil=avma;
- return gerepile(av,tetpil,rhorealform5(y,pterr,D,isqrtD,sqrtD));
- }
- else
- {
- tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
- setsigne(y[1],1);setsigne(y[3],-1);return y;
- }
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN powrealform5(x,n,D,isqrtD,sqrtD,sens,precreg)
- GEN x,D,isqrtD,sqrtD;
- long n,sens,precreg;
- {
- GEN y,p1;
- long av,tetpil,fl,toto=0;
-
- if(!n)
- {
- y=cgetg(6,17);y[1]=un;
- if(mpodd(x[2])) y[2]=(mpodd(isqrtD)) ? lcopy(isqrtD) : laddsi(-1,isqrtD);
- else {y[2]=lcopy(isqrtD);(((GEN)y[2])[lgef(isqrtD)-1])&=0xfffffffe;}
- av=avma;p1=subii(mulii(y[2],y[2]),D);tetpil=avma;
- y[3]=lpile(av,tetpil,shifti(p1,-2));y[4]=zero;affsr(1,y[5]=lgetr(precreg));
- return y;
- }
- if(n<0)
- {
- p1=cgetg(6,17);p1[1]=x[1];p1[2]=lnegi(x[2]);p1[3]=x[3];
- p1[4]=lnegi(addsi(1,x[4]));p1[5]=lshiftr(divir(un,x[5]),EXP220);
- n= -n;x=p1;
- }
- if(n==1)
- {
- tetpil=avma;
- return gerepile(av,tetpil,redrealform(x,&toto,D,isqrtD,sqrtD,sens,precreg));
- }
- for(fl=0;n>1;n>>=1)
- {
- if(n&1)
- {
- if(fl) y=comprealform5(y,x,&toto,D,isqrtD,sqrtD,sens);
- else {fl=1;y=x;}
- }
- x=sqrealform5(x,&toto,D,isqrtD,sqrtD,sens);
- }
- tetpil=avma;y=fl ? comprealform5(y,x,&toto,D,isqrtD,sqrtD,sens) : gcopy(x);
- return gerepile(av,tetpil,y);
- }
-
- #define LIMP 30000
-
- GEN lfunc(D)
- GEN D;
- {
- GEN y;
- long av=avma,tetpil,prime=0;
- byteptr p=diffptr;
-
- prime=*p++;affsr(1,y=cgetr(4));
- do
- {
- if(!*p) err(recprimer);
- y=mulsr(prime,divrs(y,prime-krogs(D,prime)));
- prime+=*p++;
- }
- while(prime<=LIMP);
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
-