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