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