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 1230123 : isvalidcoeff(GEN x)
30 : {
31 1230123 : switch (typ(x))
32 : {
33 1205597 : case t_INT: case t_REAL: case t_FRAC: return 1;
34 24512 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
35 : }
36 14 : return 0;
37 : }
38 :
39 : static void
40 267053 : checkvalidpol(GEN p, const char *f)
41 : {
42 267053 : long i,n = lg(p);
43 1448131 : for (i=2; i<n; i++)
44 1181085 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE(f, gel(p,i));
45 267046 : }
46 :
47 : /********************************************************************/
48 : /** **/
49 : /** FAST ARITHMETIC over Z[i] **/
50 : /** **/
51 : /********************************************************************/
52 :
53 : static GEN
54 16452553 : ZX_to_ZiX(GEN Pr, GEN Pi)
55 : {
56 16452553 : long i, lr = lg(Pr), li = lg(Pi), l = maxss(lr, li), m = minss(lr, li);
57 16454586 : GEN P = cgetg(l, t_POL);
58 16458639 : P[1] = Pr[1];
59 66744409 : for(i = 2; i < m; i++)
60 50286028 : gel(P,i) = signe(gel(Pi,i)) ? mkcomplex(gel(Pr,i), gel(Pi,i))
61 50286028 : : gel(Pr,i);
62 22359287 : for( ; i < lr; i++)
63 5900906 : gel(P,i) = gel(Pr, i);
64 16493252 : for( ; i < li; i++)
65 34871 : gel(P,i) = mkcomplex(gen_0, gel(Pi, i));
66 16458381 : return normalizepol_lg(P, l);
67 : }
68 :
69 : static GEN
70 99031442 : ZiX_sqr(GEN P)
71 : {
72 99031442 : pari_sp av = avma;
73 : GEN Pr2, Pi2, Qr, Qi;
74 99031442 : GEN Pr = real_i(P), Pi = imag_i(P);
75 99028448 : if (signe(Pi)==0) return gc_upto(av, ZX_sqr(Pr));
76 16518805 : if (signe(Pr)==0) return gc_upto(av, ZX_neg(ZX_sqr(Pi)));
77 16461492 : Pr2 = ZX_sqr(Pr); Pi2 = ZX_sqr(Pi);
78 16453585 : Qr = ZX_sub(Pr2, Pi2);
79 16453628 : if (degpol(Pr)==degpol(Pi))
80 10700781 : Qi = ZX_sub(ZX_sqr(ZX_add(Pr, Pi)), ZX_add(Pr2, Pi2));
81 : else
82 5757200 : Qi = ZX_shifti(ZX_mul(Pr, Pi), 1);
83 16456014 : return gc_GEN(av, ZX_to_ZiX(Qr, Qi));
84 : }
85 :
86 : static GEN
87 49503340 : graeffe(GEN p)
88 : {
89 49503340 : pari_sp av = avma;
90 : GEN p0, p1, s0, s1;
91 49503340 : long n = degpol(p);
92 :
93 49508606 : if (!n) return RgX_copy(p);
94 49508606 : RgX_even_odd(p, &p0, &p1);
95 : /* p = p0(x^2) + x p1(x^2) */
96 49517702 : s0 = ZiX_sqr(p0);
97 49525351 : s1 = ZiX_sqr(p1);
98 49524869 : return gc_upto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
99 : }
100 :
101 : GEN
102 5383 : ZX_graeffe(GEN p)
103 : {
104 5383 : pari_sp av = avma;
105 : GEN p0, p1, s0, s1;
106 5383 : long n = degpol(p);
107 :
108 5383 : if (!n) return ZX_copy(p);
109 5383 : RgX_even_odd(p, &p0, &p1);
110 : /* p = p0(x^2) + x p1(x^2) */
111 5383 : s0 = ZX_sqr(p0);
112 5383 : s1 = ZX_sqr(p1);
113 5383 : return gc_upto(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 gc_upto(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 226354504 : mydbllog2i(GEN x)
143 : {
144 : #ifdef LONG_IS_64BIT
145 194682317 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
146 : #else
147 31672187 : const double W = 1/4294967296.; /*2^-32*/
148 : #endif
149 : GEN m;
150 226354504 : long lx = lgefint(x);
151 : double l;
152 226354504 : if (lx == 2) return -pariINFINITY;
153 225537455 : m = int_MSW(x);
154 225537455 : l = (double)(ulong)*m;
155 225537455 : if (lx == 3) return log2(l);
156 70116266 : 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 70116266 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
161 : }
162 :
163 : /* return log(|x|) or -pariINFINITY */
164 : static double
165 9545425 : mydbllogr(GEN x) {
166 9545425 : if (!signe(x)) return -pariINFINITY;
167 9545425 : return M_LN2*dbllog2r(x);
168 : }
169 :
170 : /* return log2(|x|) or -pariINFINITY */
171 : static double
172 56762858 : mydbllog2r(GEN x) {
173 56762858 : if (!signe(x)) return -pariINFINITY;
174 56324413 : return dbllog2r(x);
175 : }
176 : double
177 302618439 : dbllog2(GEN z)
178 : {
179 : double x, y;
180 302618439 : switch(typ(z))
181 : {
182 226247275 : case t_INT: return mydbllog2i(z);
183 22358 : case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
184 49308588 : case t_REAL: return mydbllog2r(z);
185 27040218 : default: /*t_COMPLEX*/
186 27040218 : x = dbllog2(gel(z,1));
187 27122968 : y = dbllog2(gel(z,2));
188 27122809 : if (x == -pariINFINITY) return y;
189 26879695 : if (y == -pariINFINITY) return x;
190 26675904 : if (fabs(x-y) > 10) return maxdd(x,y);
191 26057377 : return x + 0.5*log2(1 + exp2(2*(y-x)));
192 : }
193 : }
194 : static GEN /* beware overflow */
195 6548521 : 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 40954611 : findpower(GEN p)
201 : {
202 40954611 : double x, L, mins = pariINFINITY;
203 40954611 : long n = degpol(p),i;
204 :
205 40954067 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
206 164646489 : for (i=n-1; i>=0; i--)
207 : {
208 123690892 : L += log2((double)(i+1) / (double)(n-i));
209 123690892 : x = dbllog2(gel(p,i+2));
210 123688306 : if (x != -pariINFINITY)
211 : {
212 122897768 : double s = (L - x) / (double)(n-i);
213 122897768 : if (s < mins) mins = s;
214 : }
215 : }
216 40955597 : i = (long)ceil(mins);
217 40955597 : if (i - mins > 1 - 1e-12) i--;
218 40955597 : return i;
219 : }
220 :
221 : /* returns the exponent for logmodulus(), from the Newton diagram */
222 : static long
223 5506732 : newton_polygon(GEN p, long k)
224 : {
225 5506732 : pari_sp av = avma;
226 5506732 : long n = degpol(p), i, j, h, l, *vertex = (long*)new_chunk(n+1);
227 5506712 : 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 26452883 : for (i=0; i<=n; i++) { L[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
231 5506691 : vertex[0] = 1; /* sentinel */
232 19605085 : for (i=0; i < n; i=h)
233 : {
234 : double slope;
235 14098394 : h = i+1;
236 14103263 : while (L[i] == -pariINFINITY) { vertex[h] = 1; i = h; h = i+1; }
237 14098394 : slope = L[h] - L[i];
238 37962127 : for (j = i+2; j<=n; j++) if (L[j] != -pariINFINITY)
239 : {
240 23857500 : double pij = (L[j] - L[i])/(double)(j - i);
241 23857500 : if (slope < pij) { slope = pij; h = j; }
242 : }
243 14098394 : vertex[h] = 1;
244 : }
245 6222128 : h = k; while (!vertex[h]) h++;
246 5695097 : l = k-1; while (!vertex[l]) l--;
247 5506691 : set_avma(av);
248 5506722 : 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 35504087 : myshiftrc(GEN z, long e)
254 : {
255 35504087 : if (typ(z)==t_COMPLEX)
256 : {
257 6409775 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
258 6409778 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
259 : }
260 : else
261 29094312 : if (signe(z)) shiftr_inplace(z, e);
262 35504161 : }
263 :
264 : /* return z*2^e, where z is integer or complex of integer (destroy z) */
265 : static GEN
266 135811716 : myshiftic(GEN z, long e)
267 : {
268 135811716 : if (typ(z)==t_COMPLEX)
269 : {
270 17890935 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
271 17889803 : gel(z,2) = mpshift(gel(z,2),e);
272 17886941 : return z;
273 : }
274 117920781 : return signe(z)? mpshift(z,e): gen_0;
275 : }
276 :
277 : static GEN
278 7024514 : RgX_gtofp_bit(GEN q, long bit) { return RgX_gtofp(q, nbits2prec(bit)); }
279 :
280 : static GEN
281 215009986 : mygprecrc(GEN x, long prec, long e)
282 : {
283 : GEN y;
284 215009986 : switch(typ(x))
285 : {
286 160465938 : case t_REAL:
287 160465938 : if (!signe(x)) return real_0_bit(e);
288 157192239 : return realprec(x) == prec? x: rtor(x, prec);
289 37184793 : case t_COMPLEX:
290 37184793 : y = cgetg(3,t_COMPLEX);
291 37184625 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
292 37183935 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
293 37184120 : return y;
294 17359255 : 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 63710193 : mygprec(GEN x, long bit)
302 : {
303 : long e, prec;
304 63710193 : if (bit < 0) bit = 0; /* should rarely happen */
305 63710193 : e = gexpo(x) - bit;
306 63710767 : prec = nbits2prec(bit);
307 63716582 : switch(typ(x))
308 : {
309 167987055 : case t_POL: pari_APPLY_pol_normalized(mygprecrc(gel(x,i),prec,e));
310 17251508 : default: return mygprecrc(x,prec,e);
311 : }
312 : }
313 :
314 : /* normalize a polynomial x, that is change it with coefficients in Z[i],
315 : after making product by 2^shift */
316 : static GEN
317 17587582 : pol_to_gaussint(GEN x, long shift)
318 100320380 : { pari_APPLY_pol_normalized(gtrunc2n(gel(x,i), shift)); }
319 :
320 : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
321 : static GEN
322 13137107 : eval_rel_pol(GEN p, long bit)
323 : {
324 : long i;
325 73991651 : for (i = 2; i < lg(p); i++)
326 60854504 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behavior of gexpo */
327 13137147 : return pol_to_gaussint(p, bit-gexpo(p)+1);
328 : }
329 :
330 : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
331 : static GEN
332 1609961 : homothetie(GEN p, double lrho, long bit)
333 : {
334 : GEN q, r, t, iR;
335 1609961 : long n = degpol(p), i;
336 :
337 1609957 : iR = mygprec(dblexp(-lrho),bit);
338 1609921 : q = mygprec(p, bit);
339 1609949 : r = cgetg(n+3,t_POL); r[1] = p[1];
340 1609942 : t = iR; r[n+2] = q[n+2];
341 6791795 : for (i=n-1; i>0; i--)
342 : {
343 5181881 : gel(r,i+2) = gmul(t, gel(q,i+2));
344 5181866 : t = mulrr(t, iR);
345 : }
346 1609914 : gel(r,2) = gmul(t, gel(q,2)); return r;
347 : }
348 :
349 : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) [ ~as above with R = 2^-e ]*/
350 : static void
351 9956918 : homothetie2n(GEN p, long e)
352 : {
353 9956918 : if (e)
354 : {
355 7975019 : long i,n = lg(p)-1;
356 43479378 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
357 : }
358 9957054 : }
359 :
360 : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
361 : static void
362 36494028 : homothetie_gauss(GEN p, long e, long f)
363 : {
364 36494028 : if (e || f)
365 : {
366 32601626 : long i, n = lg(p)-1;
367 168354233 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
368 : }
369 36441257 : }
370 :
371 : /* Lower bound on the modulus of the largest root z_0
372 : * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
373 : static double
374 40954124 : lower_bound(GEN p, long *k, double eps)
375 : {
376 40954124 : long n = degpol(p), i, j;
377 40954367 : pari_sp ltop = avma;
378 : GEN a, s, S, ilc;
379 : double r, R, rho;
380 :
381 40954367 : if (n < 4) { *k = n; return 0.; }
382 8235839 : S = cgetg(5,t_VEC);
383 8236549 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
384 41160853 : for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
385 : /* i = 1 split out from next loop for efficiency and initialization */
386 8235352 : s = gel(a,1);
387 8235352 : gel(S,1) = gneg(s); /* Newton sum S_i */
388 8235928 : rho = r = gtodouble(gabs(s,3));
389 8236531 : R = r / n;
390 32939140 : for (i=2; i<=4; i++)
391 : {
392 24703109 : s = gmulsg(i,gel(a,i));
393 74029339 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
394 24688045 : gel(S,i) = gneg(s); /* Newton sum S_i */
395 24697431 : r = gtodouble(gabs(s,3));
396 24702609 : if (r > 0.)
397 : {
398 24656106 : r = exp(log(r/n) / (double)i);
399 24656106 : if (r > R) R = r;
400 : }
401 : }
402 8236031 : if (R > 0. && eps < 1.2)
403 8232296 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
404 : else
405 3735 : *k = n;
406 8236031 : return gc_double(ltop, R);
407 : }
408 :
409 : /* return R such that exp(R - tau) <= rho_n(P) <= exp(R + tau)
410 : * P(0) != 0 and P non constant */
411 : static double
412 4450567 : logmax_modulus(GEN p, double tau)
413 : {
414 : GEN r, q, aux, gunr;
415 4450567 : pari_sp av, ltop = avma;
416 4450567 : long i,k,n=degpol(p),nn,bit,M,e;
417 4450568 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
418 :
419 4450568 : r = cgeti(BIGDEFAULTPREC);
420 4450552 : av = avma;
421 :
422 4450552 : eps = - 1/log(1.5*tau2); /* > 0 */
423 4450552 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
424 4450552 : gunr = real_1_bit(bit+2*n);
425 4450505 : aux = gdiv(gunr, gel(p,2+n));
426 4450467 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
427 4450343 : e = findpower(q);
428 4450489 : homothetie2n(q,e);
429 4450535 : affsi(e, r);
430 4450535 : q = pol_to_gaussint(q, bit);
431 4450170 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
432 4450170 : nn = n;
433 4450170 : for (i=0,e=0;;)
434 : { /* nn = deg(q) */
435 40956728 : rho = lower_bound(q, &k, eps);
436 40954288 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
437 40954288 : affii(shifti(addis(r,e), 1), r);
438 40911277 : if (++i == M) break;
439 :
440 36461442 : bit = (long) ((double)k * log2(1./tau2) +
441 36461442 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
442 36461442 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
443 36479319 : nn -= RgX_valrem(q, &q);
444 36483270 : q = gc_upto(av, graeffe(q));
445 36505677 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
446 36505677 : eps = -1/log(tau2); /* > 0 */
447 36505677 : e = findpower(q);
448 : }
449 4449835 : if (!signe(r)) return gc_double(ltop,0.);
450 4027424 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
451 4027844 : return gc_double(ltop, -rtodbl(r) * M_LN2); /* -log(2) sum e_i 2^-i */
452 : }
453 :
454 : static GEN
455 35887 : RgX_normalize1(GEN x)
456 : {
457 35887 : long i, n = lg(x)-1;
458 : GEN y;
459 35901 : for (i = n; i > 1; i--)
460 35894 : if (!gequal0( gel(x,i) )) break;
461 35887 : if (i == n) return x;
462 14 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
463 14 : if (i == 1) pari_err_ROOTS0("roots");
464 14 : y = cgetg(i+1, t_POL); y[1] = x[1];
465 42 : for (; i > 1; i--) gel(y,i) = gel(x,i);
466 14 : return y;
467 : }
468 :
469 : static GEN
470 27293 : polrootsbound_i(GEN P, double TAU)
471 : {
472 27293 : pari_sp av = avma;
473 : double d;
474 27293 : (void)RgX_valrem_inexact(P,&P);
475 27293 : P = RgX_normalize1(P);
476 27293 : switch(degpol(P))
477 : {
478 7 : case -1: pari_err_ROOTS0("roots");
479 140 : case 0: set_avma(av); return gen_0;
480 : }
481 27146 : d = logmax_modulus(P, TAU) + TAU;
482 : /* not dblexp: result differs on ARM emulator */
483 27146 : return gc_leaf(av, mpexp(dbltor(d)));
484 : }
485 : GEN
486 27300 : polrootsbound(GEN P, GEN tau)
487 : {
488 27300 : if (typ(P) != t_POL) pari_err_TYPE("polrootsbound",P);
489 27293 : checkvalidpol(P, "polrootsbound");
490 27293 : return polrootsbound_i(P, tau? gtodouble(tau): 0.01);
491 : }
492 :
493 : /* log of modulus of the smallest root of p, with relative error tau */
494 : static double
495 1614697 : logmin_modulus(GEN p, double tau)
496 : {
497 1614697 : pari_sp av = avma;
498 1614697 : if (gequal0(gel(p,2))) return -pariINFINITY;
499 1614694 : return gc_double(av, - logmax_modulus(RgX_recip_i(p),tau));
500 : }
501 :
502 : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
503 : static double
504 606700 : logmodulus(GEN p, long k, double tau)
505 : {
506 : GEN q;
507 606700 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
508 606700 : pari_sp av, ltop=avma;
509 606700 : double r, tau2 = tau/6;
510 :
511 606700 : bit = (long)(n * (2. + log2(3.*n/tau2)));
512 606700 : av = avma;
513 606700 : q = gprec_w(p, nbits2prec(bit));
514 606704 : q = RgX_gtofp_bit(q, bit);
515 606701 : e = newton_polygon(q,k);
516 606702 : r = (double)e;
517 606702 : homothetie2n(q,e);
518 606717 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
519 5506705 : for (i=1; i<imax; i++)
520 : {
521 4900005 : q = eval_rel_pol(q,bit);
522 4899724 : kk -= RgX_valrem(q, &q);
523 4899842 : nn = degpol(q);
524 :
525 4899834 : q = gc_upto(av, graeffe(q));
526 4900032 : e = newton_polygon(q,kk);
527 4900036 : r += e / exp2((double)i);
528 4900036 : q = RgX_gtofp_bit(q, bit);
529 4899917 : homothetie2n(q,e);
530 :
531 4899988 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
532 4899988 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
533 : }
534 606700 : return gc_double(ltop, -r * M_LN2);
535 : }
536 :
537 : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
538 : * rmin < r_k < rmax. This information helps because we may reduce precision
539 : * quicker */
540 : static double
541 606702 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
542 : {
543 : GEN q;
544 606702 : long n = degpol(p), i, imax, imax2, bit;
545 606701 : pari_sp ltop = avma, av;
546 606701 : double lrho, aux, tau2 = tau/6.;
547 :
548 606701 : aux = (lrmax - lrmin) / 2. + 4*tau2;
549 606701 : imax = (long) log2(log((double)n)/ aux);
550 606701 : if (imax <= 0) return logmodulus(p,k,tau);
551 :
552 597587 : lrho = (lrmin + lrmax) / 2;
553 597587 : av = avma;
554 597587 : bit = (long)(n*(2. + aux / M_LN2 - log2(tau2)));
555 597587 : q = homothetie(p, lrho, bit);
556 597581 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
557 597581 : if (imax > imax2) imax = imax2;
558 :
559 1584590 : for (i=0; i<imax; i++)
560 : {
561 987003 : q = eval_rel_pol(q,bit);
562 987000 : q = gc_upto(av, graeffe(q));
563 987020 : aux = 2*aux + 2*tau2;
564 987020 : tau2 *= 1.5;
565 987020 : bit = (long)(n*(2. + aux / M_LN2 - log2(1-exp(-tau2))));
566 987020 : q = RgX_gtofp_bit(q, bit);
567 : }
568 597587 : aux = exp2((double)imax);
569 597587 : return gc_double(ltop, lrho + logmodulus(q,k, aux*tau/3.) / aux);
570 : }
571 :
572 : static double
573 902604 : ind_maxlog2(GEN q)
574 : {
575 902604 : long i, k = -1;
576 902604 : double L = - pariINFINITY;
577 2222092 : for (i=0; i<=degpol(q); i++)
578 : {
579 1319483 : double d = dbllog2(gel(q,2+i));
580 1319488 : if (d > L) { L = d; k = i; }
581 : }
582 902605 : return k;
583 : }
584 :
585 : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
586 : * Assume that l <= k <= n-l */
587 : static long
588 1012373 : dual_modulus(GEN p, double lrho, double tau, long l)
589 : {
590 1012373 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
591 1012372 : double tau2 = tau * 7./8.;
592 1012372 : pari_sp av = avma;
593 : GEN q;
594 :
595 1012372 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
596 1012372 : q = homothetie(p, lrho, bit);
597 1012369 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
598 :
599 8152990 : for (i=0; i<imax; i++)
600 : {
601 7250386 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
602 7249205 : v = RgX_valrem(q, &q);
603 7249798 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
604 :
605 7249946 : nn = degpol(q); delta_k += v;
606 7249989 : if (!nn) return delta_k;
607 :
608 7140217 : q = gc_upto(av, graeffe(q));
609 7140621 : tau2 *= 7./4.;
610 7140621 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
611 : }
612 902604 : return gc_long(av, delta_k + (long)ind_maxlog2(q));
613 : }
614 :
615 : /********************************************************************/
616 : /** **/
617 : /** FACTORS THROUGH CIRCLE INTEGRATION **/
618 : /** **/
619 : /********************************************************************/
620 : /* l power of 2, W[step*j] = w_j; set f[j] = p(w_j)
621 : * if inv, w_j = exp(2IPi*j/l), else exp(-2IPi*j/l) */
622 :
623 : static void
624 7462 : fft2(GEN W, GEN p, GEN f, long step, long l)
625 : {
626 : pari_sp av;
627 : long i, s1, l1, step2;
628 :
629 7462 : if (l == 2)
630 : {
631 3766 : gel(f,0) = gadd(gel(p,0), gel(p,step));
632 3766 : gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
633 : }
634 3696 : av = avma;
635 3696 : l1 = l>>1; step2 = step<<1;
636 3696 : fft2(W,p, f, step2,l1);
637 3696 : fft2(W,p+step, f+l1,step2,l1);
638 32760 : for (i = s1 = 0; i < l1; i++, s1 += step)
639 : {
640 29064 : GEN f0 = gel(f,i);
641 29064 : GEN f1 = gmul(gel(W,s1), gel(f,i+l1));
642 29064 : gel(f,i) = gadd(f0, f1);
643 29064 : gel(f,i+l1) = gsub(f0, f1);
644 : }
645 3696 : gc_slice(av, f, l);
646 : }
647 :
648 : static void
649 14158125 : fft(GEN W, GEN p, GEN f, long step, long l, long inv)
650 : {
651 : pari_sp av;
652 : long i, s1, l1, l2, l3, step4;
653 : GEN f1, f2, f3, f02;
654 :
655 14158125 : if (l == 2)
656 : {
657 6639204 : gel(f,0) = gadd(gel(p,0), gel(p,step));
658 6638895 : gel(f,1) = gsub(gel(p,0), gel(p,step)); return;
659 : }
660 7518921 : av = avma;
661 7518921 : if (l == 4)
662 : {
663 : pari_sp av2;
664 5318285 : f1 = gadd(gel(p,0), gel(p,step<<1));
665 5317711 : f2 = gsub(gel(p,0), gel(p,step<<1));
666 5317693 : f3 = gadd(gel(p,step), gel(p,3*step));
667 5317707 : f02 = gsub(gel(p,step), gel(p,3*step));
668 5317752 : f02 = inv? mulcxI(f02): mulcxmI(f02);
669 5318156 : av2 = avma;
670 5318156 : gel(f,0) = gadd(f1, f3);
671 5317494 : gel(f,1) = gadd(f2, f02);
672 5317709 : gel(f,2) = gsub(f1, f3);
673 5317615 : gel(f,3) = gsub(f2, f02);
674 5317814 : gc_all_unsafe(av,av2,4,&gel(f,0),&gel(f,1),&gel(f,2),&gel(f,3));
675 5318454 : return;
676 : }
677 2200636 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
678 2200636 : fft(W,p, f, step4,l1,inv);
679 2201010 : fft(W,p+step, f+l1,step4,l1,inv);
680 2201031 : fft(W,p+(step<<1),f+l2,step4,l1,inv);
681 2201039 : fft(W,p+3*step, f+l3,step4,l1,inv);
682 8194420 : for (i = s1 = 0; i < l1; i++, s1 += step)
683 : {
684 5993423 : long s2 = s1 << 1, s3 = s1 + s2;
685 : GEN g02, g13, f13;
686 5993423 : f1 = gmul(gel(W,s1), gel(f,i+l1));
687 5993701 : f2 = gmul(gel(W,s2), gel(f,i+l2));
688 5993640 : f3 = gmul(gel(W,s3), gel(f,i+l3));
689 :
690 5993674 : f02 = gadd(gel(f,i),f2);
691 5993044 : g02 = gsub(gel(f,i),f2);
692 5993102 : f13 = gadd(f1,f3);
693 5993069 : g13 = gsub(f1,f3); g13 = inv? mulcxI(g13): mulcxmI(g13);
694 :
695 5993480 : gel(f,i) = gadd(f02, f13);
696 5993073 : gel(f,i+l1) = gadd(g02, g13);
697 5993112 : gel(f,i+l2) = gsub(f02, f13);
698 5993200 : gel(f,i+l3) = gsub(g02, g13);
699 : }
700 2200997 : gc_slice(av, f, l);
701 : }
702 :
703 : static GEN
704 98 : FFT_i(GEN W, GEN x)
705 : {
706 98 : long i, l = lg(W), n = lg(x), tx = typ(x), tw, pa;
707 : GEN y, z, p, pol;
708 98 : if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
709 84 : tw = RgV_type(W, &p, &pol, &pa);
710 84 : if (tx == t_POL) { x++; n--; }
711 49 : else if (!is_vec_t(tx)) pari_err_TYPE("fft",x);
712 84 : if (n > l) pari_err_DIM("fft");
713 84 : if (n < l) {
714 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
715 0 : for (i = 1; i < n; i++) gel(z,i) = gel(x,i);
716 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
717 : }
718 84 : else z = x;
719 84 : if (l == 2) return mkveccopy(gel(z,1));
720 70 : y = cgetg(l, t_VEC);
721 70 : if (tw == RgX_type_code(t_COMPLEX,t_INT) ||
722 : tw == RgX_type_code(t_COMPLEX,t_REAL))
723 0 : {
724 0 : long inv = (l >= 5 && signe(imag_i(gel(W,1+(l>>2))))==1) ? 1 : 0;
725 0 : fft(W+1, z+1, y+1, 1, l-1, inv);
726 : } else
727 70 : fft2(W+1, z+1, y+1, 1, l-1);
728 70 : return y;
729 : }
730 :
731 : GEN
732 56 : FFT(GEN W, GEN x)
733 : {
734 56 : if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
735 56 : return FFT_i(W, x);
736 : }
737 :
738 : GEN
739 56 : FFTinv(GEN W, GEN x)
740 : {
741 56 : long l = lg(W), i;
742 : GEN w;
743 56 : if (!is_vec_t(typ(W))) pari_err_TYPE("fft",W);
744 56 : if (l==1 || ((l-1) & (l-2))) pari_err_DIM("fft");
745 42 : w = cgetg(l, t_VECSMALL); /* cf stackdummy */
746 42 : gel(w,1) = gel(W,1); /* w = gconj(W), faster */
747 3773 : for (i = 2; i < l; i++) gel(w, i) = gel(W, l-i+1);
748 42 : return FFT_i(w, x);
749 : }
750 :
751 : /* returns 1 if p has only real coefficients, 0 else */
752 : static int
753 961820 : isreal(GEN p)
754 : {
755 : long i;
756 4859685 : for (i = lg(p)-1; i > 1; i--)
757 4059297 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
758 800388 : return 1;
759 : }
760 :
761 : /* x non complex */
762 : static GEN
763 778631 : abs_update_r(GEN x, double *mu) {
764 778631 : GEN y = gtofp(x, DEFAULTPREC);
765 778637 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
766 778636 : setabssign(y); return y;
767 : }
768 : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
769 : static GEN
770 8004359 : abs_update(GEN x, double *mu) {
771 : GEN y, xr, yr;
772 : double ly;
773 8004359 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
774 7239060 : xr = gel(x,1);
775 7239060 : yr = gel(x,2);
776 7239060 : if (gequal0(xr)) return abs_update_r(yr,mu);
777 7237129 : if (gequal0(yr)) return abs_update_r(xr,mu);
778 : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
779 7225989 : xr = gtofp(xr, DEFAULTPREC);
780 7226778 : yr = gtofp(yr, DEFAULTPREC);
781 7226824 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
782 7226377 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
783 7226447 : return y;
784 : }
785 :
786 : static void
787 996255 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
788 : {
789 996255 : long prec = nbits2prec(bit);
790 996255 : *Omega = grootsof1(Lmax, prec) + 1;
791 996250 : *prim = rootsof1u_cx(N, prec);
792 996259 : }
793 :
794 : static void
795 493930 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
796 : int polreal, double param, double param2)
797 : {
798 : GEN q, pc, Omega, A, RU, prim, g, TWO;
799 493930 : long n = degpol(p), bit, NN, K, i, j, Lmax;
800 493930 : pari_sp av2, av = avma;
801 :
802 493930 : bit = gexpo(p) + (long)param2+8;
803 683612 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
804 493931 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
805 493931 : K = NN/Lmax; if (K & 1) K++;
806 493931 : NN = Lmax*K;
807 493931 : if (polreal) K = K/2+1;
808 :
809 493931 : initdft(&Omega, &prim, NN, Lmax, bit);
810 493934 : q = mygprec(p,bit) + 2;
811 493931 : A = cgetg(Lmax+1,t_VEC); A++;
812 493930 : pc= cgetg(Lmax+1,t_VEC); pc++;
813 2956052 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
814 968193 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
815 :
816 493931 : *mu = pariINFINITY;
817 493931 : g = real_0_bit(-bit);
818 493930 : TWO = real2n(1, DEFAULTPREC);
819 493935 : av2 = avma;
820 493935 : RU = gen_1;
821 1737436 : for (i=0; i<K; i++)
822 : {
823 1243504 : if (i) {
824 749577 : GEN z = RU;
825 3442478 : for (j=1; j<n; j++)
826 : {
827 2692898 : gel(pc,j) = gmul(gel(q,j),z);
828 2692861 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
829 : }
830 749580 : gel(pc,n) = gmul(gel(q,n),z);
831 : }
832 :
833 1243501 : fft(Omega,pc,A,1,Lmax,1);
834 1243518 : if (polreal && i>0 && i<K-1)
835 1141896 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
836 : else
837 8105687 : for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
838 1243133 : RU = gmul(RU, prim);
839 1243503 : if (gc_needed(av,1))
840 : {
841 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
842 0 : (void)gc_all(av2,2, &g,&RU);
843 : }
844 : }
845 493932 : *gamma = mydbllog2r(divru(g,NN));
846 493931 : *LMAX = Lmax; set_avma(av);
847 493931 : }
848 :
849 : /* NN is a multiple of Lmax */
850 : static void
851 502324 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
852 : {
853 : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
854 502324 : long n = degpol(p), i, j, K;
855 : pari_sp ltop;
856 :
857 502324 : initdft(&Omega, &prim, NN, Lmax, bit);
858 502325 : RU = cgetg(n+2,t_VEC) + 1;
859 :
860 502324 : K = NN/Lmax; if (polreal) K = K/2+1;
861 502324 : q = mygprec(p,bit);
862 502324 : qd = RgX_deriv(q);
863 :
864 502318 : A = cgetg(Lmax+1,t_VEC); A++;
865 502320 : B = cgetg(Lmax+1,t_VEC); B++;
866 502321 : C = cgetg(Lmax+1,t_VEC); C++;
867 502321 : pc = cgetg(Lmax+1,t_VEC); pc++;
868 502322 : pd = cgetg(Lmax+1,t_VEC); pd++;
869 1018722 : gel(pc,0) = gel(q,2); for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
870 1521045 : gel(pd,0) = gel(qd,2); for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
871 :
872 502323 : ltop = avma;
873 502323 : W = cgetg(k+1,t_VEC);
874 502323 : U = cgetg(k+1,t_VEC);
875 1201892 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
876 :
877 502322 : gel(RU,0) = gen_1;
878 502322 : prim2 = gen_1;
879 1530185 : for (i=0; i<K; i++)
880 : {
881 1027860 : gel(RU,1) = prim2;
882 4444729 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
883 : /* RU[j] = prim^(ij)= prim2^j */
884 :
885 4444705 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
886 1027828 : fft(Omega,pd,A,1,Lmax,1);
887 5472557 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
888 1027839 : fft(Omega,pc,B,1,Lmax,1);
889 7664397 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
890 7664383 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
891 1027762 : fft(Omega,B,A,1,Lmax,1);
892 1027864 : fft(Omega,C,B,1,Lmax,1);
893 :
894 1027864 : if (polreal) /* p has real coefficients */
895 : {
896 796851 : if (i>0 && i<K-1)
897 : {
898 102741 : for (j=1; j<=k; j++)
899 : {
900 86101 : gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
901 86101 : gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
902 : }
903 : }
904 : else
905 : {
906 1830391 : for (j=1; j<=k; j++)
907 : {
908 1050193 : gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
909 1050168 : gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
910 : }
911 : }
912 : }
913 : else
914 : {
915 604896 : for (j=1; j<=k; j++)
916 : {
917 373883 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
918 373880 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
919 : }
920 : }
921 1027851 : prim2 = gmul(prim2,prim);
922 1027853 : (void)gc_all(ltop,3, &W,&U,&prim2);
923 : }
924 :
925 1201890 : for (i=1; i<=k; i++)
926 : {
927 699571 : aux=gel(W,i);
928 1098286 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
929 699572 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
930 : }
931 1201884 : for (i=0; i<k; i++)
932 : {
933 699566 : aux=gel(U,k-i);
934 1098277 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
935 699563 : gel(H,i+2) = gdivgu(aux,NN);
936 : }
937 502318 : }
938 :
939 : #define NEWTON_MAX 10
940 : static GEN
941 2460145 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
942 : {
943 2460145 : GEN H = HH, D, aux;
944 2460145 : pari_sp ltop = avma;
945 : long error, i, bit1, bit2;
946 :
947 2460145 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
948 2460119 : bit2 = bit + Sbit;
949 4496483 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
950 : {
951 2036368 : if (gc_needed(ltop,1))
952 : {
953 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
954 0 : (void)gc_all(ltop,2, &D,&H);
955 : }
956 2036368 : bit1 = -error + Sbit;
957 2036368 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
958 2036365 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
959 :
960 2036384 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
961 2036384 : H = RgX_add(mygprec(H,bit1), aux);
962 2036322 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
963 2036366 : error = gexpo(D); if (error < -bit1) error = -bit1;
964 : }
965 2460115 : if (error > -bit/2) return NULL; /* FAIL */
966 2459791 : return gc_GEN(ltop,H);
967 : }
968 :
969 : /* return 0 if fails, 1 else */
970 : static long
971 502319 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
972 : {
973 502319 : GEN f0, FF, GG, r, HH = H;
974 502319 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
975 502319 : pari_sp av = avma;
976 :
977 502319 : FF = *F; GG = RgX_divrem(p, FF, &r);
978 502326 : error = gexpo(r); if (error <= -bit) error = 1-bit;
979 502326 : normF = gexpo(FF);
980 502326 : normG = gexpo(GG);
981 502325 : enh = gexpo(H); if (enh < 0) enh = 0;
982 502325 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
983 502325 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
984 502325 : bit2 = bit + Sbit;
985 2962137 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
986 : {
987 2460143 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
988 2460143 : if (gc_needed(av,1))
989 : {
990 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
991 0 : (void)gc_all(av,4, &FF,&GG,&r,&HH);
992 : }
993 :
994 2460143 : bit1 = -error + Sbit2;
995 2460143 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
996 : 1-error, Sbit2);
997 2460145 : if (!HH) return 0; /* FAIL */
998 :
999 2459821 : bit1 = -error + Sbit;
1000 2459821 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1001 2459783 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1002 :
1003 2459816 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1004 2459816 : FF = gadd(mygprec(FF,bit1),f0);
1005 :
1006 2459781 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1007 2459781 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1008 2459819 : error = gexpo(r); if (error < -bit1) error = -bit1;
1009 : }
1010 501994 : if (error>-bit) return 0; /* FAIL */
1011 493926 : *F = FF; *G = GG; return 1;
1012 : }
1013 :
1014 : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1015 : where cd is the leading coefficient of p */
1016 : static void
1017 493934 : split_fromU(GEN p, long k, double delta, long bit,
1018 : GEN *F, GEN *G, double param, double param2)
1019 : {
1020 : GEN pp, FF, GG, H;
1021 493934 : long n = degpol(p), NN, bit2, Lmax;
1022 493934 : int polreal = isreal(p);
1023 : pari_sp ltop;
1024 : double mu, gamma;
1025 :
1026 493934 : pp = gdiv(p, gel(p,2+n));
1027 493930 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1028 :
1029 493931 : H = cgetg(k+2,t_POL); H[1] = p[1];
1030 493931 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1031 493931 : gel(FF,k+2) = gen_1;
1032 :
1033 493931 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1034 493931 : NN *= Lmax; ltop = avma;
1035 : for(;;)
1036 : {
1037 502323 : bit2 = (long)(((double)NN*delta-mu)/M_LN2) + gexpo(pp) + 8;
1038 502324 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1039 502319 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1040 8392 : NN <<= 1; set_avma(ltop);
1041 : }
1042 493934 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1043 493929 : }
1044 :
1045 : static void
1046 493933 : optimize_split(GEN p, long k, double delta, long bit,
1047 : GEN *F, GEN *G, double param, double param2)
1048 : {
1049 493933 : long n = degpol(p);
1050 : GEN FF, GG;
1051 :
1052 493933 : if (k <= n/2)
1053 383099 : split_fromU(p,k,delta,bit,F,G,param,param2);
1054 : else
1055 : {
1056 110834 : split_fromU(RgX_recip_i(p),n-k,delta,bit,&FF,&GG,param,param2);
1057 110834 : *F = RgX_recip_i(GG);
1058 110834 : *G = RgX_recip_i(FF);
1059 : }
1060 493930 : }
1061 :
1062 : /********************************************************************/
1063 : /** **/
1064 : /** SEARCH FOR SEPARATING CIRCLE **/
1065 : /** **/
1066 : /********************************************************************/
1067 :
1068 : /* return p(2^e*x) *2^(-n*e) */
1069 : static void
1070 0 : scalepol2n(GEN p, long e)
1071 : {
1072 0 : long i,n=lg(p)-1;
1073 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1074 0 : }
1075 :
1076 : /* returns p(x/R)*R^n; assume R is at the correct accuracy */
1077 : static GEN
1078 4289005 : scalepol(GEN p, GEN R, long bit)
1079 4289005 : { return RgX_rescale(mygprec(p, bit), R); }
1080 :
1081 : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1082 : static GEN
1083 1403772 : conformal_basecase(GEN p, GEN a)
1084 : {
1085 : GEN z, r, ma, ca;
1086 1403772 : long i, n = degpol(p);
1087 : pari_sp av;
1088 :
1089 1403770 : if (n <= 0) return p;
1090 1403770 : ma = gneg(a); ca = conj_i(a);
1091 1403772 : av = avma;
1092 1403772 : z = deg1pol_shallow(ca, gen_m1, 0);
1093 1403770 : r = scalarpol_shallow(gel(p,2+n), 0);
1094 3640897 : for (i=n-1; ; i--)
1095 : {
1096 3640897 : r = RgX_addmulXn_shallow(r, gmul(ma,r), 1); /* r *= (X - a) */
1097 3640857 : r = gadd(r, gmul(z, gel(p,2+i)));
1098 3640868 : if (i == 0) return gc_upto(av, r);
1099 2237107 : z = RgX_addmulXn_shallow(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1100 2237130 : if (gc_needed(av,2))
1101 : {
1102 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol (%ld/%ld)",n-i, n);
1103 0 : (void)gc_all(av,2, &r,&z);
1104 : }
1105 : }
1106 : }
1107 : static GEN
1108 1403890 : conformal_pol(GEN p, GEN a)
1109 : {
1110 1403890 : pari_sp av = avma;
1111 1403890 : long d, nR, n = degpol(p), v;
1112 : GEN Q, R, S, T;
1113 1403892 : if (n < 35) return conformal_basecase(p, a);
1114 119 : d = (n+1) >> 1; v = varn(p);
1115 119 : Q = RgX_shift_shallow(p, -d);
1116 119 : R = RgXn_red_shallow(p, d);
1117 119 : Q = conformal_pol(Q, a);
1118 119 : R = conformal_pol(R, a);
1119 119 : S = gpowgs(deg1pol_shallow(gen_1, gneg(a), v), d);
1120 119 : T = RgX_recip_i(S);
1121 119 : if (typ(a) == t_COMPLEX) T = gconj(T);
1122 119 : if (odd(d)) T = RgX_neg(T);
1123 : /* S = (X - a)^d, T = (conj(a) X - 1)^d */
1124 119 : nR = n - degpol(R) - d; /* >= 0 */
1125 119 : if (nR) T = RgX_mul(T, gpowgs(deg1pol_shallow(gconj(a), gen_m1, v), nR));
1126 119 : return gc_upto(av, RgX_add(RgX_mul(Q, S), RgX_mul(R, T)));
1127 : }
1128 :
1129 : static const double UNDEF = -100000.;
1130 :
1131 : static double
1132 493925 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1133 : {
1134 493925 : long i, n = degpol(p);
1135 : double lrho, lrmin, lrmax;
1136 493925 : if (k > 1)
1137 : {
1138 282489 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1139 207029 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1140 : }
1141 : else /* k=1 */
1142 286896 : lrmin = logmin_modulus(p,aux);
1143 493934 : radii[k] = lrmin;
1144 :
1145 493934 : if (k+1<n)
1146 : {
1147 590884 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1148 399674 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1149 : }
1150 : else /* k+1=n */
1151 94260 : lrmax = logmax_modulus(p,aux);
1152 493931 : radii[k+1] = lrmax;
1153 :
1154 493931 : lrho = radii[k];
1155 824143 : for (i=k-1; i>=1; i--)
1156 : {
1157 330212 : if (radii[i] == UNDEF || radii[i] > lrho)
1158 242313 : radii[i] = lrho;
1159 : else
1160 87899 : lrho = radii[i];
1161 : }
1162 493931 : lrho = radii[k+1];
1163 1637979 : for (i=k+1; i<=n; i++)
1164 : {
1165 1144048 : if (radii[i] == UNDEF || radii[i] < lrho)
1166 567135 : radii[i] = lrho;
1167 : else
1168 576913 : lrho = radii[i];
1169 : }
1170 493931 : *delta = (lrmax - lrmin) / 2;
1171 493931 : if (*delta > 1.) *delta = 1.;
1172 493931 : return (lrmin + lrmax) / 2;
1173 : }
1174 :
1175 : static void
1176 493931 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1177 : {
1178 493931 : double t, param = 0., param2 = 0.;
1179 : long i;
1180 2462055 : for (i=1; i<=n; i++)
1181 : {
1182 1968142 : radii[i] -= lrho;
1183 1968142 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1184 1968124 : param += t; if (t > 1.) param2 += log2(t);
1185 : }
1186 493913 : *par = param; *par2 = param2;
1187 493913 : }
1188 :
1189 : /* apply the conformal mapping then split from U */
1190 : static void
1191 467884 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1192 : double aux, GEN *F,GEN *G)
1193 : {
1194 467884 : long bit2, n = degpol(p), i;
1195 467883 : pari_sp ltop = avma, av;
1196 : GEN q, FF, GG, a, R;
1197 : double lrho, delta, param, param2;
1198 : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1199 467883 : bit2 = bit + (long)(n*3.4848775) + 1;
1200 467883 : a = sqrtr_abs( utor(3, precdbl(MEDDEFAULTPREC)) );
1201 467884 : a = divrs(a, -6);
1202 467887 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1203 :
1204 467883 : av = avma;
1205 467883 : q = conformal_pol(mygprec(p,bit2), a);
1206 2288363 : for (i=1; i<=n; i++)
1207 1820485 : if (radii[i] != UNDEF) /* update array radii */
1208 : {
1209 1540937 : pari_sp av2 = avma;
1210 1540937 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1211 : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1212 1540891 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1213 1540921 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1214 1540908 : set_avma(av2);
1215 : }
1216 467878 : lrho = logradius(radii, q,k,aux/10., &delta);
1217 467884 : update_radius(n, radii, lrho, ¶m, ¶m2);
1218 :
1219 467881 : bit2 += (long)(n * fabs(lrho)/M_LN2 + 1.);
1220 467881 : R = mygprec(dblexp(-lrho), bit2);
1221 467886 : q = scalepol(q,R,bit2);
1222 467885 : (void)gc_all(av,2, &q,&R);
1223 :
1224 467886 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1225 467883 : bit2 += n; R = invr(R);
1226 467884 : FF = scalepol(FF,R,bit2);
1227 467885 : GG = scalepol(GG,R,bit2);
1228 :
1229 467884 : a = mygprec(a,bit2);
1230 467885 : FF = conformal_pol(FF,a);
1231 467885 : GG = conformal_pol(GG,a);
1232 :
1233 467886 : a = invr(subsr(1, gnorm(a)));
1234 467884 : FF = RgX_Rg_mul(FF, powru(a,k));
1235 467887 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1236 :
1237 467886 : *F = mygprec(FF,bit+n);
1238 467886 : *G = mygprec(GG,bit+n); (void)gc_all(ltop,2, F,G);
1239 467887 : }
1240 :
1241 : /* split p, this time without scaling. returns in F and G two polynomials
1242 : * such that |p-FG|< 2^(-bit)|p| */
1243 : static void
1244 493933 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1245 : {
1246 : GEN q, FF, GG, R;
1247 : double aux, delta, param, param2;
1248 493933 : long n = degpol(p), i, j, k, bit2;
1249 : double lrmin, lrmax, lrho, *radii;
1250 :
1251 493933 : radii = (double*) stack_malloc_align((n+1) * sizeof(double), sizeof(double));
1252 :
1253 1474263 : for (i=2; i<n; i++) radii[i] = UNDEF;
1254 493933 : aux = thickness/(double)(4 * n);
1255 493933 : lrmin = logmin_modulus(p, aux);
1256 493924 : lrmax = logmax_modulus(p, aux);
1257 493926 : radii[1] = lrmin;
1258 493926 : radii[n] = lrmax;
1259 493926 : i = 1; j = n;
1260 493926 : lrho = (lrmin + lrmax) / 2;
1261 493926 : k = dual_modulus(p, lrho, aux, 1);
1262 493933 : if (5*k < n || (n < 2*k && 5*k < 4*n))
1263 77972 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1264 : else
1265 415961 : { lrmin = lrho; i=k; radii[i] = lrho; }
1266 1012379 : while (j > i+1)
1267 : {
1268 518448 : if (i+j == n+1)
1269 372201 : lrho = (lrmin + lrmax) / 2;
1270 : else
1271 : {
1272 146247 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1273 146247 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1274 116667 : else lrho = lrmin * kappa + lrmax;
1275 146247 : lrho /= 1+kappa;
1276 : }
1277 518448 : aux = (lrmax - lrmin) / (4*(j-i));
1278 518448 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1279 518446 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
1280 387238 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1281 : else
1282 131208 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1283 : }
1284 493931 : aux = lrmax - lrmin;
1285 :
1286 493931 : if (ctr)
1287 : {
1288 467885 : lrho = (lrmax + lrmin) / 2;
1289 2288407 : for (i=1; i<=n; i++)
1290 1820522 : if (radii[i] != UNDEF) radii[i] -= lrho;
1291 :
1292 467885 : bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
1293 467885 : R = mygprec(dblexp(-lrho), bit2);
1294 467882 : q = scalepol(p,R,bit2);
1295 467884 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1296 : }
1297 : else
1298 : {
1299 26046 : lrho = logradius(radii, p, k, aux/10., &delta);
1300 26047 : update_radius(n, radii, lrho, ¶m, ¶m2);
1301 :
1302 26047 : bit2 = bit + (long)(n * fabs(lrho)/M_LN2 + 1.);
1303 26047 : R = mygprec(dblexp(-lrho), bit2);
1304 26047 : q = scalepol(p,R,bit2);
1305 26047 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1306 : }
1307 493934 : bit += n;
1308 493934 : bit2 += n; R = invr(mygprec(R,bit2));
1309 493930 : *F = mygprec(scalepol(FF,R,bit2), bit);
1310 493932 : *G = mygprec(scalepol(GG,R,bit2), bit);
1311 493932 : }
1312 :
1313 : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1314 : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1315 : * where the maximum modulus of the roots of p is <=1.
1316 : * Assume sum of roots is 0. */
1317 : static void
1318 467887 : split_1(GEN p, long bit, GEN *F, GEN *G)
1319 : {
1320 467887 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1321 : GEN ctr, q, qq, FF, GG, v, gr, r, newq;
1322 : double lrmin, lrmax, lthick;
1323 467886 : const double LOG3 = 1.098613;
1324 :
1325 467886 : lrmax = logmax_modulus(p, 0.01);
1326 467886 : gr = mygprec(dblexp(-lrmax), bit2);
1327 467884 : q = scalepol(p,gr,bit2);
1328 :
1329 467886 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1330 467886 : v = cgetg(5,t_VEC);
1331 467886 : gel(v,1) = gen_2;
1332 467886 : gel(v,2) = gen_m2;
1333 467886 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1334 467886 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1335 467886 : q = mygprec(q,bit2); lthick = 0;
1336 467886 : newq = ctr = NULL; /* -Wall */
1337 467886 : imax = polreal? 3: 4;
1338 840486 : for (i=1; i<=imax; i++)
1339 : {
1340 833866 : qq = RgX_Rg_translate(q, gel(v,i));
1341 833870 : lrmin = logmin_modulus(qq,0.05);
1342 833861 : if (LOG3 > lrmin + lthick)
1343 : {
1344 821712 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1345 821716 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1346 : }
1347 833865 : if (lthick > M_LN2) break;
1348 423643 : if (polreal && i==2 && lthick > LOG3 - M_LN2) break;
1349 : }
1350 467885 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/M_LN2 + 1);
1351 467886 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1352 467886 : r = gneg(mygprec(ctr,bit2));
1353 467886 : FF = RgX_Rg_translate(FF,r);
1354 467887 : GG = RgX_Rg_translate(GG,r);
1355 :
1356 467886 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1357 467885 : *F = scalepol(FF,gr,bit2);
1358 467886 : *G = scalepol(GG,gr,bit2);
1359 467886 : }
1360 :
1361 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1362 : where the maximum modulus of the roots of p is < 0.5 */
1363 : static int
1364 468207 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1365 : {
1366 : GEN q, b;
1367 468207 : long n = degpol(p), k, bit2, eq;
1368 468208 : double aux0 = dbllog2(gel(p,n+2)); /* != -oo */
1369 468207 : double aux1 = dbllog2(gel(p,n+1)), aux;
1370 :
1371 468207 : if (aux1 == -pariINFINITY) /* p1 = 0 */
1372 9893 : aux = 0;
1373 : else
1374 : {
1375 458314 : aux = aux1 - aux0; /* log2(p1/p0) */
1376 : /* beware double overflow */
1377 458314 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
1378 458314 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1379 : }
1380 468207 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1381 468207 : q = mygprec(p,bit2);
1382 468211 : if (aux1 == -pariINFINITY) b = NULL;
1383 : else
1384 : {
1385 458318 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1386 458316 : q = RgX_Rg_translate(q,b);
1387 : }
1388 468211 : gel(q,n+1) = gen_0; eq = gexpo(q);
1389 468211 : k = 0;
1390 468765 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1391 468637 : || gequal0(gel(q,k+2)))) k++;
1392 468211 : if (k > 0)
1393 : {
1394 324 : if (k > n/2) k = n/2;
1395 324 : bit2 += k<<1;
1396 324 : *F = pol_xn(k, 0);
1397 324 : *G = RgX_shift_shallow(q, -k);
1398 : }
1399 : else
1400 : {
1401 467887 : split_1(q,bit2,F,G);
1402 467886 : bit2 = bit + gexpo(*F) + gexpo(*G) - gexpo(p) + (long)aux+1;
1403 467886 : *F = mygprec(*F,bit2);
1404 : }
1405 468210 : *G = mygprec(*G,bit2);
1406 468210 : if (b)
1407 : {
1408 458318 : GEN mb = mygprec(gneg(b), bit2);
1409 458318 : *F = RgX_Rg_translate(*F, mb);
1410 458317 : *G = RgX_Rg_translate(*G, mb);
1411 : }
1412 468210 : return 1;
1413 : }
1414 :
1415 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1416 : * Assume max_modulus(p) < 2 */
1417 : static void
1418 468207 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1419 : {
1420 : GEN FF, GG;
1421 : long n, bit2, normp;
1422 :
1423 468207 : if (split_0_2(p,bit,F,G)) return;
1424 :
1425 0 : normp = gexpo(p);
1426 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1427 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1428 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1429 0 : scalepol2n(FF,-2);
1430 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1431 0 : *F = mygprec(FF,bit2);
1432 0 : *G = mygprec(GG,bit2);
1433 : }
1434 :
1435 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1436 : static void
1437 494256 : split_0(GEN p, long bit, GEN *F, GEN *G)
1438 : {
1439 494256 : const double LOG1_9 = 0.6418539;
1440 494256 : long n = degpol(p), k = 0;
1441 : GEN q;
1442 :
1443 494256 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1444 494256 : if (k > 0)
1445 : {
1446 0 : if (k > n/2) k = n/2;
1447 0 : *F = pol_xn(k, 0);
1448 0 : *G = RgX_shift_shallow(p, -k);
1449 : }
1450 : else
1451 : {
1452 494256 : double lr = logmax_modulus(p, 0.05);
1453 494255 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1454 : else
1455 : {
1456 436758 : q = RgX_recip_i(p);
1457 436759 : lr = logmax_modulus(q,0.05);
1458 436757 : if (lr < LOG1_9)
1459 : {
1460 410710 : split_0_1(q, bit, F, G);
1461 410713 : *F = RgX_recip_i(*F);
1462 410713 : *G = RgX_recip_i(*G);
1463 : }
1464 : else
1465 26047 : split_2(p,bit,NULL, 1.2837,F,G);
1466 : }
1467 : }
1468 494257 : }
1469 :
1470 : /********************************************************************/
1471 : /** **/
1472 : /** ERROR ESTIMATE FOR THE ROOTS **/
1473 : /** **/
1474 : /********************************************************************/
1475 :
1476 : static GEN
1477 1899488 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1478 : {
1479 1899488 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1480 : long i, j;
1481 :
1482 1899488 : d = cgetg(n+1,t_VEC);
1483 12131944 : for (i=1; i<=n; i++)
1484 : {
1485 10232752 : if (i!=k)
1486 : {
1487 8333372 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1488 8333005 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1489 : }
1490 : }
1491 1899192 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1492 1899500 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1493 1899500 : eps = mulrr(rho, shatzle);
1494 1899458 : aux = shiftr(powru(rho,n), err);
1495 :
1496 5765723 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
1497 : {
1498 3866363 : GEN prod = NULL; /* 1. */
1499 3866363 : long m = n;
1500 3866363 : epsbis = mulrr(eps, dbltor(1.25));
1501 26551290 : for (i=1; i<=n; i++)
1502 : {
1503 22684656 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1504 : {
1505 18778861 : GEN dif = subrr(gel(d,i),eps);
1506 18777174 : prod = prod? mulrr(prod, dif): dif;
1507 18778128 : m--;
1508 : }
1509 : }
1510 3866634 : eps2 = prod? divrr(aux, prod): aux;
1511 3866397 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1512 3866397 : rap = divrr(eps,eps2); eps = eps2;
1513 : }
1514 1899304 : return eps;
1515 : }
1516 :
1517 : /* round a complex or real number x to an absolute value of 2^(-bit) */
1518 : static GEN
1519 4300823 : mygprec_absolute(GEN x, long bit)
1520 : {
1521 : long e;
1522 : GEN y;
1523 :
1524 4300823 : switch(typ(x))
1525 : {
1526 2953025 : case t_REAL:
1527 2953025 : e = expo(x) + bit;
1528 2953025 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1529 1217989 : case t_COMPLEX:
1530 1217989 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1531 1183424 : y = cgetg(3,t_COMPLEX);
1532 1183429 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1533 1183428 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1534 1183434 : return y;
1535 129809 : default: return x;
1536 : }
1537 : }
1538 :
1539 : static long
1540 530770 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1541 : {
1542 530770 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1543 : GEN sigma, shatzle;
1544 :
1545 530770 : err += (long)log2((double)n) + 1;
1546 530770 : if (err > -2) return 0;
1547 530770 : sigma = real2n(-err, LOWDEFAULTPREC);
1548 : /* 2 / ((s - 1)^(1/n) - 1) */
1549 530768 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1550 2430248 : for (i=1; i<=n; i++)
1551 : {
1552 1899482 : pari_sp av = avma;
1553 1899482 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1554 1899301 : long e = gexpo(x);
1555 1899357 : set_avma(av); if (e > e_max) e_max = e;
1556 1899441 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1557 : }
1558 530766 : return e_max;
1559 : }
1560 :
1561 : /********************************************************************/
1562 : /** **/
1563 : /** MAIN **/
1564 : /** **/
1565 : /********************************************************************/
1566 : static GEN
1567 1603348 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1568 :
1569 : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1570 : * returns prod (x-roots_pol[i]) */
1571 : static GEN
1572 1519277 : split_complete(GEN p, long bit, GEN roots_pol)
1573 : {
1574 1519277 : long n = degpol(p);
1575 : pari_sp ltop;
1576 : GEN p1, F, G, a, b, m1, m2;
1577 :
1578 1519275 : if (n == 1)
1579 : {
1580 446683 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1581 446684 : (void)append_clone(roots_pol,a); return p;
1582 : }
1583 1072592 : ltop = avma;
1584 1072592 : if (n == 2)
1585 : {
1586 578336 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1587 578327 : F = gsqrt(F, nbits2prec(bit));
1588 578332 : p1 = ginv(gmul2n(gel(p,4),1));
1589 578332 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1590 578335 : b = gmul(gsub(F,gel(p,3)), p1);
1591 578338 : a = append_clone(roots_pol,a);
1592 578340 : b = append_clone(roots_pol,b); set_avma(ltop);
1593 578340 : a = mygprec(a, 3*bit);
1594 578339 : b = mygprec(b, 3*bit);
1595 578339 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1596 : }
1597 494256 : split_0(p,bit,&F,&G);
1598 494257 : m1 = split_complete(F,bit,roots_pol);
1599 494256 : m2 = split_complete(G,bit,roots_pol);
1600 494256 : return gc_upto(ltop, gmul(m1,m2));
1601 : }
1602 :
1603 : static GEN
1604 6961664 : quicktofp(GEN x)
1605 : {
1606 6961664 : const long prec = DEFAULTPREC;
1607 6961664 : switch(typ(x))
1608 : {
1609 6940207 : case t_INT: return itor(x, prec);
1610 9064 : case t_REAL: return rtor(x, prec);
1611 0 : case t_FRAC: return fractor(x, prec);
1612 12395 : case t_COMPLEX: {
1613 12395 : GEN a = gel(x,1), b = gel(x,2);
1614 : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1615 12395 : if (isintzero(a)) return cxcompotor(b, prec);
1616 12353 : if (isintzero(b)) return cxcompotor(a, prec);
1617 12353 : a = cxcompotor(a, prec);
1618 12353 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1619 : }
1620 0 : default: pari_err_TYPE("quicktofp",x);
1621 : return NULL;/*LCOV_EXCL_LINE*/
1622 : }
1623 :
1624 : }
1625 :
1626 : /* bound log_2 |largest root of p| (Fujiwara's bound) */
1627 : double
1628 2265233 : fujiwara_bound(GEN p)
1629 : {
1630 2265233 : pari_sp av = avma;
1631 2265233 : long i, n = degpol(p);
1632 : GEN cc;
1633 : double loglc, Lmax;
1634 :
1635 2265231 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1636 2265231 : loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
1637 2265198 : cc = gel(p, 2);
1638 2265198 : if (gequal0(cc))
1639 783534 : Lmax = -pariINFINITY-1;
1640 : else
1641 1481690 : Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
1642 7394014 : for (i = 1; i < n; i++)
1643 : {
1644 5128824 : GEN y = gel(p,i+2);
1645 : double L;
1646 5128824 : if (gequal0(y)) continue;
1647 3214827 : L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
1648 3214826 : if (L > Lmax) Lmax = L;
1649 : }
1650 2265190 : return gc_double(av, Lmax+1);
1651 : }
1652 :
1653 : /* Fujiwara's bound, real roots. Based on the following remark: if
1654 : * p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
1655 : * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
1656 : * of q is a bound for the positive roots of p. */
1657 : double
1658 1406833 : fujiwara_bound_real(GEN p, long sign)
1659 : {
1660 1406833 : pari_sp av = avma;
1661 : GEN x;
1662 1406833 : long n = degpol(p), i, signodd, signeven;
1663 1406833 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1664 1406833 : x = shallowcopy(p);
1665 1406833 : if (gsigne(gel(x, n+2)) > 0)
1666 1406812 : { signeven = 1; signodd = sign; }
1667 : else
1668 21 : { signeven = -1; signodd = -sign; }
1669 5780124 : for (i = 0; i < n; i++)
1670 : {
1671 4373291 : if ((n - i) % 2)
1672 2490993 : { if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0; }
1673 : else
1674 1882298 : { if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0; }
1675 : }
1676 1406833 : return gc_double(av, fujiwara_bound(x));
1677 : }
1678 :
1679 : static GEN
1680 2161482 : mygprecrc_special(GEN x, long prec, long e)
1681 : {
1682 : GEN y;
1683 2161482 : switch(typ(x))
1684 : {
1685 37263 : case t_REAL:
1686 37263 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1687 36115 : return (prec > realprec(x))? rtor(x, prec): x;
1688 13678 : case t_COMPLEX:
1689 13678 : y = cgetg(3,t_COMPLEX);
1690 13678 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1691 13678 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1692 13678 : return y;
1693 2110541 : default: return x;
1694 : }
1695 : }
1696 :
1697 : /* like mygprec but keep at least the same precision as before */
1698 : static GEN
1699 530768 : mygprec_special(GEN x, long bit)
1700 : {
1701 530768 : long e = gexpo(x) - bit, prec = nbits2prec(bit);
1702 530766 : switch(typ(x))
1703 : {
1704 2664892 : case t_POL: pari_APPLY_pol_normalized(mygprecrc_special(gel(x,i),prec,e));
1705 0 : default: return mygprecrc_special(x,prec,e);
1706 : }
1707 : }
1708 :
1709 : static GEN
1710 394194 : fix_roots1(GEN R)
1711 : {
1712 394194 : long i, l = lg(R);
1713 394194 : GEN v = cgetg(l, t_VEC);
1714 1752279 : for (i=1; i < l; i++) { GEN r = gel(R,i); gel(v,i) = gcopy(r); gunclone(r); }
1715 394195 : return v;
1716 : }
1717 : static GEN
1718 530769 : fix_roots(GEN R, long h, long bit)
1719 : {
1720 : long i, j, c, n, prec;
1721 : GEN v, Z, gh;
1722 :
1723 530769 : if (h == 1) return fix_roots1(R);
1724 136575 : prec = nbits2prec(bit); Z = grootsof1(h, prec); gh = utoipos(h);
1725 136575 : n = lg(R)-1; v = cgetg(h*n + 1, t_VEC);
1726 381843 : for (c = i = 1; i <= n; i++)
1727 : {
1728 245269 : GEN s, r = gel(R,i);
1729 245269 : s = (h == 2)? gsqrt(r, prec): gsqrtn(r, gh, NULL, prec);
1730 786728 : for (j = 1; j <= h; j++) gel(v, c++) = gmul(s, gel(Z,j));
1731 245257 : gunclone(r);
1732 : }
1733 136574 : return v;
1734 : }
1735 :
1736 : static GEN
1737 529725 : all_roots(GEN p, long bit)
1738 : {
1739 529725 : long bit2, i, e, h, n = degpol(p), elc = gexpo(leading_coeff(p));
1740 529724 : GEN q, R, m, pd = RgX_deflate_max(p, &h);
1741 529724 : double fb = fujiwara_bound(pd);
1742 : pari_sp av;
1743 :
1744 529725 : if (fb < 0) fb = 0;
1745 529725 : bit2 = bit + maxss(gexpo(p), 0) + (long)ceil(log2(n / h) + 2 * fb);
1746 530768 : for (av = avma, i = 1, e = 0;; i++, set_avma(av))
1747 : {
1748 530768 : R = vectrunc_init(n+1);
1749 530768 : bit2 += e + (n << i);
1750 530768 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1751 530766 : q[1] = evalsigne(1)|evalvarn(0);
1752 530766 : m = split_complete(q, bit2, R);
1753 530769 : R = fix_roots(R, h, bit2);
1754 530769 : q = mygprec_special(pd,bit2);
1755 530768 : q[1] = evalsigne(1)|evalvarn(0);
1756 530768 : e = gexpo(RgX_sub(q, m)) - elc + (long)log2((double)n) + 1;
1757 530770 : if (e < 0)
1758 : {
1759 530770 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -oo */
1760 530770 : e = bit + a_posteriori_errors(p, R, e);
1761 530766 : if (e < 0) return R;
1762 : }
1763 1044 : if (DEBUGLEVEL)
1764 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1765 : }
1766 : }
1767 :
1768 : INLINE int
1769 931728 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1770 :
1771 : static int
1772 239746 : isexactpol(GEN p)
1773 : {
1774 239746 : long i,n = degpol(p);
1775 1162880 : for (i=0; i<=n; i++)
1776 931728 : if (!isexactscalar(gel(p,i+2))) return 0;
1777 231152 : return 1;
1778 : }
1779 :
1780 : /* p(0) != 0 [for efficiency] */
1781 : static GEN
1782 231152 : solve_exact_pol(GEN p, long bit)
1783 : {
1784 231152 : long i, j, k, m, n = degpol(p), iroots = 0;
1785 231152 : GEN ex, factors, v = zerovec(n);
1786 :
1787 231152 : factors = ZX_squff(Q_primpart(p), &ex);
1788 462304 : for (i=1; i<lg(factors); i++)
1789 : {
1790 231152 : GEN roots_fact = all_roots(gel(factors,i), bit);
1791 231152 : n = degpol(gel(factors,i));
1792 231152 : m = ex[i];
1793 922490 : for (j=1; j<=n; j++)
1794 1382676 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
1795 : }
1796 231152 : return v;
1797 : }
1798 :
1799 : /* return the roots of p with absolute error bit */
1800 : static GEN
1801 239746 : roots_com(GEN q, long bit)
1802 : {
1803 : GEN L, p;
1804 239746 : long v = RgX_valrem_inexact(q, &p);
1805 239746 : int ex = isexactpol(p);
1806 239746 : if (!ex) p = RgX_normalize1(p);
1807 239746 : if (lg(p) == 3)
1808 0 : L = cgetg(1,t_VEC); /* constant polynomial */
1809 : else
1810 239746 : L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
1811 239746 : if (v)
1812 : {
1813 3935 : GEN M, z, t = gel(q,2);
1814 : long i, x, y, l, n;
1815 :
1816 3935 : if (isrationalzero(t)) x = -bit;
1817 : else
1818 : {
1819 7 : n = gexpo(t);
1820 7 : x = n / v; l = degpol(q);
1821 35 : for (i = v; i <= l; i++)
1822 : {
1823 28 : t = gel(q,i+2);
1824 28 : if (isrationalzero(t)) continue;
1825 28 : y = (n - gexpo(t)) / i;
1826 28 : if (y < x) x = y;
1827 : }
1828 : }
1829 3935 : z = real_0_bit(x); l = v + lg(L);
1830 3935 : M = cgetg(l, t_VEC); L -= v;
1831 7933 : for (i = 1; i <= v; i++) gel(M,i) = z;
1832 11826 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
1833 3935 : L = M;
1834 : }
1835 239746 : return L;
1836 : }
1837 :
1838 : static GEN
1839 1201074 : tocomplex(GEN x, long l, long bit)
1840 : {
1841 : GEN y;
1842 1201074 : if (typ(x) == t_COMPLEX)
1843 : {
1844 1181667 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
1845 137459 : x = gel(x,2);
1846 137459 : y = cgetg(3,t_COMPLEX);
1847 137460 : gel(y,1) = real_0_bit(-bit);
1848 137459 : gel(y,2) = mygprecrc(x, l, -bit);
1849 : }
1850 : else
1851 : {
1852 19407 : y = cgetg(3,t_COMPLEX);
1853 19407 : gel(y,1) = mygprecrc(x, l, -bit);
1854 19407 : gel(y,2) = real_0_bit(-bit);
1855 : }
1856 156866 : return y;
1857 : }
1858 :
1859 : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare wrt |Im x| - |Im y|,
1860 : * then Re x - Re y, up to 2^-e absolute error */
1861 : static int
1862 2231962 : cmp_complex_appr(void *E, GEN x, GEN y)
1863 : {
1864 2231962 : long e = (long)E;
1865 : GEN z, xi, yi, xr, yr;
1866 : long sz, sxi, syi;
1867 2231962 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
1868 837858 : else { xr = x; xi = NULL; sxi = 0; }
1869 2231962 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
1870 559612 : else { yr = y; yi = NULL; syi = 0; }
1871 : /* Compare absolute values of imaginary parts */
1872 2231962 : if (!sxi)
1873 : {
1874 857247 : if (syi && expo(yi) >= e) return -1;
1875 : /* |Im x| ~ |Im y| ~ 0 */
1876 : }
1877 1374715 : else if (!syi)
1878 : {
1879 50173 : if (sxi && expo(xi) >= e) return 1;
1880 : /* |Im x| ~ |Im y| ~ 0 */
1881 : }
1882 : else
1883 : {
1884 1324542 : z = addrr_sign(xi, 1, yi, -1); sz = signe(z);
1885 1324519 : if (sz && expo(z) >= e) return (int)sz;
1886 : }
1887 : /* |Im x| ~ |Im y|, sort according to real parts */
1888 1332208 : z = subrr(xr, yr); sz = signe(z);
1889 1332199 : if (sz && expo(z) >= e) return (int)sz;
1890 : /* Re x ~ Re y. Place negative imaginary part before positive */
1891 585775 : return (int) (sxi - syi);
1892 : }
1893 :
1894 : static GEN
1895 529764 : clean_roots(GEN L, long l, long bit, long clean)
1896 : {
1897 529764 : long i, n = lg(L), ex = 5 - bit;
1898 529764 : GEN res = cgetg(n,t_COL);
1899 2429247 : for (i=1; i<n; i++)
1900 : {
1901 1899482 : GEN c = gel(L,i);
1902 1899482 : if (clean && isrealappr(c,ex))
1903 : {
1904 698415 : if (typ(c) == t_COMPLEX) c = gel(c,1);
1905 698415 : c = mygprecrc(c, l, -bit);
1906 : }
1907 : else
1908 1201067 : c = tocomplex(c, l, bit);
1909 1899483 : gel(res,i) = c;
1910 : }
1911 529765 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
1912 529762 : return res;
1913 : }
1914 :
1915 : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
1916 : static GEN
1917 239774 : roots_aux(GEN p, long l, long clean)
1918 : {
1919 239774 : pari_sp av = avma;
1920 : long bit;
1921 : GEN L;
1922 :
1923 239774 : if (typ(p) != t_POL)
1924 : {
1925 21 : if (gequal0(p)) pari_err_ROOTS0("roots");
1926 14 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
1927 7 : return cgetg(1,t_COL); /* constant polynomial */
1928 : }
1929 239753 : if (!signe(p)) pari_err_ROOTS0("roots");
1930 239753 : checkvalidpol(p,"roots");
1931 239746 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
1932 239746 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
1933 239746 : bit = prec2nbits(l);
1934 239746 : L = roots_com(p, bit);
1935 239746 : return gc_GEN(av, clean_roots(L, l, bit, clean));
1936 : }
1937 : GEN
1938 8018 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
1939 : /* clean up roots. If root is real replace it by its real part */
1940 : GEN
1941 231756 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
1942 :
1943 : /* private variant of conjvec. Allow non rational coefficients, shallow
1944 : * function. */
1945 : GEN
1946 84 : polmod_to_embed(GEN x, long prec)
1947 : {
1948 84 : GEN v, T = gel(x,1), A = gel(x,2);
1949 : long i, l;
1950 84 : if (typ(A) != t_POL || varn(A) != varn(T))
1951 : {
1952 7 : checkvalidpol(T,"polmod_to_embed");
1953 7 : return const_col(degpol(T), A);
1954 : }
1955 77 : v = cleanroots(T,prec); l = lg(v);
1956 231 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
1957 77 : return v;
1958 : }
1959 :
1960 : GEN
1961 290018 : QX_complex_roots(GEN p, long l)
1962 : {
1963 290018 : pari_sp av = avma;
1964 : long bit, v;
1965 : GEN L;
1966 :
1967 290018 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
1968 290018 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
1969 290018 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
1970 290018 : bit = prec2nbits(l);
1971 290019 : v = RgX_valrem(p, &p);
1972 290017 : L = lg(p) > 3? all_roots(Q_primpart(p), bit): cgetg(1,t_COL);
1973 290018 : if (v) L = shallowconcat(const_vec(v, real_0_bit(-bit)), L);
1974 290018 : return gc_GEN(av, clean_roots(L, l, bit, 1));
1975 : }
1976 :
1977 : /********************************************************************/
1978 : /** **/
1979 : /** REAL ROOTS OF INTEGER POLYNOMIAL **/
1980 : /** **/
1981 : /********************************************************************/
1982 :
1983 : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1)), P
1984 : * has no rational root. The inversion is implicit (we take coefficients
1985 : * backwards). */
1986 : static long
1987 5990286 : X2XP1(GEN P, GEN *Premapped)
1988 : {
1989 5990286 : const pari_sp av = avma;
1990 5990286 : GEN v = shallowcopy(P);
1991 5990374 : long i, j, nb, s, dP = degpol(P), vlim = dP+2;
1992 :
1993 34756655 : for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
1994 5989952 : s = -signe(gel(v, vlim));
1995 5989952 : vlim--; nb = 0;
1996 16667448 : for (i = 1; i < dP; i++)
1997 : {
1998 14144476 : long s2 = -signe(gel(v, 2));
1999 14144476 : int flag = (s2 == s);
2000 90211890 : for (j = 2; j < vlim; j++)
2001 : {
2002 76067389 : gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2003 76067414 : if (flag) flag = (s2 != signe(gel(v, j+1)));
2004 : }
2005 14144501 : if (s == signe(gel(v, vlim)))
2006 : {
2007 5026919 : if (++nb >= 2) return gc_long(av,2);
2008 3767745 : s = -s;
2009 : }
2010 : /* if flag is set there will be no further sign changes */
2011 12885327 : if (flag && (!Premapped || !nb)) return gc_long(av, nb);
2012 10676974 : vlim--;
2013 10676974 : if (gc_needed(av, 3))
2014 : {
2015 0 : if (DEBUGMEM>1) pari_warn(warnmem, "X2XP1, i = %ld/%ld", i, dP-1);
2016 0 : if (!Premapped) setlg(v, vlim + 2);
2017 0 : v = gc_GEN(av, v);
2018 : }
2019 : }
2020 2522972 : if (vlim >= 2 && s == signe(gel(v, vlim))) nb++;
2021 2522972 : if (Premapped && nb == 1) *Premapped = v; else set_avma(av);
2022 2522615 : return nb;
2023 : }
2024 :
2025 : static long
2026 0 : _intervalcmp(GEN x, GEN y)
2027 : {
2028 0 : if (typ(x) == t_VEC) x = gel(x, 1);
2029 0 : if (typ(y) == t_VEC) y = gel(y, 1);
2030 0 : return gcmp(x, y);
2031 : }
2032 :
2033 : static GEN
2034 11175129 : _gen_nored(void *E, GEN x) { (void)E; return x; }
2035 : static GEN
2036 24653347 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2037 : static GEN
2038 0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
2039 : static GEN
2040 4373544 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2041 : static GEN
2042 6291960 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2043 : static GEN
2044 14443110 : _gen_one(void *E) { (void)E; return gen_1; }
2045 : static GEN
2046 325985 : _gen_zero(void *E) { (void)E; return gen_0; }
2047 :
2048 : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
2049 : _mp_mul, _mp_sqr, _gen_one, _gen_zero };
2050 :
2051 : static GEN
2052 34728030 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
2053 :
2054 : /* Split the polynom P in two parts, whose coeffs have constant sign:
2055 : * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
2056 : * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
2057 : * Pprimep[i] = (i+D) Pp[i]. Return D */
2058 : static long
2059 166629 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
2060 : {
2061 166629 : long i, D, dP = degpol(P), s0 = signe(gel(P,2));
2062 : GEN Pp, Pm, Pprimep, Pprimem;
2063 512091 : for(i=1; i <= dP; i++)
2064 512093 : if (signe(gel(P, i+2)) == -s0) break;
2065 166629 : D = i;
2066 166629 : Pm = cgetg(D + 2, t_POL);
2067 166640 : Pprimem = cgetg(D + 1, t_POL);
2068 166640 : Pp = cgetg(dP-D + 3, t_POL);
2069 166640 : Pprimep = cgetg(dP-D + 3, t_POL);
2070 166641 : Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
2071 678729 : for(i=0; i < D; i++)
2072 : {
2073 512092 : GEN c = gel(P, i+2);
2074 512092 : gel(Pm, i+2) = c;
2075 512092 : if (i) gel(Pprimem, i+1) = mului(i, c);
2076 : }
2077 693303 : for(; i <= dP; i++)
2078 : {
2079 526674 : GEN c = gel(P, i+2);
2080 526674 : gel(Pp, i+2-D) = c;
2081 526674 : gel(Pprimep, i+2-D) = mului(i, c);
2082 : }
2083 166629 : *pPm = normalizepol_lg(Pm, D+2);
2084 166625 : *pPprimem = normalizepol_lg(Pprimem, D+1);
2085 166639 : *pPp = normalizepol_lg(Pp, dP-D+3);
2086 166642 : *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
2087 166643 : return dP - degpol(*pPp);
2088 : }
2089 :
2090 : static GEN
2091 5223923 : bkeval_single_power(long d, GEN V)
2092 : {
2093 5223923 : long mp = lg(V) - 2;
2094 5223923 : if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2095 5223923 : return gel(V, d+1);
2096 : }
2097 :
2098 : static GEN
2099 5223666 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
2100 : {
2101 5223666 : GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
2102 5223501 : GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
2103 5223867 : GEN xa = bkeval_single_power(D, pows);
2104 : GEN r;
2105 5223966 : if (!signe(vp)) return vm;
2106 5223966 : vp = gmul(vp, xa);
2107 5222915 : r = gadd(vp, vm);
2108 5219240 : if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
2109 342053 : return NULL;
2110 4877901 : return r;
2111 : }
2112 :
2113 : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
2114 : static GEN
2115 166643 : splitcauchy(GEN Pp, GEN Pm, long prec)
2116 : {
2117 166643 : GEN S = gel(Pp,2), A = gel(Pm,2);
2118 166643 : long i, lPm = lg(Pm), lPp = lg(Pp);
2119 509022 : for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
2120 526691 : for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
2121 166632 : return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
2122 : }
2123 :
2124 : static GEN
2125 15275 : ZX_deg1root(GEN P, long prec)
2126 : {
2127 15275 : GEN a = gel(P,3), b = gel(P,2);
2128 15275 : if (is_pm1(a))
2129 : {
2130 15275 : b = itor(b, prec); if (signe(a) > 0) togglesign(b);
2131 15275 : return b;
2132 : }
2133 0 : return rdivii(negi(b), a, prec);
2134 : }
2135 :
2136 : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
2137 : * P' has also at most one zero there */
2138 : static GEN
2139 166631 : polsolve(GEN P, long bitprec)
2140 : {
2141 : pari_sp av;
2142 : GEN Pp, Pm, Pprimep, Pprimem, Pprime, Pprime2, ra, rb, rc, Pc;
2143 166631 : long dP = degpol(P), prec = nbits2prec(bitprec);
2144 : long expoold, iter, D, rt, s0, bitaddprec, cprec, PREC;
2145 :
2146 166634 : if (dP == 1) return ZX_deg1root(P, prec);
2147 166634 : Pprime = ZX_deriv(P);
2148 166630 : Pprime2 = ZX_deriv(Pprime);
2149 166629 : bitaddprec = 1 + 2*expu(dP); PREC = prec + nbits2prec(bitaddprec);
2150 166628 : D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
2151 166643 : s0 = signe(gel(P, 2));
2152 166643 : rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
2153 166643 : rb = splitcauchy(Pp, Pm, DEFAULTPREC);
2154 166631 : for (cprec = DEFAULTPREC, expoold = LONG_MAX;;)
2155 0 : {
2156 166631 : GEN pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2157 166636 : Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
2158 166637 : if (!Pc) { cprec += EXTRAPREC64; rb = rtor(rb, cprec); continue; }
2159 166637 : if (signe(Pc) != s0) break;
2160 0 : shiftr_inplace(rb,1);
2161 : }
2162 166637 : for (iter = 0, ra = NULL;;)
2163 1817286 : {
2164 : GEN wdth;
2165 1983923 : iter++;
2166 1983923 : if (ra)
2167 907236 : rc = shiftr(addrr(ra, rb), -1);
2168 : else
2169 1076687 : rc = shiftr(rb, -1);
2170 : for(;;)
2171 0 : {
2172 1984055 : GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2173 1983861 : Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
2174 1983646 : if (Pc) break;
2175 0 : cprec += EXTRAPREC64;
2176 0 : rc = rtor(rc, cprec);
2177 : }
2178 1983646 : if (signe(Pc) == s0)
2179 596687 : ra = rc;
2180 : else
2181 1386959 : rb = rc;
2182 1983646 : if (!ra) continue;
2183 1073568 : wdth = subrr(rb, ra);
2184 1073765 : if (!(iter % 8))
2185 : {
2186 167812 : GEN m1 = poleval(Pprime, ra), M2;
2187 167816 : if (signe(m1) == s0) continue;
2188 166661 : M2 = poleval(Pprime2, rb);
2189 166663 : if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2190 163505 : break;
2191 : }
2192 905953 : else if (gexpo(wdth) <= -bitprec)
2193 3165 : break;
2194 : }
2195 166670 : rc = rb; av = avma;
2196 1371297 : for(;; rc = gc_leaf(av, rc))
2197 1371482 : {
2198 : long exponew;
2199 1538152 : GEN Ppc, dist, rcold = rc;
2200 1538152 : GEN pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2201 1537956 : Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
2202 1537630 : if (Ppc) Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
2203 1537758 : if (!Ppc || !Pc)
2204 : {
2205 342057 : if (cprec >= PREC)
2206 44287 : cprec += EXTRAPREC64;
2207 : else
2208 297770 : cprec = minss(2*cprec, PREC);
2209 342059 : rc = rtor(rc, cprec); continue; /* backtrack one step */
2210 : }
2211 1195701 : dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
2212 1195894 : rc = subrr(rc, dist);
2213 1195303 : if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2214 : {
2215 0 : if (cprec >= PREC) break;
2216 0 : cprec = minss(2*cprec, PREC);
2217 0 : rc = rtor(rcold, cprec); continue; /* backtrack one step */
2218 : }
2219 1195854 : if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2220 975416 : exponew = expo(dist);
2221 975416 : if (exponew < -bitprec - 1)
2222 : {
2223 232390 : if (cprec >= PREC) break;
2224 65756 : cprec = minss(2*cprec, PREC);
2225 65761 : rc = rtor(rc, cprec); continue;
2226 : }
2227 743026 : if (exponew > expoold - 2)
2228 : {
2229 53814 : if (cprec >= PREC) break;
2230 53814 : expoold = LONG_MAX;
2231 53814 : cprec = minss(2*cprec, PREC);
2232 53815 : rc = rtor(rc, cprec); continue;
2233 : }
2234 689212 : expoold = exponew;
2235 : }
2236 166634 : return rtor(rc, prec);
2237 : }
2238 :
2239 : /* Return primpart(P(x / 2)) */
2240 : static GEN
2241 2233749 : ZX_rescale2prim(GEN P)
2242 : {
2243 2233749 : long i, l = lg(P), v, n;
2244 : GEN Q;
2245 2233749 : if (l==2) return pol_0(varn(P));
2246 2233749 : Q = cgetg(l,t_POL); v = vali(gel(P,l-1));
2247 10777203 : for (i = l-2, n = 1; v > n && i >= 2; i--, n++)
2248 8543391 : v = minss(v, vali(gel(P,i)) + n);
2249 2233812 : gel(Q,l-1) = v? shifti(gel(P,l-1), -v): gel(P,l-1);
2250 12720891 : for (i = l-2, n = 1-v; i >= 2; i--, n++)
2251 10487161 : gel(Q,i) = shifti(gel(P,i), n);
2252 2233730 : Q[1] = P[1]; return Q;
2253 : }
2254 :
2255 : /* assume Q0 has no rational root */
2256 : static GEN
2257 1125974 : usp(GEN Q0, long flag, long bitprec)
2258 : {
2259 1125974 : const pari_sp av = avma;
2260 : GEN Qremapped, Q, c, Lc, Lk, sol;
2261 1125974 : GEN *pQremapped = flag == 1? &Qremapped: NULL;
2262 1125974 : const long prec = nbits2prec(bitprec), deg = degpol(Q0);
2263 1125974 : long listsize = 64, nbr = 0, nb_todo, ind, indf, i, k, nb;
2264 :
2265 1125974 : sol = zerocol(deg);
2266 1126011 : Lc = zerovec(listsize);
2267 1126032 : Lk = cgetg(listsize+1, t_VECSMALL);
2268 1126030 : k = Lk[1] = 0;
2269 1126030 : ind = 1; indf = 2;
2270 1126030 : Q = Q0;
2271 1126030 : c = gen_0;
2272 1126030 : nb_todo = 1;
2273 7116114 : while (nb_todo)
2274 : {
2275 5990116 : GEN nc = gel(Lc, ind);
2276 : pari_sp av2;
2277 5990116 : if (Lk[ind] == k + 1)
2278 : {
2279 2233751 : Q = Q0 = ZX_rescale2prim(Q0);
2280 2233737 : c = gen_0;
2281 : }
2282 5990102 : if (!equalii(nc, c)) Q = ZX_Z_translate(Q, subii(nc, c));
2283 5990156 : av2 = avma;
2284 5990156 : k = Lk[ind];
2285 5990156 : ind++;
2286 5990156 : c = nc;
2287 5990156 : nb_todo--;
2288 5990156 : nb = X2XP1(Q, pQremapped);
2289 :
2290 5989860 : if (nb == 1)
2291 : { /* exactly one root */
2292 1911155 : GEN s = gen_0;
2293 1911155 : if (flag == 0)
2294 : {
2295 0 : s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
2296 0 : s = gc_GEN(av2, s);
2297 : }
2298 1911155 : else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
2299 : {
2300 166631 : s = polsolve(*pQremapped, bitprec+1);
2301 166639 : s = addir(c, divrr(s, addsr(1, s)));
2302 166634 : shiftr_inplace(s, -k);
2303 166633 : if (realprec(s) != prec) s = rtor(s, prec);
2304 166641 : s = gc_upto(av2, s);
2305 : }
2306 1744524 : else set_avma(av2);
2307 1911181 : gel(sol, ++nbr) = s;
2308 : }
2309 4078705 : else if (nb)
2310 : { /* unknown, add two nodes to refine */
2311 2432146 : if (indf + 2 > listsize)
2312 : {
2313 1788 : if (ind>1)
2314 : {
2315 5297 : for (i = ind; i < indf; i++)
2316 : {
2317 3509 : gel(Lc, i-ind+1) = gel(Lc, i);
2318 3509 : Lk[i-ind+1] = Lk[i];
2319 : }
2320 1788 : indf -= ind-1;
2321 1788 : ind = 1;
2322 : }
2323 1788 : if (indf + 2 > listsize)
2324 : {
2325 0 : listsize *= 2;
2326 0 : Lc = vec_lengthen(Lc, listsize);
2327 0 : Lk = vecsmall_lengthen(Lk, listsize);
2328 : }
2329 112711 : for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2330 : }
2331 2432146 : gel(Lc, indf) = nc = shifti(c, 1);
2332 2432146 : gel(Lc, indf + 1) = addiu(nc, 1);
2333 2432154 : Lk[indf] = Lk[indf + 1] = k + 1;
2334 2432154 : indf += 2;
2335 2432154 : nb_todo += 2;
2336 : }
2337 5989894 : if (gc_needed(av, 2))
2338 : {
2339 0 : (void)gc_all(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
2340 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
2341 : }
2342 : }
2343 1125998 : setlg(sol, nbr+1);
2344 1126000 : return gc_GEN(av, sol);
2345 : }
2346 :
2347 : static GEN
2348 14 : ZX_Uspensky_equal_yes(GEN a, long flag, long bit)
2349 : {
2350 14 : if (flag == 2) return gen_1;
2351 7 : if (flag == 1 && typ(a) != t_REAL)
2352 : {
2353 7 : if (typ(a) == t_INT && !signe(a))
2354 0 : a = real_0_bit(bit);
2355 : else
2356 7 : a = gtofp(a, nbits2prec(bit));
2357 : }
2358 7 : return mkcol(a);
2359 : }
2360 : static GEN
2361 21 : ZX_Uspensky_no(long flag)
2362 21 : { return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
2363 : /* ZX_Uspensky(P, [a,a], flag) */
2364 : static GEN
2365 28 : ZX_Uspensky_equal(GEN P, GEN a, long flag, long bit)
2366 : {
2367 28 : if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
2368 14 : return ZX_Uspensky_equal_yes(a, flag, bit);
2369 : else
2370 14 : return ZX_Uspensky_no(flag);
2371 : }
2372 : static int
2373 3378 : sol_ok(GEN r, GEN a, GEN b) { return gcmp(a, r) <= 0 && gcmp(r, b) <= 0; }
2374 :
2375 : /* P a ZX without real double roots; better if primitive and squarefree but
2376 : * caller should ensure that. If flag & 4 assume that P has no rational root
2377 : * (modest speedup) */
2378 : GEN
2379 1313844 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
2380 : {
2381 1313844 : pari_sp av = avma;
2382 : GEN a, b, res, sol;
2383 : double fb;
2384 : long l, nbz, deg;
2385 :
2386 1313844 : if (ab)
2387 : {
2388 1206263 : if (typ(ab) == t_VEC)
2389 : {
2390 1178570 : if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
2391 1178570 : a = gel(ab, 1);
2392 1178570 : b = gel(ab, 2);
2393 : }
2394 : else
2395 : {
2396 27693 : a = ab;
2397 27693 : b = mkoo();
2398 : }
2399 : }
2400 : else
2401 : {
2402 107581 : a = mkmoo();
2403 107581 : b = mkoo();
2404 : }
2405 1313844 : if (flag & 4)
2406 : {
2407 129320 : if (gcmp(a, b) >= 0) { set_avma(av); return ZX_Uspensky_no(flag); }
2408 129320 : flag &= ~4;
2409 129320 : sol = cgetg(1, t_COL);
2410 : }
2411 : else
2412 : {
2413 1184524 : switch (gcmp(a, b))
2414 : {
2415 7 : case 1: set_avma(av); return ZX_Uspensky_no(flag);
2416 28 : case 0: return gc_GEN(av, ZX_Uspensky_equal(P, a, flag, bitprec));
2417 : }
2418 1184489 : sol = nfrootsQ(P);
2419 : }
2420 1313809 : nbz = 0; l = lg(sol);
2421 1313809 : if (l > 1)
2422 : {
2423 : long i, j;
2424 2706 : P = RgX_div(P, roots_to_pol(sol, varn(P)));
2425 2706 : if (!RgV_is_ZV(sol)) P = Q_primpart(P);
2426 6084 : for (i = j = 1; i < l; i++)
2427 3378 : if (sol_ok(gel(sol,i), a, b)) gel(sol,j++) = gel(sol,i);
2428 2706 : setlg(sol, j);
2429 2706 : if (flag == 2) { nbz = j-1; sol = utoi(nbz); }
2430 2559 : else if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
2431 : }
2432 1311103 : else if (flag == 2) sol = gen_0;
2433 1313809 : deg = degpol(P);
2434 1313810 : if (deg == 0) return gc_GEN(av, sol);
2435 1311863 : if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
2436 : {
2437 28 : fb = fujiwara_bound_real(P, -1);
2438 28 : if (fb <= -pariINFINITY) a = gen_0;
2439 21 : else if (fb < 0) a = gen_m1;
2440 21 : else a = negi(int2n((long)ceil(fb)));
2441 : }
2442 1311863 : if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
2443 : {
2444 21 : fb = fujiwara_bound_real(P, 1);
2445 21 : if (fb <= -pariINFINITY) b = gen_0;
2446 21 : else if (fb < 0) b = gen_1;
2447 7 : else b = int2n((long)ceil(fb));
2448 : }
2449 1311862 : if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
2450 : {
2451 : GEN d, ad, bd, diff;
2452 : long i;
2453 : /* can occur if one of a,b was initially a t_INFINITY */
2454 12582 : if (gequal(a,b)) return gc_GEN(av, sol);
2455 12575 : d = lcmii(Q_denom(a), Q_denom(b));
2456 12575 : if (is_pm1(d)) { d = NULL; ad = a; bd = b; }
2457 : else
2458 14 : { P = ZX_rescale(P, d); ad = gmul(a, d); bd = gmul(b, d); }
2459 12575 : diff = subii(bd, ad);
2460 12575 : P = ZX_affine(P, diff, ad);
2461 12575 : res = usp(P, flag, bitprec);
2462 12575 : if (flag <= 1)
2463 : {
2464 34176 : for (i = 1; i < lg(res); i++)
2465 : {
2466 21916 : GEN z = gmul(diff, gel(res, i));
2467 21916 : if (typ(z) == t_VEC)
2468 : {
2469 0 : gel(z, 1) = gadd(ad, gel(z, 1));
2470 0 : gel(z, 2) = gadd(ad, gel(z, 2));
2471 : }
2472 : else
2473 21916 : z = gadd(ad, z);
2474 21916 : if (d) z = gdiv(z, d);
2475 21916 : gel(res, i) = z;
2476 : }
2477 12260 : sol = shallowconcat(sol, res);
2478 : }
2479 : else
2480 315 : nbz += lg(res) - 1;
2481 : }
2482 1311855 : if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(P, 1)) > -pariINFINITY)
2483 : {
2484 1023074 : long bp = maxss((long)ceil(fb), 0);
2485 1023074 : res = usp(ZX_unscale2n(P, bp), flag, bitprec);
2486 1023083 : if (flag <= 1)
2487 71241 : sol = shallowconcat(sol, gmul2n(res, bp));
2488 : else
2489 951842 : nbz += lg(res)-1;
2490 : }
2491 1311863 : if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(P,-1)) > -pariINFINITY)
2492 : {
2493 90372 : long i, bm = maxss((long)ceil(fb), 0);
2494 90372 : res = usp(ZX_unscale2n(ZX_z_unscale(P, -1), bm), flag, bitprec);
2495 90373 : if (flag <= 1)
2496 : {
2497 75627 : for (i = 1; i < lg(res); i++)
2498 : {
2499 47406 : GEN z = gneg(gmul2n(gel(res, i), bm));
2500 47406 : if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
2501 47406 : gel(res, i) = z;
2502 : }
2503 28221 : sol = shallowconcat(res, sol);
2504 : }
2505 : else
2506 62152 : nbz += lg(res)-1;
2507 : }
2508 1311862 : if (flag >= 2) return utoi(nbz);
2509 83585 : if (flag)
2510 83585 : sol = sort(sol);
2511 : else
2512 0 : sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
2513 83584 : return gc_upto(av, sol);
2514 : }
2515 :
2516 : /* x a scalar */
2517 : static GEN
2518 42 : rootsdeg0(GEN x)
2519 : {
2520 42 : if (!is_real_t(typ(x))) pari_err_TYPE("realroots",x);
2521 35 : if (gequal0(x)) pari_err_ROOTS0("realroots");
2522 14 : return cgetg(1,t_COL); /* constant polynomial */
2523 : }
2524 : static void
2525 2354852 : checkbound(GEN a)
2526 : {
2527 2354852 : switch(typ(a))
2528 : {
2529 2354852 : case t_INT: case t_FRAC: case t_INFINITY: break;
2530 0 : default: pari_err_TYPE("polrealroots", a);
2531 : }
2532 2354852 : }
2533 : static GEN
2534 1178861 : check_ab(GEN ab)
2535 : {
2536 : GEN a, b;
2537 1178861 : if (!ab) return NULL;
2538 1177426 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
2539 1177426 : a = gel(ab,1); checkbound(a);
2540 1177426 : b = gel(ab,2); checkbound(b);
2541 1177427 : if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
2542 448 : typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
2543 1177427 : return ab;
2544 : }
2545 : /* e^(1/h) assuming the h-th root is real, beware that sqrtnr assumes e >= 0 */
2546 : static GEN
2547 22710 : _sqrtnr(GEN e, long h)
2548 : {
2549 : long s;
2550 : GEN r;
2551 22710 : if (h == 2) return sqrtr(e);
2552 14 : s = signe(e); setsigne(e, 1); /* e < 0 is possible, implies h is odd */
2553 14 : r = sqrtnr(e, h); if (s < 0) setsigne(r, -1);
2554 14 : return r;
2555 : }
2556 : GEN
2557 50627 : realroots(GEN P, GEN ab, long prec)
2558 : {
2559 50627 : pari_sp av = avma;
2560 50627 : GEN sol = NULL, fa, ex;
2561 : long i, j, v, l;
2562 :
2563 50627 : ab = check_ab(ab);
2564 50627 : if (typ(P) != t_POL) return rootsdeg0(P);
2565 50606 : if (!RgX_is_ZX(P)) P = RgX_rescale_to_int(P);
2566 50606 : switch(degpol(P))
2567 : {
2568 14 : case -1: return rootsdeg0(gen_0);
2569 7 : case 0: return rootsdeg0(gel(P,2));
2570 : }
2571 50585 : v = ZX_valrem(Q_primpart(P), &P);
2572 50585 : fa = ZX_squff(P, &ex); l = lg(fa); sol = cgetg(l + 1, t_VEC);
2573 102711 : for (i = 1; i < l; i++)
2574 : {
2575 52126 : GEN Pi = gel(fa, i), soli, soli2;
2576 : long n, h;
2577 52126 : if (ab) h = 1; else Pi = ZX_deflate_max(Pi, &h);
2578 52126 : soli = ZX_Uspensky(Pi, odd(h)? ab: gen_0, 1, prec2nbits(prec));
2579 52126 : n = lg(soli); soli2 = odd(h)? NULL: cgetg(n, t_COL);
2580 119122 : for (j = 1; j < n; j++)
2581 : {
2582 66996 : GEN r = gel(soli, j); /* != 0 */
2583 66996 : if (typ(r) != t_REAL) gel(soli, j) = r = gtofp(r, prec);
2584 66996 : if (h > 1)
2585 : {
2586 77 : gel(soli, j) = r = _sqrtnr(r, h);
2587 77 : if (soli2) gel(soli2, j) = negr(r);
2588 : }
2589 : }
2590 52126 : if (soli2) soli = shallowconcat(soli, soli2);
2591 52126 : if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
2592 52126 : gel(sol, i) = soli;
2593 : }
2594 50585 : if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
2595 84 : gel(sol, i++) = const_col(v, real_0(prec));
2596 50585 : setlg(sol, i); if (i == 1) retgc_const(av, cgetg(1, t_COL));
2597 50571 : return gc_upto(av, sort(shallowconcat1(sol)));
2598 : }
2599 : GEN
2600 48624 : ZX_realroots_irred(GEN P, long prec)
2601 : {
2602 48624 : long dP = degpol(P), j, n, h;
2603 : GEN sol, sol2;
2604 : pari_sp av;
2605 48624 : if (dP == 1) retmkvec(ZX_deg1root(P, prec));
2606 45316 : av = avma; P = ZX_deflate_max(P, &h);
2607 45316 : if (h == dP)
2608 : {
2609 11967 : GEN r = _sqrtnr(ZX_deg1root(P, prec), h);
2610 11967 : return gc_GEN(av, odd(h)? mkvec(r): mkvec2(negr(r), r));
2611 : }
2612 33349 : sol = ZX_Uspensky(P, odd(h)? NULL: gen_0, 1 | 4, prec2nbits(prec));
2613 33349 : n = lg(sol); sol2 = odd(h)? NULL: cgetg(n, t_COL);
2614 133754 : for (j = 1; j < n; j++)
2615 : {
2616 100403 : GEN r = gel(sol, j);
2617 100403 : if (typ(r) != t_REAL) gel(sol, j) = r = gtofp(r, prec);
2618 100403 : if (h > 1)
2619 : {
2620 10666 : gel(sol, j) = r = _sqrtnr(r, h);
2621 10666 : if (sol2) gel(sol2, j) = negr(r);
2622 : }
2623 : }
2624 33351 : if (sol2) sol = shallowconcat(sol, sol2);
2625 33351 : return gc_upto(av, sort(sol));
2626 : }
2627 :
2628 : static long
2629 123155 : ZX_sturm_i(GEN P, long flag)
2630 : {
2631 : pari_sp av;
2632 123155 : long h, r, dP = degpol(P);
2633 123154 : if (dP == 1) return 1;
2634 119812 : av = avma; P = ZX_deflate_max(P, &h);
2635 119813 : if (h == dP)
2636 : { /* now deg P = 1 */
2637 18274 : if (odd(h))
2638 665 : r = 1;
2639 : else
2640 17609 : r = (signe(gel(P,2)) != signe(gel(P,3)))? 2: 0;
2641 18274 : return gc_long(av, r);
2642 : }
2643 101539 : if (odd(h))
2644 78767 : r = itou(ZX_Uspensky(P, NULL, flag, 0));
2645 : else
2646 22773 : r = 2*itou(ZX_Uspensky(P, gen_0, flag, 0));
2647 101540 : return gc_long(av,r);
2648 : }
2649 : /* P nonconstant, squarefree ZX */
2650 : long
2651 1128234 : ZX_sturmpart(GEN P, GEN ab)
2652 : {
2653 1128234 : pari_sp av = avma;
2654 1128234 : if (!check_ab(ab)) return ZX_sturm(P);
2655 1126829 : return gc_long(av, itou(ZX_Uspensky(P, ab, 2, 0)));
2656 : }
2657 : /* P nonconstant, squarefree ZX */
2658 : long
2659 6395 : ZX_sturm(GEN P) { return ZX_sturm_i(P, 2); }
2660 : /* P irreducible ZX */
2661 : long
2662 116760 : ZX_sturm_irred(GEN P) { return ZX_sturm_i(P, 2 + 4); }
|
