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