1 /* Copyright (C) 2000 The PARI group.
2
3 This file is part of the PARI/GP package.
4
5 PARI/GP is free software; you can redistribute it and/or modify it under the
6 terms of the GNU General Public License as published by the Free Software
7 Foundation; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 /*******************************************************************/
16 /* */
17 /* ROOTS OF COMPLEX POLYNOMIALS */
18 /* (original code contributed by Xavier Gourdon, INRIA RR 1852) */
19 /* */
20 /*******************************************************************/
21 #include "pari.h"
22 #include "paripriv.h"
23
24 static const double pariINFINITY = 1./0.;
25
26 static long
isvalidcoeff(GEN x)27 isvalidcoeff(GEN x)
28 {
29 switch (typ(x))
30 {
31 case t_INT: case t_REAL: case t_FRAC: return 1;
32 case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
33 }
34 return 0;
35 }
36
37 static void
checkvalidpol(GEN p,const char * f)38 checkvalidpol(GEN p, const char *f)
39 {
40 long i,n = lg(p);
41 for (i=2; i<n; i++)
42 if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
43 }
44
45 /********************************************************************/
46 /** **/
47 /** FAST ARITHMETIC over Z[i] **/
48 /** **/
49 /********************************************************************/
50
51 static GEN
ZX_to_ZiX(GEN Pr,GEN Pi)52 ZX_to_ZiX(GEN Pr, GEN Pi)
53 {
54 long i, lr = lg(Pr), li = lg(Pi), l = maxss(lr, li), m = minss(lr, li);
55 GEN P = cgetg(l, t_POL);
56 for(i = 2; i < m; i++)
57 gel(P,i) = signe(gel(Pi,i)) ? mkcomplex(gel(Pr,i), gel(Pi,i))
58 : gel(Pr,i);
59 for( ; i < lr; i++)
60 gel(P,i) = gel(Pr, i);
61 for( ; i < li; i++)
62 gel(P,i) = mkcomplex(gen_0, gel(Pi, i));
63 return P;
64 }
65
66 static GEN
ZiX_sqr(GEN P)67 ZiX_sqr(GEN P)
68 {
69 pari_sp av = avma;
70 GEN Pr2, Pi2, Qr, Qi;
71 GEN Pr = real_i(P), Pi = imag_i(P);
72 if (signe(Pi)==0) return gerepileupto(av, ZX_sqr(Pr));
73 if (signe(Pr)==0) return gerepileupto(av, ZX_neg(ZX_sqr(Pi)));
74 Pr2 = ZX_sqr(Pr); Pi2 = ZX_sqr(Pi);
75 Qr = ZX_sub(Pr2, Pi2);
76 if (degpol(Pr)==degpol(Pi))
77 Qi = ZX_sub(ZX_sqr(ZX_add(Pr, Pi)), ZX_add(Pr2, Pi2));
78 else
79 Qi = ZX_shifti(ZX_mul(Pr, Pi), 1);
80 return gerepilecopy(av, ZX_to_ZiX(Qr, Qi));
81 }
82
83 static GEN
graeffe(GEN p)84 graeffe(GEN p)
85 {
86 pari_sp av = avma;
87 GEN p0, p1, s0, s1;
88 long n = degpol(p);
89
90 if (!n) return RgX_copy(p);
91 RgX_even_odd(p, &p0, &p1);
92 /* p = p0(x^2) + x p1(x^2) */
93 s0 = ZiX_sqr(p0);
94 s1 = ZiX_sqr(p1);
95 return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
96 }
97
98 GEN
ZX_graeffe(GEN p)99 ZX_graeffe(GEN p)
100 {
101 pari_sp av = avma;
102 GEN p0, p1, s0, s1;
103 long n = degpol(p);
104
105 if (!n) return ZX_copy(p);
106 RgX_even_odd(p, &p0, &p1);
107 /* p = p0(x^2) + x p1(x^2) */
108 s0 = ZX_sqr(p0);
109 s1 = ZX_sqr(p1);
110 return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
111 }
112 GEN
polgraeffe(GEN p)113 polgraeffe(GEN p)
114 {
115 pari_sp av = avma;
116 GEN p0, p1, s0, s1;
117 long n = degpol(p);
118
119 if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
120 n = degpol(p);
121 if (!n) return gcopy(p);
122 RgX_even_odd(p, &p0, &p1);
123 /* p = p0(x^2) + x p1(x^2) */
124 s0 = RgX_sqr(p0);
125 s1 = RgX_sqr(p1);
126 return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
127 }
128
129 /********************************************************************/
130 /** **/
131 /** MODULUS OF ROOTS **/
132 /** **/
133 /********************************************************************/
134
135 /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
136 * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
137 * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
138 static double
mydbllog2i(GEN x)139 mydbllog2i(GEN x)
140 {
141 #ifdef LONG_IS_64BIT
142 const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
143 #else
144 const double W = 1/4294967296.; /*2^-32*/
145 #endif
146 GEN m;
147 long lx = lgefint(x);
148 double l;
149 if (lx == 2) return -pariINFINITY;
150 m = int_MSW(x);
151 l = (double)(ulong)*m;
152 if (lx == 3) return log2(l);
153 l += ((double)(ulong)*int_precW(m)) * W;
154 /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
155 * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
156 * exponent BIL(lx-3) causes 1ulp further round-off error */
157 return log2(l) + (double)(BITS_IN_LONG*(lx-3));
158 }
159
160 /* return log(|x|) or -pariINFINITY */
161 static double
mydbllogr(GEN x)162 mydbllogr(GEN x) {
163 if (!signe(x)) return -pariINFINITY;
164 return M_LN2*dbllog2r(x);
165 }
166
167 /* return log2(|x|) or -pariINFINITY */
168 static double
mydbllog2r(GEN x)169 mydbllog2r(GEN x) {
170 if (!signe(x)) return -pariINFINITY;
171 return dbllog2r(x);
172 }
173 double
dbllog2(GEN z)174 dbllog2(GEN z)
175 {
176 double x, y;
177 switch(typ(z))
178 {
179 case t_INT: return mydbllog2i(z);
180 case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
181 case t_REAL: return mydbllog2r(z);
182 default: /*t_COMPLEX*/
183 x = dbllog2(gel(z,1));
184 y = dbllog2(gel(z,2));
185 if (x == -pariINFINITY) return y;
186 if (y == -pariINFINITY) return x;
187 if (fabs(x-y) > 10) return maxdd(x,y);
188 return x + 0.5*log2(1 + exp2(2*(y-x)));
189 }
190 }
191 static GEN /* beware overflow */
dblexp(double x)192 dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
193
194 /* find s such that A_h <= 2^s <= 2 A_i for one h and all i < n = deg(p),
195 * with A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and p = sum p_i X^i */
196 static long
findpower(GEN p)197 findpower(GEN p)
198 {
199 double x, L, mins = pariINFINITY;
200 long n = degpol(p),i;
201
202 L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
203 for (i=n-1; i>=0; i--)
204 {
205 L += log2((double)(i+1) / (double)(n-i));
206 x = dbllog2(gel(p,i+2));
207 if (x != -pariINFINITY)
208 {
209 double s = (L - x) / (double)(n-i);
210 if (s < mins) mins = s;
211 }
212 }
213 i = (long)ceil(mins);
214 if (i - mins > 1 - 1e-12) i--;
215 return i;
216 }
217
218 /* returns the exponent for logmodulus(), from the Newton diagram */
219 static long
newton_polygon(GEN p,long k)220 newton_polygon(GEN p, long k)
221 {
222 pari_sp av = avma;
223 double *logcoef, slope;
224 long n = degpol(p), i, j, h, l, *vertex;
225
226 logcoef = (double*)stack_malloc_align((n+1)*sizeof(double), sizeof(double));
227 vertex = (long*)new_chunk(n+1);
228
229 /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
230 for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
231 vertex[0] = 1; /* sentinel */
232 for (i=0; i < n; i=h)
233 {
234 slope = logcoef[i+1]-logcoef[i];
235 for (j = h = i+1; j<=n; j++)
236 {
237 double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
238 if (slope < pij) { slope = pij; h = j; }
239 }
240 vertex[h] = 1;
241 }
242 h = k; while (!vertex[h]) h++;
243 l = k-1; while (!vertex[l]) l--;
244 set_avma(av);
245 return (long)floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
246 }
247
248 /* change z into z*2^e, where z is real or complex of real */
249 static void
myshiftrc(GEN z,long e)250 myshiftrc(GEN z, long e)
251 {
252 if (typ(z)==t_COMPLEX)
253 {
254 if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
255 if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
256 }
257 else
258 if (signe(z)) shiftr_inplace(z, e);
259 }
260
261 /* return z*2^e, where z is integer or complex of integer (destroy z) */
262 static GEN
myshiftic(GEN z,long e)263 myshiftic(GEN z, long e)
264 {
265 if (typ(z)==t_COMPLEX)
266 {
267 gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
268 gel(z,2) = mpshift(gel(z,2),e);
269 return z;
270 }
271 return signe(z)? mpshift(z,e): gen_0;
272 }
273
274 static GEN
RgX_gtofp_bit(GEN q,long bit)275 RgX_gtofp_bit(GEN q, long bit) { return RgX_gtofp(q, nbits2prec(bit)); }
276
277 static GEN
mygprecrc(GEN x,long prec,long e)278 mygprecrc(GEN x, long prec, long e)
279 {
280 GEN y;
281 switch(typ(x))
282 {
283 case t_REAL:
284 if (!signe(x)) return real_0_bit(e);
285 return realprec(x) == prec? x: rtor(x, prec);
286 case t_COMPLEX:
287 y = cgetg(3,t_COMPLEX);
288 gel(y,1) = mygprecrc(gel(x,1),prec,e);
289 gel(y,2) = mygprecrc(gel(x,2),prec,e);
290 return y;
291 default: return x;
292 }
293 }
294
295 /* gprec behaves badly with the zero for polynomials.
296 The second parameter in mygprec is the precision in base 2 */
297 static GEN
mygprec(GEN x,long bit)298 mygprec(GEN x, long bit)
299 {
300 long lx, i, e, prec;
301 GEN y;
302
303 if (bit < 0) bit = 0; /* should rarely happen */
304 e = gexpo(x) - bit;
305 prec = nbits2prec(bit);
306 switch(typ(x))
307 {
308 case t_POL:
309 y = cgetg_copy(x, &lx); y[1] = x[1];
310 for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
311 break;
312
313 default: y = mygprecrc(x,prec,e);
314 }
315 return y;
316 }
317
318 /* normalize a polynomial p, that is change it with coefficients in Z[i],
319 after making product by 2^shift */
320 static GEN
pol_to_gaussint(GEN p,long shift)321 pol_to_gaussint(GEN p, long shift)
322 {
323 long i, l = lg(p);
324 GEN q = cgetg(l, t_POL); q[1] = p[1];
325 for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
326 return q;
327 }
328
329 /* returns a polynomial q in Z[i][x] keeping bit bits of p */
330 static GEN
eval_rel_pol(GEN p,long bit)331 eval_rel_pol(GEN p, long bit)
332 {
333 long i;
334 for (i = 2; i < lg(p); i++)
335 if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
336 return pol_to_gaussint(p, bit-gexpo(p)+1);
337 }
338
339 /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
340 static GEN
homothetie(GEN p,double lrho,long bit)341 homothetie(GEN p, double lrho, long bit)
342 {
343 GEN q, r, t, iR;
344 long n = degpol(p), i;
345
346 iR = mygprec(dblexp(-lrho),bit);
347 q = mygprec(p, bit);
348 r = cgetg(n+3,t_POL); r[1] = p[1];
349 t = iR; r[n+2] = q[n+2];
350 for (i=n-1; i>0; i--)
351 {
352 gel(r,i+2) = gmul(t, gel(q,i+2));
353 t = mulrr(t, iR);
354 }
355 gel(r,2) = gmul(t, gel(q,2)); return r;
356 }
357
358 /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) [ ~as above with R = 2^-e ]*/
359 static void
homothetie2n(GEN p,long e)360 homothetie2n(GEN p, long e)
361 {
362 if (e)
363 {
364 long i,n = lg(p)-1;
365 for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
366 }
367 }
368
369 /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
370 static void
homothetie_gauss(GEN p,long e,long f)371 homothetie_gauss(GEN p, long e, long f)
372 {
373 if (e || f)
374 {
375 long i, n = lg(p)-1;
376 for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
377 }
378 }
379
380 /* Lower bound on the modulus of the largest root z_0
381 * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
382 static double
lower_bound(GEN p,long * k,double eps)383 lower_bound(GEN p, long *k, double eps)
384 {
385 long n = degpol(p), i, j;
386 pari_sp ltop = avma;
387 GEN a, s, S, ilc;
388 double r, R, rho;
389
390 if (n < 4) { *k = n; return 0.; }
391 S = cgetg(5,t_VEC);
392 a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
393 for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
394 /* i = 1 split out from next loop for efficiency and initialization */
395 s = gel(a,1);
396 gel(S,1) = gneg(s); /* Newton sum S_i */
397 rho = r = gtodouble(gabs(s,3));
398 R = r / n;
399 for (i=2; i<=4; i++)
400 {
401 s = gmulsg(i,gel(a,i));
402 for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
403 gel(S,i) = gneg(s); /* Newton sum S_i */
404 r = gtodouble(gabs(s,3));
405 if (r > 0.)
406 {
407 r = exp(log(r/n) / (double)i);
408 if (r > R) R = r;
409 }
410 }
411 if (R > 0. && eps < 1.2)
412 *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
413 else
414 *k = n;
415 return gc_double(ltop, R);
416 }
417
418 /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
419 * P(0) != 0 and P non constant */
420 static double
logmax_modulus(GEN p,double tau)421 logmax_modulus(GEN p, double tau)
422 {
423 GEN r, q, aux, gunr;
424 pari_sp av, ltop = avma;
425 long i,k,n=degpol(p),nn,bit,M,e;
426 double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
427
428 r = cgeti(BIGDEFAULTPREC);
429 av = avma;
430
431 eps = - 1/log(1.5*tau2); /* > 0 */
432 bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
433 gunr = real_1_bit(bit+2*n);
434 aux = gdiv(gunr, gel(p,2+n));
435 q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
436 e = findpower(q);
437 homothetie2n(q,e);
438 affsi(e, r);
439 q = pol_to_gaussint(q, bit);
440 M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
441 nn = n;
442 for (i=0,e=0;;)
443 { /* nn = deg(q) */
444 rho = lower_bound(q, &k, eps);
445 if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
446 affii(shifti(addis(r,e), 1), r);
447 if (++i == M) break;
448
449 bit = (long) ((double)k * log2(1./tau2) +
450 (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
451 homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
452 nn -= RgX_valrem(q, &q);
453 q = gerepileupto(av, graeffe(q));
454 tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
455 eps = -1/log(tau2); /* > 0 */
456 e = findpower(q);
457 }
458 if (!signe(r)) return gc_double(ltop,0.);
459 r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
460 return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
461 }
462
463 static GEN
RgX_normalize1(GEN x)464 RgX_normalize1(GEN x)
465 {
466 long i, n = lg(x)-1;
467 GEN y;
468 for (i = n; i > 1; i--)
469 if (!gequal0( gel(x,i) )) break;
470 if (i == n) return x;
471 pari_warn(warner,"normalizing a polynomial with 0 leading term");
472 if (i == 1) pari_err_ROOTS0("roots");
473 y = cgetg(i+1, t_POL); y[1] = x[1];
474 for (; i > 1; i--) gel(y,i) = gel(x,i);
475 return y;
476 }
477
478 static GEN
polrootsbound_i(GEN P,double TAU)479 polrootsbound_i(GEN P, double TAU)
480 {
481 pari_sp av = avma;
482 double d;
483 (void)RgX_valrem_inexact(P,&P);
484 P = RgX_normalize1(P);
485 switch(degpol(P))
486 {
487 case -1: pari_err_ROOTS0("roots");
488 case 0: set_avma(av); return gen_0;
489 }
490 d = logmax_modulus(P, TAU) + TAU;
491 /* not dblexp: result differs on ARM emulator */
492 return gerepileuptoleaf(av, mpexp(dbltor(d)));
493 }
494 GEN
polrootsbound(GEN P,GEN tau)495 polrootsbound(GEN P, GEN tau)
496 {
497 if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
498 checkvalidpol(P, "polrootsbound");
499 return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
500 }
501
502 /* log of modulus of the smallest root of p, with relative error tau */
503 static double
logmin_modulus(GEN p,double tau)504 logmin_modulus(GEN p, double tau)
505 {
506 pari_sp av = avma;
507 if (gequal0(gel(p,2))) return -pariINFINITY;
508 return gc_double(av, - logmax_modulus(RgX_recip_shallow(p),tau));
509 }
510
511 /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
512 static double
logmodulus(GEN p,long k,double tau)513 logmodulus(GEN p, long k, double tau)
514 {
515 GEN q;
516 long i, kk = k, imax, n = degpol(p), nn, bit, e;
517 pari_sp av, ltop=avma;
518 double r, tau2 = tau/6;
519
520 bit = (long)(n * (2. + log2(3.*n/tau2)));
521 av = avma;
522 q = gprec_w(p, nbits2prec(bit));
523 q = RgX_gtofp_bit(q, bit);
524 e = newton_polygon(q,k);
525 r = (double)e;
526 homothetie2n(q,e);
527 imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
528 for (i=1; i<imax; i++)
529 {
530 q = eval_rel_pol(q,bit);
531 kk -= RgX_valrem(q, &q);
532 nn = degpol(q);
533
534 q = gerepileupto(av, graeffe(q));
535 e = newton_polygon(q,kk);
536 r += e / exp2((double)i);
537 q = RgX_gtofp_bit(q, bit);
538 homothetie2n(q,e);
539
540 tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
541 bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
542 }
543 return gc_double(ltop, -r * M_LN2);
544 }
545
546 /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
547 * rmin < r_k < rmax. This information helps because we may reduce precision
548 * quicker */
549 static double
logpre_modulus(GEN p,long k,double tau,double lrmin,double lrmax)550 logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
551 {
552 GEN q;
553 long n = degpol(p), i, imax, imax2, bit;
554 pari_sp ltop = avma, av;
555 double lrho, aux, tau2 = tau/6.;
556
557 aux = (lrmax - lrmin) / 2. + 4*tau2;
558 imax = (long) log2(log((double)n)/ aux);
559 if (imax <= 0) return logmodulus(p,k,tau);
560
561 lrho = (lrmin + lrmax) / 2;
562 av = avma;
563 bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
564 q = homothetie(p, lrho, bit);
565 imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
566 if (imax > imax2) imax = imax2;
567
568 for (i=0; i<imax; i++)
569 {
570 q = eval_rel_pol(q,bit);
571 q = gerepileupto(av, graeffe(q));
572 aux = 2*aux + 2*tau2;
573 tau2 *= 1.5;
574 bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
575 q = RgX_gtofp_bit(q, bit);
576 }
577 aux = exp2((double)imax);
578 return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
579 }
580
581 static double
ind_maxlog2(GEN q)582 ind_maxlog2(GEN q)
583 {
584 long i, k = -1;
585 double L = - pariINFINITY;
586 for (i=0; i<=degpol(q); i++)
587 {
588 double d = dbllog2(gel(q,2+i));
589 if (d > L) { L = d; k = i; }
590 }
591 return k;
592 }
593
594 /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
595 * Assume that l <= k <= n-l */
596 static long
dual_modulus(GEN p,double lrho,double tau,long l)597 dual_modulus(GEN p, double lrho, double tau, long l)
598 {
599 long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
600 double tau2 = tau * 7./8.;
601 pari_sp av = avma;
602 GEN q;
603
604 bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
605 q = homothetie(p, lrho, bit);
606 imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
607
608 for (i=0; i<imax; i++)
609 {
610 q = eval_rel_pol(q,bit); v2 = n - degpol(q);
611 v = RgX_valrem(q, &q);
612 ll -= maxss(v, v2); if (ll < 0) ll = 0;
613
614 nn = degpol(q); delta_k += v;
615 if (!nn) return delta_k;
616
617 q = gerepileupto(av, graeffe(q));
618 tau2 *= 7./4.;
619 bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
620 }
621 return gc_long(av, delta_k + (long)ind_maxlog2(q));
622 }
623
624 /********************************************************************/
625 /** **/
626 /** FACTORS THROUGH CIRCLE INTEGRATION **/
627 /** **/
628 /********************************************************************/
629 /* l power of 2, W[step*j] = w_j; set f[j] = p(w_j)
630 * if inv, w_j = exp(2IPi*j/l), else exp(-2IPi*j/l) */
631
632 static void
fft2(GEN W,GEN p,GEN f,long step,long l)633 fft2(GEN W, GEN p, GEN f, long step, long l)
634 {
635 pari_sp av;
636 long i, s1, l1, step2;
637
638 if (l == 2)
639 {
640 gel(f,0) = gadd(gel(p,0), gel(p,step));
641 gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
642 }
643 av = avma;
644 l1 = l>>1; step2 = step<<1;
645 fft2(W,p, f, step2,l1);
646 fft2(W,p+step, f+l1,step2,l1);
647 for (i = s1 = 0; i < l1; i++, s1 += step)
648 {
649 GEN f0 = gel(f,i);
650 GEN f1 = gmul(gel(W,s1), gel(f,i+l1));
651 gel(f,i) = gadd(f0, f1);
652 gel(f,i+l1) = gsub(f0, f1);
653 }
654 gerepilecoeffs(av, f, l);
655 }
656
657 static void
fft(GEN W,GEN p,GEN f,long step,long l,long inv)658 fft(GEN W, GEN p, GEN f, long step, long l, long inv)
659 {
660 pari_sp av;
661 long i, s1, l1, l2, l3, step4;
662 GEN f1, f2, f3, f02;
663
664 if (l == 2)
665 {
666 gel(f,0) = gadd(gel(p,0), gel(p,step));
667 gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
668 }
669 av = avma;
670 if (l == 4)
671 {
672 pari_sp av2;
673 f1 = gadd(gel(p,0), gel(p,step<<1));
674 f2 = gsub(gel(p,0), gel(p,step<<1));
675 f3 = gadd(gel(p,step), gel(p,3*step));
676 f02 = gsub(gel(p,step), gel(p,3*step));
677 f02 = inv? mulcxI(f02): mulcxmI(f02);
678 av2 = avma;
679 gel(f,0) = gadd(f1, f3);
680 gel(f,1) = gadd(f2, f02);
681 gel(f,2) = gsub(f1, f3);
682 gel(f,3) = gsub(f2, f02);
683 gerepileallsp(av,av2,4,&gel(f,0),&gel(f,1),&gel(f,2),&gel(f,3));
684 return;
685 }
686 l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
687 fft(W,p, f, step4,l1,inv);
688 fft(W,p+step, f+l1,step4,l1,inv);
689 fft(W,p+(step<<1),f+l2,step4,l1,inv);
690 fft(W,p+3*step, f+l3,step4,l1,inv);
691 for (i = s1 = 0; i < l1; i++, s1 += step)
692 {
693 long s2 = s1 << 1, s3 = s1 + s2;
694 GEN g02, g13, f13;
695 f1 = gmul(gel(W,s1), gel(f,i+l1));
696 f2 = gmul(gel(W,s2), gel(f,i+l2));
697 f3 = gmul(gel(W,s3), gel(f,i+l3));
698
699 f02 = gadd(gel(f,i),f2);
700 g02 = gsub(gel(f,i),f2);
701 f13 = gadd(f1,f3);
702 g13 = gsub(f1,f3); g13 = inv? mulcxI(g13): mulcxmI(g13);
703
704 gel(f,i) = gadd(f02, f13);
705 gel(f,i+l1) = gadd(g02, g13);
706 gel(f,i+l2) = gsub(f02, f13);
707 gel(f,i+l3) = gsub(g02, g13);
708 }
709 gerepilecoeffs(av, f, l);
710 }
711
712 #define code(t1,t2) ((t1 << 6) | t2)
713
714 static GEN
FFT_i(GEN W,GEN x)715 FFT_i(GEN W, GEN x)
716 {
717 long i, l = lg(W), n = lg(x), tx = typ(x), tw, pa;
718 GEN y, z, p, pol;
719 if ((l-1) & (l-2)) pari_err_DIM("fft");
720 tw = RgV_type(W, &p, &pol, &pa);
721 if (tx == t_POL) { x++; n--; }
722 else if (!is_vec_t(tx)) pari_err_TYPE("fft",x);
723 if (n > l) pari_err_DIM("fft");
724
725 if (n < l) {
726 z = cgetg(l, t_VECSMALL); /* cf stackdummy */
727 for (i = 1; i < n; i++) gel(z,i) = gel(x,i);
728 for ( ; i < l; i++) gel(z,i) = gen_0;
729 }
730 else z = x;
731 y = cgetg(l, t_VEC);
732 if (tw==code(t_COMPLEX,t_INT) || tw==code(t_COMPLEX,t_REAL))
733 {
734 long inv = (l >= 5 && signe(imag_i(gel(W,1+(l>>2))))==1) ? 1 : 0;
735 fft(W+1, z+1, y+1, 1, l-1, inv);
736 } else
737 fft2(W+1, z+1, y+1, 1, l-1);
738 return y;
739 }
740
741 #undef code
742
743 GEN
FFT(GEN W,GEN x)744 FFT(GEN W, GEN x)
745 {
746 if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
747 return FFT_i(W, x);
748 }
749
750 GEN
FFTinv(GEN W,GEN x)751 FFTinv(GEN W, GEN x)
752 {
753 long l = lg(W), i;
754 GEN w;
755 if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
756 w = cgetg(l, t_VECSMALL); /* cf stackdummy */
757 gel(w,1) = gel(W,1); /* w = gconj(W), faster */
758 for (i = 2; i < l; i++) gel(w, i) = gel(W, l-i+1);
759 return FFT_i(w, x);
760 }
761
762 /* returns 1 if p has only real coefficients, 0 else */
763 static int
isreal(GEN p)764 isreal(GEN p)
765 {
766 long i;
767 for (i = lg(p)-1; i > 1; i--)
768 if (typ(gel(p,i)) == t_COMPLEX) return 0;
769 return 1;
770 }
771
772 /* x non complex */
773 static GEN
abs_update_r(GEN x,double * mu)774 abs_update_r(GEN x, double *mu) {
775 GEN y = gtofp(x, DEFAULTPREC);
776 double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
777 setabssign(y); return y;
778 }
779 /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
780 static GEN
abs_update(GEN x,double * mu)781 abs_update(GEN x, double *mu) {
782 GEN y, xr, yr;
783 double ly;
784 if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
785 xr = gel(x,1);
786 yr = gel(x,2);
787 if (gequal0(xr)) return abs_update_r(yr,mu);
788 if (gequal0(yr)) return abs_update_r(xr,mu);
789 /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
790 xr = gtofp(xr, DEFAULTPREC);
791 yr = gtofp(yr, DEFAULTPREC);
792 y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
793 ly = mydbllogr(y); if (ly < *mu) *mu = ly;
794 return y;
795 }
796
797 static void
initdft(GEN * Omega,GEN * prim,long N,long Lmax,long bit)798 initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
799 {
800 long prec = nbits2prec(bit);
801 *Omega = grootsof1(Lmax, prec) + 1;
802 *prim = rootsof1u_cx(N, prec);
803 }
804
805 static void
parameters(GEN p,long * LMAX,double * mu,double * gamma,int polreal,double param,double param2)806 parameters(GEN p, long *LMAX, double *mu, double *gamma,
807 int polreal, double param, double param2)
808 {
809 GEN q, pc, Omega, A, RU, prim, g, TWO;
810 long n = degpol(p), bit, NN, K, i, j, Lmax;
811 pari_sp av2, av = avma;
812
813 bit = gexpo(p) + (long)param2+8;
814 Lmax = 4; while (Lmax <= n) Lmax <<= 1;
815 NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
816 K = NN/Lmax; if (K & 1) K++;
817 NN = Lmax*K;
818 if (polreal) K = K/2+1;
819
820 initdft(&Omega, &prim, NN, Lmax, bit);
821 q = mygprec(p,bit) + 2;
822 A = cgetg(Lmax+1,t_VEC); A++;
823 pc= cgetg(Lmax+1,t_VEC); pc++;
824 for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
825 for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
826
827 *mu = pariINFINITY;
828 g = real_0_bit(-bit);
829 TWO = real2n(1, DEFAULTPREC);
830 av2 = avma;
831 RU = gen_1;
832 for (i=0; i<K; i++)
833 {
834 if (i) {
835 GEN z = RU;
836 for (j=1; j<n; j++)
837 {
838 gel(pc,j) = gmul(gel(q,j),z);
839 z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
840 }
841 gel(pc,n) = gmul(gel(q,n),z);
842 }
843
844 fft(Omega,pc,A,1,Lmax,1);
845 if (polreal && i>0 && i<K-1)
846 for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
847 else
848 for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
849 RU = gmul(RU, prim);
850 if (gc_needed(av,1))
851 {
852 if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
853 gerepileall(av2,2, &g,&RU);
854 }
855 }
856 *gamma = mydbllog2r(divru(g,NN));
857 *LMAX = Lmax; set_avma(av);
858 }
859
860 /* NN is a multiple of Lmax */
861 static void
dft(GEN p,long k,long NN,long Lmax,long bit,GEN F,GEN H,long polreal)862 dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
863 {
864 GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
865 long n = degpol(p), i, j, K;
866 pari_sp ltop;
867
868 initdft(&Omega, &prim, NN, Lmax, bit);
869 RU = cgetg(n+2,t_VEC) + 1;
870
871 K = NN/Lmax; if (polreal) K = K/2+1;
872 q = mygprec(p,bit);
873 qd = RgX_deriv(q);
874
875 A = cgetg(Lmax+1,t_VEC); A++;
876 B = cgetg(Lmax+1,t_VEC); B++;
877 C = cgetg(Lmax+1,t_VEC); C++;
878 pc = cgetg(Lmax+1,t_VEC); pc++;
879 pd = cgetg(Lmax+1,t_VEC); pd++;
880 gel(pc,0) = gel(q,2); for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
881 gel(pd,0) = gel(qd,2); for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
882
883 ltop = avma;
884 W = cgetg(k+1,t_VEC);
885 U = cgetg(k+1,t_VEC);
886 for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
887
888 gel(RU,0) = gen_1;
889 prim2 = gen_1;
890 for (i=0; i<K; i++)
891 {
892 gel(RU,1) = prim2;
893 for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
894 /* RU[j] = prim^(ij)= prim2^j */
895
896 for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
897 fft(Omega,pd,A,1,Lmax,1);
898 for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
899 fft(Omega,pc,B,1,Lmax,1);
900 for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
901 for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
902 fft(Omega,B,A,1,Lmax,1);
903 fft(Omega,C,B,1,Lmax,1);
904
905 if (polreal) /* p has real coefficients */
906 {
907 if (i>0 && i<K-1)
908 {
909 for (j=1; j<=k; j++)
910 {
911 gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
912 gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
913 }
914 }
915 else
916 {
917 for (j=1; j<=k; j++)
918 {
919 gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
920 gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
921 }
922 }
923 }
924 else
925 {
926 for (j=1; j<=k; j++)
927 {
928 gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
929 gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
930 }
931 }
932 prim2 = gmul(prim2,prim);
933 gerepileall(ltop,3, &W,&U,&prim2);
934 }
935
936 for (i=1; i<=k; i++)
937 {
938 aux=gel(W,i);
939 for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
940 gel(F,k+2-i) = gdivgs(aux,-i*NN);
941 }
942 for (i=0; i<k; i++)
943 {
944 aux=gel(U,k-i);
945 for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
946 gel(H,i+2) = gdivgs(aux,NN);
947 }
948 }
949
950 #define NEWTON_MAX 10
951 static GEN
refine_H(GEN F,GEN G,GEN HH,long bit,long Sbit)952 refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
953 {
954 GEN H = HH, D, aux;
955 pari_sp ltop = avma;
956 long error, i, bit1, bit2;
957
958 D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
959 bit2 = bit + Sbit;
960 for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
961 {
962 if (gc_needed(ltop,1))
963 {
964 if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
965 gerepileall(ltop,2, &D,&H);
966 }
967 bit1 = -error + Sbit;
968 aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
969 aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
970
971 bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
972 H = RgX_add(mygprec(H,bit1), aux);
973 D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
974 error = gexpo(D); if (error < -bit1) error = -bit1;
975 }
976 if (error > -bit/2) return NULL; /* FAIL */
977 return gerepilecopy(ltop,H);
978 }
979
980 /* return 0 if fails, 1 else */
981 static long
refine_F(GEN p,GEN * F,GEN * G,GEN H,long bit,double gamma)982 refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
983 {
984 GEN f0, FF, GG, r, HH = H;
985 long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
986 pari_sp av = avma;
987
988 FF = *F; GG = RgX_divrem(p, FF, &r);
989 error = gexpo(r); if (error <= -bit) error = 1-bit;
990 normF = gexpo(FF);
991 normG = gexpo(GG);
992 enh = gexpo(H); if (enh < 0) enh = 0;
993 Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
994 Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
995 bit2 = bit + Sbit;
996 for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
997 {
998 if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
999 if (gc_needed(av,1))
1000 {
1001 if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1002 gerepileall(av,4, &FF,&GG,&r,&HH);
1003 }
1004
1005 bit1 = -error + Sbit2;
1006 HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1007 1-error, Sbit2);
1008 if (!HH) return 0; /* FAIL */
1009
1010 bit1 = -error + Sbit;
1011 r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1012 f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1013
1014 bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1015 FF = gadd(mygprec(FF,bit1),f0);
1016
1017 bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1018 GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1019 error = gexpo(r); if (error < -bit1) error = -bit1;
1020 }
1021 if (error>-bit) return 0; /* FAIL */
1022 *F = FF; *G = GG; return 1;
1023 }
1024
1025 /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1026 where cd is the leading coefficient of p */
1027 static void
split_fromU(GEN p,long k,double delta,long bit,GEN * F,GEN * G,double param,double param2)1028 split_fromU(GEN p, long k, double delta, long bit,
1029 GEN *F, GEN *G, double param, double param2)
1030 {
1031 GEN pp, FF, GG, H;
1032 long n = degpol(p), NN, bit2, Lmax;
1033 int polreal = isreal(p);
1034 pari_sp ltop;
1035 double mu, gamma;
1036
1037 pp = gdiv(p, gel(p,2+n));
1038 parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1039
1040 H = cgetg(k+2,t_POL); H[1] = p[1];
1041 FF = cgetg(k+3,t_POL); FF[1]= p[1];
1042 gel(FF,k+2) = gen_1;
1043
1044 NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1045 NN *= Lmax; ltop = avma;
1046 for(;;)
1047 {
1048 bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
1049 dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1050 if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1051 NN <<= 1; set_avma(ltop);
1052 }
1053 *G = gmul(GG,gel(p,2+n)); *F = FF;
1054 }
1055
1056 static void
optimize_split(GEN p,long k,double delta,long bit,GEN * F,GEN * G,double param,double param2)1057 optimize_split(GEN p, long k, double delta, long bit,
1058 GEN *F, GEN *G, double param, double param2)
1059 {
1060 long n = degpol(p);
1061 GEN FF, GG;
1062
1063 if (k <= n/2)
1064 split_fromU(p,k,delta,bit,F,G,param,param2);
1065 else
1066 {
1067 split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1068 *F = RgX_recip_shallow(GG);
1069 *G = RgX_recip_shallow(FF);
1070 }
1071 }
1072
1073 /********************************************************************/
1074 /** **/
1075 /** SEARCH FOR SEPARATING CIRCLE **/
1076 /** **/
1077 /********************************************************************/
1078
1079 /* return p(2^e*x) *2^(-n*e) */
1080 static void
scalepol2n(GEN p,long e)1081 scalepol2n(GEN p, long e)
1082 {
1083 long i,n=lg(p)-1;
1084 for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1085 }
1086
1087 /* returns p(x/R)*R^n; assume R is at the correct accuracy */
1088 static GEN
scalepol(GEN p,GEN R,long bit)1089 scalepol(GEN p, GEN R, long bit)
1090 { return RgX_rescale(mygprec(p, bit), R); }
1091
1092 /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1093 static GEN
conformal_basecase(GEN p,GEN a)1094 conformal_basecase(GEN p, GEN a)
1095 {
1096 GEN z, r, ma, ca;
1097 long i, n = degpol(p);
1098 pari_sp av;
1099
1100 if (n <= 0) return p;
1101 ma = gneg(a); ca = conj_i(a);
1102 av = avma;
1103 z = deg1pol_shallow(ca, gen_m1, 0);
1104 r = scalarpol_shallow(gel(p,2+n), 0);
1105 for (i=n-1; ; i--)
1106 {
1107 r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
1108 r = gadd(r, gmul(z, gel(p,2+i)));
1109 if (i == 0) return gerepileupto(av, r);
1110 z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1111 if (gc_needed(av,2))
1112 {
1113 if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol (%ld/%ld)",n-i, n);
1114 gerepileall(av,2, &r,&z);
1115 }
1116 }
1117 }
1118 static GEN
conformal_pol(GEN p,GEN a)1119 conformal_pol(GEN p, GEN a)
1120 {
1121 pari_sp av = avma;
1122 long d, nR, n = degpol(p), v;
1123 GEN Q, R, S, T;
1124 if (n < 35) return conformal_basecase(p, a);
1125 d = (n+1) >> 1; v = varn(p);
1126 Q = RgX_shift_shallow(p, -d);
1127 R = RgXn_red_shallow(p, d);
1128 Q = conformal_pol(Q, a);
1129 R = conformal_pol(R, a);
1130 S = gpowgs(deg1pol_shallow(gen_1, gneg(a), v), d);
1131 T = RgX_recip_shallow(S);
1132 if (typ(a) == t_COMPLEX) T = gconj(T);
1133 if (odd(d)) T = RgX_neg(T);
1134 /* S = (X - a)^d, T = (conj(a) X - 1)^d */
1135 nR = n - degpol(R) - d; /* >= 0 */
1136 if (nR) T = RgX_mul(T, gpowgs(deg1pol_shallow(gconj(a), gen_m1, v), nR));
1137 return gerepileupto(av, RgX_add(RgX_mul(Q, S), RgX_mul(R, T)));
1138 }
1139
1140 static const double UNDEF = -100000.;
1141
1142 static double
logradius(double * radii,GEN p,long k,double aux,double * delta)1143 logradius(double *radii, GEN p, long k, double aux, double *delta)
1144 {
1145 long i, n = degpol(p);
1146 double lrho, lrmin, lrmax;
1147 if (k > 1)
1148 {
1149 i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1150 lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1151 }
1152 else /* k=1 */
1153 lrmin = logmin_modulus(p,aux);
1154 radii[k] = lrmin;
1155
1156 if (k+1<n)
1157 {
1158 i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1159 lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1160 }
1161 else /* k+1=n */
1162 lrmax = logmax_modulus(p,aux);
1163 radii[k+1] = lrmax;
1164
1165 lrho = radii[k];
1166 for (i=k-1; i>=1; i--)
1167 {
1168 if (radii[i] == UNDEF || radii[i] > lrho)
1169 radii[i] = lrho;
1170 else
1171 lrho = radii[i];
1172 }
1173 lrho = radii[k+1];
1174 for (i=k+1; i<=n; i++)
1175 {
1176 if (radii[i] == UNDEF || radii[i] < lrho)
1177 radii[i] = lrho;
1178 else
1179 lrho = radii[i];
1180 }
1181 *delta = (lrmax - lrmin) / 2;
1182 if (*delta > 1.) *delta = 1.;
1183 return (lrmin + lrmax) / 2;
1184 }
1185
1186 static void
update_radius(long n,double * radii,double lrho,double * par,double * par2)1187 update_radius(long n, double *radii, double lrho, double *par, double *par2)
1188 {
1189 double t, param = 0., param2 = 0.;
1190 long i;
1191 for (i=1; i<=n; i++)
1192 {
1193 radii[i] -= lrho;
1194 t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1195 param += t; if (t > 1.) param2 += log2(t);
1196 }
1197 *par = param; *par2 = param2;
1198 }
1199
1200 /* apply the conformal mapping then split from U */
1201 static void
conformal_mapping(double * radii,GEN ctr,GEN p,long k,long bit,double aux,GEN * F,GEN * G)1202 conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1203 double aux, GEN *F,GEN *G)
1204 {
1205 long bit2, n = degpol(p), i;
1206 pari_sp ltop = avma, av;
1207 GEN q, FF, GG, a, R;
1208 double lrho, delta, param, param2;
1209 /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1210 bit2 = bit + (long)(n*3.4848775) + 1;
1211 a = sqrtr_abs( utor(3, 2*MEDDEFAULTPREC - 2) );
1212 a = divrs(a, -6);
1213 a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1214
1215 av = avma;
1216 q = conformal_pol(mygprec(p,bit2), a);
1217 for (i=1; i<=n; i++)
1218 if (radii[i] != UNDEF) /* update array radii */
1219 {
1220 pari_sp av2 = avma;
1221 GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1222 /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1223 t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1224 radii[i] = mydbllogr(addsr(1,t)) / 2;
1225 set_avma(av2);
1226 }
1227 lrho = logradius(radii, q,k,aux/10., &delta);
1228 update_radius(n, radii, lrho, ¶m, ¶m2);
1229
1230 bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
1231 R = mygprec(dblexp(-lrho), bit2);
1232 q = scalepol(q,R,bit2);
1233 gerepileall(av,2, &q,&R);
1234
1235 optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1236 bit2 += n; R = invr(R);
1237 FF = scalepol(FF,R,bit2);
1238 GG = scalepol(GG,R,bit2);
1239
1240 a = mygprec(a,bit2);
1241 FF = conformal_pol(FF,a);
1242 GG = conformal_pol(GG,a);
1243
1244 a = invr(subsr(1, gnorm(a)));
1245 FF = RgX_Rg_mul(FF, powru(a,k));
1246 GG = RgX_Rg_mul(GG, powru(a,n-k));
1247
1248 *F = mygprec(FF,bit+n);
1249 *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1250 }
1251
1252 /* split p, this time without scaling. returns in F and G two polynomials
1253 * such that |p-FG|< 2^(-bit)|p| */
1254 static void
split_2(GEN p,long bit,GEN ctr,double thickness,GEN * F,GEN * G)1255 split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1256 {
1257 GEN q, FF, GG, R;
1258 double aux, delta, param, param2;
1259 long n = degpol(p), i, j, k, bit2;
1260 double lrmin, lrmax, lrho, *radii;
1261
1262 radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
1263
1264 for (i=2; i<n; i++) radii[i] = UNDEF;
1265 aux = thickness/(double)(4 * n);
1266 lrmin = logmin_modulus(p, aux);
1267 lrmax = logmax_modulus(p, aux);
1268 radii[1] = lrmin;
1269 radii[n] = lrmax;
1270 i = 1; j = n;
1271 lrho = (lrmin + lrmax) / 2;
1272 k = dual_modulus(p, lrho, aux, 1);
1273 if (5*k < n || (n < 2*k && 5*k < 4*n))
1274 { lrmax = lrho; j=k+1; radii[j] = lrho; }
1275 else
1276 { lrmin = lrho; i=k; radii[i] = lrho; }
1277 while (j > i+1)
1278 {
1279 if (i+j == n+1)
1280 lrho = (lrmin + lrmax) / 2;
1281 else
1282 {
1283 double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1284 if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1285 else lrho = lrmin * kappa + lrmax;
1286 lrho /= 1+kappa;
1287 }
1288 aux = (lrmax - lrmin) / (4*(j-i));
1289 k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1290 if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
1291 { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1292 else
1293 { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1294 }
1295 aux = lrmax - lrmin;
1296
1297 if (ctr)
1298 {
1299 lrho = (lrmax + lrmin) / 2;
1300 for (i=1; i<=n; i++)
1301 if (radii[i] != UNDEF) radii[i] -= lrho;
1302
1303 bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
1304 R = mygprec(dblexp(-lrho), bit2);
1305 q = scalepol(p,R,bit2);
1306 conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1307 }
1308 else
1309 {
1310 lrho = logradius(radii, p, k, aux/10., &delta);
1311 update_radius(n, radii, lrho, ¶m, ¶m2);
1312
1313 bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
1314 R = mygprec(dblexp(-lrho), bit2);
1315 q = scalepol(p,R,bit2);
1316 optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1317 }
1318 bit += n;
1319 bit2 += n; R = invr(mygprec(R,bit2));
1320 *F = mygprec(scalepol(FF,R,bit2), bit);
1321 *G = mygprec(scalepol(GG,R,bit2), bit);
1322 }
1323
1324 /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1325 /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1326 * where the maximum modulus of the roots of p is <=1.
1327 * Assume sum of roots is 0. */
1328 static void
split_1(GEN p,long bit,GEN * F,GEN * G)1329 split_1(GEN p, long bit, GEN *F, GEN *G)
1330 {
1331 long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1332 GEN ctr, q, qq, FF, GG, v, gr, r, newq;
1333 double lrmin, lrmax, lthick;
1334 const double LOG3 = 1.098613;
1335
1336 lrmax = logmax_modulus(p, 0.01);
1337 gr = mygprec(dblexp(-lrmax), bit2);
1338 q = scalepol(p,gr,bit2);
1339
1340 bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1341 v = cgetg(5,t_VEC);
1342 gel(v,1) = gen_2;
1343 gel(v,2) = gen_m2;
1344 gel(v,3) = mkcomplex(gen_0, gel(v,1));
1345 gel(v,4) = mkcomplex(gen_0, gel(v,2));
1346 q = mygprec(q,bit2); lthick = 0;
1347 newq = ctr = NULL; /* -Wall */
1348 imax = polreal? 3: 4;
1349 for (i=1; i<=imax; i++)
1350 {
1351 qq = RgX_translate(q, gel(v,i));
1352 lrmin = logmin_modulus(qq,0.05);
1353 if (LOG3 > lrmin + lthick)
1354 {
1355 double lquo = logmax_modulus(qq,0.05) - lrmin;
1356 if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1357 }
1358 if (lthick > M_LN2) break;
1359 if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
1360 }
1361 bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
1362 split_2(newq, bit2, ctr, lthick, &FF, &GG);
1363 r = gneg(mygprec(ctr,bit2));
1364 FF = RgX_translate(FF,r);
1365 GG = RgX_translate(GG,r);
1366
1367 gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1368 *F = scalepol(FF,gr,bit2);
1369 *G = scalepol(GG,gr,bit2);
1370 }
1371
1372 /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1373 where the maximum modulus of the roots of p is < 0.5 */
1374 static int
split_0_2(GEN p,long bit,GEN * F,GEN * G)1375 split_0_2(GEN p, long bit, GEN *F, GEN *G)
1376 {
1377 GEN q, b;
1378 long n = degpol(p), k, bit2, eq;
1379 double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
1380 double aux1 = dbllog2(gel(p,n+1)), aux;
1381
1382 if (aux1 == -pariINFINITY) /* p1 = 0 */
1383 aux = 0;
1384 else
1385 {
1386 aux = aux1 - aux0; /* log2(p1/p0) */
1387 /* beware double overflow */
1388 if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
1389 aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1390 }
1391 bit2 = bit+1 + (long)(log2((double)n) + aux);
1392 q = mygprec(p,bit2);
1393 if (aux1 == -pariINFINITY) b = NULL;
1394 else
1395 {
1396 b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1397 q = RgX_translate(q,b);
1398 }
1399 gel(q,n+1) = gen_0; eq = gexpo(q);
1400 k = 0;
1401 while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1402 || gequal0(gel(q,k+2)))) k++;
1403 if (k > 0)
1404 {
1405 if (k > n/2) k = n/2;
1406 bit2 += k<<1;
1407 *F = pol_xn(k, 0);
1408 *G = RgX_shift_shallow(q, -k);
1409 }
1410 else
1411 {
1412 split_1(q,bit2,F,G);
1413 bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
1414 *F = mygprec(*F,bit2);
1415 }
1416 *G = mygprec(*G,bit2);
1417 if (b)
1418 {
1419 GEN mb = mygprec(gneg(b), bit2);
1420 *F = RgX_translate(*F, mb);
1421 *G = RgX_translate(*G, mb);
1422 }
1423 return 1;
1424 }
1425
1426 /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1427 * Assume max_modulus(p) < 2 */
1428 static void
split_0_1(GEN p,long bit,GEN * F,GEN * G)1429 split_0_1(GEN p, long bit, GEN *F, GEN *G)
1430 {
1431 GEN FF, GG;
1432 long n, bit2, normp;
1433
1434 if (split_0_2(p,bit,F,G)) return;
1435
1436 normp = gexpo(p);
1437 scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1438 n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1439 split_1(mygprec(p,bit2), bit2,&FF,&GG);
1440 scalepol2n(FF,-2);
1441 scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1442 *F = mygprec(FF,bit2);
1443 *G = mygprec(GG,bit2);
1444 }
1445
1446 /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1447 static void
split_0(GEN p,long bit,GEN * F,GEN * G)1448 split_0(GEN p, long bit, GEN *F, GEN *G)
1449 {
1450 const double LOG1_9 = 0.6418539;
1451 long n = degpol(p), k = 0;
1452 GEN q;
1453
1454 while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1455 if (k > 0)
1456 {
1457 if (k > n/2) k = n/2;
1458 *F = pol_xn(k, 0);
1459 *G = RgX_shift_shallow(p, -k);
1460 }
1461 else
1462 {
1463 double lr = logmax_modulus(p, 0.05);
1464 if (lr < LOG1_9) split_0_1(p, bit, F, G);
1465 else
1466 {
1467 q = RgX_recip_shallow(p);
1468 lr = logmax_modulus(q,0.05);
1469 if (lr < LOG1_9)
1470 {
1471 split_0_1(q, bit, F, G);
1472 *F = RgX_recip_shallow(*F);
1473 *G = RgX_recip_shallow(*G);
1474 }
1475 else
1476 split_2(p,bit,NULL, 1.2837,F,G);
1477 }
1478 }
1479 }
1480
1481 /********************************************************************/
1482 /** **/
1483 /** ERROR ESTIMATE FOR THE ROOTS **/
1484 /** **/
1485 /********************************************************************/
1486
1487 static GEN
root_error(long n,long k,GEN roots_pol,long err,GEN shatzle)1488 root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1489 {
1490 GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1491 long i, j;
1492
1493 d = cgetg(n+1,t_VEC);
1494 for (i=1; i<=n; i++)
1495 {
1496 if (i!=k)
1497 {
1498 aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1499 gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1500 }
1501 }
1502 rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1503 if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1504 eps = mulrr(rho, shatzle);
1505 aux = shiftr(powru(rho,n), err);
1506
1507 for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
1508 {
1509 GEN prod = NULL; /* 1. */
1510 long m = n;
1511 epsbis = mulrr(eps, dbltor(1.25));
1512 for (i=1; i<=n; i++)
1513 {
1514 if (i != k && cmprr(gel(d,i),epsbis) > 0)
1515 {
1516 GEN dif = subrr(gel(d,i),eps);
1517 prod = prod? mulrr(prod, dif): dif;
1518 m--;
1519 }
1520 }
1521 eps2 = prod? divrr(aux, prod): aux;
1522 if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1523 rap = divrr(eps,eps2); eps = eps2;
1524 }
1525 return eps;
1526 }
1527
1528 /* round a complex or real number x to an absolute value of 2^(-bit) */
1529 static GEN
mygprec_absolute(GEN x,long bit)1530 mygprec_absolute(GEN x, long bit)
1531 {
1532 long e;
1533 GEN y;
1534
1535 switch(typ(x))
1536 {
1537 case t_REAL:
1538 e = expo(x) + bit;
1539 return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1540 case t_COMPLEX:
1541 if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1542 y = cgetg(3,t_COMPLEX);
1543 gel(y,1) = mygprec_absolute(gel(x,1),bit);
1544 gel(y,2) = mygprec_absolute(gel(x,2),bit);
1545 return y;
1546 default: return x;
1547 }
1548 }
1549
1550 static long
a_posteriori_errors(GEN p,GEN roots_pol,long err)1551 a_posteriori_errors(GEN p, GEN roots_pol, long err)
1552 {
1553 long i, n = degpol(p), e_max = -(long)EXPOBITS;
1554 GEN sigma, shatzle;
1555
1556 err += (long)log2((double)n) + 1;
1557 if (err > -2) return 0;
1558 sigma = real2n(-err, LOWDEFAULTPREC);
1559 /* 2 / ((s - 1)^(1/n) - 1) */
1560 shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1561 for (i=1; i<=n; i++)
1562 {
1563 pari_sp av = avma;
1564 GEN x = root_error(n,i,roots_pol,err,shatzle);
1565 long e = gexpo(x);
1566 set_avma(av); if (e > e_max) e_max = e;
1567 gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1568 }
1569 return e_max;
1570 }
1571
1572 /********************************************************************/
1573 /** **/
1574 /** MAIN **/
1575 /** **/
1576 /********************************************************************/
1577 static GEN
append_clone(GEN r,GEN a)1578 append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1579
1580 /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1581 * returns prod (x-roots_pol[i]) */
1582 static GEN
split_complete(GEN p,long bit,GEN roots_pol)1583 split_complete(GEN p, long bit, GEN roots_pol)
1584 {
1585 long n = degpol(p);
1586 pari_sp ltop;
1587 GEN p1, F, G, a, b, m1, m2;
1588
1589 if (n == 1)
1590 {
1591 a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1592 (void)append_clone(roots_pol,a); return p;
1593 }
1594 ltop = avma;
1595 if (n == 2)
1596 {
1597 F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1598 F = gsqrt(F, nbits2prec(bit));
1599 p1 = ginv(gmul2n(gel(p,4),1));
1600 a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1601 b = gmul(gsub(F,gel(p,3)), p1);
1602 a = append_clone(roots_pol,a);
1603 b = append_clone(roots_pol,b); set_avma(ltop);
1604 a = mygprec(a, 3*bit);
1605 b = mygprec(b, 3*bit);
1606 return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1607 }
1608 split_0(p,bit,&F,&G);
1609 m1 = split_complete(F,bit,roots_pol);
1610 m2 = split_complete(G,bit,roots_pol);
1611 return gerepileupto(ltop, gmul(m1,m2));
1612 }
1613
1614 static GEN
quicktofp(GEN x)1615 quicktofp(GEN x)
1616 {
1617 const long prec = DEFAULTPREC;
1618 switch(typ(x))
1619 {
1620 case t_INT: return itor(x, prec);
1621 case t_REAL: return rtor(x, prec);
1622 case t_FRAC: return fractor(x, prec);
1623 case t_COMPLEX: {
1624 GEN a = gel(x,1), b = gel(x,2);
1625 /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1626 if (isintzero(a)) return cxcompotor(b, prec);
1627 if (isintzero(b)) return cxcompotor(a, prec);
1628 a = cxcompotor(a, prec);
1629 b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1630 }
1631 default: pari_err_TYPE("quicktofp",x);
1632 return NULL;/*LCOV_EXCL_LINE*/
1633 }
1634
1635 }
1636
1637 /* bound log_2 |largest root of p| (Fujiwara's bound) */
1638 double
fujiwara_bound(GEN p)1639 fujiwara_bound(GEN p)
1640 {
1641 pari_sp av = avma;
1642 long i, n = degpol(p);
1643 GEN cc;
1644 double loglc, Lmax;
1645
1646 if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1647 loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
1648 cc = gel(p, 2);
1649 if (gequal0(cc))
1650 Lmax = -pariINFINITY-1;
1651 else
1652 Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
1653 for (i = 1; i < n; i++)
1654 {
1655 GEN y = gel(p,i+2);
1656 double L;
1657 if (gequal0(y)) continue;
1658 L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
1659 if (L > Lmax) Lmax = L;
1660 }
1661 return gc_double(av, Lmax+1);
1662 }
1663
1664 /* Fujiwara's bound, real roots. Based on the following remark: if
1665 * p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
1666 * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
1667 * of q is a bound for the positive roots of p. */
1668 double
fujiwara_bound_real(GEN p,long sign)1669 fujiwara_bound_real(GEN p, long sign)
1670 {
1671 pari_sp av = avma;
1672 GEN x;
1673 long n = degpol(p), i, signodd, signeven;
1674 if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1675 x = shallowcopy(p);
1676 if (gsigne(gel(x, n+2)) > 0)
1677 { signeven = 1; signodd = sign; }
1678 else
1679 { signeven = -1; signodd = -sign; }
1680 for (i = 0; i < n; i++)
1681 {
1682 if ((n - i) % 2)
1683 { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
1684 else
1685 { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
1686 }
1687 return gc_double(av, fujiwara_bound(x));
1688 }
1689
1690 static GEN
mygprecrc_special(GEN x,long prec,long e)1691 mygprecrc_special(GEN x, long prec, long e)
1692 {
1693 GEN y;
1694 switch(typ(x))
1695 {
1696 case t_REAL:
1697 if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1698 return (prec > realprec(x))? rtor(x, prec): x;
1699 case t_COMPLEX:
1700 y = cgetg(3,t_COMPLEX);
1701 gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1702 gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1703 return y;
1704 default: return x;
1705 }
1706 }
1707
1708 /* like mygprec but keep at least the same precision as before */
1709 static GEN
mygprec_special(GEN x,long bit)1710 mygprec_special(GEN x, long bit)
1711 {
1712 long lx, i, e, prec;
1713 GEN y;
1714
1715 if (bit < 0) bit = 0; /* should not happen */
1716 e = gexpo(x) - bit;
1717 prec = nbits2prec(bit);
1718 switch(typ(x))
1719 {
1720 case t_POL:
1721 y = cgetg_copy(x, &lx); y[1] = x[1];
1722 for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1723 break;
1724
1725 default: y = mygprecrc_special(x,prec,e);
1726 }
1727 return y;
1728 }
1729
1730 static GEN
fix_roots1(GEN r)1731 fix_roots1(GEN r)
1732 {
1733 long i, l = lg(r);
1734 GEN allr = cgetg(l, t_VEC);
1735 for (i=1; i<l; i++)
1736 {
1737 GEN t = gel(r,i);
1738 gel(allr,i) = gcopy(t);
1739 gunclone(t);
1740 }
1741 return allr;
1742 }
1743 static GEN
fix_roots(GEN r,long h,long bit)1744 fix_roots(GEN r, long h, long bit)
1745 {
1746 long i, j, k, l, prec;
1747 GEN allr, ro1;
1748 if (h == 1) return fix_roots1(r);
1749 prec = nbits2prec(bit);
1750 ro1 = grootsof1(h, prec) + 1;
1751 l = lg(r)-1;
1752 allr = cgetg(h*l+1, t_VEC);
1753 for (k=1,i=1; i<=l; i++)
1754 {
1755 GEN p2, p1 = gel(r,i);
1756 p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1757 for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1758 gunclone(p1);
1759 }
1760 return allr;
1761 }
1762
1763 static GEN
all_roots(GEN p,long bit)1764 all_roots(GEN p, long bit)
1765 {
1766 GEN pd, q, roots_pol, m;
1767 long bit2, i, e, h, elc, n = degpol(p);
1768 double fb;
1769 pari_sp av;
1770
1771 pd = RgX_deflate_max(p, &h); elc = gexpo(leading_coeff(pd));
1772 fb = fujiwara_bound(pd);
1773 e = (fb < 0)? 0: (long)(2 * fb);
1774 bit2 = bit + maxss(0, gexpo(pd)) + (long)log2(n/h)+1+e;
1775 e = 0;
1776 for (av=avma,i=1;; i++,set_avma(av))
1777 {
1778 roots_pol = vectrunc_init(n+1);
1779 bit2 += e + (n << i);
1780 q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1781 q[1] = evalsigne(1)|evalvarn(0);
1782 m = split_complete(q,bit2,roots_pol);
1783 roots_pol = fix_roots(roots_pol, h, bit2);
1784 q = mygprec_special(pd,bit2);
1785 q[1] = evalsigne(1)|evalvarn(0);
1786 e = gexpo(RgX_sub(q, m)) - elc + (long)log2((double)n) + 1;
1787 if (e < -2*bit2) e = -2*bit2; /* avoid e = -oo */
1788 if (e < 0)
1789 {
1790 e = bit + a_posteriori_errors(p,roots_pol,e);
1791 if (e < 0) return roots_pol;
1792 }
1793 if (DEBUGLEVEL > 7)
1794 err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1795 }
1796 }
1797
1798 INLINE int
isexactscalar(GEN x)1799 isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1800
1801 static int
isexactpol(GEN p)1802 isexactpol(GEN p)
1803 {
1804 long i,n = degpol(p);
1805 for (i=0; i<=n; i++)
1806 if (!isexactscalar(gel(p,i+2))) return 0;
1807 return 1;
1808 }
1809
1810 static GEN
solve_exact_pol(GEN p,long bit)1811 solve_exact_pol(GEN p, long bit)
1812 {
1813 long i, j, k, m, n = degpol(p), iroots = 0;
1814 GEN ex, factors, v = zerovec(n);
1815
1816 factors = ZX_squff(Q_primpart(p), &ex);
1817 for (i=1; i<lg(factors); i++)
1818 {
1819 GEN roots_fact = all_roots(gel(factors,i), bit);
1820 n = degpol(gel(factors,i));
1821 m = ex[i];
1822 for (j=1; j<=n; j++)
1823 for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
1824 }
1825 return v;
1826 }
1827
1828 /* return the roots of p with absolute error bit */
1829 static GEN
roots_com(GEN q,long bit)1830 roots_com(GEN q, long bit)
1831 {
1832 GEN L, p;
1833 long v = RgX_valrem_inexact(q, &p);
1834 int ex = isexactpol(p);
1835 if (!ex) p = RgX_normalize1(p);
1836 if (lg(p) == 3)
1837 L = cgetg(1,t_VEC); /* constant polynomial */
1838 else
1839 L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
1840 if (v)
1841 {
1842 GEN M, z, t = gel(q,2);
1843 long i, x, y, l, n;
1844
1845 if (isrationalzero(t)) x = -bit;
1846 else
1847 {
1848 n = gexpo(t);
1849 x = n / v; l = degpol(q);
1850 for (i = v; i <= l; i++)
1851 {
1852 t = gel(q,i+2);
1853 if (isrationalzero(t)) continue;
1854 y = (n - gexpo(t)) / i;
1855 if (y < x) x = y;
1856 }
1857 }
1858 z = real_0_bit(x); l = v + lg(L);
1859 M = cgetg(l, t_VEC); L -= v;
1860 for (i = 1; i <= v; i++) gel(M,i) = z;
1861 for ( ; i < l; i++) gel(M,i) = gel(L,i);
1862 L = M;
1863 }
1864 return L;
1865 }
1866
1867 static GEN
tocomplex(GEN x,long l,long bit)1868 tocomplex(GEN x, long l, long bit)
1869 {
1870 GEN y;
1871 if (typ(x) == t_COMPLEX)
1872 {
1873 if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
1874 x = gel(x,2);
1875 y = cgetg(3,t_COMPLEX);
1876 gel(y,1) = real_0_bit(-bit);
1877 gel(y,2) = mygprecrc(x, l, -bit);
1878 }
1879 else
1880 {
1881 y = cgetg(3,t_COMPLEX);
1882 gel(y,1) = mygprecrc(x, l, -bit);
1883 gel(y,2) = real_0_bit(-bit);
1884 }
1885 return y;
1886 }
1887
1888 /* x,y are t_COMPLEX of t_REALs or t_REAL, compare wrt |Im x| - |Im y|,
1889 * then Re x - Re y, up to 2^-e absolute error */
1890 static int
cmp_complex_appr(void * E,GEN x,GEN y)1891 cmp_complex_appr(void *E, GEN x, GEN y)
1892 {
1893 long e = (long)E;
1894 GEN z, xi, yi, xr, yr;
1895 long sz, sxi, syi;
1896 if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
1897 else { xr = x; xi = NULL; sxi = 0; }
1898 if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
1899 else { yr = y; yi = NULL; syi = 0; }
1900 /* Compare absolute values of imaginary parts */
1901 if (!sxi)
1902 {
1903 if (syi && expo(yi) >= e) return -1;
1904 /* |Im x| ~ |Im y| ~ 0 */
1905 }
1906 else if (!syi)
1907 {
1908 if (sxi && expo(xi) >= e) return 1;
1909 /* |Im x| ~ |Im y| ~ 0 */
1910 }
1911 else
1912 {
1913 z = addrr_sign(xi, 1, yi, -1); sz = signe(z);
1914 if (sz && expo(z) >= e) return (int)sz;
1915 }
1916 /* |Im x| ~ |Im y|, sort according to real parts */
1917 z = subrr(xr, yr); sz = signe(z);
1918 if (sz && expo(z) >= e) return (int)sz;
1919 /* Re x ~ Re y. Place negative imaginary part before positive */
1920 return (int) (sxi - syi);
1921 }
1922
1923 static GEN
clean_roots(GEN L,long l,long bit,long clean)1924 clean_roots(GEN L, long l, long bit, long clean)
1925 {
1926 long i, n = lg(L), ex = 5 - bit;
1927 GEN res = cgetg(n,t_COL);
1928 for (i=1; i<n; i++)
1929 {
1930 GEN c = gel(L,i);
1931 if (clean && isrealappr(c,ex))
1932 {
1933 if (typ(c) == t_COMPLEX) c = gel(c,1);
1934 c = mygprecrc(c, l, -bit);
1935 }
1936 else
1937 c = tocomplex(c, l, bit);
1938 gel(res,i) = c;
1939 }
1940 gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
1941 return res;
1942 }
1943
1944 /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
1945 static GEN
roots_aux(GEN p,long l,long clean)1946 roots_aux(GEN p, long l, long clean)
1947 {
1948 pari_sp av = avma;
1949 long bit;
1950 GEN L;
1951
1952 if (typ(p) != t_POL)
1953 {
1954 if (gequal0(p)) pari_err_ROOTS0("roots");
1955 if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
1956 return cgetg(1,t_COL); /* constant polynomial */
1957 }
1958 if (!signe(p)) pari_err_ROOTS0("roots");
1959 checkvalidpol(p,"roots");
1960 if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
1961 if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
1962 bit = prec2nbits(l);
1963 L = roots_com(p, bit);
1964 return gerepilecopy(av, clean_roots(L, l, bit, clean));
1965 }
1966 GEN
roots(GEN p,long l)1967 roots(GEN p, long l) { return roots_aux(p,l, 0); }
1968 /* clean up roots. If root is real replace it by its real part */
1969 GEN
cleanroots(GEN p,long l)1970 cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
1971
1972 /* private variant of conjvec. Allow non rational coefficients, shallow
1973 * function. */
1974 GEN
polmod_to_embed(GEN x,long prec)1975 polmod_to_embed(GEN x, long prec)
1976 {
1977 GEN v, T = gel(x,1), A = gel(x,2);
1978 long i, l;
1979 if (typ(A) != t_POL || varn(A) != varn(T))
1980 {
1981 checkvalidpol(T,"polmod_to_embed");
1982 return const_col(degpol(T), A);
1983 }
1984 v = cleanroots(T,prec); l = lg(v);
1985 for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
1986 return v;
1987 }
1988
1989 GEN
QX_complex_roots(GEN p,long l)1990 QX_complex_roots(GEN p, long l)
1991 {
1992 pari_sp av = avma;
1993 long bit, v;
1994 GEN L;
1995
1996 if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
1997 if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
1998 if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
1999 bit = prec2nbits(l);
2000 v = RgX_valrem(p, &p);
2001 L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
2002 if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
2003 return gerepilecopy(av, clean_roots(L, l, bit, 1));
2004 }
2005
2006 /********************************************************************/
2007 /** **/
2008 /** REAL ROOTS OF INTEGER POLYNOMIAL **/
2009 /** **/
2010 /********************************************************************/
2011
2012 /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1)), P
2013 * has no rational root. The inversion is implicit (we take coefficients
2014 * backwards). */
2015 static long
X2XP1(GEN P,GEN * Premapped)2016 X2XP1(GEN P, GEN *Premapped)
2017 {
2018 const pari_sp av = avma;
2019 GEN v = shallowcopy(P);
2020 long i, j, nb, s, dP = degpol(P), vlim = dP+2;
2021
2022 for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2023 s = -signe(gel(v, vlim));
2024 vlim--; nb = 0;
2025 for (i = 1; i < dP; i++)
2026 {
2027 long s2 = -signe(gel(v, 2));
2028 int flag = (s2 == s);
2029 for (j = 2; j < vlim; j++)
2030 {
2031 gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2032 if (flag) flag = (s2 != signe(gel(v, j+1)));
2033 }
2034 if (s == signe(gel(v, vlim)))
2035 {
2036 if (++nb >= 2) return gc_long(av,2);
2037 s = -s;
2038 }
2039 /* if flag is set there will be no further sign changes */
2040 if (flag && (!Premapped || !nb)) return gc_long(av, nb);
2041 vlim--;
2042 if (gc_needed(av, 3))
2043 {
2044 if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, dP-1);
2045 if (!Premapped) setlg(v, vlim + 2);
2046 v = gerepilecopy(av, v);
2047 }
2048 }
2049 if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
2050 if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
2051 return nb;
2052 }
2053
2054 static long
_intervalcmp(GEN x,GEN y)2055 _intervalcmp(GEN x, GEN y)
2056 {
2057 if (typ(x) == t_VEC) x = gel(x, 1);
2058 if (typ(y) == t_VEC) y = gel(y, 1);
2059 return gcmp(x, y);
2060 }
2061
2062 static GEN
_gen_nored(void * E,GEN x)2063 _gen_nored(void *E, GEN x) { (void)E; return x; }
2064 static GEN
_mp_add(void * E,GEN x,GEN y)2065 _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2066 static GEN
_mp_sub(void * E,GEN x,GEN y)2067 _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
2068 static GEN
_mp_mul(void * E,GEN x,GEN y)2069 _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2070 static GEN
_mp_sqr(void * E,GEN x)2071 _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2072 static GEN
_gen_one(void * E)2073 _gen_one(void *E) { (void)E; return gen_1; }
2074 static GEN
_gen_zero(void * E)2075 _gen_zero(void *E) { (void)E; return gen_0; }
2076
2077 static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
2078 _mp_mul, _mp_sqr, _gen_one, _gen_zero };
2079
2080 static GEN
_mp_cmul(void * E,GEN P,long a,GEN x)2081 _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
2082
2083 /* Split the polynom P in two parts, whose coeffs have constant sign:
2084 * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
2085 * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
2086 * Pprimep[i] = (i+D) Pp[i]. Return D */
2087 static long
split_pols(GEN P,GEN * pPp,GEN * pPm,GEN * pPprimep,GEN * pPprimem)2088 split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
2089 {
2090 long i, D, dP = degpol(P), s0 = signe(gel(P,2));
2091 GEN Pp, Pm, Pprimep, Pprimem;
2092 for(i=1; i <= dP; i++)
2093 if (signe(gel(P, i+2)) == -s0) break;
2094 D = i;
2095 Pm = cgetg(D + 2, t_POL);
2096 Pprimem = cgetg(D + 1, t_POL);
2097 Pp = cgetg(dP-D + 3, t_POL);
2098 Pprimep = cgetg(dP-D + 3, t_POL);
2099 Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
2100 for(i=0; i < D; i++)
2101 {
2102 GEN c = gel(P, i+2);
2103 gel(Pm, i+2) = c;
2104 if (i) gel(Pprimem, i+1) = mului(i, c);
2105 }
2106 for(; i <= dP; i++)
2107 {
2108 GEN c = gel(P, i+2);
2109 gel(Pp, i+2-D) = c;
2110 gel(Pprimep, i+2-D) = mului(i, c);
2111 }
2112 *pPm = normalizepol_lg(Pm, D+2);
2113 *pPprimem = normalizepol_lg(Pprimem, D+1);
2114 *pPp = normalizepol_lg(Pp, dP-D+3);
2115 *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
2116 return dP - degpol(*pPp);
2117 }
2118
2119 static GEN
bkeval_single_power(long d,GEN V)2120 bkeval_single_power(long d, GEN V)
2121 {
2122 long mp = lg(V) - 2;
2123 if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2124 return gel(V, d+1);
2125 }
2126
2127 static GEN
splitpoleval(GEN Pp,GEN Pm,GEN pows,long D,long bitprec)2128 splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
2129 {
2130 GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
2131 GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
2132 GEN xa = bkeval_single_power(D, pows);
2133 GEN r;
2134 if (!signe(vp)) return vm;
2135 vp = gmul(vp, xa);
2136 r = gadd(vp, vm);
2137 if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
2138 return NULL;
2139 return r;
2140 }
2141
2142 /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
2143 static GEN
splitcauchy(GEN Pp,GEN Pm,long prec)2144 splitcauchy(GEN Pp, GEN Pm, long prec)
2145 {
2146 GEN S = gel(Pp,2), A = gel(Pm,2);
2147 long i, lPm = lg(Pm), lPp = lg(Pp);
2148 for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
2149 for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
2150 return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
2151 }
2152
2153 static GEN
ZX_deg1root(GEN P,long prec)2154 ZX_deg1root(GEN P, long prec)
2155 {
2156 GEN a = gel(P,3), b = gel(P,2);
2157 if (is_pm1(a))
2158 {
2159 b = itor(b, prec); if (signe(a) > 0) togglesign(b);
2160 return b;
2161 }
2162 return rdivii(negi(b), a, prec);
2163 }
2164
2165 /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
2166 * P' has also at most one zero there */
2167 static GEN
polsolve(GEN P,long bitprec)2168 polsolve(GEN P, long bitprec)
2169 {
2170 pari_sp av;
2171 GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
2172 long dP = degpol(P), prec = nbits2prec(bitprec);
2173 long expoold, iter, D, rt, s0, bitaddprec, cprec, PREC;
2174
2175 if (dP == 1) return ZX_deg1root(P, prec);
2176 Pprime = ZX_deriv(P);
2177 Pprime2 = ZX_deriv(Pprime);
2178 bitaddprec = 1 + 2*expu(dP); PREC = prec + nbits2prec(bitaddprec);
2179 D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
2180 s0 = signe(gel(P, 2));
2181 rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
2182 rb = splitcauchy(Pp, Pm, DEFAULTPREC);
2183 for (cprec = DEFAULTPREC, expoold = LONG_MAX;;)
2184 {
2185 GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2186 Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
2187 if (!Pc) { cprec += EXTRAPRECWORD; rb = rtor(rb, cprec); continue; }
2188 if (signe(Pc) != s0) break;
2189 shiftr_inplace(rb,1);
2190 }
2191 for (iter = 0, ra = NULL;;)
2192 {
2193 GEN wdth;
2194 iter++;
2195 if (ra)
2196 rc = shiftr(addrr(ra, rb), -1);
2197 else
2198 rc = shiftr(rb, -1);
2199 for(;;)
2200 {
2201 GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2202 Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
2203 if (Pc) break;
2204 cprec += EXTRAPRECWORD;
2205 rc = rtor(rc, cprec);
2206 }
2207 if (signe(Pc) == s0)
2208 ra = rc;
2209 else
2210 rb = rc;
2211 if (!ra) continue;
2212 wdth = subrr(rb, ra);
2213 if (!(iter % 8))
2214 {
2215 GEN m1 = poleval(Pprime, ra), M2;
2216 if (signe(m1) == s0) continue;
2217 M2 = poleval(Pprime2, rb);
2218 if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2219 break;
2220 }
2221 else if (gexpo(wdth) <= -bitprec)
2222 break;
2223 }
2224 rc = rb; av = avma;
2225 for(;; rc = gerepileuptoleaf(av, rc))
2226 {
2227 long exponew;
2228 GEN Ppc, dist, rcold = rc;
2229 GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2230 Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
2231 if (Ppc) Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
2232 if (!Ppc || !Pc)
2233 {
2234 if (cprec >= PREC)
2235 cprec += EXTRAPRECWORD;
2236 else
2237 cprec = minss(2*cprec, PREC);
2238 rc = rtor(rc, cprec); continue; /* backtrack one step */
2239 }
2240 dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
2241 rc = subrr(rc, dist);
2242 if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2243 {
2244 if (cprec >= PREC) break;
2245 cprec = minss(2*cprec, PREC);
2246 rc = rtor(rcold, cprec); continue; /* backtrack one step */
2247 }
2248 if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2249 exponew = expo(dist);
2250 if (exponew < -bitprec - 1)
2251 {
2252 if (cprec >= PREC) break;
2253 cprec = minss(2*cprec, PREC);
2254 rc = rtor(rc, cprec); continue;
2255 }
2256 if (exponew > expoold - 2)
2257 {
2258 if (cprec >= PREC) break;
2259 expoold = LONG_MAX;
2260 cprec = minss(2*cprec, PREC);
2261 rc = rtor(rc, cprec); continue;
2262 }
2263 expoold = exponew;
2264 }
2265 return rtor(rc, prec);
2266 }
2267
2268 /* Return primpart(P(x / 2)) */
2269 static GEN
ZX_rescale2prim(GEN P)2270 ZX_rescale2prim(GEN P)
2271 {
2272 long i, l = lg(P), v, n;
2273 GEN Q;
2274 if (l==2) return pol_0(varn(P));
2275 Q = cgetg(l,t_POL); v = vali(gel(P,l-1));
2276 for (i = l-2, n = 1; v > n && i >= 2; i--, n++)
2277 v = minss(v, vali(gel(P,i)) + n);
2278 gel(Q,l-1) = v? shifti(gel(P,l-1), -v): gel(P,l-1);
2279 for (i = l-2, n = 1-v; i >= 2; i--, n++)
2280 gel(Q,i) = shifti(gel(P,i), n);
2281 Q[1] = P[1]; return Q;
2282 }
2283
2284 /* assume Q0 has no rational root */
2285 static GEN
usp(GEN Q0,long flag,long bitprec)2286 usp(GEN Q0, long flag, long bitprec)
2287 {
2288 const pari_sp av = avma;
2289 GEN Qremapped, Q, c, Lc, Lk, sol;
2290 GEN *pQremapped = flag == 1? &Qremapped: NULL;
2291 const long prec = nbits2prec(bitprec), deg = degpol(Q0);
2292 long listsize = 64, nbr = 0, nb_todo, ind, indf, i, k, nb;
2293
2294 sol = zerocol(deg);
2295 Lc = zerovec(listsize);
2296 Lk = cgetg(listsize+1, t_VECSMALL);
2297 k = Lk[1] = 0;
2298 ind = 1; indf = 2;
2299 Q = Q0;
2300 c = gen_0;
2301 nb_todo = 1;
2302 while (nb_todo)
2303 {
2304 GEN nc = gel(Lc, ind);
2305 pari_sp av2;
2306 if (Lk[ind] == k + 1)
2307 {
2308 Q = Q0 = ZX_rescale2prim(Q0);
2309 c = gen_0;
2310 }
2311 if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
2312 av2 = avma;
2313 k = Lk[ind];
2314 ind++;
2315 c = nc;
2316 nb_todo--;
2317 nb = X2XP1(Q, pQremapped);
2318
2319 if (nb == 1)
2320 { /* exactly one root */
2321 GEN s = gen_0;
2322 if (flag == 0)
2323 {
2324 s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
2325 s = gerepilecopy(av2, s);
2326 }
2327 else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
2328 {
2329 s = polsolve(*pQremapped, bitprec+1);
2330 s = addir(c, divrr(s, addsr(1, s)));
2331 shiftr_inplace(s, -k);
2332 if (realprec(s) != prec) s = rtor(s, prec);
2333 s = gerepileupto(av2, s);
2334 }
2335 else set_avma(av2);
2336 gel(sol, ++nbr) = s;
2337 }
2338 else if (nb)
2339 { /* unknown, add two nodes to refine */
2340 if (indf + 2 > listsize)
2341 {
2342 if (ind>1)
2343 {
2344 for (i = ind; i < indf; i++)
2345 {
2346 gel(Lc, i-ind+1) = gel(Lc, i);
2347 Lk[i-ind+1] = Lk[i];
2348 }
2349 indf -= ind-1;
2350 ind = 1;
2351 }
2352 if (indf + 2 > listsize)
2353 {
2354 listsize *= 2;
2355 Lc = vec_lengthen(Lc, listsize);
2356 Lk = vecsmall_lengthen(Lk, listsize);
2357 }
2358 for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2359 }
2360 gel(Lc, indf) = nc = shifti(c, 1);
2361 gel(Lc, indf + 1) = addiu(nc, 1);
2362 Lk[indf] = Lk[indf + 1] = k + 1;
2363 indf += 2;
2364 nb_todo += 2;
2365 }
2366 if (gc_needed(av, 2))
2367 {
2368 gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
2369 if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
2370 }
2371 }
2372 setlg(sol, nbr+1);
2373 return gerepilecopy(av, sol);
2374 }
2375
2376 static GEN
ZX_Uspensky_equal_yes(GEN a,long flag,long bit)2377 ZX_Uspensky_equal_yes(GEN a, long flag, long bit)
2378 {
2379 if (flag == 2) return gen_1;
2380 if (flag == 1 && typ(a) != t_REAL)
2381 {
2382 if (typ(a) == t_INT && !signe(a))
2383 a = real_0_bit(bit);
2384 else
2385 a = gtofp(a, nbits2prec(bit));
2386 }
2387 return mkcol(a);
2388 }
2389 static GEN
ZX_Uspensky_no(long flag)2390 ZX_Uspensky_no(long flag)
2391 { return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
2392 /* ZX_Uspensky(P, [a,a], flag) */
2393 static GEN
ZX_Uspensky_equal(GEN P,GEN a,long flag,long bit)2394 ZX_Uspensky_equal(GEN P, GEN a, long flag, long bit)
2395 {
2396 if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
2397 return ZX_Uspensky_equal_yes(a, flag, bit);
2398 else
2399 return ZX_Uspensky_no(flag);
2400 }
2401 static int
sol_ok(GEN r,GEN a,GEN b)2402 sol_ok(GEN r, GEN a, GEN b) { return gcmp(a, r) <= 0 && gcmp(r, b) <= 0; }
2403
2404 /* P a ZX without real double roots; better if primitive and squarefree but
2405 * caller should ensure that. If flag & 4 assume that P has no rational root
2406 * (modest speedup) */
2407 GEN
ZX_Uspensky(GEN P,GEN ab,long flag,long bitprec)2408 ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
2409 {
2410 pari_sp av = avma;
2411 GEN a, b, res, sol;
2412 double fb;
2413 long l, nbz, deg;
2414
2415 if (ab)
2416 {
2417 if (typ(ab) == t_VEC)
2418 {
2419 if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
2420 a = gel(ab, 1);
2421 b = gel(ab, 2);
2422 }
2423 else
2424 {
2425 a = ab;
2426 b = mkoo();
2427 }
2428 }
2429 else
2430 {
2431 a = mkmoo();
2432 b = mkoo();
2433 }
2434 if (flag & 4)
2435 {
2436 if (gcmp(a, b) >= 0) { set_avma(av); return ZX_Uspensky_no(flag); }
2437 flag &= ~4;
2438 sol = cgetg(1, t_COL);
2439 }
2440 else
2441 {
2442 switch (gcmp(a, b))
2443 {
2444 case 1: set_avma(av); return ZX_Uspensky_no(flag);
2445 case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag, bitprec));
2446 }
2447 sol = nfrootsQ(P);
2448 }
2449 nbz = 0; l = lg(sol);
2450 if (l > 1)
2451 {
2452 long i, j;
2453 P = RgX_div(P, roots_to_pol(sol, varn(P)));
2454 if (!RgV_is_ZV(sol)) P = Q_primpart(P);
2455 for (i = j = 1; i < l; i++)
2456 if (sol_ok(gel(sol,i), a, b)) gel(sol,j++) = gel(sol,i);
2457 setlg(sol, j);
2458 if (flag == 2) { nbz = j-1; sol = utoi(nbz); }
2459 else if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
2460 }
2461 deg = degpol(P);
2462 if (deg == 0) return gerepilecopy(av, sol);
2463 if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
2464 {
2465 fb = fujiwara_bound_real(P, -1);
2466 if (fb <= -pariINFINITY) a = gen_0;
2467 else if (fb < 0) a = gen_m1;
2468 else a = negi(int2n((long)ceil(fb)));
2469 }
2470 if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
2471 {
2472 fb = fujiwara_bound_real(P, 1);
2473 if (fb <= -pariINFINITY) b = gen_0;
2474 else if (fb < 0) b = gen_1;
2475 else b = int2n((long)ceil(fb));
2476 }
2477 if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
2478 {
2479 GEN d, ad, bd, diff;
2480 long i;
2481 /* can occur if one of a,b was initially a t_INFINITY */
2482 if (gequal(a,b)) return gerepilecopy(av, sol);
2483 d = lcmii(Q_denom(a), Q_denom(b));
2484 if (is_pm1(d)) { d = NULL; ad = a; bd = b; }
2485 else
2486 { P = ZX_rescale(P, d); ad = gmul(a, d); bd = gmul(b, d); }
2487 diff = subii(bd, ad);
2488 P = ZX_unscale(ZX_translate(P, ad), diff);
2489 res = usp(P, flag, bitprec);
2490 if (flag <= 1)
2491 {
2492 for (i = 1; i < lg(res); i++)
2493 {
2494 GEN z = gmul(diff, gel(res, i));
2495 if (typ(z) == t_VEC)
2496 {
2497 gel(z, 1) = gadd(ad, gel(z, 1));
2498 gel(z, 2) = gadd(ad, gel(z, 2));
2499 }
2500 else
2501 z = gadd(ad, z);
2502 if (d) z = gdiv(z, d);
2503 gel(res, i) = z;
2504 }
2505 sol = shallowconcat(sol, res);
2506 }
2507 else
2508 nbz += lg(res) - 1;
2509 }
2510 if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(P, 1)) > -pariINFINITY)
2511 {
2512 long bp = maxss((long)ceil(fb), 0);
2513 res = usp(ZX_unscale2n(P, bp), flag, bitprec);
2514 if (flag <= 1)
2515 sol = shallowconcat(sol, gmul2n(res, bp));
2516 else
2517 nbz += lg(res)-1;
2518 }
2519 if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(P,-1)) > -pariINFINITY)
2520 {
2521 long i, bm = maxss((long)ceil(fb), 0);
2522 res = usp(ZX_unscale2n(ZX_z_unscale(P, -1), bm), flag, bitprec);
2523 if (flag <= 1)
2524 {
2525 for (i = 1; i < lg(res); i++)
2526 {
2527 GEN z = gneg(gmul2n(gel(res, i), bm));
2528 if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
2529 gel(res, i) = z;
2530 }
2531 sol = shallowconcat(res, sol);
2532 }
2533 else
2534 nbz += lg(res)-1;
2535 }
2536 if (flag >= 2) return utoi(nbz);
2537 if (flag)
2538 sol = sort(sol);
2539 else
2540 sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
2541 return gerepileupto(av, sol);
2542 }
2543
2544 /* x a scalar */
2545 static GEN
rootsdeg0(GEN x)2546 rootsdeg0(GEN x)
2547 {
2548 if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
2549 if (gequal0(x)) pari_err_ROOTS0("realroots");
2550 return cgetg(1,t_COL); /* constant polynomial */
2551 }
2552 static void
checkbound(GEN a)2553 checkbound(GEN a)
2554 {
2555 switch(typ(a))
2556 {
2557 case t_INT: case t_FRAC: case t_INFINITY: break;
2558 default: pari_err_TYPE("polrealroots", a);
2559 }
2560 }
2561 static GEN
check_ab(GEN ab)2562 check_ab(GEN ab)
2563 {
2564 GEN a, b;
2565 if (!ab) return NULL;
2566 if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
2567 a = gel(ab,1); checkbound(a);
2568 b = gel(ab,2); checkbound(b);
2569 if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
2570 typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
2571 return ab;
2572 }
2573 /* e^(1/h) assuming the h-th root is real, beware that sqrtnr assumes e >= 0 */
2574 static GEN
_sqrtnr(GEN e,long h)2575 _sqrtnr(GEN e, long h)
2576 {
2577 long s;
2578 GEN r;
2579 if (h == 2) return sqrtr(e);
2580 s = signe(e); setsigne(e, 1); /* e < 0 is possible, implies h is odd */
2581 r = sqrtnr(e, h); if (s < 0) setsigne(r, -1);
2582 return r;
2583 }
2584 GEN
realroots(GEN P,GEN ab,long prec)2585 realroots(GEN P, GEN ab, long prec)
2586 {
2587 pari_sp av = avma;
2588 GEN sol = NULL, fa, ex;
2589 long i, j, v, l;
2590
2591 ab = check_ab(ab);
2592 if (typ(P) != t_POL) return rootsdeg0(P);
2593 switch(degpol(P))
2594 {
2595 case -1: return rootsdeg0(gen_0);
2596 case 0: return rootsdeg0(gel(P,2));
2597 }
2598 if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
2599 v = ZX_valrem(Q_primpart(P), &P);
2600 fa = ZX_squff(P, &ex); l = lg(fa); sol = cgetg(l + 1, t_VEC);
2601 for (i = 1; i < l; i++)
2602 {
2603 GEN Pi = gel(fa, i), soli, soli2;
2604 long n, h;
2605 if (ab) h = 1; else Pi = ZX_deflate_max(Pi, &h);
2606 soli = ZX_Uspensky(Pi, odd(h)? ab: gen_0, 1, prec2nbits(prec));
2607 n = lg(soli); soli2 = odd(h)? NULL: cgetg(n, t_COL);
2608 for (j = 1; j < n; j++)
2609 {
2610 GEN r = gel(soli, j); /* != 0 */
2611 if (typ(r) != t_REAL) gel(soli, j) = r = gtofp(r, prec);
2612 if (h > 1)
2613 {
2614 gel(soli, j) = r = _sqrtnr(r, h);
2615 if (soli2) gel(soli2, j) = negr(r);
2616 }
2617 }
2618 if (soli2) soli = shallowconcat(soli, soli2);
2619 if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
2620 gel(sol, i) = soli;
2621 }
2622 if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
2623 gel(sol, i++) = const_col(v, real_0(prec));
2624 setlg(sol, i); if (i == 1) { set_avma(av); return cgetg(1,t_COL); }
2625 return gerepileupto(av, sort(shallowconcat1(sol)));
2626 }
2627 GEN
ZX_realroots_irred(GEN P,long prec)2628 ZX_realroots_irred(GEN P, long prec)
2629 {
2630 long dP = degpol(P), j, n, h;
2631 GEN sol, sol2;
2632 pari_sp av;
2633 if (dP == 1) retmkvec(ZX_deg1root(P, prec));
2634 av = avma; P = ZX_deflate_max(P, &h);
2635 if (h == dP)
2636 {
2637 GEN r = _sqrtnr(ZX_deg1root(P, prec), h);
2638 return gerepilecopy(av, odd(h)? mkvec(r): mkvec2(negr(r), r));
2639 }
2640 sol = ZX_Uspensky(P, odd(h)? NULL: gen_0, 1 | 4, prec2nbits(prec));
2641 n = lg(sol); sol2 = odd(h)? NULL: cgetg(n, t_COL);
2642 for (j = 1; j < n; j++)
2643 {
2644 GEN r = gel(sol, j);
2645 if (typ(r) != t_REAL) gel(sol, j) = r = gtofp(r, prec);
2646 if (h > 1)
2647 {
2648 gel(sol, j) = r = _sqrtnr(r, h);
2649 if (sol2) gel(sol2, j) = negr(r);
2650 }
2651 }
2652 if (sol2) sol = shallowconcat(sol, sol2);
2653 return gerepileupto(av, sort(sol));
2654 }
2655
2656 static long
ZX_sturm_i(GEN P,long flag)2657 ZX_sturm_i(GEN P, long flag)
2658 {
2659 pari_sp av;
2660 long h, r, dP = degpol(P);
2661 if (dP == 1) return 1;
2662 av = avma; P = ZX_deflate_max(P, &h);
2663 if (h == dP)
2664 { /* now deg P = 1 */
2665 if (odd(h))
2666 r = 1;
2667 else
2668 r = (signe(gel(P,2)) != signe(gel(P,3)))? 2: 0;
2669 return gc_long(av, r);
2670 }
2671 if (odd(h))
2672 r = itou(ZX_Uspensky(P, NULL, flag, 0));
2673 else
2674 r = 2*itou(ZX_Uspensky(P, gen_0, flag, 0));
2675 return gc_long(av,r);
2676 }
2677 /* P nonconstant, squarefree ZX */
2678 long
ZX_sturmpart(GEN P,GEN ab)2679 ZX_sturmpart(GEN P, GEN ab)
2680 {
2681 pari_sp av = avma;
2682 if (!check_ab(ab)) return ZX_sturm(P);
2683 return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
2684 }
2685 /* P nonconstant, squarefree ZX */
2686 long
ZX_sturm(GEN P)2687 ZX_sturm(GEN P) { return ZX_sturm_i(P, 2); }
2688 /* P irreducible ZX */
2689 long
ZX_sturm_irred(GEN P)2690 ZX_sturm_irred(GEN P) { return ZX_sturm_i(P, 2 + 4); }
2691