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