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