| Karim Belabas on Sun, 22 Mar 1998 04:15:16 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| patch, was Re: polroots algorithm |
Hum, I wrote:
> (I might have overdone the garbage collecting simplification though, I'll
> have to recheck that on _huge_ examples, but that doesn't have high
> priority)
After a first check, I increased the priority of this...
If you are dealing with polynomials of huge degree (> 60, say), you should
apply this patch (if not, you don't need it). Otherwise, the stack overflows
much too easily...
Karim.
========================= patch 15 (2.0.7.alpha) ============================
*** src/basemath/rootpol.c.orig Sat Mar 21 04:31:12 1998
--- src/basemath/rootpol.c Sun Mar 22 03:55:40 1998
***************
*** 238,261 ****
graeffe(GEN p)
{
GEN p0,p1,s0,s1,ss1;
! long n=lgef(p)-3,n0,n1,ns1,i,ltop=avma,lbot,auxi;
! if (n==0) return p;
n0=(n>>1); n1=((n-1)>>1);
! auxi=evalsigne(1)+evalvarn(varn(p));
! p0=cgetg(n0+3,t_POL); p0[1]=auxi+evallgef(n0+3);
! p1=cgetg(n1+3,t_POL); p1[1]=auxi+evallgef(n1+3);
for (i=0; i<=n0; i++) p0[i+2]=p[2+(i<<1)];
for (i=0; i<=n1; i++) p1[i+2]=p[3+(i<<1)];
s0=cook_square(p0);
! s1=cook_square(p1);
! ns1=lgef(s1)-3;
! ss1=cgetg(ns1+4,t_POL); /* will contain x*s1 */
! ss1[1]=auxi+evallgef(ns1+4);
ss1[2]=zero;
! for (i=0; i<=ns1; i++) ss1[3+i]=s1[2+i];
! lbot=avma; return gerepile(ltop,lbot,gsub(s0,ss1));
}
/********************************************************************/
--- 238,260 ----
graeffe(GEN p)
{
GEN p0,p1,s0,s1,ss1;
! long n=lgef(p)-3,n0,n1,i,auxi,ns1;
! if (n==0) return gcopy(p);
n0=(n>>1); n1=((n-1)>>1);
! auxi=evalsigne(1)|evalvarn(varn(p));
! p0=cgetg(n0+3,t_POL); p0[1]=auxi|evallgef(n0+3);
! p1=cgetg(n1+3,t_POL); p1[1]=auxi|evallgef(n1+3);
for (i=0; i<=n0; i++) p0[i+2]=p[2+(i<<1)];
for (i=0; i<=n1; i++) p1[i+2]=p[3+(i<<1)];
s0=cook_square(p0);
! s1=cook_square(p1); ns1 = lgef(s1)-3;
! ss1 = cgetg(ns1+4, t_POL);
! ss1[1] = auxi | evallgef(ns1+4);
ss1[2]=zero;
! for (i=0; i<=ns1; i++) ss1[3+i]=s1[2+i]; /* now ss1 contains x * s1 */
! ss1 = gneg(ss1); return gadd(s0,ss1);
}
/********************************************************************/
***************
*** 822,828 ****
for (j=0; j<=nn; j++) q[2+j]=q[2+valuat+j];
}
set_karasquare_limit(gexpo(q));
! q=graeffe(q); tau2=1.5*tau2; eps=1/log(1./tau2);
}
}
avma=ltop; return gpui(dbltor(2.),dbltor(r),DEFAULTPREC);
--- 821,828 ----
for (j=0; j<=nn; j++) q[2+j]=q[2+valuat+j];
}
set_karasquare_limit(gexpo(q));
! q = gerepileupto(ltop, graeffe(q));
! tau2=1.5*tau2; eps=1/log(1./tau2);
}
}
avma=ltop; return gpui(dbltor(2.),dbltor(r),DEFAULTPREC);
***************
*** 839,851 ****
modulus(GEN p, long k, double tau)
{
GEN q,new_gunr;
! long i,j,kk=k,imax,n=lgef(p)-3,nn,nnn,valuat,ltop=avma,bitprec,decprec,e;
double tau2,r;
tau2=tau/6; nn=n;
bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));
decprec=(long) ((double) bitprec * L2SL10)+1;
new_gunr=gprec(gunr,decprec);
q=gprec(p,decprec);
q=gmul(new_gunr,q);
e=polygone_newton(q,k);
--- 839,852 ----
modulus(GEN p, long k, double tau)
{
GEN q,new_gunr;
! long i,j,kk=k,imax,n=lgef(p)-3,nn,nnn,valuat,av,ltop=avma,bitprec,decprec,e;
double tau2,r;
tau2=tau/6; nn=n;
bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));
decprec=(long) ((double) bitprec * L2SL10)+1;
new_gunr=gprec(gunr,decprec);
+ av = avma;
q=gprec(p,decprec);
q=gmul(new_gunr,q);
e=polygone_newton(q,k);
***************
*** 866,872 ****
nn=nnn-valuat;
set_karasquare_limit(bitprec);
! q=graeffe(q);
e=polygone_newton(q,kk);
r=r+e/exp2((double)i);
q=gmul(new_gunr,q);
--- 867,873 ----
nn=nnn-valuat;
set_karasquare_limit(bitprec);
! q = gerepileupto(av, graeffe(q));
e=polygone_newton(q,kk);
r=r+e/exp2((double)i);
q=gmul(new_gunr,q);
***************
*** 902,908 ****
{
q=eval_rel_pol(q,bitprec);
set_karasquare_limit(bitprec);
! q=graeffe(q);
aux=aux*aux*exp(2*tau2);
tau2=1.5*tau2;
bitprec= (long) ((double) n*(2.+log2(aux)+log2(1./(1-exp(-tau2)))));
--- 903,909 ----
{
q=eval_rel_pol(q,bitprec);
set_karasquare_limit(bitprec);
! q = gerepileupto(ltop, graeffe(q));
aux=aux*aux*exp(2*tau2);
tau2=1.5*tau2;
bitprec= (long) ((double) n*(2.+log2(aux)+log2(1./(1-exp(-tau2)))));
***************
*** 963,969 ****
if (nn==0) return delta_k;
set_karasquare_limit(bitprec);
! q=graeffe(q); tau2=tau2*7./4.;
}
k=-1; logmax=- (double) pariINFINITY;
for (i=0; i<=lgef(q)-3; i++)
--- 964,971 ----
if (nn==0) return delta_k;
set_karasquare_limit(bitprec);
! q = gerepileupto(ltop, graeffe(q));
! tau2=tau2*7./4.;
}
k=-1; logmax=- (double) pariINFINITY;
for (i=0; i<=lgef(q)-3; i++)
***************
*** 1363,1375 ****
Lmax=4; while (Lmax<=n) Lmax=(Lmax<<1);
parameters(pp,&mu,&gamma,polreal,param,param2);
! FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(k+3);
FF[k+2]=un;
- H=cgetg(k+2,t_POL); H[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(k+2);
NN=(long) (0.5/delta); NN+=(NN%2); if (NN<2) NN=2;
NN=NN*Lmax; ltop=avma;
-
for(;;)
{
bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8;
--- 1365,1376 ----
Lmax=4; while (Lmax<=n) Lmax=(Lmax<<1);
parameters(pp,&mu,&gamma,polreal,param,param2);
! H =cgetg(k+2,t_POL); H[1] =evalsigne(1) | evalvarn(varn(p)) | evallgef(k+2);
! FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
FF[k+2]=un;
NN=(long) (0.5/delta); NN+=(NN%2); if (NN<2) NN=2;
NN=NN*Lmax; ltop=avma;
for(;;)
{
bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8;
***************
*** 1389,1395 ****
if (k<=n/2) split_fromU(p,k,delta,bitprec,F,G,param,param2);
else
! { /* we start from the reciprocal of p */
q=cgetg(n+3,t_POL); q[1]=p[1];
for (i=0; i<=n; i++) q[i+2]=p[n-i+2];
split_fromU(q,n-k,delta,bitprec,&FF,&GG,param,param2);
--- 1390,1396 ----
if (k<=n/2) split_fromU(p,k,delta,bitprec,F,G,param,param2);
else
! { /* we start from the reciprocal of p */
q=cgetg(n+3,t_POL); q[1]=p[1];
for (i=0; i<=n; i++) q[i+2]=p[n-i+2];
split_fromU(q,n-k,delta,bitprec,&FF,&GG,param,param2);
***************
*** 1431,1445 ****
static GEN
shiftpol(GEN p, GEN b)
{
! long ltop=avma,i,lbot=avma;
! GEN r,q=gzero;
for (i=lgef(p)-1; i>=2; i--)
{
! r=gmul(q,b);
! lbot=avma; q=gadd((GEN)p[i],r);
}
! return gerepile(ltop,lbot,q);
}
/* return (aX-1)^n * p[ (X-a) / (aX-1) ] */
--- 1432,1446 ----
static GEN
shiftpol(GEN p, GEN b)
{
! long av = avma,i, limit = (av+bot)>>1;
! GEN q = gzero;
for (i=lgef(p)-1; i>=2; i--)
{
! q = gadd((GEN)p[i], gmul(q,b));
! if (low_stack(limit, (av+bot)>>1)) q = gerepileupto(av,q);
}
! return gerepileupto(av,q);
}
/* return (aX-1)^n * p[ (X-a) / (aX-1) ] */
***************
*** 1449,1459 ****
GEN r,pui,num,aux;
long n=lgef(p)-3, i;
! pui=cgetg(4,t_POL); pui[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(4);
! pui[2]=(long) mygprec(gneg(gunr),bitprec); pui[3]=lconj(a);
! aux=pui;
! num=cgetg(4,t_POL); num[1]=pui[1];
! num[2]=lneg(a); num[3]=(long) mygprec(gunr,bitprec);
r=(GEN)p[2+n];
for (i=n-1; ; i--)
{
--- 1450,1463 ----
GEN r,pui,num,aux;
long n=lgef(p)-3, i;
! aux = pui = cgetg(4,t_POL);
! pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4);
! pui[2] = (long) mygprec(gneg(gunr),bitprec);
! pui[3] = lconj(a);
! num = cgetg(4,t_POL);
! num[1] = pui[1];
! num[2] = lneg(a);
! num[3] = (long) mygprec(gunr,bitprec);
r=(GEN)p[2+n];
for (i=n-1; ; i--)
{
***************
*** 1467,1474 ****
static void
conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
{
! long bitprec2,bitprec3,n=lgef(p)-3,decprec,i;
! GEN q,FF,GG,a,R;
double rmin,rmax,rho,delta,aux2,param,param2,r1,r2;
bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
--- 1471,1478 ----
static void
conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
{
! long bitprec2,bitprec3,n=lgef(p)-3,decprec,i,ltop = avma, av;
! GEN q,FF,GG,a,R, *gptr[2];
double rmin,rmax,rho,delta,aux2,param,param2,r1,r2;
bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
***************
*** 1476,1483 ****
a=gmul(mygprec(a,bitprec2),mygprec(globalu,bitprec2));
a=gdivgs(a,-6); /* a=-globalu/2/sqrt(3) */
! q=mygprec(p,bitprec2);
! q=conformal_pol(q,a,bitprec2);
for (i=1; i<=n; i++)
if (radius[i]!=0.) /* updating array radius */
{
--- 1480,1487 ----
a=gmul(mygprec(a,bitprec2),mygprec(globalu,bitprec2));
a=gdivgs(a,-6); /* a=-globalu/2/sqrt(3) */
! av = avma; q=mygprec(p,bitprec2);
! q = conformal_pol(q,a,bitprec2);
for (i=1; i<=n; i++)
if (radius[i]!=0.) /* updating array radius */
{
***************
*** 1514,1519 ****
--- 1518,1524 ----
R=mygprec(dbltor(1/rho),bitprec2);
q=scalepol(q,R,bitprec2);
+ gptr[0]=&q; gptr[1]=&R; gerepilemany(av,gptr,2);
aux2=radius[k];
for (i=k-1; i>=1; i--)
***************
*** 1536,1546 ****
if (aux2<1.) param2-=log2(aux2);
}
optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
! bitprec2=bitprec2+n; R=ginv(R);
! FF=scalepol(FF,R,bitprec2); GG=scalepol(GG,R,bitprec2);
a=mygprec(a,bitprec2);
! FF=conformal_pol(FF,a,bitprec2); GG=conformal_pol(GG,a,bitprec2);
a=ginv(gsub(gun,gmul(a,gconj(a))));
a=glog(a,(long) (bitprec2 * L2SL10)+1);
--- 1541,1553 ----
if (aux2<1.) param2-=log2(aux2);
}
optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
! bitprec2=bitprec2+n; R = ginv(R);
! FF=scalepol(FF,R,bitprec2);
! GG=scalepol(GG,R,bitprec2);
a=mygprec(a,bitprec2);
! FF=conformal_pol(FF,a,bitprec2);
! GG=conformal_pol(GG,a,bitprec2);
a=ginv(gsub(gun,gmul(a,gconj(a))));
a=glog(a,(long) (bitprec2 * L2SL10)+1);
***************
*** 1549,1557 ****
GG=gmul(GG,gexp(gmulgs(a,n-k),decprec));
*F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
}
! /* split p, this time with no scaling . returns in F and G two polynomials
such that |p-FG|< 2^(-bitprec)|p| */
static void
split_2(GEN p, long bitprec, double thickness, GEN *F, GEN *G)
--- 1556,1565 ----
GG=gmul(GG,gexp(gmulgs(a,n-k),decprec));
*F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
+ gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2);
}
! /* split p, this time with no scaling. returns in F and G two polynomials
such that |p-FG|< 2^(-bitprec)|p| */
static void
split_2(GEN p, long bitprec, double thickness, GEN *F, GEN *G)
***************
*** 1699,1706 ****
rmin=min_modulus(qq,0.05);
if (3./rmin > thickness)
{
! rmax=max_modulus(qq,0.05);
! quo = rmax/rmin;
if ((float)quo > (float)thickness)
{
thickness=quo; newq=qq; globalu=(GEN)v[i];
--- 1707,1713 ----
rmin=min_modulus(qq,0.05);
if (3./rmin > thickness)
{
! rmax=max_modulus(qq,0.05); quo = rmax/rmin;
if ((float)quo > (float)thickness)
{
thickness=quo; newq=qq; globalu=(GEN)v[i];
***************
*** 1798,1808 ****
if (k>0)
{
if (k>n/2) k=n/2;
! FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(k+3);
for (i=0; i<k; i++) FF[i+2]=lcopy(gzero);
FF[k+2]=(long) mygprec(gunr,bitprec);
GG=cgetg(n-k+3,t_POL);
! GG[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(n-k+3);
for (i=0; i<=n-k; i++) GG[i+2]=p[i+k+2];
}
else
--- 1805,1816 ----
if (k>0)
{
if (k>n/2) k=n/2;
! FF=cgetg(k+3,t_POL);
! FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
for (i=0; i<k; i++) FF[i+2]=lcopy(gzero);
FF[k+2]=(long) mygprec(gunr,bitprec);
GG=cgetg(n-k+3,t_POL);
! GG[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(n-k+3);
for (i=0; i<=n-k; i++) GG[i+2]=p[i+k+2];
}
else
***************
*** 1939,1945 ****
{
long n=lgef(p)-3,decprec,ltop,lbot;
GEN F,G,a,b,m1,m2,m;
! GEN *gptr[3];
if (n==1)
{
--- 1947,1953 ----
{
long n=lgef(p)-3,decprec,ltop,lbot;
GEN F,G,a,b,m1,m2,m;
! GEN *gptr[2];
if (n==1)
{
--
Karim Belabas e-mail:
Max-Planck-Institut fuer Mathematik karim@mpim-bonn.mpg.de
Gottfried-Claren-Str. 26 tel:
53225 Bonn (Germany) (00 49 228) 402-245