version 1.1, 2001/10/02 11:17:05 |
version 1.2, 2002/09/11 07:26:52 |
Line 22 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 22 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
#include "pari.h" |
#include "pari.h" |
|
|
extern GEN polrecip_i(GEN x); |
extern GEN polrecip_i(GEN x); |
extern GEN poldeflate(GEN x0, long *m); |
|
extern GEN roots_to_pol(GEN a, long v); |
|
#define pariINFINITY 100000 |
#define pariINFINITY 100000 |
#define NEWTON_MAX 10 |
#define NEWTON_MAX 10 |
|
|
Line 65 quickmulcc(GEN x, GEN y) |
|
Line 63 quickmulcc(GEN x, GEN y) |
|
} |
} |
if (ty==t_COMPLEX) |
if (ty==t_COMPLEX) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN p1,p2; |
GEN p1,p2; |
|
|
z=cgetg(3,t_COMPLEX); av=avma; |
z=cgetg(3,t_COMPLEX); av=avma; |
|
|
mysquare(GEN p) |
mysquare(GEN p) |
{ |
{ |
GEN s,aux1,aux2; |
GEN s,aux1,aux2; |
long i,j,n=degpol(p),nn,ltop,lbot; |
long i, j, n=degpol(p), nn; |
|
gpmem_t ltop, lbot; |
|
|
if (n==-1) return gcopy(p); |
if (n==-1) return gcopy(p); |
nn=n<<1; s=cgetg(nn+3,t_POL); |
nn=n<<1; s=cgetg(nn+3,t_POL); |
Line 142 karasquare(GEN p) |
|
Line 141 karasquare(GEN p) |
|
{ |
{ |
GEN p1,s0,s1,s2,aux; |
GEN p1,s0,s1,s2,aux; |
long n=degpol(p),n0,n1,i,var,nn0; |
long n=degpol(p),n0,n1,i,var,nn0; |
ulong ltop; |
gpmem_t ltop; |
|
|
if (n<=KARASQUARE_LIMIT) return mysquare(p); |
if (n<=KARASQUARE_LIMIT) return mysquare(p); |
ltop=avma; |
ltop=avma; |
Line 174 cook_square(GEN p) |
|
Line 173 cook_square(GEN p) |
|
{ |
{ |
GEN p0,p1,p2,p3,q,aux0,aux1,r,aux,plus,moins; |
GEN p0,p1,p2,p3,q,aux0,aux1,r,aux,plus,moins; |
long n=degpol(p),n0,n3,i,j,var; |
long n=degpol(p),n0,n3,i,j,var; |
ulong ltop = avma; |
gpmem_t ltop = avma; |
|
|
if (n<=COOK_SQUARE_LIMIT) return karasquare(p); |
if (n<=COOK_SQUARE_LIMIT) return karasquare(p); |
|
|
|
|
return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG; |
return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG; |
} |
} |
|
|
static double |
double |
mylog2(GEN z) |
mylog2(GEN z) |
{ |
{ |
double x,y; |
double x,y; |
|
|
return x+0.5*log2( 1 + exp2(2*(y-x))); |
return x+0.5*log2( 1 + exp2(2*(y-x))); |
} |
} |
|
|
|
/* find s such that A_h <= 2^s <= 2 A_i for one h and all i < n = deg(p), |
|
* with A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and p = sum p_i X^i */ |
static long |
static long |
findpower(GEN p) |
findpower(GEN p) |
{ |
{ |
double x, logbinomial,pente,pentemax=-pariINFINITY; |
double x, logbinomial, mins = pariINFINITY; |
long n=degpol(p),i; |
long n = degpol(p),i; |
|
|
logbinomial = mylog2((GEN) p[n+2]); |
logbinomial = mylog2((GEN)p[n+2]); /* log2(lc * binom(n,i)) */ |
for (i=n-1; i>=0; i--) |
for (i=n-1; i>=0; i--) |
{ |
{ |
logbinomial += log2((double) (i+1) / (double) (n-i)); |
logbinomial += log2((double) (i+1) / (double) (n-i)); |
x = mylog2((GEN) p[2+i])-logbinomial; |
x = mylog2((GEN)p[i+2]); |
if (x>-pariINFINITY) |
if (x != -pariINFINITY) |
{ |
{ |
pente = x/ (double) (n-i); |
double s = (logbinomial - x) / (double) (n-i); |
if (pente > pentemax) pentemax = pente; |
if (s < mins) mins = s; |
} |
} |
} |
} |
return (long) -floor(pentemax); |
i = (long)ceil(mins); |
|
if (i - mins > 1 - 1e-12) i--; |
|
return i; |
} |
} |
|
|
/* returns the exponent for the procedure modulus, from the newton diagram */ |
/* returns the exponent for the procedure modulus, from the newton diagram */ |
static long |
static long |
polygone_newton(GEN p, long k) |
newton_polygon(GEN p, long k) |
{ |
{ |
double *logcoef,pente; |
double *logcoef,pente; |
long n=degpol(p),i,j,h,l,*sommet,pentelong; |
long n=degpol(p),i,j,h,l,*vertex,pentelong; |
|
|
logcoef=(double*) gpmalloc((n+1)*sizeof(double)); |
logcoef = (double*) gpmalloc((n+1)*sizeof(double)); |
sommet=(long*) gpmalloc((n+1)*sizeof(long)); |
vertex = (long*) gpmalloc((n+1)*sizeof(long)); |
|
|
/* sommet[i]=1 si i est un sommet, =0 sinon */ |
/* vertex[i]=1 if i a vertex of convex hull, 0 otherwise */ |
for (i=0; i<=n; i++) { logcoef[i]=mylog2((GEN)p[2+i]); sommet[i]=0; } |
for (i=0; i<=n; i++) { logcoef[i] = mylog2((GEN)p[2+i]); vertex[i] = 0; } |
sommet[0]=1; i=0; |
vertex[0]=1; |
while (i<n) |
for (i=0; i < n; i=h) |
{ |
{ |
pente=logcoef[i+1]-logcoef[i]; |
pente = logcoef[i+1]-logcoef[i]; |
h=i+1; |
h = i+1; |
for (j=i+1; j<=n; j++) |
for (j=i+1; j<=n; j++) |
{ |
if (pente < (logcoef[j]-logcoef[i])/(double)(j-i)) |
if (pente<(logcoef[j]-logcoef[i])/(double)(j-i)) |
|
{ |
{ |
h=j; |
h = j; |
pente=(logcoef[j]-logcoef[i])/(double)(j-i); |
pente = (logcoef[j]-logcoef[i])/(double)(j-i); |
} |
} |
} |
vertex[h] = 1; |
i=h; |
|
sommet[h]=1; |
|
} |
} |
h=k; while (!sommet[h]) h++; |
h = k; while (!vertex[h]) h++; |
l=k-1; while (!sommet[l]) l--; |
l = k-1; while (!vertex[l]) l--; |
pentelong=(long) floor((logcoef[h]-logcoef[l])/(double)(h-l)+0.5); |
pentelong = (long) floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5); |
free(logcoef); free(sommet); return pentelong; |
free(logcoef); free(vertex); return pentelong; |
} |
} |
|
|
/* change z into z*2^e, where z is real or complex of real */ |
/* change z into z*2^e, where z is real or complex of real */ |
Line 444 myshiftic(GEN z, long e) |
|
Line 444 myshiftic(GEN z, long e) |
|
static GEN |
static GEN |
myrealun(long bitprec) |
myrealun(long bitprec) |
{ |
{ |
GEN x; |
|
if (bitprec < 0) bitprec = 0; |
if (bitprec < 0) bitprec = 0; |
x = cgetr(bitprec/BITS_IN_LONG+3); |
return realun(bitprec/BITS_IN_LONG+3); |
affsr(1,x); return x; |
|
} |
} |
|
|
static GEN |
static GEN |
Line 547 homothetie(GEN p, GEN R, long bitprec) |
|
Line 545 homothetie(GEN p, GEN R, long bitprec) |
|
return r; |
return r; |
} |
} |
|
|
/* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) */ |
/* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) [ ~as above with R = 2^-e ]*/ |
static void |
static void |
homothetie2n(GEN p, long e) |
homothetie2n(GEN p, long e) |
{ |
{ |
Line 569 homothetie_gauss(GEN p, long e,long f) |
|
Line 567 homothetie_gauss(GEN p, long e,long f) |
|
} |
} |
} |
} |
|
|
static long |
|
valuation(GEN p) |
|
{ |
|
long j=0,n=degpol(p); |
|
|
|
while ((j<=n) && isexactzero((GEN)p[j+2])) j++; |
|
return j; |
|
} |
|
|
|
/* provides usually a good lower bound on the largest modulus of the roots, |
/* provides usually a good lower bound on the largest modulus of the roots, |
puts in k an upper bound of the number of roots near the largest root |
puts in k an upper bound of the number of roots near the largest root |
at a distance eps */ |
at a distance eps */ |
static double |
static double |
lower_bound(GEN p, long *k, double eps) |
lower_bound(GEN p, long *k, double eps) |
{ |
{ |
long n=degpol(p),i,j,ltop=avma; |
long n=degpol(p), i, j; |
GEN a,s,icd; |
gpmem_t ltop=avma; |
|
GEN a,s,S,icd; |
double r,*rho; |
double r,*rho; |
|
|
if (n<4) { *k=n; return 0.; } |
if (n<4) { *k=n; return 0.; } |
a=cgetg(6,t_POL); s=cgetg(6,t_POL); |
S = cgetg(5,t_VEC); |
|
a = cgetg(5,t_VEC); |
rho=(double *) gpmalloc(4*sizeof(double)); |
rho=(double *) gpmalloc(4*sizeof(double)); |
icd = gdiv(realun(DEFAULTPREC), (GEN) p[n+2]); |
icd = gdiv(realun(DEFAULTPREC), (GEN) p[n+2]); |
for (i=1; i<=4; i++) a[i+1]=lmul(icd,(GEN)p[n+2-i]); |
for (i=1; i<=4; i++) a[i] = lmul(icd,(GEN)p[n+2-i]); |
for (i=1; i<=4; i++) |
for (i=1; i<=4; i++) |
{ |
{ |
s[i+1]=lmulsg(i,(GEN)a[i+1]); |
s = gmulsg(i,(GEN)a[i]); |
for (j=1; j<i; j++) |
for (j=1; j<i; j++) |
s[i+1]=ladd((GEN)s[i+1],gmul((GEN)s[j+1],(GEN)a[i+1-j])); |
s = gadd(s, gmul((GEN)S[j],(GEN)a[i-j])); |
s[i+1]=lneg((GEN)s[i+1]); |
S[i] = lneg(s); /* Newton sum S_i */ |
r=gtodouble(gabs((GEN) s[i+1],3)); |
r = gtodouble(gabs(s,3)); |
if (r<=0.) /* should not be strictly negative */ |
if (r == 0.) |
rho[i-1]=0.; |
rho[i-1] = 0.; |
else |
else |
rho[i-1]=exp(log(r/(double)n)/(double) i); |
rho[i-1] = exp(log(r/(double)n)/(double) i); |
} |
} |
r=0.; |
r = 0.; |
for (i=0; i<4; i++) if (r<rho[i]) r=rho[i]; |
for (i=0; i<4; i++) if (r<rho[i]) r=rho[i]; |
if (r>0. && eps<1.2) |
if (r>0. && eps<1.2) |
*k=(long) floor((n*rho[0]/r+n)/(1+exp(-eps)*cos(eps))); |
*k = (long) floor((n*rho[0]/r+n)/(1+exp(-eps)*cos(eps))); |
else |
else |
*k=n; |
*k = n; |
free(rho); avma=ltop; return r; |
free(rho); avma = ltop; return r; |
} |
} |
|
|
/* returns the maximum of the modulus of p with a relative error tau */ |
/* returns the maximum of the modulus of p with a relative error tau */ |
static GEN |
GEN |
max_modulus(GEN p, double tau) |
max_modulus(GEN p, double tau) |
{ |
{ |
GEN q,aux,gunr; |
GEN r,q,aux,gunr; |
long i,j,k,valuat,n=degpol(p),nn,ltop=avma,bitprec,imax,e; |
gpmem_t av, ltop = avma; |
double r,rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.; |
long i,k,n=degpol(p),nn,bitprec,M,e; |
|
double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.; |
|
|
|
r = cgeti(BIGDEFAULTPREC); |
|
av = avma; |
|
|
eps = - 1/log(1.5*tau2); /* > 0 */ |
eps = - 1/log(1.5*tau2); /* > 0 */ |
bitprec=(long) ((double) n*log2(1./tau2)+3*log2((double) n))+1; |
bitprec = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1; |
gunr=myrealun(bitprec+2*n); |
gunr = myrealun(bitprec+2*n); |
aux=gdiv(gunr,(GEN) p[2+n]); |
aux = gdiv(gunr,(GEN) p[2+n]); |
q=gmul(aux,p); q[2+n]=lcopy(gunr); |
q = gmul(aux,p); q[2+n] = (long)gunr; |
k=nn=n; |
e = findpower(q); |
e=findpower(q); homothetie2n(q,e); r=-(double) e; |
homothetie2n(q,e); |
q=mygprec(q,bitprec+(n<<1)); |
affsi(e, r); |
|
q = mygprec(q,bitprec+(n<<1)); |
pol_to_gaussint(q,bitprec); |
pol_to_gaussint(q,bitprec); |
imax=(long) ((log(log(4.*n)/(2*tau2))) / log(2.)) + 2; |
M = (long) (log2( log(4.*n) / (2*tau2) )) + 2; |
|
nn = n; |
for (i=0,e=0;;) |
for (i=0,e=0;;) |
{ |
{ /* nn = deg(q) */ |
rho=lower_bound(q,&k,eps); |
rho = lower_bound(q,&k,eps); |
if (rho>exp2(-(double) e)) e = (long) -floor(log2(rho)); |
if (rho > exp2(-(double) e)) e = (long) -floor(log2(rho)); |
r -= e / exp2((double)i); |
affii(shifti(addis(r,e), 1), r); |
if (++i == imax) { |
if (++i == M) break; |
avma=ltop; |
|
return gpui(dbltor(2.),dbltor(r),DEFAULTPREC); |
|
} |
|
|
|
if (k<nn) |
bitprec = (long) ((double) k* log2(1./tau2) + |
bitprec=(long) ((double) k* log2(1./tau2)+ |
(double) (nn-k)*log2(1./eps) + 3*log2((double) nn)) + 1; |
(double) (nn-k)*log2(1./eps)+ |
|
3*log2((double) nn))+1; |
|
else |
|
bitprec=(long) ((double) nn* log2(1./tau2)+ |
|
3.*log2((double) nn))+1; |
|
homothetie_gauss(q,e,bitprec-(long)floor(mylog2((GEN) q[2+nn])+0.5)); |
homothetie_gauss(q,e,bitprec-(long)floor(mylog2((GEN) q[2+nn])+0.5)); |
valuat=valuation(q); |
nn -= polvaluation(q, &q); |
if (valuat>0) |
|
{ |
|
nn-=valuat; setlgef(q,nn+3); |
|
for (j=0; j<=nn; j++) q[2+j]=q[2+valuat+j]; |
|
} |
|
set_karasquare_limit(gexpo(q)); |
set_karasquare_limit(gexpo(q)); |
q = gerepileupto(ltop, graeffe(q)); |
q = gerepileupto(av, graeffe(q)); |
tau2=1.5*tau2; eps=1/log(1./tau2); |
tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5; |
|
eps = -1/log(tau2); /* > 0 */ |
e = findpower(q); |
e = findpower(q); |
} |
} |
|
if (!signe(r)) { avma = ltop; return realun(DEFAULTPREC); } |
|
r = itor(r, DEFAULTPREC); |
|
setexpo(r, expo(r) - M); |
|
rho = rtodbl(r); |
|
/* rho = sum e_i 2^-i */ |
|
avma = ltop; |
|
return mpexp(dbltor(-rho * LOG2)); /* 2^-r */ |
} |
} |
|
|
/* return the k-th modulus (in ascending order) of p, rel. error tau*/ |
/* return the k-th modulus (in ascending order) of p, rel. error tau*/ |
|
|
modulus(GEN p, long k, double tau) |
modulus(GEN p, long k, double tau) |
{ |
{ |
GEN q,gunr; |
GEN q,gunr; |
long i,j,kk=k,imax,n=degpol(p),nn,nnn,valuat,av,ltop=avma,bitprec,decprec,e; |
long i, kk=k, imax, n=degpol(p), nn, bitprec, decprec, e; |
|
gpmem_t av, ltop=avma; |
double tau2,r; |
double tau2,r; |
|
|
tau2=tau/6; nn=n; |
tau2 = tau/6; |
bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2))); |
bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2))); |
decprec=(long) ((double) bitprec * L2SL10)+1; |
decprec=(long) ((double) bitprec * L2SL10)+1; |
gunr=myrealun(bitprec); |
gunr=myrealun(bitprec); |
av = avma; |
av = avma; |
q=gprec(p,decprec); |
q = gprec(p,decprec); |
q=gmul(gunr,q); |
q = gmul(gunr,q); |
e=polygone_newton(q,k); |
e = newton_polygon(q,k); |
homothetie2n(q,e); |
homothetie2n(q,e); |
r=(double) e; |
r = (double) e; |
imax=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1; |
imax=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1; |
for (i=1; i<imax; i++) |
for (i=1; i<imax; i++) |
{ |
{ |
q=eval_rel_pol(q,bitprec); |
q = eval_rel_pol(q,bitprec); |
|
kk -= polvaluation(q, &q); |
|
nn = degpol(q); |
|
|
nnn=degpol(q); valuat=valuation(q); |
|
if (valuat>0) |
|
{ |
|
kk-=valuat; |
|
for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j]; |
|
setlgef(q,nnn-valuat+3); |
|
} |
|
nn=nnn-valuat; |
|
|
|
set_karasquare_limit(bitprec); |
set_karasquare_limit(bitprec); |
q = gerepileupto(av, graeffe(q)); |
q = gerepileupto(av, graeffe(q)); |
e=polygone_newton(q,kk); |
e = newton_polygon(q,kk); |
r += e / exp2((double)i); |
r += e / exp2((double)i); |
q=gmul(gunr,q); |
q=gmul(gunr,q); |
homothetie2n(q,e); |
homothetie2n(q,e); |
|
|
tau2=1.5*tau2; if (tau2>1.) tau2=1.; |
tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.; |
bitprec= 1+(long) ((double) nn*(2.+log2(3.*(double) nn)+log2(1./tau2))); |
bitprec= 1+(long) ((double) nn*(2.+log2(3.*(double) nn)+log2(1./tau2))); |
} |
} |
avma=ltop; return mpexp(dbltor(-r * LOG2)); |
avma=ltop; return mpexp(dbltor(-r * LOG2)); |
|
|
pre_modulus(GEN p, long k, double tau, GEN rmin, GEN rmax) |
pre_modulus(GEN p, long k, double tau, GEN rmin, GEN rmax) |
{ |
{ |
GEN R, q, aux; |
GEN R, q, aux; |
long n=degpol(p),i,imax,imax2,bitprec,ltop=avma, av; |
long n=degpol(p), i, imax, imax2, bitprec; |
|
gpmem_t ltop=avma, av; |
double tau2, aux2; |
double tau2, aux2; |
|
|
tau2=tau/6.; |
tau2=tau/6.; |
Line 750 pre_modulus(GEN p, long k, double tau, GEN rmin, GEN r |
|
Line 737 pre_modulus(GEN p, long k, double tau, GEN rmin, GEN r |
|
static GEN |
static GEN |
min_modulus(GEN p, double tau) |
min_modulus(GEN p, double tau) |
{ |
{ |
long av=avma; |
gpmem_t av=avma; |
GEN r; |
GEN r; |
|
|
if (isexactzero((GEN)p[2])) return realzero(3); |
if (isexactzero((GEN)p[2])) return realzero(3); |
|
|
dual_modulus(GEN p, GEN R, double tau, long l) |
dual_modulus(GEN p, GEN R, double tau, long l) |
{ |
{ |
GEN q; |
GEN q; |
long i,j,imax,k,delta_k=0,n=degpol(p),nn,nnn,valuat,ltop=avma,bitprec,ll=l; |
long i, imax, k, delta_k=0, n=degpol(p), nn, nnn, valuat, bitprec, ll=l; |
|
gpmem_t ltop=avma; |
double logmax,aux,tau2; |
double logmax,aux,tau2; |
|
|
tau2=7.*tau/8.; |
tau2=7.*tau/8.; |
Line 778 dual_modulus(GEN p, GEN R, double tau, long l) |
|
Line 766 dual_modulus(GEN p, GEN R, double tau, long l) |
|
bitprec=6*nn-5*ll+(long) ((double) nn*(log2(1/tau2)+8.*tau2/7.)); |
bitprec=6*nn-5*ll+(long) ((double) nn*(log2(1/tau2)+8.*tau2/7.)); |
|
|
q=eval_rel_pol(q,bitprec); |
q=eval_rel_pol(q,bitprec); |
nnn=degpol(q); valuat=valuation(q); |
nnn = degpol(q); valuat = polvaluation(q, &q); |
if (valuat>0) |
|
{ |
|
delta_k+=valuat; |
|
for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j]; |
|
setlgef(q,nnn-valuat+3); |
|
} |
|
ll= (-valuat<nnn-n)? ll-valuat: ll+nnn-n; /* min(ll-valuat,ll+nnn-n) */ |
ll= (-valuat<nnn-n)? ll-valuat: ll+nnn-n; /* min(ll-valuat,ll+nnn-n) */ |
if (ll<0) ll=0; |
if (ll<0) ll=0; |
|
|
nn=nnn-valuat; |
nn = nnn-valuat; delta_k += valuat; |
if (nn==0) return delta_k; |
if (nn==0) return delta_k; |
|
|
set_karasquare_limit(bitprec); |
set_karasquare_limit(bitprec); |
q = gerepileupto(ltop, graeffe(q)); |
q = gerepileupto(ltop, graeffe(q)); |
tau2=tau2*7./4.; |
tau2=tau2*7./4.; |
} |
} |
k=-1; logmax=- (double) pariINFINITY; |
k = -1; logmax = - (double)pariINFINITY; |
for (i=0; i<=degpol(q); i++) |
for (i=0; i<=degpol(q); i++) |
{ |
{ |
aux=mylog2((GEN)q[2+i]); |
aux=mylog2((GEN)q[2+i]); |
|
|
static void |
static void |
fft(GEN Omega, GEN p, GEN f, long Step, long l) |
fft(GEN Omega, GEN p, GEN f, long Step, long l) |
{ |
{ |
ulong ltop; |
gpmem_t ltop; |
long i,l1,l2,l3,rap=Lmax/l,rapi,Step4; |
long i,l1,l2,l3,rap=Lmax/l,rapi,Step4; |
GEN f1,f2,f3,f02,f13,g02,g13,ff; |
GEN f1,f2,f3,f02,f13,g02,g13,ff; |
|
|
Line 885 fft(GEN Omega, GEN p, GEN f, long Step, long l) |
|
Line 867 fft(GEN Omega, GEN p, GEN f, long Step, long l) |
|
for (i=0; i<l; i++) f[i]=ff[i+1]; |
for (i=0; i<l; i++) f[i]=ff[i+1]; |
} |
} |
|
|
extern void mpsincos(GEN x, GEN *s, GEN *c); |
|
|
|
/* return exp(ix), x a t_REAL */ |
|
static GEN |
|
exp_i(GEN x) |
|
{ |
|
GEN v; |
|
|
|
if (!signe(x)) return realun(lg(x)); /* should not happen */ |
|
v = cgetg(3,t_COMPLEX); |
|
mpsincos(x, (GEN*)(v+2), (GEN*)(v+1)); |
|
return v; |
|
} |
|
|
|
/* e(1/N) */ |
/* e(1/N) */ |
static GEN |
static GEN |
RUgen(long N, long bitprec) |
RUgen(long N, long bitprec) |
Line 907 RUgen(long N, long bitprec) |
|
Line 875 RUgen(long N, long bitprec) |
|
if (N == 2) return mpneg(realun(bitprec)); |
if (N == 2) return mpneg(realun(bitprec)); |
if (N == 4) return gi; |
if (N == 4) return gi; |
pi2 = gmul2n(mppi(bitprec/BITS_IN_LONG+3), 1); |
pi2 = gmul2n(mppi(bitprec/BITS_IN_LONG+3), 1); |
return exp_i(gdivgs(pi2,N)); |
return exp_Ir(gdivgs(pi2,N)); |
} |
} |
|
|
/* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */ |
/* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */ |
Line 947 initRUgen(long N, long bitprec) |
|
Line 915 initRUgen(long N, long bitprec) |
|
} |
} |
|
|
/* returns 1 if p has only real coefficients, 0 else */ |
/* returns 1 if p has only real coefficients, 0 else */ |
static long |
static int |
isreal(GEN p) |
isreal(GEN p) |
{ |
{ |
long n=degpol(p),i=0; |
long n=degpol(p),i=0; |
|
|
|
|
static void |
static void |
parameters(GEN p, double *mu, double *gamma, |
parameters(GEN p, double *mu, double *gamma, |
long polreal, double param, double param2) |
int polreal, double param, double param2) |
{ |
{ |
GEN q,pc,Omega,coef,RU,prim,aux,aux0,ggamma,gx,mygpi; |
GEN q,pc,Omega,coef,RU,prim,aux,aux0,ggamma,gx,mygpi; |
long ltop=avma,limite=stack_lim(ltop,1),n=degpol(p),bitprec,NN,K,i,j,ltop2; |
long n=degpol(p), bitprec, NN, K, i, j; |
|
gpmem_t ltop2, ltop=avma, limite=stack_lim(ltop, 1); |
double lx; |
double lx; |
|
|
bitprec=gexpo(p)+(long)param2+8; |
bitprec=gexpo(p)+(long)param2+8; |
Line 969 parameters(GEN p, double *mu, double *gamma, |
|
Line 938 parameters(GEN p, double *mu, double *gamma, |
|
K=NN/Lmax; if (K%2==1) K++; NN=Lmax*K; |
K=NN/Lmax; if (K%2==1) K++; NN=Lmax*K; |
mygpi=mppi(bitprec/BITS_IN_LONG+3); |
mygpi=mppi(bitprec/BITS_IN_LONG+3); |
aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */ |
aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */ |
prim = exp_i(aux); |
prim = exp_Ir(aux); |
aux = gmulbyi(aux); |
aux = gmulbyi(aux); |
RU = myrealun(bitprec); |
RU = myrealun(bitprec); |
|
|
Line 1026 parameters(GEN p, double *mu, double *gamma, |
|
Line 995 parameters(GEN p, double *mu, double *gamma, |
|
static void |
static void |
dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H, long polreal) |
dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H, long polreal) |
{ |
{ |
GEN Omega,q,qd,pc,pdc,alpha,beta,gamma,RU,aux,U,W,mygpi,prim,prim2; |
GEN Omega,q,qd,pc,pdc,alpha,beta,gamma,RU,aux,U,W,mygpi,prim,prim2,*gptr[3]; |
long limite,n=degpol(p),i,j,K,ltop; |
long n=degpol(p),i,j,K; |
|
gpmem_t ltop; |
|
|
mygpi=mppi(bitprec/BITS_IN_LONG+3); |
mygpi=mppi(bitprec/BITS_IN_LONG+3); |
aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */ |
aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */ |
prim = exp_i(aux); |
prim = exp_Ir(aux); |
aux = gmulbyi(aux); |
aux = gmulbyi(aux); |
prim2 = myrealun(bitprec); |
prim2 = myrealun(bitprec); |
RU=cgetg(n+2,t_VEC); RU++; |
RU=cgetg(n+2,t_VEC); RU++; |
Line 1051 dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H |
|
Line 1021 dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H |
|
|
|
if (polreal) K=K/2+1; |
if (polreal) K=K/2+1; |
|
|
ltop=avma; limite = stack_lim(ltop,1); |
ltop=avma; |
W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC); |
W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC); |
for (i=1; i<=k; i++) W[i]=U[i]=zero; |
for (i=1; i<=k; i++) W[i]=U[i]=zero; |
|
|
Line 1102 dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H |
|
Line 1072 dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H |
|
} |
} |
} |
} |
prim2=gmul(prim2,prim); |
prim2=gmul(prim2,prim); |
if (low_stack(limite, stack_lim(ltop,1))) |
gptr[0]=&W; gptr[1]=&U; gptr[2]=&prim2; |
{ |
gerepilemany(ltop,gptr,3); |
GEN *gptr[3]; |
|
if(DEBUGMEM>1) err(warnmem,"dft"); |
|
gptr[0]=&W; gptr[1]=&U; gptr[2]=&prim2; |
|
gerepilemany(ltop,gptr,3); |
|
} |
|
} |
} |
|
|
for (i=1; i<=k; i++) |
for (i=1; i<=k; i++) |
|
|
refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec) |
refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec) |
{ |
{ |
GEN H=HH,D,aux; |
GEN H=HH,D,aux; |
ulong ltop=avma, limite=stack_lim(ltop,1); |
gpmem_t ltop=avma, limite=stack_lim(ltop, 1); |
long error=0,i,bitprec1,bitprec2; |
long error=0,i,bitprec1,bitprec2; |
|
|
D=gsub(gun,gres(gmul(HH,G),F)); error=gexpo(D); |
D=gsub(gun,gres(gmul(HH,G),F)); error=gexpo(D); |
Line 1139 refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif |
|
Line 1104 refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif |
|
{ |
{ |
if (low_stack(limite, stack_lim(ltop,1))) |
if (low_stack(limite, stack_lim(ltop,1))) |
{ |
{ |
GEN *gptr[2]; |
GEN *gptr[2]; gptr[0]=&D; gptr[1]=&H; |
if(DEBUGMEM>1) err(warnmem,"refine_H"); |
if(DEBUGMEM>1) err(warnmem,"refine_H"); |
gptr[0]=&D; gptr[1]=&H; gerepilemany(ltop,gptr,2); |
gerepilemany(ltop,gptr,2); |
} |
} |
bitprec1=-error+shiftbitprec; |
bitprec1=-error+shiftbitprec; |
aux=gmul(mygprec(H,bitprec1),mygprec(D,bitprec1)); |
aux=gmul(mygprec(H,bitprec1),mygprec(D,bitprec1)); |
Line 1154 refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif |
|
Line 1119 refine_H(GEN F, GEN G, GEN HH, long bitprec, long shif |
|
D=gsub(gun,gres(gmul(H,G),F)); |
D=gsub(gun,gres(gmul(H,G),F)); |
error=gexpo(D); if (error<-bitprec1) error=-bitprec1; |
error=gexpo(D); if (error<-bitprec1) error=-bitprec1; |
} |
} |
if (error<=-bitprec/2) return gerepilecopy(ltop,H); |
if (error>-bitprec/2) return NULL; /* FAIL */ |
avma=ltop; return gzero; /* procedure failed */ |
return gerepilecopy(ltop,H); |
} |
} |
|
|
/* return 0 if fails, 1 else */ |
/* return 0 if fails, 1 else */ |
|
|
refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, double gamma) |
refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, double gamma) |
{ |
{ |
GEN pp,FF,GG,r,HH,f0; |
GEN pp,FF,GG,r,HH,f0; |
long error,i,bitprec1=0,bitprec2,ltop=avma,shiftbitprec; |
long error, i, bitprec1=0, bitprec2, shiftbitprec; |
long shiftbitprec2,n=degpol(p),enh,normF,normG,limite=stack_lim(ltop,1); |
long shiftbitprec2, n=degpol(p), enh, normF, normG; |
|
gpmem_t ltop=avma, limite=stack_lim(ltop, 1); |
|
|
FF=*F; HH=H; |
FF=*F; HH=H; |
GG=poldivres(p,*F,&r); |
GG=poldivres(p,*F,&r); |
Line 1184 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d |
|
Line 1150 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d |
|
} |
} |
if (low_stack(limite, stack_lim(ltop,1))) |
if (low_stack(limite, stack_lim(ltop,1))) |
{ |
{ |
GEN *gptr[4]; |
GEN *gptr[4]; gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH; |
if(DEBUGMEM>1) err(warnmem,"refine_F"); |
if(DEBUGMEM>1) err(warnmem,"refine_F"); |
gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH; |
|
gerepilemany(ltop,gptr,4); |
gerepilemany(ltop,gptr,4); |
} |
} |
|
|
bitprec1=-error+shiftbitprec2; |
bitprec1=-error+shiftbitprec2; |
HH=refine_H(mygprec(FF,bitprec1),mygprec(GG,bitprec1), |
HH=refine_H(mygprec(FF,bitprec1),mygprec(GG,bitprec1), |
mygprec(HH,bitprec1),1-error,shiftbitprec2); |
mygprec(HH,bitprec1),1-error,shiftbitprec2); |
if (HH==gzero) return 0; /* procedure failed */ |
if (!HH) return 0; /* procedure failed */ |
|
|
bitprec1=-error+shiftbitprec; |
bitprec1=-error+shiftbitprec; |
r=gmul(mygprec(HH,bitprec1),mygprec(r,bitprec1)); |
r=gmul(mygprec(HH,bitprec1),mygprec(r,bitprec1)); |
Line 1210 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d |
|
Line 1175 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, d |
|
GG=poldivres(pp,mygprec(FF,bitprec1),&r); |
GG=poldivres(pp,mygprec(FF,bitprec1),&r); |
error=gexpo(r); if (error<-bitprec1) error=-bitprec1; |
error=gexpo(r); if (error<-bitprec1) error=-bitprec1; |
} |
} |
if (error<=-bitprec) |
if (error>-bitprec) return 0; /* FAIL */ |
{ |
*F=FF; *G=GG; return 1; |
*F=FF; *G=GG; |
|
return 1; /* procedure succeeds */ |
|
} |
|
return 0; /* procedure failed */ |
|
} |
} |
|
|
/* returns F and G from the unit circle U such that |p-FG|<2^(-bitprec) |cd|, |
/* returns F and G from the unit circle U such that |p-FG|<2^(-bitprec) |cd|, |
Line 1225 split_fromU(GEN p, long k, double delta, long bitprec, |
|
Line 1186 split_fromU(GEN p, long k, double delta, long bitprec, |
|
GEN *F, GEN *G, double param, double param2) |
GEN *F, GEN *G, double param, double param2) |
{ |
{ |
GEN pp,FF,GG,H; |
GEN pp,FF,GG,H; |
long n=degpol(p),NN,bitprec2, |
long n=degpol(p),NN,bitprec2; |
ltop=avma,polreal=isreal(p); |
int polreal=isreal(p); |
|
gpmem_t ltop; |
double mu,gamma; |
double mu,gamma; |
|
|
pp=gdiv(p,(GEN)p[2+n]); |
pp=gdiv(p,(GEN)p[2+n]); |
Line 1244 split_fromU(GEN p, long k, double delta, long bitprec, |
|
Line 1206 split_fromU(GEN p, long k, double delta, long bitprec, |
|
bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8; |
bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8; |
dft(pp,k,NN,bitprec2,FF,H,polreal); |
dft(pp,k,NN,bitprec2,FF,H,polreal); |
if (refine_F(pp,&FF,&GG,H,bitprec,gamma)) break; |
if (refine_F(pp,&FF,&GG,H,bitprec,gamma)) break; |
NN=(NN<<1); avma=ltop; |
NN <<= 1; avma=ltop; |
} |
} |
*G=gmul(GG,(GEN)p[2+n]); *F=FF; |
*G=gmul(GG,(GEN)p[2+n]); *F=FF; |
} |
} |
Line 1301 extern GEN addshiftpol(GEN x, GEN y, long d); |
|
Line 1263 extern GEN addshiftpol(GEN x, GEN y, long d); |
|
static GEN |
static GEN |
shiftpol(GEN p, GEN b) |
shiftpol(GEN p, GEN b) |
{ |
{ |
long av = avma,i, limit = stack_lim(av,1); |
long i; |
|
gpmem_t av = avma, limit = stack_lim(av, 1); |
GEN q = gzero; |
GEN q = gzero; |
|
|
if (gcmp0(b)) return p; |
if (gcmp0(b)) return p; |
Line 1326 conformal_pol(GEN p, GEN a, long bitprec) |
|
Line 1289 conformal_pol(GEN p, GEN a, long bitprec) |
|
{ |
{ |
GEN r,pui,num,aux, unr = myrealun(bitprec); |
GEN r,pui,num,aux, unr = myrealun(bitprec); |
long n=degpol(p), i; |
long n=degpol(p), i; |
ulong av, limit; |
gpmem_t av, limit; |
|
|
aux = pui = cgetg(4,t_POL); |
aux = pui = cgetg(4,t_POL); |
pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4); |
pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4); |
Line 1388 compute_radius(GEN* radii, GEN p, long k, double aux, |
|
Line 1351 compute_radius(GEN* radii, GEN p, long k, double aux, |
|
{ |
{ |
if (!signe(radii[i]) || cmprr(radii[i], p1) < 0) |
if (!signe(radii[i]) || cmprr(radii[i], p1) < 0) |
affrr(p1, radii[i]); |
affrr(p1, radii[i]); |
else |
else |
p1 = radii[i]; |
p1 = radii[i]; |
} |
} |
*delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1)); |
*delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1)); |
|
|
conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, long bitprec, |
conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, long bitprec, |
double aux, GEN *F,GEN *G) |
double aux, GEN *F,GEN *G) |
{ |
{ |
long bitprec2,n=degpol(p),decprec,i,ltop = avma, av; |
long bitprec2, n=degpol(p), decprec, i; |
|
gpmem_t ltop = avma, av; |
GEN q,FF,GG,a,R, *gptr[2]; |
GEN q,FF,GG,a,R, *gptr[2]; |
GEN rho,invrho; |
GEN rho,invrho; |
double delta,param,param2; |
double delta,param,param2; |
Line 1432 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
Line 1396 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
if (signe(radii[i])) /* updating array radii */ |
if (signe(radii[i])) /* updating array radii */ |
{ |
{ |
long a = avma; |
gpmem_t a = avma; |
GEN p1 = gsqr(radii[i]); |
GEN p1 = gsqr(radii[i]); |
/* 2(r^2 - 1) / (r^2 - 3(r-1)) */ |
/* 2(r^2 - 1) / (r^2 - 3(r-1)) */ |
p1 = divrr(gmul2n((subrs(p1,1)),1), |
p1 = divrr(gmul2n((subrs(p1,1)),1), |
Line 1443 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
Line 1407 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
|
|
rho = compute_radius(radii, q,k,aux/10., &delta); |
rho = compute_radius(radii, q,k,aux/10., &delta); |
invrho = update_radius(radii, rho, ¶m, ¶m2); |
invrho = update_radius(radii, rho, ¶m, ¶m2); |
|
|
bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.); |
bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.); |
R = mygprec(invrho,bitprec2); |
R = mygprec(invrho,bitprec2); |
q = scalepol(q,R,bitprec2); |
q = scalepol(q,R,bitprec2); |
Line 1470 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
Line 1434 conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, |
|
gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2); |
gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2); |
} |
} |
|
|
|
static GEN |
|
myrealzero(void) |
|
{ |
|
GEN x = cgetr(3); |
|
x[1] = evalexpo(-bit_accuracy(3)); |
|
return x; |
|
} |
|
|
/* split p, this time with no scaling. returns in F and G two polynomials |
/* split p, this time with no scaling. returns in F and G two polynomials |
such that |p-FG|< 2^(-bitprec)|p| */ |
such that |p-FG|< 2^(-bitprec)|p| */ |
static void |
static void |
Line 1480 split_2(GEN p, long bitprec, GEN ctr, double thickness |
|
Line 1452 split_2(GEN p, long bitprec, GEN ctr, double thickness |
|
long n=degpol(p),i,j,k,bitprec2; |
long n=degpol(p),i,j,k,bitprec2; |
GEN q,FF,GG,R; |
GEN q,FF,GG,R; |
GEN *radii = (GEN*) cgetg(n+1, t_VEC); |
GEN *radii = (GEN*) cgetg(n+1, t_VEC); |
for (i=2; i<n; i++) radii[i]=realzero(3); |
for (i=2; i<n; i++) radii[i] = myrealzero(); |
aux = thickness/(double) n/4.; |
aux = thickness/(double) n/4.; |
radii[1] = rmin = min_modulus(p, aux); |
radii[1] = rmin = min_modulus(p, aux); |
radii[n] = rmax = max_modulus(p, aux); |
radii[n] = rmax = max_modulus(p, aux); |
Line 1752 mygprec_absolute(GEN x, long bitprec) |
|
Line 1724 mygprec_absolute(GEN x, long bitprec) |
|
{ |
{ |
case t_REAL: |
case t_REAL: |
e=gexpo(x); |
e=gexpo(x); |
if (e<-bitprec || !signe(x)) { y=dbltor(0.); setexpo(y,-bitprec); } |
if (e<-bitprec || !signe(x)) |
else y=mygprec(x,bitprec+e); |
y=realzero_bit(-bitprec); |
|
else |
|
y=mygprec(x,bitprec+e); |
break; |
break; |
case t_COMPLEX: |
case t_COMPLEX: |
if (gexpo((GEN)x[2])<-bitprec) |
if (gexpo((GEN)x[2])<-bitprec) |
Line 1781 a_posteriori_errors(GEN p, GEN roots_pol, long err) |
|
Line 1755 a_posteriori_errors(GEN p, GEN roots_pol, long err) |
|
setexpo(sigma, err + (long)log2((double)n) + 1); |
setexpo(sigma, err + (long)log2((double)n) + 1); |
overn=dbltor(1./n); |
overn=dbltor(1./n); |
shatzle=gdiv(gpui(sigma,overn,0), |
shatzle=gdiv(gpui(sigma,overn,0), |
gsub(gpui(gsub(gun,sigma),overn,0), |
gsub(gpui(gsub(realun(DEFAULTPREC),sigma),overn,0), |
gpui(sigma,overn,0))); |
gpui(sigma,overn,0))); |
shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY; |
shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY; |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
Line 1799 a_posteriori_errors(GEN p, GEN roots_pol, long err) |
|
Line 1773 a_posteriori_errors(GEN p, GEN roots_pol, long err) |
|
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
static GEN |
static GEN |
append_root(GEN roots_pol, GEN a) |
append_clone(GEN r, GEN a) { a = gclone(a); appendL(r, a); return a; } |
{ |
|
long l = lg(roots_pol); |
|
setlg(roots_pol, l+1); return (GEN)(roots_pol[l] = lclone(a)); |
|
} |
|
|
|
/* put roots in placeholder roots_pol so that |P-L_1...L_n|<2^(-bitprec)|P| |
/* put roots in placeholder roots_pol so that |P-L_1...L_n|<2^(-bitprec)|P| |
* and returns prod (x-roots_pol[i]) for i=1..degree(p) */ |
* and returns prod (x-roots_pol[i]) for i=1..degree(p) */ |
static GEN |
static GEN |
split_complete(GEN p, long bitprec, GEN roots_pol) |
split_complete(GEN p, long bitprec, GEN roots_pol) |
{ |
{ |
long n=degpol(p),decprec,ltop; |
long n=degpol(p), decprec; |
|
gpmem_t ltop; |
GEN p1,F,G,a,b,m1,m2,m; |
GEN p1,F,G,a,b,m1,m2,m; |
|
|
if (n==1) |
if (n==1) |
{ |
{ |
a=gneg_i(gdiv((GEN)p[2],(GEN)p[3])); |
a=gneg_i(gdiv((GEN)p[2],(GEN)p[3])); |
append_root(roots_pol,a); return p; |
(void)append_clone(roots_pol,a); return p; |
} |
} |
ltop = avma; |
ltop = avma; |
if (n==2) |
if (n==2) |
Line 1827 split_complete(GEN p, long bitprec, GEN roots_pol) |
|
Line 1798 split_complete(GEN p, long bitprec, GEN roots_pol) |
|
p1 = gmul2n((GEN)p[4],1); |
p1 = gmul2n((GEN)p[4],1); |
a = gneg_i(gdiv(gadd(F,(GEN)p[3]), p1)); |
a = gneg_i(gdiv(gadd(F,(GEN)p[3]), p1)); |
b = gdiv(gsub(F,(GEN)p[3]), p1); |
b = gdiv(gsub(F,(GEN)p[3]), p1); |
a = append_root(roots_pol,a); |
a = append_clone(roots_pol,a); |
b = append_root(roots_pol,b); avma = ltop; |
b = append_clone(roots_pol,b); avma = ltop; |
m=gmul(gsub(polx[varn(p)],mygprec(a,3*bitprec)), |
m=gmul(gsub(polx[varn(p)],mygprec(a,3*bitprec)), |
gsub(polx[varn(p)],mygprec(b,3*bitprec))); |
gsub(polx[varn(p)],mygprec(b,3*bitprec))); |
return gmul(m,(GEN)p[4]); |
return gmul(m,(GEN)p[4]); |
Line 1840 split_complete(GEN p, long bitprec, GEN roots_pol) |
|
Line 1811 split_complete(GEN p, long bitprec, GEN roots_pol) |
|
} |
} |
|
|
/* compute a bound on the maximum modulus of roots of p */ |
/* compute a bound on the maximum modulus of roots of p */ |
static GEN |
GEN |
cauchy_bound(GEN p) |
cauchy_bound(GEN p) |
{ |
{ |
long i,n=degpol(p); |
long i, n = degpol(p), prec = DEFAULTPREC; |
GEN x=gzero,y,lc; |
GEN lc, y, q = gmul(p, realun(prec)), x = gzero; |
|
|
lc=gabs((GEN)p[n+2],DEFAULTPREC); /* leading coefficient */ |
if (n <= 0) err(constpoler,"cauchy_bound"); |
lc=gdiv(dbltor(1.),lc); |
|
|
lc = gabs((GEN)q[n+2],prec); /* leading coefficient */ |
|
lc = ginv(lc); |
for (i=0; i<n; i++) |
for (i=0; i<n; i++) |
{ |
{ |
y=gmul(gabs((GEN) p[i+2],DEFAULTPREC),lc); |
y = (GEN)q[i+2]; if (gcmp0(y)) continue; |
y=gpui(y,dbltor(1./(n-i)),DEFAULTPREC); |
y = gmul(gabs(y,prec), lc); |
if (gcmp(y,x) > 0) x=y; |
y = divrs(mplog(y), n-i); |
|
if (gcmp(y,x) > 0) x = y; |
} |
} |
return x; |
return gexp(x, prec); |
} |
} |
|
|
static GEN |
static GEN |
Line 1925 fix_roots(GEN r, GEN *m, long h, long bitprec) |
|
Line 1899 fix_roots(GEN r, GEN *m, long h, long bitprec) |
|
static GEN |
static GEN |
all_roots(GEN p, long bitprec) |
all_roots(GEN p, long bitprec) |
{ |
{ |
GEN pd,q,roots_pol,m; |
GEN lc,pd,q,roots_pol,m; |
long bitprec0, bitprec2,n=degpol(p),i,e,h; |
long bitprec0, bitprec2,n=degpol(p),i,e,h; |
ulong av; |
gpmem_t av; |
|
|
#if 0 |
pd = poldeflate(p, &h); lc = leading_term(pd); |
pd = poldeflate(p, &h); |
|
#else |
|
pd = p; h = 1; |
|
#endif |
|
e = 2*gexpo(cauchy_bound(pd)); if (e<0) e=0; |
e = 2*gexpo(cauchy_bound(pd)); if (e<0) e=0; |
bitprec0=bitprec + gexpo(pd) - gexpo(leading_term(pd)) + (long)log2(n/h)+1+e; |
bitprec0=bitprec + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e; |
|
bitprec2 = bitprec0; e = 0; |
for (av=avma,i=1;; i++,avma=av) |
for (av=avma,i=1;; i++,avma=av) |
{ |
{ |
roots_pol = cgetg(n+1,t_VEC); setlg(roots_pol,1); |
roots_pol = cget1(n+1,t_VEC); |
bitprec2 = bitprec0 + (1<<i)*n; |
bitprec2 += e + (1<<i)*n; |
q = gmul(myrealun(bitprec2), mygprec(pd,bitprec2)); |
q = gmul(myrealun(bitprec2), mygprec(pd,bitprec2)); |
m = split_complete(q,bitprec2,roots_pol); |
m = split_complete(q,bitprec2,roots_pol); |
roots_pol = fix_roots(roots_pol, &m, h, bitprec2); |
roots_pol = fix_roots(roots_pol, &m, h, bitprec2); |
|
q = mygprec_special(p,bitprec2); lc = leading_term(q); |
|
if (h > 1) m = gmul(m,lc); |
|
|
e = gexpo(gsub(mygprec_special(p,bitprec2), m)) |
e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1; |
- gexpo(leading_term(q)) + (long)log2((double)n) + 1; |
|
if (e<-2*bitprec2) e=-2*bitprec2; /* to avoid e=-pariINFINITY */ |
if (e<-2*bitprec2) e=-2*bitprec2; /* to avoid e=-pariINFINITY */ |
if (e < 0) |
if (e < 0) |
{ |
{ |
e = a_posteriori_errors(q,roots_pol,e); |
e = bitprec + a_posteriori_errors(p,roots_pol,e); |
if (e < -bitprec) return roots_pol; |
if (e < 0) return roots_pol; |
} |
} |
if (DEBUGLEVEL > 7) |
if (DEBUGLEVEL > 7) |
fprintferr("all_roots: restarting, i = %ld, e = %ld\n", i,e); |
fprintferr("all_roots: restarting, i = %ld, e = %ld\n", i,e); |
Line 2079 isrealappr(GEN x, long e) |
|
Line 2051 isrealappr(GEN x, long e) |
|
static int |
static int |
isconj(GEN x, GEN y, long e) |
isconj(GEN x, GEN y, long e) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long i= (gexpo( gsub((GEN)x[1],(GEN)y[1]) ) < e |
long i= (gexpo( gsub((GEN)x[1],(GEN)y[1]) ) < e |
&& gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e); |
&& gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e); |
avma = av; return i; |
avma = av; return i; |
Line 2091 isconj(GEN x, GEN y, long e) |
|
Line 2063 isconj(GEN x, GEN y, long e) |
|
GEN |
GEN |
roots(GEN p, long l) |
roots(GEN p, long l) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long n,i,k,s,t,e; |
long n,i,k,s,t,e; |
GEN c,L,p1,res,rea,com; |
GEN c,L,p1,res,rea,com; |
|
|