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