Annotation of OpenXM_contrib/pari/src/basemath/rootpol.c, Revision 1.1.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