OpenXM_contrib/pari-2.2/src/basemath/rootpol.c - diff

Return to rootpol.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/rootpol.c between version 1.1 and 1.2

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,

#include "pari.h"

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 NEWTON_MAX 10

Line 65 quickmulcc(GEN x, GEN y)

Line 63 quickmulcc(GEN x, GEN y)

}

if (ty==t_COMPLEX)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN p1,p2;

z=cgetg(3,t_COMPLEX); av=avma;

Line 97 static GEN

Line 95 static GEN

mysquare(GEN p)

{

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);

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;

long n=degpol(p),n0,n1,i,var,nn0;

ulong ltop;

gpmem_t ltop;

if (n<=KARASQUARE_LIMIT) return mysquare(p);

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;

long n=degpol(p),n0,n3,i,j,var;

ulong ltop = avma;

gpmem_t ltop = avma;

if (n<=COOK_SQUARE_LIMIT) return karasquare(p);

Line 347 log2ir(GEN x)

Line 346 log2ir(GEN x)

return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG;

}

static double

double

mylog2(GEN z)

{

double x,y;

Line 360 mylog2(GEN z)

Line 359 mylog2(GEN z)

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

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--)

{

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 */

static long

polygone_newton(GEN p, long k)

newton_polygon(GEN p, long k)

{

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++)

{

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 */

Line 444 myshiftic(GEN z, long e)

static GEN

myrealun(long bitprec)

{

GEN x;

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

Line 547 homothetie(GEN p, GEN R, long bitprec)

Line 545 homothetie(GEN p, GEN R, long bitprec)

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

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,

puts in k an upper bound of the number of roots near the largest root

at a distance eps */

static double

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;

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));

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++)

{

s[i+1]=lmulsg(i,(GEN)a[i+1]);

s = gmulsg(i,(GEN)a[i]);

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

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];

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

*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 */

static GEN

GEN

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 */

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);

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;;)

{

{ /* 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));

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));

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);

}

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*/

Line 668 static GEN

Line 660 static GEN

modulus(GEN p, long k, double tau)

{

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;

tau2=tau/6; nn=n;

tau2 = tau/6;

bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));

decprec=(long) ((double) bitprec * L2SL10)+1;

gunr=myrealun(bitprec);

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);

r=(double) e;

r = (double) e;

imax=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1;

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);

q = gerepileupto(av, graeffe(q));

e=polygone_newton(q,kk);

e = newton_polygon(q,kk);

r += e / exp2((double)i);

q=gmul(gunr,q);

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)));

}

avma=ltop; return mpexp(dbltor(-r * LOG2));

Line 715 static GEN

Line 701 static GEN

pre_modulus(GEN p, long k, double tau, GEN rmin, GEN rmax)

{

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;

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

min_modulus(GEN p, double tau)

{

long av=avma;

gpmem_t av=avma;

GEN r;

if (isexactzero((GEN)p[2])) return realzero(3);

Line 764 static long

Line 751 static long

dual_modulus(GEN p, GEN R, double tau, long l)

{

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;

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.));

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) */

if (ll<0) ll=0;

nn=nnn-valuat;

nn = nnn-valuat; delta_k += valuat;

if (nn==0) return delta_k;

set_karasquare_limit(bitprec);

q = gerepileupto(ltop, graeffe(q));

tau2=tau2*7./4.;

}

k=-1; logmax=- (double) pariINFINITY;

k = -1; logmax = - (double)pariINFINITY;

for (i=0; i<=degpol(q); i++)

{

aux=mylog2((GEN)q[2+i]);

Line 831 gmulbyi(GEN z)

Line 813 gmulbyi(GEN z)

static void

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;

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];

}

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) */

static GEN

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 == 4) return gi;

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 */

Line 947 initRUgen(long N, long bitprec)

Line 915 initRUgen(long N, long bitprec)

}

/* returns 1 if p has only real coefficients, 0 else */

static long

static int

isreal(GEN p)

{

long n=degpol(p),i=0;

Line 958 isreal(GEN p)

Line 926 isreal(GEN p)

static void

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;

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;

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;

mygpi=mppi(bitprec/BITS_IN_LONG+3);

aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */

prim = exp_i(aux);

prim = exp_Ir(aux);

aux = gmulbyi(aux);

RU = myrealun(bitprec);

Line 1026 parameters(GEN p, double *mu, double *gamma,

Line 995 parameters(GEN p, double *mu, double *gamma,

static void

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);

aux = gdivgs(mygpi,NN/2); /* 2 Pi/NN */

prim = exp_i(aux);

prim = exp_Ir(aux);

aux = gmulbyi(aux);

prim2 = myrealun(bitprec);

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;

ltop=avma; limite = stack_lim(ltop,1);

ltop=avma;

W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC);

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);

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++)

Line 1129 static GEN

Line 1094 static GEN

refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec)

{

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;

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)))

{

GEN *gptr[2];

GEN *gptr[2]; gptr[0]=&D; gptr[1]=&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;

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));

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 */

Line 1163 static long

Line 1128 static long

refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, double gamma)

{

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;

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)))

{

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");

gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH;

gerepilemany(ltop,gptr,4);

}

bitprec1=-error+shiftbitprec2;

HH=refine_H(mygprec(FF,bitprec1),mygprec(GG,bitprec1),

mygprec(HH,bitprec1),1-error,shiftbitprec2);

if (HH==gzero) return 0; /* procedure failed */

if (!HH) return 0; /* procedure failed */

bitprec1=-error+shiftbitprec;

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);

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|,

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 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;

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;

dft(pp,k,NN,bitprec2,FF,H,polreal);

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;

}

Line 1301 extern GEN addshiftpol(GEN x, GEN y, long d);

Line 1263 extern GEN addshiftpol(GEN x, GEN y, long d);

static GEN

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;

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);

long n=degpol(p), i;

ulong av, limit;

gpmem_t av, limit;

aux = pui = cgetg(4,t_POL);

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)

affrr(p1, radii[i]);

else

p1 = radii[i];

}

*delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1));

Line 1417 static void

Line 1380 static void

conformal_mapping(GEN *radii, GEN ctr, GEN p, long k, long bitprec,

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 rho,invrho;

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++)

if (signe(radii[i])) /* updating array radii */

{

long a = avma;

gpmem_t a = avma;

GEN p1 = gsqr(radii[i]);

/* 2(r^2 - 1) / (r^2 - 3(r-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);

invrho = update_radius(radii, rho, &param, &param2);

bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.);

R = mygprec(invrho,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);

}

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

such that |p-FG|< 2^(-bitprec)|p| */

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;

GEN q,FF,GG,R;

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.;

radii[1] = rmin = min_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:

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;

case t_COMPLEX:

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);

overn=dbltor(1./n);

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)));

shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY;

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

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|

* and returns prod (x-roots_pol[i]) for i=1..degree(p) */

static GEN

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;

if (n==1)

{

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;

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);

a = gneg_i(gdiv(gadd(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)),

gsub(polx[varn(p)],mygprec(b,3*bitprec)));

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 */

static GEN

GEN

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++)

{

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

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

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;

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;

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)

{

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));

m = split_complete(q,bitprec2,roots_pol);

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 < 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)

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

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

&& gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e);

avma = av; return i;

Line 2091 isconj(GEN x, GEN y, long e)

Line 2063 isconj(GEN x, GEN y, long e)

GEN

roots(GEN p, long l)

{

ulong av = avma;

gpmem_t av = avma;

long n,i,k,s,t,e;

GEN c,L,p1,res,rea,com;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>