Annotation of OpenXM_contrib/pari/src/basemath/rootpol.c, Revision 1.1
1.1 ! maekawa 1: /*******************************************************************/
! 2: /* */
! 3: /* ROOTS OF COMPLEX POLYNOMIALS */
! 4: /* (written by Xavier Gourdon) */
! 5: /* */
! 6: /*******************************************************************/
! 7: /* $Id: rootpol.c,v 1.2 1999/09/23 17:50:56 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: GEN polrecip_i(GEN x);
! 11: #define pariINFINITY 100000
! 12: #define NEWTON_MAX 10
! 13:
! 14: static GEN gunr, globalu;
! 15: static long KARASQUARE_LIMIT, COOK_SQUARE_LIMIT, Lmax, step4;
! 16: static double *radius;
! 17:
! 18: /********************************************************************/
! 19: /** **/
! 20: /** ARITHMETIQUE RAPIDE **/
! 21: /** **/
! 22: /********************************************************************/
! 23:
! 24: /* fast product of x,y which must be integer or complex of integer */
! 25: static GEN
! 26: quickmulcc(GEN x, GEN y)
! 27: {
! 28: long tx=typ(x),ty=typ(y);
! 29: GEN z;
! 30:
! 31: if (tx==t_INT)
! 32: {
! 33: if (ty==t_INT) return mulii(x,y);
! 34: if (ty==t_COMPLEX)
! 35: {
! 36: z=cgetg(3,t_COMPLEX);
! 37: z[1]=(long) mulii(x,(GEN) y[1]);
! 38: z[2]=(long) mulii(x,(GEN) y[2]);
! 39: return z;
! 40: }
! 41: }
! 42:
! 43: if (tx==t_COMPLEX)
! 44: {
! 45: if (ty==t_INT)
! 46: {
! 47: z=cgetg(3,t_COMPLEX);
! 48: z[1]=(long) mulii((GEN)x[1],y);
! 49: z[2]=(long) mulii((GEN)x[2],y);
! 50: return z;
! 51: }
! 52: if (ty==t_COMPLEX)
! 53: {
! 54: long av,tetpil;
! 55: GEN p1,p2;
! 56:
! 57: z=cgetg(3,t_COMPLEX); av=avma;
! 58: p1=mulii((GEN)x[1],(GEN)y[1]); p2=mulii((GEN)x[2],(GEN)y[2]);
! 59: x=addii((GEN)x[1],(GEN)x[2]); y=addii((GEN)y[1],(GEN)y[2]);
! 60: y=mulii(x,y); x=addii(p1,p2);
! 61: tetpil=avma; z[1]=lsubii(p1,p2); z[2]=lsubii(y,x);
! 62: gerepilemanyvec(av,tetpil,z+1,2);
! 63: return z;
! 64: }
! 65: }
! 66: err(talker,"bug in quickmulcc");
! 67: return NULL; /* not reached */
! 68: }
! 69:
! 70: static void
! 71: set_karasquare_limit(long bitprec)
! 72: {
! 73: if (bitprec<600) { KARASQUARE_LIMIT=8; COOK_SQUARE_LIMIT=400; return; }
! 74: if (bitprec<2000) { KARASQUARE_LIMIT=4; COOK_SQUARE_LIMIT=200; return; }
! 75: if (bitprec<3000) { KARASQUARE_LIMIT=4; COOK_SQUARE_LIMIT=125; return; }
! 76: if (bitprec<5000) { KARASQUARE_LIMIT=2; COOK_SQUARE_LIMIT=75; return; }
! 77: KARASQUARE_LIMIT=1; COOK_SQUARE_LIMIT=50;
! 78: }
! 79:
! 80: /* the pari library does not have specific procedure for the square of
! 81: polynomials. This one is twice faster than gmul */
! 82: static GEN
! 83: mysquare(GEN p)
! 84: {
! 85: GEN s,aux1,aux2;
! 86: long i,j,n=lgef(p)-3,nn,ltop,lbot;
! 87:
! 88: if (n==-1) return gcopy(p);
! 89: nn=n<<1; s=cgetg(nn+3,t_POL);
! 90: s[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(nn+3);
! 91: for (i=0; i<=n; i++)
! 92: {
! 93: aux1=gzero; ltop=avma;
! 94: for (j=0; j<(i+1)>>1; j++)
! 95: {
! 96: aux2=quickmulcc((GEN) p[j+2],(GEN)p[i-j+2]);
! 97: aux1=gadd(aux1,aux2);
! 98: }
! 99: if (i%2==1) { lbot=avma; s[i+2]=lpile(ltop,lbot,gshift(aux1,1)); }
! 100: else
! 101: {
! 102: aux1=gshift(aux1,1);
! 103: aux2=quickmulcc((GEN)p[2+(i>>1)],(GEN)p[2+(i>>1)]);
! 104: lbot=avma; s[i+2]=lpile(ltop,lbot,gadd(aux1,aux2));
! 105: }
! 106: }
! 107: for (i=n+1; i<=nn; i++)
! 108: {
! 109: aux1=gzero; ltop=avma;
! 110: for (j=i-n; j<(i+1)>>1; j++)
! 111: {
! 112: aux2=quickmulcc((GEN)p[j+2],(GEN)p[i-j+2]);
! 113: aux1=gadd(aux1,aux2);
! 114: }
! 115: if (i%2==1) { lbot=avma; s[i+2]=lpile(ltop,lbot,gshift(aux1,1)); }
! 116: else
! 117: {
! 118: aux1=gshift(aux1,1);
! 119: aux2=quickmulcc((GEN)p[2+(i>>1)],(GEN)p[2+(i>>1)]);
! 120: lbot=avma; s[i+2]=lpile(ltop,lbot,gadd(aux1,aux2));
! 121: }
! 122: }
! 123: return s;
! 124: }
! 125:
! 126: static GEN
! 127: karasquare(GEN p)
! 128: {
! 129: GEN p1,s0,s1,s2,aux;
! 130: long n=lgef(p)-3,n0,n1,i,var,nn0,ltop,lbot;
! 131:
! 132: if (n<=KARASQUARE_LIMIT) return mysquare(p);
! 133: ltop=avma;
! 134: var=evalsigne(1)+evalvarn(varn(p)); n0=n>>1; n1=n-n0-1;
! 135: setlgef(p,n0+3); /* hack to have the first half of p */
! 136: s0=karasquare(p);
! 137: p1=cgetg(n1+3,t_POL); p1[1]=var+evallgef(n1+3);
! 138: for (i=2; i<=n1+2; i++) p1[i]=p[1+i+n0];
! 139: s2=karasquare(p1);
! 140: s1=karasquare(gadd(p,p1));
! 141: s1=gsub(s1,gadd(s0,s2));
! 142: nn0=n0<<1;
! 143: aux=cgetg((n<<1)+3,t_POL); aux[1]=var+evallgef(2*n+3);
! 144: var=lgef(s0);
! 145: for (i=2; i<var; i++) aux[i]=s0[i];
! 146: for ( ; i<=nn0+2; i++) aux[i]=zero;
! 147: var=lgef(s2);
! 148: for (i=2; i<var; i++) aux[2+i+nn0]=s2[i];
! 149: for (i=var-2; i<=(n1<<1); i++) aux[4+i+nn0]=zero;
! 150: aux[3+nn0]=zero;
! 151: for (i=3; i<=lgef(s1); i++)
! 152: aux[i+n0]=ladd((GEN) aux[i+n0],(GEN) s1[i-1]);
! 153: setlgef(p,n+3); /* recover all the polynomial p */
! 154: lbot=avma; return gerepile(ltop,lbot,gcopy(aux));
! 155: }
! 156:
! 157: static GEN
! 158: cook_square(GEN p)
! 159: {
! 160: GEN p0,p1,p2,p3,q,aux0,aux1,r,aux,plus,moins;
! 161: long n=lgef(p)-3,n0,n3,i,j,ltop=avma,lbot,var;
! 162:
! 163: if (n<=COOK_SQUARE_LIMIT) return karasquare(p);
! 164:
! 165: n0=(n+1)/4; n3=n+1-3*n0;
! 166: p0=cgetg(n0+2,t_POL); p1=cgetg(n0+2,t_POL); p2=cgetg(n0+2,t_POL);
! 167: p3=cgetg(n3+2,t_POL);
! 168: var=evalsigne(1)|evalvarn(varn(p));
! 169: p0[1]=p1[1]=p2[1]=var|evallgef(n0+2); p3[1]=var|evallgef(n3+2);
! 170:
! 171: for (i=0; i<n0; i++)
! 172: {
! 173: p0[i+2]=p[i+2]; p1[i+2]=p[i+n0+2]; p2[i+2]=p[i+2*n0+2];
! 174: }
! 175: for (i=0; i<n3; i++) p3[i+2]=p[i+3*n0+2];
! 176:
! 177: q=cgetg(8,t_VEC); q=q+4;
! 178:
! 179: q[0]=(long) p0;
! 180: aux0=gadd(p0,p2); aux1=gadd(p1,p3);
! 181: q[-1]=lsub(aux0,aux1); q[1]=ladd(aux0,aux1);
! 182: aux0=gadd(p0,gmulgs(p2,4)); aux1=gmulgs(gadd(p1,gmulgs(p3,4)),2);
! 183: q[-2]=lsub(aux0,aux1); q[2]=ladd(aux0,aux1);
! 184: aux0=gadd(p0,gmulgs(p2,9)); aux1=gmulgs(gadd(p1,gmulgs(p3,9)),3);
! 185: q[-3]=lsub(aux0,aux1); q[3]=ladd(aux0,aux1);
! 186: for (i=-3; i<=3; i++) q[i]=(long) cook_square((GEN)q[i]);
! 187: r=new_chunk(7);
! 188: plus=cgetg(4,t_VEC); moins=cgetg(4,t_VEC);
! 189: for (i=1; i<=3; i++)
! 190: {
! 191: plus[i]=ladd((GEN)q[-i],(GEN)q[i]);
! 192: moins[i]=lsub((GEN)q[-i],(GEN)q[i]);
! 193: }
! 194: r[0]=q[0];
! 195: r[1]=ldivgs(
! 196: gsub(
! 197: gsub(gmulgs((GEN)moins[2],9),(GEN)moins[3]),
! 198: gmulgs((GEN)moins[1],45)),
! 199: 60);
! 200: r[2]=ldivgs(
! 201: gadd(
! 202: gadd(gmulgs((GEN)plus[1],270),gmulgs((GEN)q[0],-490)),
! 203: gadd(gmulgs((GEN)plus[2],-27),gmulgs((GEN)plus[3],2))),
! 204: 360);
! 205: r[3]=ldivgs(
! 206: gadd(
! 207: gadd(gmulgs((GEN)moins[1],13),gmulgs((GEN)moins[2],-8)),
! 208: (GEN)moins[3]),
! 209: 48);
! 210: r[4]=ldivgs(
! 211: gadd(
! 212: gadd(gmulgs((GEN)q[0],56),gmulgs((GEN)plus[1],-39)),
! 213: gsub(gmulgs((GEN)plus[2],12),(GEN)plus[3])),
! 214: 144);
! 215: r[5]=ldivgs(
! 216: gsub(
! 217: gadd(gmulgs((GEN)moins[1],-5),gmulgs((GEN)moins[2],4)),
! 218: (GEN)moins[3]),
! 219: 240);
! 220: r[6]=ldivgs(
! 221: gadd(
! 222: gadd(gmulgs((GEN)q[0],-20),gmulgs((GEN)plus[1],15)),
! 223: gadd(gmulgs((GEN)plus[2],-6),(GEN)plus[3])),
! 224: 720);
! 225: q=cgetg(2*n+3,t_POL); q[1]=var|evallgef(2*n+3);
! 226: for (i=0; i<=2*n; i++) q[i+2]=zero;
! 227: for (i=0; i<=6; i++)
! 228: {
! 229: aux=(GEN) r[i];
! 230: for (j=0; j<=lgef(aux)-3; j++)
! 231: q[n0*i+j+2]=ladd((GEN)q[n0*i+j+2],(GEN)aux[j+2]);
! 232: }
! 233: lbot=avma; return gerepile(ltop,lbot,gcopy(q));
! 234: }
! 235:
! 236: static GEN
! 237: graeffe(GEN p)
! 238: {
! 239: GEN p0,p1,s0,s1,ss1;
! 240: long n=lgef(p)-3,n0,n1,i,auxi,ns1;
! 241:
! 242: if (n==0) return gcopy(p);
! 243: n0=n>>1; n1=(n-1)>>1;
! 244: auxi=evalsigne(1)|evalvarn(varn(p));
! 245: p0=cgetg(n0+3,t_POL); p0[1]=auxi|evallgef(n0+3);
! 246: p1=cgetg(n1+3,t_POL); p1[1]=auxi|evallgef(n1+3);
! 247: for (i=0; i<=n0; i++) p0[i+2]=p[2+(i<<1)];
! 248: for (i=0; i<=n1; i++) p1[i+2]=p[3+(i<<1)];
! 249:
! 250: s0=cook_square(p0);
! 251: s1=cook_square(p1); ns1 = lgef(s1)-3;
! 252: ss1 = cgetg(ns1+4, t_POL);
! 253: ss1[1] = auxi | evallgef(ns1+4);
! 254: ss1[2]=zero;
! 255: for (i=0; i<=ns1; i++) ss1[3+i]=lneg((GEN)s1[2+i]);
! 256: /* now ss1 contains -x * s1 */
! 257: return gadd(s0,ss1);
! 258: }
! 259:
! 260: /********************************************************************/
! 261: /** **/
! 262: /** FACTORISATION SQUAREFREE AVEC LE GCD MODULAIRE **/
! 263: /** **/
! 264: /********************************************************************/
! 265:
! 266: /* return a n x 2 matrix:
! 267: * col 1 contains the i's such that A_i non constant
! 268: * col 2 the A_i's, s.t. pol = A_i1^i1.A_i2^i2...A_in^in.
! 269: * if pol is constant return [;]
! 270: */
! 271: GEN
! 272: square_free_factorization(GEN pol)
! 273: {
! 274: long deg,i,j,m;
! 275: GEN p1,x,t1,v1,t,v,A;
! 276:
! 277: if (typ(pol)!=t_POL) err(typeer,"square_free_factorization");
! 278: deg=lgef(pol)-3; if (deg<1) return cgetg(1,t_MAT);
! 279: p1 = content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);
! 280:
! 281: x=cgetg(3,t_MAT);
! 282: if (deg > 1)
! 283: {
! 284: t1 = modulargcd(pol,derivpol(pol));
! 285: if (isscalar(t1)) deg = 1;
! 286: }
! 287: if (deg==1)
! 288: {
! 289: x[1]=lgetg(2,t_COL); p1=(GEN)x[1]; p1[1]=un;
! 290: x[2]=lgetg(2,t_COL); p1=(GEN)x[2]; p1[1]=(long)pol; return x;
! 291: }
! 292: A=new_chunk(deg+1); v1=gdivexact(pol,t1); v=v1; i=0;
! 293: while (lgef(v)>3)
! 294: {
! 295: v=modulargcd(t1,v1); i++;
! 296: A[i]=(long)gdivexact(v1,v);
! 297: t=gdivexact(t1,v); v1=v; t1=t;
! 298: }
! 299: m=1; x[1]=lgetg(deg+1,t_COL); x[2]=lgetg(deg+1,t_COL);
! 300: for (j=1; j<=i; j++)
! 301: if (isnonscalar(A[j]))
! 302: {
! 303: p1=(GEN)x[1]; p1[m] = lstoi(j);
! 304: p1=(GEN)x[2]; p1[m] = A[j];
! 305: m++;
! 306: }
! 307: setlg(x[1],m); setlg(x[2],m); return x;
! 308: }
! 309:
! 310: /********************************************************************/
! 311: /** **/
! 312: /** CALCUL DU MODULE DES RACINES **/
! 313: /** **/
! 314: /********************************************************************/
! 315:
! 316: static double
! 317: log2ir(GEN x)
! 318: {
! 319: double l;
! 320:
! 321: if (signe(x)==0) return (double) -pariINFINITY;
! 322: if (typ(x)==t_INT)
! 323: {
! 324: if (lgef(x)==3) return (double) log2( (double)(ulong) x[2]);
! 325: l=(double)(ulong) x[2]+
! 326: (double)(ulong) x[3] / exp2((double) BITS_IN_LONG);
! 327: return log2(l)+ (double) BITS_IN_LONG * (lgef(x)-3.);
! 328: }
! 329: /* else x is real */
! 330: return 1.+ (double) expo(x)+log2( (double)(ulong) x[2]) - (double) BITS_IN_LONG;
! 331: }
! 332:
! 333: static double
! 334: mylog2(GEN z)
! 335: {
! 336: double x,y;
! 337:
! 338: if (typ(z)!=t_COMPLEX) return log2ir(z);
! 339:
! 340: x=log2ir((GEN) z[1]); y=log2ir((GEN) z[2]);
! 341: if (fabs(x-y)>10) return (x>y)? x: y;
! 342: return x+0.5*log2( 1 + exp2(2*(y-x)));
! 343: }
! 344:
! 345: static long
! 346: findpower(GEN p)
! 347: {
! 348: double x, logbinomial,pente,pentemax=-pariINFINITY;
! 349: long n=lgef(p)-3,i;
! 350:
! 351: logbinomial=mylog2((GEN) p[n+2]);
! 352: for (i=n-1; i>=0; i--)
! 353: {
! 354: logbinomial=logbinomial+log2((double) (i+1) / (double) (n-i));
! 355: x=mylog2((GEN) p[2+i])-logbinomial;
! 356: if (x>-pariINFINITY)
! 357: {
! 358: pente=x/ (double) (n-i);
! 359: if (pente>pentemax) pentemax=pente;
! 360: }
! 361: }
! 362: return (long) -floor(pentemax);
! 363: }
! 364:
! 365: /* returns the exponent for the procedure modulus, from the newton diagram */
! 366: static long
! 367: polygone_newton(GEN p, long k)
! 368: {
! 369: double *logcoef,pente;
! 370: long n=lgef(p)-3,i,j,h,l,*sommet,pentelong;
! 371:
! 372: logcoef=(double*) gpmalloc((n+1)*sizeof(double));
! 373: sommet=(long*) gpmalloc((n+1)*sizeof(long));
! 374:
! 375: /* sommet[i]=1 si i est un sommet, =0 sinon */
! 376: for (i=0; i<=n; i++) { logcoef[i]=mylog2((GEN)p[2+i]); sommet[i]=0; }
! 377: sommet[0]=1; i=0;
! 378: while (i<n)
! 379: {
! 380: pente=logcoef[i+1]-logcoef[i];
! 381: h=i+1;
! 382: for (j=i+1; j<=n; j++)
! 383: {
! 384: if (pente<(logcoef[j]-logcoef[i])/(double)(j-i))
! 385: {
! 386: h=j;
! 387: pente=(logcoef[j]-logcoef[i])/(double)(j-i);
! 388: }
! 389: }
! 390: i=h;
! 391: sommet[h]=1;
! 392: }
! 393: h=k; while (!sommet[h]) h++;
! 394: l=k-1; while (!sommet[l]) l--;
! 395: pentelong=(long) floor((logcoef[h]-logcoef[l])/(double)(h-l)+0.5);
! 396: free(logcoef); free(sommet); return pentelong;
! 397: }
! 398:
! 399: /* change z into z*2^e, where z is real or complex of real */
! 400: static void
! 401: myshiftrc(GEN z, long e)
! 402: {
! 403: if (typ(z)==t_COMPLEX)
! 404: {
! 405: if (signe(z[1])!=0) setexpo(z[1], expo(z[1])+e);
! 406: if (signe(z[2])!=0) setexpo(z[2], expo(z[2])+e);
! 407: }
! 408: else
! 409: if (signe(z)!=0) setexpo(z,expo(z)+e);
! 410: }
! 411:
! 412: /* return z*2^e, where z is integer or complex of integer (destroy z) */
! 413: static GEN
! 414: myshiftic(GEN z, long e)
! 415: {
! 416: if (typ(z)==t_COMPLEX)
! 417: {
! 418: z[1]=signe(z[1])? lmpshift((GEN) z[1],e): zero;
! 419: z[2]=lmpshift((GEN) z[2],e);
! 420: return z;
! 421: }
! 422: return signe(z)? mpshift(z,e): gzero;
! 423: }
! 424:
! 425: static GEN
! 426: mygprecrc(GEN x, long bitprec, long e)
! 427: {
! 428: long tx=typ(x);
! 429: GEN y;
! 430:
! 431: if (bitprec<0) bitprec=0; /* should rarely happen */
! 432: switch(tx)
! 433: {
! 434: case t_REAL:
! 435: y=cgetr(bitprec/BITS_IN_LONG+3); affrr(x,y);
! 436: if (!signe(x)) setexpo(y,-bitprec+e);
! 437: break;
! 438: case t_COMPLEX:
! 439: y=cgetg(3,t_COMPLEX);
! 440: y[1]=(long) mygprecrc((GEN)x[1],bitprec,e);
! 441: y[2]=(long) mygprecrc((GEN)x[2],bitprec,e);
! 442: break;
! 443: default: y=gcopy(x);
! 444: }
! 445: return y;
! 446: }
! 447:
! 448: /* gprec behaves badly with the zero for polynomials.
! 449: The second parameter in mygprec is the precision in base 2 */
! 450: static GEN
! 451: mygprec(GEN x, long bitprec)
! 452: {
! 453: long tx=typ(x),lx,i,e;
! 454: GEN y;
! 455:
! 456: switch(tx)
! 457: {
! 458: case t_POL:
! 459: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1]; e=gexpo(x);
! 460: for (i=2; i<lx; i++) y[i]=(long) mygprecrc((GEN)x[i],bitprec,e);
! 461: break;
! 462:
! 463: default: y=mygprecrc(x,bitprec,0);
! 464: }
! 465: return y;
! 466: }
! 467:
! 468: /* the round fonction has a bug in pari. Thus I create mygfloor, using gfloor
! 469: which has no bug (destroy z)*/
! 470: static GEN
! 471: mygfloor(GEN z)
! 472: {
! 473: if (typ(z)!=t_COMPLEX) return gfloor(z);
! 474: z[1]=lfloor((GEN)z[1]); z[2]=lfloor((GEN)z[2]); return z;
! 475: }
! 476:
! 477: /* returns a polynomial q in (Z[i])[x] keeping bitprec bits of p */
! 478: static GEN
! 479: eval_rel_pol(GEN p,long bitprec)
! 480: {
! 481: long e=gexpo(p),n=lgef(p),i,shift;
! 482: GEN q = gprec(p,(long) ((double) bitprec * L2SL10)+2);
! 483:
! 484: shift=bitprec-e+1;
! 485: for (i=2; i<n; i++)
! 486: q[i]=(long) mygfloor(myshiftic((GEN)q[i],shift));
! 487: return q;
! 488: }
! 489:
! 490: /* normalize a polynomial p, that is change it with coefficients in Z[i],
! 491: after making product by 2^bitprec */
! 492: static void
! 493: pol_to_gaussint(GEN p, long bitprec)
! 494: {
! 495: long i,n=lgef(p);
! 496: for (i=2; i<n; i++)
! 497: {
! 498: myshiftrc((GEN) p[i],bitprec);
! 499: p[i]=(long) mygfloor((GEN) p[i]);
! 500: }
! 501: }
! 502:
! 503: /* returns p(R*x)/R^n (in R or R[i]), n=deg(p), bitprec bits of precision */
! 504: static GEN
! 505: homothetie(GEN p, double R, long bitprec)
! 506: {
! 507: GEN q,r,gR,aux;
! 508: long n=lgef(p)-3,i;
! 509:
! 510: gR=mygprec(dbltor(1/R),bitprec);
! 511: q=mygprec(p,bitprec);
! 512: r=cgetg(n+3,t_POL); r[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(n+3);
! 513: aux=gR; r[n+2]=q[n+2];
! 514: for (i=n-1; i>=0; i--)
! 515: {
! 516: r[i+2]=lmul(aux,(GEN)q[i+2]);
! 517: aux=gmul(aux,gR);
! 518: }
! 519: return r;
! 520: }
! 521:
! 522: /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) */
! 523: static void
! 524: homothetie2n(GEN p, long e)
! 525: {
! 526: if (e)
! 527: {
! 528: long i,n=lgef(p)-1;
! 529: for (i=2; i<=n; i++) myshiftrc((GEN) p[i],(n-i)*e);
! 530: }
! 531: }
! 532:
! 533: /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
! 534: static void
! 535: homothetie_gauss(GEN p, long e,long f)
! 536: {
! 537: if (e||f)
! 538: {
! 539: long i, n=lgef(p)-1;
! 540: for (i=2; i<=n; i++) p[i]=(long) myshiftic((GEN) p[i],f+(n-i)*e);
! 541: }
! 542: }
! 543:
! 544: static long
! 545: valuation(GEN p)
! 546: {
! 547: long j=0,n=lgef(p)-3;
! 548:
! 549: while ((j<=n) && isexactzero((GEN)p[j+2])) j++;
! 550: return j;
! 551: }
! 552:
! 553: /* provides usually a good lower bound on the largest modulus of the roots,
! 554: puts in k an upper bound of the number of roots near the largest root
! 555: at a distance eps */
! 556: static double
! 557: lower_bound(GEN p, long *k, double eps)
! 558: {
! 559: long n=lgef(p)-3,i,j,ltop=avma;
! 560: GEN a,s,icd;
! 561: double r,*rho;
! 562:
! 563: if (n<4) { *k=n; return 0.; }
! 564: a=cgetg(6,t_POL); s=cgetg(6,t_POL);
! 565: rho=(double *) gpmalloc(4*sizeof(double));
! 566: icd=gdiv(gunr,(GEN) p[n+2]);
! 567: for (i=1; i<=4; i++) a[i+1]=lmul(icd,(GEN)p[n+2-i]);
! 568: for (i=1; i<=4; i++)
! 569: {
! 570: s[i+1]=lmulsg(i,(GEN)a[i+1]);
! 571: for (j=1; j<i; j++)
! 572: s[i+1]=ladd((GEN)s[i+1],gmul((GEN)s[j+1],(GEN)a[i+1-j]));
! 573: s[i+1]=lneg((GEN)s[i+1]);
! 574: r=gtodouble(gabs((GEN) s[i+1],3));
! 575: if (r<=0.) /* should not be strictly negative */
! 576: rho[i-1]=0.;
! 577: else
! 578: rho[i-1]=exp(log(r/(double)n)/(double) i);
! 579: }
! 580: r=0.;
! 581: for (i=0; i<4; i++) if (r<rho[i]) r=rho[i];
! 582: if (r>0. && eps<1.2)
! 583: *k=(long) floor((n*rho[0]/r+n)/(1+exp(-eps)*cos(eps)));
! 584: else
! 585: *k=n;
! 586: free(rho); avma=ltop; return r;
! 587: }
! 588:
! 589: /* returns the maximum of the modulus of p with a relative error tau */
! 590: static GEN
! 591: gmax_modulus(GEN p, double tau)
! 592: {
! 593: GEN q,aux,new_gunr;
! 594: long i,j,k,valuat,n=lgef(p)-3,nn,ltop=avma,bitprec,imax,e;
! 595: double r=0.,rho,tau2,eps;
! 596:
! 597: if (tau>3.0) tau=3.0; /* fix PZ 7.2.98: ensures eps is positive */
! 598: eps=1/log(4./tau); tau2=tau/6.;
! 599: bitprec=(long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
! 600: new_gunr=mygprec(gunr,bitprec+2*n);
! 601: aux=gdiv(new_gunr,(GEN) p[2+n]);
! 602: q=gmul(aux,p); q[2+n]=lcopy(new_gunr);
! 603: k=nn=n;
! 604: e=findpower(q); homothetie2n(q,e); r=-(double) e;
! 605: q=mygprec(q,bitprec+(n<<1));
! 606: pol_to_gaussint(q,bitprec);
! 607: imax=(long) ((log(log(4.*n))+log(3./tau))/log(2.))+2;
! 608: for (i=0,e=0;;)
! 609: {
! 610: rho=lower_bound(q,&k,eps);
! 611: if (rho>exp2(-(double) e)) e = (long) -floor(log2(rho));
! 612: r = r-(double) e/ exp2( (double) i);
! 613: if (++i == imax)
! 614: {
! 615: avma=ltop;
! 616: return gpui(dbltor(2.),dbltor(r),DEFAULTPREC);
! 617: }
! 618:
! 619: if (k<nn)
! 620: bitprec=(long) ((double) k* log2(1./tau2)+
! 621: (double) (nn-k)*log2(1./eps)+
! 622: 3*log2((double) nn))+1;
! 623: else
! 624: bitprec=(long) ((double) nn* log2(1./tau2)+
! 625: 3.*log2((double) nn))+1;
! 626: homothetie_gauss(q,e,bitprec-(long)floor(mylog2((GEN) q[2+nn])+0.5));
! 627: valuat=valuation(q);
! 628: if (valuat>0)
! 629: {
! 630: nn-=valuat; setlgef(q,nn+3);
! 631: for (j=0; j<=nn; j++) q[2+j]=q[2+valuat+j];
! 632: }
! 633: set_karasquare_limit(gexpo(q));
! 634: q = gerepileupto(ltop, graeffe(q));
! 635: tau2=1.5*tau2; eps=1/log(1./tau2);
! 636: e = findpower(q);
! 637: }
! 638: }
! 639:
! 640: static double
! 641: max_modulus(GEN p, double tau)
! 642: {
! 643: return gtodouble(gmax_modulus(p,tau));
! 644: }
! 645:
! 646: /* return the k-th modulus (in ascending order) of p, rel. error tau*/
! 647: static double
! 648: modulus(GEN p, long k, double tau)
! 649: {
! 650: GEN q,new_gunr;
! 651: long i,j,kk=k,imax,n=lgef(p)-3,nn,nnn,valuat,av,ltop=avma,bitprec,decprec,e;
! 652: double tau2,r;
! 653:
! 654: tau2=tau/6; nn=n;
! 655: bitprec= (long) ((double) n*(2.+log2(3.*(double) n)+log2(1./tau2)));
! 656: decprec=(long) ((double) bitprec * L2SL10)+1;
! 657: new_gunr=gprec(gunr,decprec);
! 658: av = avma;
! 659: q=gprec(p,decprec);
! 660: q=gmul(new_gunr,q);
! 661: e=polygone_newton(q,k);
! 662: homothetie2n(q,e);
! 663: r=(double) e;
! 664: imax=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1;
! 665: for (i=1; i<imax; i++)
! 666: {
! 667: q=eval_rel_pol(q,bitprec);
! 668:
! 669: nnn=lgef(q)-3; valuat=valuation(q);
! 670: if (valuat>0)
! 671: {
! 672: kk-=valuat;
! 673: for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j];
! 674: setlgef(q,nnn-valuat+3);
! 675: }
! 676: nn=nnn-valuat;
! 677:
! 678: set_karasquare_limit(bitprec);
! 679: q = gerepileupto(av, graeffe(q));
! 680: e=polygone_newton(q,kk);
! 681: r=r+e/exp2((double)i);
! 682: q=gmul(new_gunr,q);
! 683: homothetie2n(q,e);
! 684:
! 685: tau2=1.5*tau2; if (tau2>1.) tau2=1.;
! 686: bitprec= 1+(long) ((double) nn*(2.+log2(3.*(double) nn)+log2(1./tau2)));
! 687: }
! 688: avma=ltop; return exp2(-r);
! 689: }
! 690:
! 691: /* return the k-th modulus r_k of p, rel. error tau, knowing that
! 692: rmin < r_k < rmax. This helps because the information enable us to use
! 693: less precision... quicker ! */
! 694: static double
! 695: pre_modulus(GEN p, long k, double tau, double rmin, double rmax)
! 696: {
! 697: GEN q;
! 698: long n=lgef(p)-3,i,imax,imax2,bitprec,ltop=avma;
! 699: double aux,tau2,R;
! 700:
! 701: tau2=tau/6.; aux=sqrt(rmax/rmin)*exp(4*tau2);
! 702: imax=(long) log2(log((double)n)/log(aux));
! 703: if (imax<=0) return modulus(p,k,tau);
! 704:
! 705: R=sqrt(rmin*rmax);
! 706: bitprec=(long) ((double) n*(2.+log2(aux)+log2(1./tau2)));
! 707: q=homothetie(p,R,bitprec);
! 708: imax2=(long) ((log2(3./tau)+log2(log(4.*(double) n)) ))+1;
! 709: if (imax>imax2) imax=imax2;
! 710:
! 711: for (i=0; i<imax; i++)
! 712: {
! 713: q=eval_rel_pol(q,bitprec);
! 714: set_karasquare_limit(bitprec);
! 715: q = gerepileupto(ltop, graeffe(q));
! 716: aux=aux*aux*exp(2*tau2);
! 717: tau2=1.5*tau2;
! 718: bitprec= (long) ((double) n*(2.+log2(aux)+log2(1./(1-exp(-tau2)))));
! 719: q=gmul(mygprec(gunr,bitprec),q);
! 720: }
! 721:
! 722: aux=modulus(q,k,exp2((double)imax)*tau/3.);
! 723: R=R*exp(log(aux)*exp2(-(double)imax));
! 724: avma=ltop; return R;
! 725: }
! 726:
! 727: /* returns the minimum of the modulus of p with a relative error tau */
! 728: static double
! 729: min_modulus(GEN p, double tau)
! 730: {
! 731: long ltop=avma;
! 732: double r;
! 733:
! 734: if (isexactzero((GEN)p[2])) return 0.;
! 735: r=1./max_modulus(polrecip_i(p),tau); avma=ltop; return r;
! 736: }
! 737:
! 738: /* returns k such that r_k e^(-tau) < R < r_{ k+1 } e^tau.
! 739: l is such that you know in advance that l<= k <= n-l */
! 740: static long
! 741: dual_modulus(GEN p, double R, double tau, long l)
! 742: {
! 743: GEN q;
! 744: long i,j,imax,k,delta_k=0,n=lgef(p)-3,nn,nnn,valuat,ltop=avma,bitprec,ll=l;
! 745: double logmax,aux,tau2;
! 746:
! 747: tau2=7.*tau/8.;
! 748: bitprec=6*n-5*l+(long) ((double) n*(log2(1/tau2)+8.*tau2/7.));
! 749: q=homothetie(p,R,bitprec);
! 750: nn=n;
! 751: imax=(long)(log(log(2.*(double)n)/tau2)/log(7./4.)+1);
! 752:
! 753: for (i=0; i<imax; i++)
! 754: {
! 755: bitprec=6*nn-5*ll+(long) ((double) nn*(log2(1/tau2)+8.*tau2/7.));
! 756:
! 757: q=eval_rel_pol(q,bitprec);
! 758: nnn=lgef(q)-3; valuat=valuation(q);
! 759: if (valuat>0)
! 760: {
! 761: delta_k+=valuat;
! 762: for (j=0; j<=nnn-valuat; j++) q[2+j]=q[2+valuat+j];
! 763: setlgef(q,nnn-valuat+3);
! 764: }
! 765: ll= (-valuat<nnn-n)? ll-valuat: ll+nnn-n; /* min(ll-valuat,ll+nnn-n) */
! 766: if (ll<0) ll=0;
! 767:
! 768: nn=nnn-valuat;
! 769: if (nn==0) return delta_k;
! 770:
! 771: set_karasquare_limit(bitprec);
! 772: q = gerepileupto(ltop, graeffe(q));
! 773: tau2=tau2*7./4.;
! 774: }
! 775: k=-1; logmax=- (double) pariINFINITY;
! 776: for (i=0; i<=lgef(q)-3; i++)
! 777: {
! 778: aux=mylog2((GEN)q[2+i]);
! 779: if (aux>logmax) { logmax=aux; k=i; }
! 780: }
! 781: avma=ltop; return delta_k+k;
! 782: }
! 783:
! 784: /********************************************************************/
! 785: /** **/
! 786: /** CALCUL D'UN FACTEUR PAR INTEGRATION SUR LE CERCLE **/
! 787: /** **/
! 788: /********************************************************************/
! 789:
! 790: static GEN
! 791: gmulbyi(GEN z)
! 792: {
! 793: GEN aux = cgetg(3,t_COMPLEX);
! 794:
! 795: if (typ(z)!=t_COMPLEX)
! 796: {
! 797: aux[1]=zero;
! 798: aux[2]=(long) z;
! 799: }
! 800: else
! 801: {
! 802: aux[1]=lneg((GEN)z[2]);
! 803: aux[2]=z[1];
! 804: }
! 805: return aux;
! 806: }
! 807:
! 808: static void
! 809: fft(GEN Omega, GEN p, GEN f, long step, long l)
! 810: {
! 811: long i,l1,l2,l3,rap=Lmax/l,rapi,step4,ltop,lbot;
! 812: GEN f1,f2,f3,f02,f13,g02,g13,ff;
! 813:
! 814: if (l==2)
! 815: {
! 816: f[0]=ladd((GEN)p[0],(GEN)p[step]);
! 817: f[1]=lsub((GEN)p[0],(GEN)p[step]); return;
! 818: }
! 819: if (l==4)
! 820: {
! 821: f1=gadd((GEN)p[0],(GEN)p[(step<<1)]);
! 822: f3=gadd((GEN)p[step],(GEN)p[3*step]);
! 823: f[0]=ladd(f1,f3);
! 824: f[2]=lsub(f1,f3);
! 825:
! 826: f2=gsub((GEN)p[0],(GEN)p[(step<<1)]);
! 827: f02=gsub((GEN)p[step],(GEN)p[3*step]);
! 828: f02=gmulbyi(f02);
! 829: f[1]=ladd(f2,f02);
! 830: f[3]=lsub(f2,f02);
! 831: return;
! 832: }
! 833:
! 834: l1=(l>>2); l2=(l1<<1); l3=l1+l2; step4=(step<<2);
! 835:
! 836: ltop=avma;
! 837: fft(Omega,p,f,step4,l1);
! 838: fft(Omega,p+step,f+l1,step4,l1);
! 839: fft(Omega,p+(step<<1),f+l2,step4,l1);
! 840: fft(Omega,p+3*step,f+l3,step4,l1);
! 841:
! 842: ff=cgetg(l+1,t_VEC);
! 843: for (i=0; i<l1; i++)
! 844: {
! 845: rapi=rap*i;
! 846: f1=gmul((GEN)Omega[rapi],(GEN) f[i+l1]);
! 847: f2=gmul((GEN)Omega[(rapi<<1)],(GEN) f[i+l2]);
! 848: f3=gmul((GEN)Omega[3*rapi],(GEN) f[i+l3]);
! 849:
! 850: f02=gadd((GEN)f[i],f2);
! 851: g02=gsub((GEN)f[i],f2);
! 852: f13=gadd(f1,f3);
! 853: g13=gmulbyi(gsub(f1,f3));
! 854:
! 855: ff[i+1]=ladd(f02,f13);
! 856: ff[i+l1+1]=ladd(g02,g13);
! 857: ff[i+l2+1]=lsub(f02,f13);
! 858: ff[i+l3+1]=lsub(g02,g13);
! 859: }
! 860: lbot=avma; ff=gerepile(ltop,lbot,gcopy(ff));
! 861: for (i=0; i<l; i++) f[i]=ff[i+1];
! 862: }
! 863:
! 864: /* returns a vector RU which contains exp(2*i*k*Pi/NN), k=0..NN-1 */
! 865: static GEN
! 866: initRU(long NN, long bitprec)
! 867: {
! 868: GEN prim,aux,RU,mygpi;
! 869: long i,N2=(NN>>1),N4=(NN>>2),N8=(NN>>3),decprec;
! 870:
! 871: RU=cgetg(NN+1,t_VEC); RU++;
! 872: mygpi=mppi(bitprec/BITS_IN_LONG+3);
! 873: aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
! 874: decprec=(long) ((double) bitprec * L2SL10)+1;
! 875: prim=gexp(aux,decprec);
! 876: RU[0]=lprec(gunr,decprec);
! 877:
! 878: for (i=1; i<=N8; i++) RU[i]=lmul(prim,(GEN)RU[i-1]);
! 879: for (i=1; i<N8; i++)
! 880: {
! 881: aux=cgetg(3,t_COMPLEX);
! 882: aux[1]=((GEN)RU[i])[2]; aux[2]=((GEN)RU[i])[1];
! 883: RU[N4-i]=(long) aux;
! 884: }
! 885: for (i=0; i<N4; i++) RU[i+N4]=(long)gmulbyi((GEN)RU[i]);
! 886: for (i=0; i<N2; i++) RU[i+N2]=lneg((GEN)RU[i]);
! 887: return RU;
! 888: }
! 889:
! 890: /* returns 1 if p has only real coefficients, 0 else */
! 891: static long
! 892: isreal(GEN p)
! 893: {
! 894: long n=lgef(p)-3,i=0;
! 895:
! 896: while (i<=n && typ(p[i+2])!=t_COMPLEX) i++;
! 897: return (i>n);
! 898: }
! 899:
! 900: static void
! 901: parameters(GEN p, double *mu, double *gamma,
! 902: long polreal, double param, double param2)
! 903: {
! 904: GEN q,pc,Omega,coef,RU,prim,aux,ggamma,gx,mygpi;
! 905: long ltop=avma,limite=stack_lim(ltop,1),n=lgef(p)-3,bitprec,NN,K,i,j,ltop2;
! 906: double lx;
! 907:
! 908: bitprec=gexpo(p)+(long)param2+8;
! 909: NN=(long) (param*3.14)+1; if (NN<Lmax) NN=Lmax;
! 910: K=NN/Lmax; if (K%2==1) K++; NN=Lmax*K;
! 911: mygpi=mppi(bitprec/BITS_IN_LONG+3);
! 912: aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
! 913: prim=gexp(aux,(long) ((double) bitprec * L2SL10)+1);
! 914: RU=mygprec(gunr,bitprec);
! 915:
! 916: Omega=initRU(Lmax,bitprec);
! 917:
! 918: q=mygprec(p,bitprec);
! 919: pc=cgetg(Lmax+1,t_VEC); pc++;
! 920: for (i=n+1; i<Lmax; i++) pc[i]=zero;
! 921:
! 922: coef=cgetg(Lmax+1,t_VEC); coef++;
! 923: *mu=(double)pariINFINITY; *gamma=0.;
! 924: ggamma=gmul(gzero,gunr);
! 925: if (polreal) K=K/2+1;
! 926: ltop2=avma;
! 927: for (i=0; i<K; i++)
! 928: {
! 929: aux=mygprec(gunr,bitprec);
! 930: for (j=0; j<=n; j++)
! 931: {
! 932: pc[j]=lmul((GEN)q[j+2],aux);
! 933: aux=gmul(aux,RU); /* RU = prim^i, aux=prim^(ij) */
! 934: }
! 935:
! 936: fft(Omega,pc,coef,1,Lmax);
! 937: for (j=0; j<Lmax; j++)
! 938: {
! 939: aux=gprec((GEN)coef[j],DEFAULTPREC);
! 940: gx=gabs(aux,DEFAULTPREC);
! 941: lx=gtodouble(mplog(gx));
! 942: if (lx<*mu) *mu=lx;
! 943: if (polreal && (i>0 && i<K-1))
! 944: {
! 945: gx=gdiv(gdeux,gx);
! 946: ggamma=gadd(ggamma,gx);
! 947: }
! 948: else
! 949: ggamma=gadd(ggamma,ginv(gx));
! 950: }
! 951: RU=gmul(RU,prim);
! 952: if (low_stack(limite, stack_lim(ltop,1)))
! 953: {
! 954: GEN *gptr[2];
! 955: if(DEBUGMEM>1) err(warnmem,"parameters");
! 956: gptr[0]=&ggamma; gptr[1]=&RU; gerepilemany(ltop2,gptr,2);
! 957: }
! 958: }
! 959: ggamma=gdivgs(ggamma,NN);
! 960: *gamma=gtodouble(glog(ggamma,DEFAULTPREC))/log(2.);
! 961: avma=ltop;
! 962: }
! 963:
! 964: /* NN is a multiple of Lmax */
! 965: static void
! 966: dft(GEN p, long k, long NN, long bitprec, GEN F, GEN H, long polreal)
! 967: {
! 968: GEN Omega,q,qd,pc,pdc,alpha,beta,gamma,RU,aux,U,W,mygpi,prim,prim2;
! 969: long limite,n=lgef(p)-3,i,j,K,ltop;
! 970:
! 971: mygpi=mppi(bitprec/BITS_IN_LONG+3);
! 972: aux=gmul(gi,gdivgs(mygpi,NN/2)); /* 2i Pi/NN */
! 973: prim=gexp(aux,(long) ((double) bitprec * L2SL10)+1);
! 974: prim2=mygprec(gunr,bitprec);
! 975: RU=cgetg(n+2,t_VEC); RU++;
! 976:
! 977: Omega=initRU(Lmax,bitprec);
! 978: K=NN/Lmax; q=mygprec(p,bitprec);
! 979: qd=derivpol(q);
! 980:
! 981: pc=cgetg(Lmax+1,t_VEC); pc++;
! 982: for (i=n+1; i<Lmax; i++) pc[i]=zero;
! 983: pdc=cgetg(Lmax+1,t_VEC); pdc++;
! 984: for (i=n; i<Lmax; i++) pdc[i]=zero;
! 985:
! 986: alpha=cgetg(Lmax+1,t_VEC); alpha++;
! 987: beta=cgetg(Lmax+1,t_VEC); beta++;
! 988: gamma=cgetg(Lmax+1,t_VEC); gamma++;
! 989:
! 990: if (polreal) K=K/2+1;
! 991:
! 992: ltop=avma; limite = stack_lim(ltop,1);
! 993: W=cgetg(k+1,t_VEC); U=cgetg(k+1,t_VEC);
! 994: for (i=1; i<=k; i++) W[i]=U[i]=zero;
! 995:
! 996: for (i=0; i<K; i++)
! 997: {
! 998: RU[0]=(long) gun;
! 999: for (j=1; j<=n; j++) RU[j]=lmul((GEN)RU[j-1],prim2);
! 1000: /* RU[j]=prim^{ ij }=prim2^j */
! 1001:
! 1002: for (j=0; j<n; j++) pdc[j]=lmul((GEN)qd[j+2],(GEN)RU[j]);
! 1003: fft(Omega,pdc,alpha,1,Lmax);
! 1004: for (j=0; j<=n; j++) pc[j]=lmul((GEN)q[j+2],(GEN)RU[j]);
! 1005: fft(Omega,pc,beta,1,Lmax);
! 1006: for (j=0; j<Lmax; j++) gamma[j]=linv((GEN)beta[j]);
! 1007: for (j=0; j<Lmax; j++) beta[j]=lmul((GEN)alpha[j],(GEN)gamma[j]);
! 1008: fft(Omega,beta,alpha,1,Lmax);
! 1009: fft(Omega,gamma,beta,1,Lmax);
! 1010:
! 1011: if (polreal) /* p has real coefficients */
! 1012: {
! 1013: if (i>0 && i<K-1)
! 1014: {
! 1015: for (j=1; j<=k; j++)
! 1016: {
! 1017: aux=gmul((GEN)alpha[j+1],(GEN)RU[j+1]);
! 1018: W[j]=ladd((GEN)W[j],gshift(greal(aux),1));
! 1019: aux=gmul((GEN)beta[j],(GEN)RU[j]);
! 1020: U[j]=ladd((GEN)U[j],gshift(greal(aux),1));
! 1021: }
! 1022: }
! 1023: else
! 1024: {
! 1025: for (j=1; j<=k; j++)
! 1026: {
! 1027: aux=gmul((GEN)alpha[j+1],(GEN)RU[j+1]);
! 1028: W[j]=ladd((GEN)W[j],greal(aux));
! 1029: aux=gmul((GEN)beta[j],(GEN)RU[j]);
! 1030: U[j]=ladd((GEN)U[j],greal(aux));
! 1031: }
! 1032: }
! 1033: }
! 1034: else
! 1035: {
! 1036: for (j=1; j<=k; j++)
! 1037: {
! 1038: W[j]=ladd((GEN)W[j],gmul((GEN)alpha[j+1],(GEN)RU[j+1]));
! 1039: U[j]=ladd((GEN)U[j],gmul((GEN)beta[j],(GEN)RU[j]));
! 1040: }
! 1041: }
! 1042: prim2=gmul(prim2,prim);
! 1043: if (low_stack(limite, stack_lim(ltop,1)))
! 1044: {
! 1045: GEN *gptr[3];
! 1046: if(DEBUGMEM>1) err(warnmem,"dft");
! 1047: gptr[0]=&W; gptr[1]=&U; gptr[2]=&prim2;
! 1048: gerepilemany(ltop,gptr,3);
! 1049: }
! 1050: }
! 1051:
! 1052: for (i=1; i<=k; i++)
! 1053: {
! 1054: aux=(GEN)W[i];
! 1055: for (j=1; j<i; j++) aux=gadd(aux,gmul((GEN)W[i-j],(GEN)F[k+2-j]));
! 1056: F[k+2-i] = ldivgs(aux,-i*NN);
! 1057: }
! 1058: for (i=0; i<k; i++)
! 1059: {
! 1060: aux=(GEN)U[k-i];
! 1061: for (j=1+i; j<k; j++) aux=gadd(aux,gmul((GEN)F[2+j],(GEN)U[j-i]));
! 1062: H[i+2] = ldivgs(aux,NN);
! 1063: }
! 1064: }
! 1065:
! 1066: static GEN
! 1067: refine_H(GEN F, GEN G, GEN HH, long bitprec, long shiftbitprec)
! 1068: {
! 1069: GEN H=HH,D,aux;
! 1070: long ltop=avma, limite=stack_lim(ltop,1),error=0,i,bitprec1,bitprec2,lbot;
! 1071:
! 1072: D=gsub(gun,gres(gmul(HH,G),F)); error=gexpo(D);
! 1073: bitprec2=bitprec+shiftbitprec;
! 1074:
! 1075: for (i=0; (error>-bitprec && i<NEWTON_MAX) && error<=0; i++)
! 1076: {
! 1077: if (low_stack(limite, stack_lim(ltop,1)))
! 1078: {
! 1079: GEN *gptr[2];
! 1080: if(DEBUGMEM>1) err(warnmem,"refine_H");
! 1081: gptr[0]=&D; gptr[1]=&H; gerepilemany(ltop,gptr,2);
! 1082: }
! 1083: bitprec1=-error+shiftbitprec;
! 1084: aux=gmul(mygprec(H,bitprec1),mygprec(D,bitprec1));
! 1085: aux=mygprec(aux,bitprec1);
! 1086: aux=gres(aux,mygprec(F,bitprec1));
! 1087:
! 1088: bitprec1=-error*2+shiftbitprec;
! 1089: if (bitprec1>bitprec2) bitprec1=bitprec2;
! 1090: H=gadd(mygprec(H,bitprec1),aux);
! 1091: D=gsub(gun,gres(gmul(H,G),F));
! 1092: error=gexpo(D); if (error<-bitprec1) error=-bitprec1;
! 1093: }
! 1094: if (error<=-bitprec/2) { lbot=avma; return gerepile(ltop,lbot,gcopy(H)); }
! 1095: avma=ltop; return gzero; /* procedure failed */
! 1096: }
! 1097:
! 1098: /* return 0 if fails, 1 else */
! 1099: static long
! 1100: refine_F(GEN p, GEN *F, GEN *G, GEN H, long bitprec, double gamma)
! 1101: {
! 1102: GEN pp,FF,GG,r,HH,f0;
! 1103: long error,i,bitprec1=0,bitprec2,ltop=avma,shiftbitprec;
! 1104: long shiftbitprec2,n=lgef(p)-3,enh,normF,normG,limite=stack_lim(ltop,1);
! 1105:
! 1106: FF=*F; HH=H;
! 1107: GG=poldivres(p,*F,&r);
! 1108: normF=gexpo(FF);
! 1109: normG=gexpo(GG);
! 1110: enh=gexpo(H); if (enh<0) enh=0;
! 1111: shiftbitprec=normF+2*normG+enh+(long) (4.*log2((double)n)+gamma) +1;
! 1112: shiftbitprec2=enh+2*(normF+normG)+(long) (2.*gamma+5.*log2((double)n))+1;
! 1113: bitprec2=bitprec+shiftbitprec;
! 1114: error=gexpo(r);
! 1115: if (error<-bitprec) error=1-bitprec;
! 1116: for (i=0; (error>-bitprec && i<NEWTON_MAX) && error<=0; i++)
! 1117: {
! 1118: if ((bitprec1==bitprec2) && (i>=2))
! 1119: {
! 1120: shiftbitprec+=n; shiftbitprec2+=n; bitprec2+=n;
! 1121: }
! 1122: if (low_stack(limite, stack_lim(ltop,1)))
! 1123: {
! 1124: GEN *gptr[4];
! 1125: if(DEBUGMEM>1) err(warnmem,"refine_F");
! 1126: gptr[0]=&FF; gptr[1]=&GG; gptr[2]=&r; gptr[3]=&HH;
! 1127: gerepilemany(ltop,gptr,4);
! 1128: }
! 1129:
! 1130: bitprec1=-error+shiftbitprec2;
! 1131: HH=refine_H(mygprec(FF,bitprec1),mygprec(GG,bitprec1),
! 1132: mygprec(HH,bitprec1),1-error,shiftbitprec2);
! 1133: if (HH==gzero) return 0; /* procedure failed */
! 1134:
! 1135: bitprec1=-error+shiftbitprec;
! 1136: r=gmul(mygprec(HH,bitprec1),mygprec(r,bitprec1));
! 1137: r=mygprec(r,bitprec1);
! 1138: f0=gres(r,mygprec(FF,bitprec1));
! 1139:
! 1140: bitprec1=-2*error+shiftbitprec;
! 1141: if (bitprec1>bitprec2) bitprec1=bitprec2;
! 1142: FF=gadd(mygprec(FF,bitprec1),f0);
! 1143:
! 1144: bitprec1=-3*error+shiftbitprec;
! 1145: if (bitprec1>bitprec2) bitprec1=bitprec2;
! 1146: pp=mygprec(p,bitprec1);
! 1147: GG=poldivres(pp,mygprec(FF,bitprec1),&r);
! 1148: error=gexpo(r); if (error<-bitprec1) error=-bitprec1;
! 1149: }
! 1150: if (error<=-bitprec)
! 1151: {
! 1152: *F=FF; *G=GG;
! 1153: return 1; /* procedure succeeds */
! 1154: }
! 1155: return 0; /* procedure failed */
! 1156: }
! 1157:
! 1158: /* returns F and G from the unit circle U such that |p-FG|<2^(-bitprec) |cd|,
! 1159: where cd is the leading coefficient of p */
! 1160: static void
! 1161: split_fromU(GEN p, long k, double delta, long bitprec,
! 1162: GEN *F, GEN *G, double param, double param2)
! 1163: {
! 1164: GEN pp,FF,GG,H;
! 1165: long n=lgef(p)-3,NN,bitprec2,
! 1166: ltop=avma,polreal=isreal(p);
! 1167: double mu,gamma;
! 1168:
! 1169: pp=gdiv(p,(GEN)p[2+n]);
! 1170: Lmax=4; while (Lmax<=n) Lmax=(Lmax<<1);
! 1171: parameters(pp,&mu,&gamma,polreal,param,param2);
! 1172:
! 1173: H =cgetg(k+2,t_POL); H[1] =evalsigne(1) | evalvarn(varn(p)) | evallgef(k+2);
! 1174: FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
! 1175: FF[k+2]=un;
! 1176:
! 1177: NN=(long) (0.5/delta); NN+=(NN%2); if (NN<2) NN=2;
! 1178: NN=NN*Lmax; ltop=avma;
! 1179: for(;;)
! 1180: {
! 1181: bitprec2=(long) (((double) NN*delta-mu)/log(2.))+gexpo(pp)+8;
! 1182: dft(pp,k,NN,bitprec2,FF,H,polreal);
! 1183: if (refine_F(pp,&FF,&GG,H,bitprec,gamma)) break;
! 1184: NN=(NN<<1); avma=ltop;
! 1185: }
! 1186: *G=gmul(GG,(GEN)p[2+n]); *F=FF;
! 1187: }
! 1188:
! 1189: static void
! 1190: optimize_split(GEN p, long k, double delta, long bitprec,
! 1191: GEN *F, GEN *G, double param, double param2)
! 1192: {
! 1193: long n=lgef(p)-3;
! 1194: GEN FF,GG;
! 1195:
! 1196: if (k<=n/2)
! 1197: split_fromU(p,k,delta,bitprec,F,G,param,param2);
! 1198: else
! 1199: { /* start from the reciprocal of p */
! 1200: split_fromU(polrecip_i(p),n-k,delta,bitprec,&FF,&GG,param,param2);
! 1201: *F=polrecip(GG); *G=polrecip(FF);
! 1202: }
! 1203: }
! 1204:
! 1205: /********************************************************************/
! 1206: /** **/
! 1207: /** RECHERCHE DU CERCLE DE SEPARATION **/
! 1208: /** **/
! 1209: /********************************************************************/
! 1210:
! 1211: /* return p(2^e*x) *2^(-n*e) */
! 1212: static void
! 1213: scalepol2n(GEN p, long e)
! 1214: {
! 1215: long i,n=lgef(p)-1;
! 1216: for (i=2; i<=n; i++) p[i]=lmul2n((GEN)p[i],(i-n)*e);
! 1217: }
! 1218:
! 1219: /* returns p(x/R)*R^n */
! 1220: static GEN
! 1221: scalepol(GEN p, GEN R, long bitprec)
! 1222: {
! 1223: GEN q,aux,gR;
! 1224: long i;
! 1225:
! 1226: aux = gR = mygprec(R,bitprec); q = mygprec(p,bitprec);
! 1227: for (i=lgef(p)-2; i>=2; i--)
! 1228: {
! 1229: q[i]=lmul(aux,(GEN)q[i]);
! 1230: aux = gmul(aux,gR);
! 1231: }
! 1232: return q;
! 1233: }
! 1234:
! 1235: /* returns q(x)=p(b), where b is (usually) a polynomial */
! 1236: static GEN
! 1237: shiftpol(GEN p, GEN b)
! 1238: {
! 1239: long av = avma,i, limit = stack_lim(av,1);
! 1240: GEN q = gzero;
! 1241:
! 1242: for (i=lgef(p)-1; i>=2; i--)
! 1243: {
! 1244: q = gadd((GEN)p[i], gmul(q,b));
! 1245: if (low_stack(limit, stack_lim(av,1))) q = gerepileupto(av,q);
! 1246: }
! 1247: return gerepileupto(av,q);
! 1248: }
! 1249:
! 1250: /* return (aX-1)^n * p[ (X-a) / (aX-1) ] */
! 1251: static GEN
! 1252: conformal_pol(GEN p, GEN a, long bitprec)
! 1253: {
! 1254: GEN r,pui,num,aux;
! 1255: long n=lgef(p)-3, i;
! 1256:
! 1257: aux = pui = cgetg(4,t_POL);
! 1258: pui[1] = evalsigne(1) | evalvarn(varn(p)) | evallgef(4);
! 1259: pui[2] = (long) mygprec(gneg_i(gunr),bitprec);
! 1260: pui[3] = lconj(a);
! 1261: num = cgetg(4,t_POL);
! 1262: num[1] = pui[1];
! 1263: num[2] = lneg(a);
! 1264: num[3] = (long) mygprec(gunr,bitprec);
! 1265: r=(GEN)p[2+n];
! 1266: for (i=n-1; ; i--)
! 1267: {
! 1268: r=gadd(gmul(r,num),gmul(aux,(GEN) p[2+i]));
! 1269: if (i==0) return r;
! 1270: aux=gmul(pui,aux);
! 1271: }
! 1272: }
! 1273:
! 1274: /* apply the conformal mapping then split from U */
! 1275: static void
! 1276: conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
! 1277: {
! 1278: long bitprec2,bitprec3,n=lgef(p)-3,decprec,i,ltop = avma, av;
! 1279: GEN q,FF,GG,a,R, *gptr[2];
! 1280: double rmin,rmax,rho,delta,aux2,param,param2,r1,r2;
! 1281:
! 1282: bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
! 1283: a=gsqrt(stoi(3),10);
! 1284: a=gmul(mygprec(a,bitprec2),mygprec(globalu,bitprec2));
! 1285: a=gdivgs(a,-6); /* a=-globalu/2/sqrt(3) */
! 1286:
! 1287: av = avma; q=mygprec(p,bitprec2);
! 1288: q = conformal_pol(q,a,bitprec2);
! 1289: for (i=1; i<=n; i++)
! 1290: if (radius[i]!=0.) /* updating array radius */
! 1291: {
! 1292: aux2=radius[i]*radius[i];
! 1293: aux2=2.*(aux2-1)/(aux2-3.*radius[i]+3.);
! 1294: radius[i]=sqrt(1+aux2);
! 1295: }
! 1296: if (k>1)
! 1297: {
! 1298: i=k-1; while (i>0 && radius[i]==0.) i--;
! 1299: r1=radius[i]; r2=radius[k];
! 1300: rmin=pre_modulus(q,k,aux/10.,r1,r2);
! 1301: }
! 1302: else
! 1303: rmin=min_modulus(q,aux/10.);
! 1304: radius[k]=rmin;
! 1305:
! 1306: if (k+1<n)
! 1307: {
! 1308: i=k+2; while (i<=n && radius[i]==0.) i++;
! 1309: r2=radius[i]; r1=radius[k+1];
! 1310: rmax=pre_modulus(q,k+1,aux/10.,r1,r2);
! 1311: }
! 1312: else /* k+1=n */
! 1313: rmax=max_modulus(q,aux/10.);
! 1314: radius[k+1]=rmax;
! 1315:
! 1316: rho=sqrt(rmin*rmax); delta=0.5*log(rmax/rmin);
! 1317: if (delta>1.) delta=1.;
! 1318:
! 1319: if (rho<1.) bitprec3=(long) ((double)n*log2(1./rho))+1;
! 1320: else bitprec3=(long) ((double)n*log2(rho))+1;
! 1321: bitprec2=bitprec2+bitprec3;
! 1322:
! 1323: R=mygprec(dbltor(1/rho),bitprec2);
! 1324: q=scalepol(q,R,bitprec2);
! 1325: gptr[0]=&q; gptr[1]=&R; gerepilemany(av,gptr,2);
! 1326:
! 1327: aux2=radius[k];
! 1328: for (i=k-1; i>=1; i--)
! 1329: {
! 1330: if (radius[i]==0.) radius[i]=aux2;
! 1331: else aux2=radius[i];
! 1332: }
! 1333: aux2=radius[k+1];
! 1334: for (i=k+1; i<=n; i++)
! 1335: {
! 1336: if (radius[i]==0.) radius[i]=aux2;
! 1337: else aux2=radius[i];
! 1338: }
! 1339: param=0.; param2=0.;
! 1340: for (i=1; i<=n; i++)
! 1341: {
! 1342: radius[i]=radius[i]/rho;
! 1343: aux2=fabs(1-radius[i]);
! 1344: param+=1./aux2;
! 1345: if (aux2<1.) param2-=log2(aux2);
! 1346: }
! 1347: optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
! 1348: bitprec2=bitprec2+n; R = ginv(R);
! 1349: FF=scalepol(FF,R,bitprec2);
! 1350: GG=scalepol(GG,R,bitprec2);
! 1351:
! 1352: a=mygprec(a,bitprec2);
! 1353: FF=conformal_pol(FF,a,bitprec2);
! 1354: GG=conformal_pol(GG,a,bitprec2);
! 1355: a=ginv(gsub(gun,gmul(a,gconj(a))));
! 1356: a=glog(a,(long) (bitprec2 * L2SL10)+1);
! 1357:
! 1358: decprec=(long) ((bitprec+n) * L2SL10)+1;
! 1359: FF=gmul(FF,gexp(gmulgs(a,k),decprec));
! 1360: GG=gmul(GG,gexp(gmulgs(a,n-k),decprec));
! 1361:
! 1362: *F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
! 1363: gptr[0]=F; gptr[1]=G; gerepilemany(ltop,gptr,2);
! 1364: }
! 1365:
! 1366: /* split p, this time with no scaling. returns in F and G two polynomials
! 1367: such that |p-FG|< 2^(-bitprec)|p| */
! 1368: static void
! 1369: split_2(GEN p, long bitprec, double thickness, GEN *F, GEN *G)
! 1370: {
! 1371: double rmin,rmax,rho,kappa,aux,delta,param,param2,r1,r2;
! 1372: long n=lgef(p)-3,i,j,k,disti,diste,bitprec2;
! 1373: GEN q,FF,GG,R;
! 1374:
! 1375: radius=(double *) gpmalloc((n+1)*sizeof(double));
! 1376: for (i=2; i<n; i++) radius[i]=0.;
! 1377: rmin=min_modulus(p,thickness/(double) n/4.);
! 1378: rmax=max_modulus(p,thickness/(double) n/4.);
! 1379: i=1; j=n; radius[1]=rmin; radius[n]=rmax;
! 1380: rho=sqrt(rmin*rmax);
! 1381: k=dual_modulus(p,rho,thickness/(double) n/4.,1);
! 1382: if (k<n/5. || (n/2.<k && k<(4*n)/5.)) { rmax=rho; j=k+1; radius[j]=rho; }
! 1383: else { rmin=rho; i=k; radius[i]=rho; }
! 1384: while (j>i+1)
! 1385: {
! 1386: disti= (i<n-j)? i: n-j; /* min(i,n-j) */
! 1387: diste= (j<n-i)? j: n-i;
! 1388: kappa=1.-log(1.+(double)disti)/log(1.+(double)diste);
! 1389: if (i+j<n+1)
! 1390: rho=exp( (log(rmin)+log(rmax)*(1+kappa))/(2+kappa));
! 1391: else if (i+j==n+1) rho=sqrt(rmin*rmax);
! 1392: else
! 1393: rho=exp( (log(rmin)*(1+kappa)+log(rmax))/(2+kappa));
! 1394: /* use log(rmax) - log(rmin) since rmax / rmin can overflow */
! 1395: aux=(log(rmax)-log(rmin))/(j-i)/4.;
! 1396:
! 1397: k=(i<n+1-j)? i: n+1-j; /* min(i,n+1-j) */
! 1398: k=dual_modulus(p,rho,aux,k);
! 1399: if (k-i < j-k-1) { rmax=rho; j=k+1; radius[j]=rho*exp(-aux); }
! 1400: else
! 1401: {
! 1402: if (k-i > j-k-1) { rmin=rho; i=k; radius[i]=rho*exp(aux); }
! 1403: else
! 1404: {
! 1405: if (2*k>n) { rmax=rho; j=k+1; radius[j]=rho*exp(-aux); }
! 1406: else { rmin=rho; i=k; radius[i]=rho*exp(aux); }
! 1407: }
! 1408: }
! 1409: }
! 1410: aux=log(rmax)-log(rmin);
! 1411:
! 1412: if (step4==0)
! 1413: {
! 1414: if (k>1)
! 1415: {
! 1416: i=k-1; while ((i>0) && (radius[i]==0.)) i--;
! 1417: r1=radius[i]; r2=radius[k];
! 1418: rmin=pre_modulus(p,k,aux/10.,r1,r2);
! 1419: }
! 1420: else /* k=1 */
! 1421: rmin=min_modulus(p,aux/10.);
! 1422: radius[k]=rmin;
! 1423:
! 1424: if (k+1<n)
! 1425: {
! 1426: i=k+2; while ((i<=n) && (radius[i]==0.)) i++;
! 1427: r2=radius[i]; r1=radius[k+1];
! 1428: rmax=pre_modulus(p,k+1,aux/10.,r1,r2);
! 1429: }
! 1430: else /* k+1=n */
! 1431: rmax=max_modulus(p,aux/10.);
! 1432: radius[k+1]=rmax;
! 1433: rho=sqrt(rmin*rmax); delta=0.5*(log(rmax)-log(rmin));
! 1434:
! 1435: aux=radius[k];
! 1436: for (i=k-1; i>=1; i--)
! 1437: {
! 1438: if (radius[i]==0.) radius[i]=aux;
! 1439: else aux=radius[i];
! 1440: }
! 1441: aux=radius[k+1];
! 1442: for (i=k+1; i<=n; i++)
! 1443: {
! 1444: if (radius[i]==0.) radius[i]=aux;
! 1445: else aux=radius[i];
! 1446: }
! 1447: param=0.; param2=0.;
! 1448: for (i=1; i<=n; i++)
! 1449: {
! 1450: radius[i]=radius[i]/rho;
! 1451: aux=fabs(1.-radius[i]);
! 1452: param+=1./aux;
! 1453: if (aux<1) param2-=log2(aux);
! 1454: }
! 1455: if (rho<1.) bitprec2=(long) ((double)n*log2(1./rho))+1;
! 1456: else bitprec2=(long) ((double)n*log2(rho))+1;
! 1457: bitprec2 += bitprec;
! 1458:
! 1459: R=mygprec(dbltor(1./rho),bitprec2);
! 1460: q=scalepol(p,R,bitprec2);
! 1461: optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
! 1462: bitprec2 += n; R=ginv(R);
! 1463: }
! 1464: else
! 1465: {
! 1466: rho=sqrt(rmax*rmin);
! 1467: if (rho<1.) bitprec2=(long) ((double)n*log2(1./rho))+1;
! 1468: else bitprec2=(long) ((double)n*log2(rho))+1;
! 1469: bitprec2=bitprec+bitprec2;
! 1470:
! 1471: R=mygprec(dbltor(1/rho),bitprec2);
! 1472: q=scalepol(p,R,bitprec2);
! 1473: for (i=1; i<=n; i++)
! 1474: if (radius[i]!=0.) radius[i]=radius[i]/rho;
! 1475:
! 1476: conformal_mapping(q, k, bitprec2, aux, &FF, &GG);
! 1477: bitprec2 += n; R = ginv(mygprec(R,bitprec2));
! 1478: }
! 1479: free(radius);
! 1480: FF=scalepol(FF,R,bitprec2); GG=scalepol(GG,R,bitprec2);
! 1481: *F=mygprec(FF,bitprec+n); *G=mygprec(GG,bitprec+n);
! 1482: }
! 1483:
! 1484: /* procedure corresponding to steps 5,6,.. page 44 in the RR n. 1852 */
! 1485: /* put in F and G two polynomial such that |p-FG|<2^(-bitprec)|p|
! 1486: where the maximum modulus of the roots of p is <=1 and the sum of roots
! 1487: is zero */
! 1488:
! 1489: static void
! 1490: split_1(GEN p, long bitprec, GEN *F, GEN *G)
! 1491: {
! 1492: long bitprec2,bitprec3,i,imax,n=lgef(p)-3, polreal = isreal(p);
! 1493: double rmax,rmin,thickness,quo;
! 1494: GEN q,qq,newq,FF,GG,v,gr,r;
! 1495:
! 1496: r=gmax_modulus(p,0.01);
! 1497: bitprec2=bitprec+n;
! 1498: gr=mygprec(ginv(r),bitprec2);
! 1499: q=scalepol(p,gr,bitprec2);
! 1500:
! 1501: bitprec2=bitprec+gexpo(q)-gexpo(p);
! 1502:
! 1503: bitprec3=bitprec2+(long)((double)n*2.*log2(3.)+1);
! 1504: v=cgetg(5,t_VEC); v++;
! 1505: v[0]=lmulgs(mygprec(gunr,bitprec3),2); v[1]=lneg((GEN)v[0]);
! 1506: v[2]=lmul((GEN)v[0],gi); v[3]=lneg((GEN)v[2]);
! 1507: imax = polreal? 3: 4;
! 1508: q=mygprec(q,bitprec3); thickness=1.;
! 1509: for (i=0; i<imax; i++)
! 1510: {
! 1511: qq=shiftpol(q,gadd(polx[varn(p)],(GEN)v[i]));
! 1512: rmin=min_modulus(qq,0.05);
! 1513: if (3./rmin > thickness)
! 1514: {
! 1515: rmax=max_modulus(qq,0.05); quo = rmax/rmin;
! 1516: if ((float)quo > (float)thickness)
! 1517: {
! 1518: thickness=quo; newq=qq; globalu=(GEN)v[i];
! 1519: }
! 1520: }
! 1521: if (thickness>2.) break;
! 1522: if (polreal && (i==1 && thickness>1.5)) break;
! 1523: }
! 1524: bitprec3=bitprec2+(long)((double) n*log2(3.)+1)+gexpo(newq)-gexpo(q);
! 1525: split_2(newq,bitprec3,log(thickness),&FF,&GG);
! 1526: globalu=gsub(polx[varn(p)],mygprec(globalu,bitprec3));
! 1527: FF=shiftpol(FF,globalu); GG=shiftpol(GG,globalu);
! 1528:
! 1529: gr=ginv(gr);
! 1530: bitprec2=bitprec2+gexpo(FF)+gexpo(GG)-gexpo(q);
! 1531: *F=scalepol(FF,gr,bitprec2); *G=scalepol(GG,gr,bitprec2);
! 1532: }
! 1533:
! 1534: /* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P|,
! 1535: where the maximum modulus of the roots of p is <0.5 */
! 1536: static void
! 1537: split_0_2(GEN p, long bitprec, GEN *F, GEN *G)
! 1538: {
! 1539: GEN q,b,FF,GG;
! 1540: long n=lgef(p)-3,k,bitprec2,i;
! 1541: double auxnorm;
! 1542:
! 1543: step4=1;
! 1544: auxnorm=(double) n*log2(1+exp2(mylog2((GEN)p[n+1])
! 1545: -mylog2((GEN)p[n+2]))/(double)n);
! 1546: bitprec2=bitprec+1+(long) (log2((double)n)+auxnorm);
! 1547:
! 1548: q=mygprec(p,bitprec2);
! 1549: b=gdivgs(gdiv((GEN)q[n+1],(GEN)q[n+2]),-n);
! 1550: q=shiftpol(q,gadd(polx[varn(p)],b));
! 1551:
! 1552: k=0; while
! 1553: (gexpo((GEN)q[k+2])<-(bitprec2+2*(n-k)+gexpo(q))
! 1554: || gcmp0((GEN)q[k+2])) k++;
! 1555: if (k>0)
! 1556: {
! 1557: if (k>n/2) k=n/2;
! 1558: bitprec2+=(k<<1);
! 1559: FF=cgetg(k+3,t_POL); FF[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(k+3);
! 1560: for (i=0; i<k; i++) FF[i+2]=zero;
! 1561: FF[k+2]=(long) mygprec(gunr,bitprec2);
! 1562: GG=cgetg(n-k+3,t_POL); GG[1]=evalsigne(1)+evalvarn(varn(p))+evallgef(n-k+3);
! 1563: for (i=0; i<=n-k; i++) GG[i+2]=q[i+k+2];
! 1564: b=gsub(polx[varn(p)],mygprec(b,bitprec2));
! 1565: }
! 1566: else
! 1567: {
! 1568: split_1(q,bitprec2,&FF,&GG);
! 1569: bitprec2=bitprec+gexpo(FF)+gexpo(GG)-gexpo(p)+(long)auxnorm+1;
! 1570: b=gsub(polx[varn(p)],mygprec(b,bitprec2));
! 1571: FF=mygprec(FF,bitprec2);
! 1572: }
! 1573: GG=mygprec(GG,bitprec2);
! 1574: *F=shiftpol(FF,b); *G=shiftpol(GG,b);
! 1575: }
! 1576:
! 1577: /* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P|,
! 1578: where the maximum modulus of the roots of p is <2 */
! 1579: static void
! 1580: split_0_1(GEN p, long bitprec, GEN *F, GEN *G)
! 1581: {
! 1582: GEN q,FF,GG;
! 1583: long n=lgef(p)-3,bitprec2,normp,pow;
! 1584: double aux;
! 1585:
! 1586: normp=gexpo(p);
! 1587: aux=exp2(mylog2((GEN)p[n+1])-mylog2((GEN)p[n+2]))/(double)n;
! 1588: pow=2;
! 1589: if (aux<=2.5){ split_0_2(p,bitprec,F,G); return; }
! 1590:
! 1591: if (aux<=1.) pow=1;
! 1592: scalepol2n(p,pow); /* p <- 4^(-n) p(4*x) */
! 1593: bitprec2=bitprec+pow*n+gexpo(p)-normp;
! 1594: q=mygprec(p,bitprec2);
! 1595: split_1(q,bitprec2,&FF,&GG);
! 1596: scalepol2n(FF,-pow); scalepol2n(GG,-pow);
! 1597: bitprec2=bitprec+gexpo(FF)+gexpo(GG)-normp;
! 1598: *F=mygprec(FF,bitprec2); *G=mygprec(GG,bitprec2);
! 1599: }
! 1600:
! 1601: /* put in F and G two polynomials such that |P-FG|<2^(-bitprec)|P| */
! 1602: static void
! 1603: split_0(GEN p, long bitprec, GEN *F, GEN *G)
! 1604: {
! 1605: GEN FF,GG,q,R;
! 1606: long n=lgef(p)-3,k=0,i;
! 1607:
! 1608: while (gexpo((GEN)p[k+2])<-bitprec) k++;
! 1609: if (k>0)
! 1610: {
! 1611: if (k>n/2) k=n/2;
! 1612: FF=cgetg(k+3,t_POL);
! 1613: FF[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(k+3);
! 1614: for (i=0; i<k; i++) FF[i+2]=lcopy(gzero);
! 1615: FF[k+2]=(long) mygprec(gunr,bitprec);
! 1616: GG=cgetg(n-k+3,t_POL);
! 1617: GG[1]=evalsigne(1) | evalvarn(varn(p)) | evallgef(n-k+3);
! 1618: for (i=0; i<=n-k; i++) GG[i+2]=p[i+k+2];
! 1619: }
! 1620: else
! 1621: {
! 1622: R = gmax_modulus(p,0.05);
! 1623: if (gexpo(R)<1 && gtodouble(R)<1.9) split_0_1(p,bitprec,&FF,&GG);
! 1624: else
! 1625: {
! 1626: q=cgetg(n+3,t_POL); q[1]=p[1];
! 1627: for (i=0; i<=n; i++) q[i+2]=p[n-i+2]; /* p^* with copy of ptr */
! 1628: R = gmax_modulus(q,0.05);
! 1629: if (gexpo(R)<1 && gtodouble(R)<1.9)
! 1630: {
! 1631: split_0_1(q,bitprec,&FF,&GG);
! 1632: FF=polrecip(FF); GG=polrecip(GG);
! 1633: }
! 1634: else
! 1635: {
! 1636: step4=0;
! 1637: split_2(p,bitprec,1.2837,&FF,&GG);
! 1638: }
! 1639: }
! 1640: }
! 1641: *F=FF; *G=GG;
! 1642: }
! 1643:
! 1644: /********************************************************************/
! 1645: /** **/
! 1646: /** CALCUL A POSTERIORI DE L'ERREUR ABSOLUE SUR LES RACINES **/
! 1647: /** **/
! 1648: /********************************************************************/
! 1649:
! 1650: static GEN
! 1651: root_error(long n, long k, GEN roots_pol, GEN sigma, GEN shatzle)
! 1652: {
! 1653: GEN rho,d,eps,epsbis,eps2,prod,aux,rap=NULL;
! 1654: long i,j,m;
! 1655:
! 1656: d=cgetg(n+1,t_VEC);
! 1657: for (i=1; i<=n; i++)
! 1658: {
! 1659: if (i!=k)
! 1660: {
! 1661: aux=gsub((GEN)roots_pol[i],(GEN)roots_pol[k]);
! 1662: d[i]=(long) gabs(mygprec(aux,31),4);
! 1663: }
! 1664: }
! 1665: rho=gabs(mygprec((GEN)roots_pol[k],31),4);
! 1666: if (gcmp(rho,dbltor(1.))==-1) rho=gun;
! 1667: eps=gmul(rho,shatzle);
! 1668: aux=gmul(gpui(rho,stoi(n),4),sigma);
! 1669:
! 1670: for (j=1; j<=2 || (j<=5 && gcmp(rap,dbltor(1.2))==1); j++)
! 1671: {
! 1672: m=n; prod=gun;
! 1673: epsbis=gdivgs(gmulgs(eps,5),4);
! 1674: for (i=1; i<=n; i++)
! 1675: {
! 1676: if (i!=k && gcmp((GEN)d[i],epsbis)==1)
! 1677: {
! 1678: m--;
! 1679: prod=gmul(prod,gsub((GEN)d[i],eps));
! 1680: }
! 1681: }
! 1682: eps2=gdiv(gmul2n(aux,2*m-2),prod);
! 1683: eps2=gpui(eps2,dbltor(1./m),4);
! 1684: rap=gdiv(eps,eps2); eps=eps2;
! 1685: }
! 1686: return eps;
! 1687: }
! 1688:
! 1689: /* round a complex or real number x to an absolute value of 2^(-e) */
! 1690: static GEN
! 1691: mygprec_absolute(GEN x, long bitprec)
! 1692: {
! 1693: long tx=typ(x),e;
! 1694: GEN y;
! 1695:
! 1696: switch(tx)
! 1697: {
! 1698: case t_REAL:
! 1699: e=gexpo(x);
! 1700: if (e<-bitprec || !signe(x)) { y=dbltor(0.); setexpo(y,-bitprec); }
! 1701: else y=mygprec(x,bitprec+e);
! 1702: break;
! 1703: case t_COMPLEX:
! 1704: if (gexpo((GEN)x[2])<-bitprec)
! 1705: y=mygprec_absolute((GEN)x[1],bitprec);
! 1706: else
! 1707: {
! 1708: y=cgetg(3,t_COMPLEX);
! 1709: y[1]=(long) mygprec_absolute((GEN)x[1],bitprec);
! 1710: y[2]=(long) mygprec_absolute((GEN)x[2],bitprec);
! 1711: }
! 1712: break;
! 1713:
! 1714: default: y=mygprec(x,bitprec);
! 1715: }
! 1716: return y;
! 1717: }
! 1718:
! 1719: static long
! 1720: a_posteriori_errors(GEN p, GEN roots_pol, long err)
! 1721: {
! 1722: GEN sigma,overn,shatzle,x;
! 1723: long i,n=lgef(p)-3,e,e_max;
! 1724:
! 1725: sigma=mygprec(dbltor(1.),10);
! 1726: setexpo(sigma,err+(long) log2((double) n)+1);
! 1727: overn=dbltor(1./n);
! 1728: shatzle=gdiv(gpui(sigma,overn,0),
! 1729: gsub(gpui(gsub(gun,sigma),overn,0),
! 1730: gpui(sigma,overn,0)));
! 1731: shatzle=gmul2n(shatzle,1); e_max=-pariINFINITY;
! 1732: for (i=1; i<=n; i++)
! 1733: {
! 1734: x=root_error(n,i,roots_pol,sigma,shatzle);
! 1735: e=gexpo(x); if (e>e_max) e_max=e;
! 1736: roots_pol[i]=(long) mygprec_absolute((GEN)roots_pol[i],-e);
! 1737: }
! 1738: return e_max;
! 1739: }
! 1740:
! 1741: /********************************************************************/
! 1742: /** **/
! 1743: /** MAIN **/
! 1744: /** **/
! 1745: /********************************************************************/
! 1746:
! 1747: /* compute roots in roots_pol so that |P-L_1...L_n|<2^(-bitprec)|P| , and
! 1748: returns prod (x-roots_pol[i]) for i=1..degree(p) */
! 1749: static GEN
! 1750: split_complete(GEN p, long bitprec, GEN *roots_pol, long *iroots)
! 1751: {
! 1752: long n=lgef(p)-3,decprec,ltop,lbot;
! 1753: GEN F,G,a,b,m1,m2,m;
! 1754: GEN *gptr[2];
! 1755:
! 1756: if (n==1)
! 1757: {
! 1758: a=gneg_i(gdiv((GEN)p[2],(GEN)p[3]));
! 1759: (*iroots)++; (*roots_pol)[*iroots]=(long) a;
! 1760: return p;
! 1761: }
! 1762: ltop=avma;
! 1763: if (n==2)
! 1764: {
! 1765: F=gsub(gsqr((GEN)p[3]),gmul2n(gmul((GEN)p[2],(GEN)p[4]),2));
! 1766: decprec=(long) ((double) bitprec * L2SL10)+1;
! 1767: F=gsqrt(F,decprec);
! 1768: a=gneg_i(gdiv(gadd((GEN)p[3],F),gmul2n((GEN)p[4],1)));
! 1769: b=gdiv(gsub(F,(GEN)p[3]),gmul2n((GEN)p[4],1));
! 1770:
! 1771: gptr[0]=&a; gptr[1]=&b;
! 1772: gerepilemany(ltop,gptr,2);
! 1773: (*iroots)++; (*roots_pol)[*iroots]=(long) a;
! 1774: (*iroots)++; (*roots_pol)[*iroots]=(long) b;
! 1775: m=gmul(gsub(polx[varn(p)],mygprec(a,3*bitprec)),
! 1776: gsub(polx[varn(p)],mygprec(b,3*bitprec)));
! 1777: return gmul(m,(GEN)p[4]);
! 1778: }
! 1779: split_0(p,bitprec,&F,&G);
! 1780: m1=split_complete(F,bitprec,roots_pol,iroots);
! 1781: m2=split_complete(G,bitprec,roots_pol,iroots);
! 1782: lbot=avma;
! 1783: m=gmul(m1,m2); *roots_pol=gcopy(*roots_pol);
! 1784: gptr[0]=roots_pol; gptr[1]=&m;
! 1785: gerepilemanysp(ltop,lbot,gptr,2); return m;
! 1786: }
! 1787:
! 1788: /* compute a bound on the maximum modulus of roots of p */
! 1789: static GEN
! 1790: cauchy_bound(GEN p)
! 1791: {
! 1792: long i,n=lgef(p)-3;
! 1793: GEN x=gzero,y,lc;
! 1794:
! 1795: lc=gabs((GEN)p[n+2],DEFAULTPREC); /* leading coefficient */
! 1796: lc=gdiv(dbltor(1.),lc);
! 1797: for (i=0; i<n; i++)
! 1798: {
! 1799: y=gmul(gabs((GEN) p[i+2],DEFAULTPREC),lc);
! 1800: y=gpui(y,dbltor(1./(n-i)),DEFAULTPREC);
! 1801: if (gcmp(y,x) > 0) x=y;
! 1802: }
! 1803: return x;
! 1804: }
! 1805:
! 1806: static GEN
! 1807: mygprecrc_special(GEN x, long bitprec, long e)
! 1808: {
! 1809: long tx=typ(x),lx,ex;
! 1810: GEN y;
! 1811:
! 1812: if (bitprec<=0) bitprec=0; /* should not happen */
! 1813: switch(tx)
! 1814: {
! 1815: case t_REAL:
! 1816: lx=bitprec/BITS_IN_LONG+3;
! 1817: if (lx<lg(x)) lx=lg(x);
! 1818: y=cgetr(lx); affrr(x,y); ex=-bitprec+e;
! 1819: if (!signe(x) && expo(x)>ex) setexpo(y,ex);
! 1820: break;
! 1821: case t_COMPLEX:
! 1822: y=cgetg(3,t_COMPLEX);
! 1823: y[1]=(long) mygprecrc_special((GEN)x[1],bitprec,e);
! 1824: y[2]=(long) mygprecrc_special((GEN)x[2],bitprec,e);
! 1825: break;
! 1826: default: y=gcopy(x);
! 1827: }
! 1828: return y;
! 1829: }
! 1830:
! 1831: /* like mygprec but keep at least the same precision as before */
! 1832: static GEN
! 1833: mygprec_special(GEN x, long bitprec)
! 1834: {
! 1835: long tx=typ(x),lx,i,e;
! 1836: GEN y;
! 1837:
! 1838: switch(tx)
! 1839: {
! 1840: case t_POL:
! 1841: lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1]; e=gexpo(x);
! 1842: for (i=2; i<lx; i++) y[i]=(long) mygprecrc_special((GEN)x[i],bitprec,e);
! 1843: break;
! 1844:
! 1845: default: y=mygprecrc_special(x,bitprec,0);
! 1846: }
! 1847: return y;
! 1848: }
! 1849:
! 1850: static GEN
! 1851: all_roots(GEN p, long bitprec)
! 1852: {
! 1853: GEN q,roots_pol,m;
! 1854: long bitprec2,n=lgef(p)-3,iroots,i,j,e,av,tetpil;
! 1855:
! 1856: roots_pol=cgetg(n+1,t_VEC); av=avma;
! 1857: for (i=1; i<=n; i++) roots_pol[i]=zero;
! 1858: for (j=1; j<=10; j++)
! 1859: {
! 1860: iroots=0; e = 2*gexpo(cauchy_bound(p)); if (e<0) e=0;
! 1861: bitprec2=bitprec+(1<<j)*n+gexpo(p)-gexpo((GEN)p[2+n])
! 1862: +(long) log2(n)+1+e;
! 1863: q=gmul(mygprec(gunr,bitprec2),mygprec(p,bitprec2));
! 1864: m=split_complete(q,bitprec2,&roots_pol,&iroots);
! 1865: e = gexpo(gsub(mygprec_special(p,bitprec2),m))
! 1866: - gexpo((GEN)q[2+n])+(long) log2((double)n)+1;
! 1867: if (e<-2*bitprec2) e=-2*bitprec2; /* to avoid e=-pariINFINITY */
! 1868: if (a_posteriori_errors(q,roots_pol,e) < -bitprec) return roots_pol;
! 1869:
! 1870: tetpil=avma; roots_pol=gerepile(av,tetpil,gcopy(roots_pol));
! 1871: }
! 1872: err(bugparier,"all_roots");
! 1873: return NULL; /* not reached */
! 1874: }
! 1875:
! 1876: /* true if x is an exact scalar, that is integer or rational */
! 1877: static int
! 1878: isexactscalar(GEN x)
! 1879: {
! 1880: long tx=typ(x);
! 1881: return (tx==t_INT || is_frac_t(tx));
! 1882: }
! 1883:
! 1884: static int
! 1885: isexactpol(GEN p)
! 1886: {
! 1887: long i,n=lgef(p)-3;
! 1888:
! 1889: for (i=0; i<=n; i++)
! 1890: if (isexactscalar((GEN)p[i+2])==0) return 0;
! 1891: return 1;
! 1892: }
! 1893:
! 1894: static long
! 1895: isvalidcoeff(GEN x)
! 1896: {
! 1897: long tx=typ(x);
! 1898:
! 1899: switch(tx)
! 1900: {
! 1901: case t_INT: case t_REAL: case t_FRAC: case t_FRACN: return 1;
! 1902: case t_COMPLEX:
! 1903: if (isvalidcoeff((GEN)x[1]) && isvalidcoeff((GEN)x[2])) return 1;
! 1904: }
! 1905: return 0;
! 1906: }
! 1907:
! 1908: static long
! 1909: isvalidpol(GEN p)
! 1910: {
! 1911: long i,n = lgef(p);
! 1912: for (i=2; i<n; i++)
! 1913: if (!isvalidcoeff((GEN)p[i])) return 0;
! 1914: return 1;
! 1915: }
! 1916:
! 1917: static GEN
! 1918: solve_exact_pol(GEN p, long bitprec)
! 1919: {
! 1920: GEN S,ex,factors,roots_pol,roots_fact;
! 1921: long i,j,k,m,n,iroots;
! 1922:
! 1923: n=lgef(p)-3;
! 1924:
! 1925: iroots=0;
! 1926: roots_pol=cgetg(n+1,t_VEC); for (i=1; i<=n; i++) roots_pol[i]=zero;
! 1927:
! 1928: S=square_free_factorization(p);
! 1929: ex=(GEN) S[1]; factors=(GEN) S[2];
! 1930: for (i=1; i<lg(factors); i++)
! 1931: {
! 1932: roots_fact=all_roots((GEN)factors[i],bitprec);
! 1933: n=lgef(factors[i])-3; m=itos((GEN)ex[i]);
! 1934: for (j=1; j<=n; j++)
! 1935: for (k=1; k<=m; k++) roots_pol[++iroots] = roots_fact[j];
! 1936: }
! 1937: return roots_pol;
! 1938: }
! 1939:
! 1940: /* return the roots of p with absolute error bitprec */
! 1941: static GEN
! 1942: roots_com(GEN p, long l)
! 1943: {
! 1944: long bitprec;
! 1945:
! 1946: if (gcmp0(p)) err(zeropoler,"roots");
! 1947: if (typ(p)!=t_POL)
! 1948: {
! 1949: if (!isvalidcoeff(p)) err(typeer,"roots");
! 1950: return cgetg(1,t_VEC); /* constant polynomial */
! 1951: }
! 1952: if (!isvalidpol(p)) err(talker,"invalid coefficients in roots");
! 1953: if (lgef(p) == 3) return cgetg(1,t_VEC); /* constant polynomial */
! 1954: if (l<3) l=3;
! 1955: bitprec=bit_accuracy(l); gunr=realun(l);
! 1956: return isexactpol(p)? solve_exact_pol(p,bitprec): all_roots(p,bitprec);
! 1957: }
! 1958:
! 1959: static GEN
! 1960: tocomplex(GEN x, long l)
! 1961: {
! 1962: GEN y=cgetg(3,t_COMPLEX);
! 1963:
! 1964: y[1]=lgetr(l);
! 1965: if (typ(x) == t_COMPLEX)
! 1966: { y[2]=lgetr(l); gaffect(x,y); }
! 1967: else
! 1968: { gaffect(x,(GEN)y[1]); y[2]=(long)realzero(l); }
! 1969: return y;
! 1970: }
! 1971:
! 1972: /* Check if x is approximately real with precision e */
! 1973: int
! 1974: isrealappr(GEN x, long e)
! 1975: {
! 1976: long tx=typ(x),lx,i;
! 1977: switch(tx)
! 1978: {
! 1979: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 1980: return 1;
! 1981: case t_COMPLEX:
! 1982: return (gexpo((GEN)x[2]) < e);
! 1983: case t_QUAD:
! 1984: err(impl,"isrealappr for type t_QUAD");
! 1985: case t_POL: case t_SER: case t_RFRAC: case t_RFRACN:
! 1986: case t_VEC: case t_COL: case t_MAT:
! 1987: lx = (tx==t_POL)?lgef(x): lg(x);
! 1988: for (i=lontyp[tx]; i<lx; i++)
! 1989: if (! isrealappr((GEN)x[i],e)) return 0;
! 1990: return 1;
! 1991: default: err(typeer,"isrealappr"); return 0;
! 1992: }
! 1993: }
! 1994:
! 1995: /* x,y sont de type t_COMPLEX */
! 1996: static int
! 1997: isconj(GEN x, GEN y, long e)
! 1998: {
! 1999: return (gexpo( gsub((GEN)x[1],(GEN)y[1]) ) < e
! 2000: && gexpo( gadd((GEN)x[2],(GEN)y[2]) ) < e);
! 2001: }
! 2002:
! 2003: /* returns the vector of roots of p, with guaranteed absolute error
! 2004: * BASE^(-(l-2)), where BASE=2^{ BITS_IN_LONG }.
! 2005: */
! 2006: GEN
! 2007: roots(GEN p, long l)
! 2008: {
! 2009: long av,av1,tetpil,n,j,k,s, e = 5 - bit_accuracy(l);
! 2010: GEN p1,p2,p3,p22,res;
! 2011:
! 2012: av=avma; p1=roots_com(p,l); n=lg(p1);
! 2013: if (isrealappr(p,e))
! 2014: {
! 2015: p3=cgetg(n,t_COL); s=0;
! 2016: for (j=1; j<n; j++)
! 2017: {
! 2018: p2=(GEN)p1[j];
! 2019: if (typ(p2) != t_COMPLEX) { p3[++s]=(long)p2; p1[j]=zero; }
! 2020: else if (isrealappr(p2,e)) { p3[++s]=p2[1]; p1[j]=zero; }
! 2021: }
! 2022: setlg(p3,s+1); p2=sort(p3); setlg(p3,n);
! 2023: tetpil=avma; res=cgetg(n,t_COL);
! 2024: for (j=1; j<=s; j++) res[j]=(long)tocomplex((GEN)p2[j],l);
! 2025: for (j=1; j<n; j++)
! 2026: {
! 2027: p2=(GEN)p1[j];
! 2028: if (typ(p2) == t_COMPLEX)
! 2029: {
! 2030: res[++s]=(long)tocomplex(p2,l);
! 2031: av1=avma;
! 2032: for (k=j+1; k<n; k++)
! 2033: {
! 2034: p22=(GEN)p1[k];
! 2035: if (typ(p22) == t_COMPLEX && isconj(p2,p22,e))
! 2036: {
! 2037: avma=av1; res[++s]=(long)tocomplex(p22,l);
! 2038: p1[k]=zero; break;
! 2039: }
! 2040: }
! 2041: if (k==n) err(bugparier,"roots (conjugates)");
! 2042: }
! 2043: }
! 2044: return gerepile(av,tetpil,res);
! 2045: }
! 2046: tetpil=avma; res=cgetg(n,t_COL);
! 2047: for (j=1; j<n; j++) res[j]=(long)tocomplex((GEN)p1[j],l);
! 2048: return gerepile(av,tetpil,res);
! 2049: }
! 2050:
! 2051: GEN
! 2052: roots0(GEN p, long flag,long l)
! 2053: {
! 2054: switch(flag)
! 2055: {
! 2056: case 0: return roots(p,l);
! 2057: case 1: return rootsold(p,l);
! 2058: default: err(flagerr,"polroots");
! 2059: }
! 2060: return NULL; /* not reached */
! 2061: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rootpol.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath