Branch data 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 = 100000.;
24 : :
25 : : /********************************************************************/
26 : : /** **/
27 : : /** FAST ARITHMETIC over Z[i] **/
28 : : /** **/
29 : : /********************************************************************/
30 : : static THREAD long KARASQUARE_LIMIT, COOKSQUARE_LIMIT;
31 : :
32 : : /* fast sum of x,y: t_INT or t_COMPLEX(t_INT) */
33 : : static GEN
34 : 9390951 : addCC(GEN x, GEN y)
35 : : {
36 : : GEN z;
37 : :
38 [ + + ]: 9390951 : if (typ(x) == t_INT)
39 : : {
40 [ + + ]: 8314356 : if (typ(y) == t_INT) return addii(x,y);
41 : : /* ty == t_COMPLEX */
42 : 1036 : z = cgetg(3,t_COMPLEX);
43 : 1036 : gel(z,1) = addii(x, gel(y,1));
44 : 1036 : gel(z,2) = icopy(gel(y,2)); return z;
45 : : }
46 : : /* tx == t_COMPLEX */
47 : 1076595 : z = cgetg(3,t_COMPLEX);
48 [ + + ]: 1076595 : if (typ(y) == t_INT)
49 : : {
50 : 14667 : gel(z,1) = addii(gel(x,1),y);
51 : 14667 : gel(z,2) = icopy(gel(x,2)); return z;
52 : : }
53 : : /* ty == t_COMPLEX */
54 : 1061928 : gel(z,1) = addii(gel(x,1),gel(y,1));
55 : 9390951 : gel(z,2) = addii(gel(x,2),gel(y,2)); return z;
56 : : }
57 : : /* fast product of x,y: t_INT or t_COMPLEX(t_INT) */
58 : : static GEN
59 : 17763117 : mulCC(GEN x, GEN y)
60 : : {
61 : : GEN z;
62 : :
63 [ + + ]: 17763117 : if (typ(x) == t_INT)
64 : : {
65 [ + + ]: 15172030 : if (typ(y) == t_INT) return mulii(x,y);
66 : : /* ty == t_COMPLEX */
67 : 3744 : z = cgetg(3,t_COMPLEX);
68 : 3744 : gel(z,1) = mulii(x, gel(y,1));
69 : 3744 : gel(z,2) = mulii(x, gel(y,2)); return z;
70 : : }
71 : : /* tx == t_COMPLEX */
72 : 2591087 : z = cgetg(3,t_COMPLEX);
73 [ + + ]: 2591087 : if (typ(y) == t_INT)
74 : : {
75 : 747248 : gel(z,1) = mulii(gel(x,1),y);
76 : 747248 : gel(z,2) = mulii(gel(x,2),y); return z;
77 : : }
78 : : /* ty == t_COMPLEX */
79 : : {
80 : 1843839 : pari_sp av = avma, tetpil;
81 : : GEN p1, p2;
82 : :
83 : 1843839 : p1 = mulii(gel(x,1),gel(y,1));
84 : 1843839 : p2 = mulii(gel(x,2),gel(y,2));
85 : 1843839 : y = mulii(addii(gel(x,1),gel(x,2)),
86 : 3687678 : addii(gel(y,1),gel(y,2)));
87 : 1843839 : x = addii(p1,p2); tetpil = avma;
88 : 1843839 : gel(z,1) = subii(p1,p2);
89 : 1843839 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
90 : 17763117 : return z;
91 : : }
92 : : }
93 : : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
94 : : static GEN
95 : 14605028 : sqrCC(GEN x)
96 : : {
97 : : GEN z;
98 : :
99 [ + + ]: 14605028 : if (typ(x) == t_INT) return sqri(x);
100 : : /* tx == t_COMPLEX */
101 : 2506560 : z = cgetg(3,t_COMPLEX);
102 : : {
103 : 2506560 : pari_sp av = avma, tetpil;
104 : : GEN y, p1, p2;
105 : :
106 : 2506560 : p1 = sqri(gel(x,1));
107 : 2506560 : p2 = sqri(gel(x,2));
108 : 2506560 : y = sqri(addii(gel(x,1),gel(x,2)));
109 : 2506560 : x = addii(p1,p2); tetpil = avma;
110 : 2506560 : gel(z,1) = subii(p1,p2);
111 : 2506560 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
112 : 14605028 : return z;
113 : : }
114 : : }
115 : :
116 : : static void
117 : 2428606 : set_karasquare_limit(long bit)
118 : : {
119 [ + + ]: 2428606 : if (bit<600) { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
120 [ + - ]: 630 : else if (bit<2000) { KARASQUARE_LIMIT=4; COOKSQUARE_LIMIT=200; }
121 [ # # ]: 0 : else if (bit<3000) { KARASQUARE_LIMIT=4; COOKSQUARE_LIMIT=125; }
122 [ # # ]: 0 : else if (bit<5000) { KARASQUARE_LIMIT=2; COOKSQUARE_LIMIT= 75; }
123 : 0 : else { KARASQUARE_LIMIT=1; COOKSQUARE_LIMIT= 50; }
124 : 2428606 : }
125 : :
126 : : /* assume lP > 0, lP = lgpol(P) */
127 : : static GEN
128 : 5057169 : CX_square_spec(GEN P, long lP)
129 : : {
130 : : GEN s, t;
131 : 5057169 : long i, j, l, nn, n = lP - 1;
132 : : pari_sp av;
133 : :
134 : 5057169 : nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
135 : 5057169 : gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
136 [ + + ]: 13429335 : for (i=1; i<=n; i++)
137 : : {
138 : 8372166 : av = avma; l = (i+1)>>1;
139 : 8372166 : t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
140 [ + + ]: 11599755 : for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
141 : 8372166 : t = gmul2n(t,1);
142 [ + + ]: 8372166 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
143 : 8372166 : gel(s,i+2) = gerepileupto(av, t);
144 : : }
145 : 5057169 : gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
146 [ + + ]: 9547859 : for ( ; i<nn; i++)
147 : : {
148 : 4490690 : av = avma; l = (i+1)>>1;
149 : 4490690 : t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
150 [ + + ]: 6163362 : for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
151 : 4490690 : t = gmul2n(t,1);
152 [ + + ]: 4490690 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
153 : 4490690 : gel(s,i+2) = gerepileupto(av, t);
154 : : }
155 : 5057169 : return normalizepol_lg(s, nn+3);
156 : : }
157 : : /* not stack clean */
158 : : static GEN
159 : 100437 : RgX_addspec(GEN x, long nx, GEN y, long ny)
160 : : {
161 : : GEN z, t;
162 : : long i;
163 [ + + ]: 100437 : if (nx == ny) {
164 : 53781 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
165 [ + + ]: 404779 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
166 : 53781 : return normalizepol_lg(z, nx+2);
167 : : }
168 [ + - ]: 46656 : if (ny < nx) {
169 : 46656 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
170 [ + + ]: 321166 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
171 [ + + ]: 93312 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
172 : 46656 : return normalizepol_lg(z, nx+2);
173 : : } else {
174 : 0 : z = cgetg(ny+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
175 [ # # ]: 0 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
176 [ # # ]: 0 : for ( ; i < ny; i++) gel(t,i) = gel(y,i);
177 : 100437 : return normalizepol_lg(z, ny+2);
178 : : }
179 : : }
180 : : /* nx = lgpol(x) */
181 : : static GEN
182 : 0 : RgX_s_mulspec(GEN x, long nx, long s)
183 : : {
184 : : GEN z, t;
185 : : long i;
186 [ # # ][ # # ]: 0 : if (!s || !nx) return pol_0(0);
187 : 0 : z = cgetg(nx+2, t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z + 2;
188 [ # # ]: 0 : for (i=0; i < nx; i++) gel(t,i) = gmulgs(gel(x,i), s);
189 : 0 : return z;
190 : : }
191 : : /* nx = lgpol(x), return x << s. Inefficient if s = 0... */
192 : : static GEN
193 : 0 : RgX_shiftspec(GEN x, long nx, long s)
194 : : {
195 : : GEN z, t;
196 : : long i;
197 [ # # ]: 0 : if (!nx) return pol_0(0);
198 : 0 : z = cgetg(nx+2, t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z + 2;
199 [ # # ]: 0 : for (i=0; i < nx; i++) gel(t,i) = gmul2n(gel(x,i), s);
200 : 0 : return z;
201 : : }
202 : :
203 : : /* spec function. nP = lgpol(P) */
204 : : static GEN
205 : 5158523 : karasquare(GEN P, long nP)
206 : : {
207 : : GEN Q, s0, s1, s2, a, t;
208 : 5158523 : long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
209 : : pari_sp av;
210 : :
211 [ + + ][ + + ]: 5158523 : if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
212 : 100437 : av = avma;
213 : 100437 : n0 = (n>>1) + 1; n1 = nP - n0;
214 : 100437 : s0 = karasquare(P, n0); Q = P + n0;
215 : 100437 : s2 = karasquare(Q, n1);
216 : 100437 : s1 = RgX_addspec(P, n0, Q, n1);
217 : 100437 : s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
218 : 100437 : N = (n<<1) + 1;
219 : 100437 : a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
220 : 100437 : t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
221 [ + + ]: 1303350 : for (i=0; i < l; i++) gel(t,i) = gel(s0,i);
222 [ + + ]: 241852 : for ( ; i < N0; i++) gel(t,i) = gen_0;
223 : 100437 : t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
224 [ + + ]: 1170910 : for (i=0; i < l; i++) gel(t,i) = gel(s2,i);
225 [ + + ]: 180543 : for ( ; i < N1; i++) gel(t,i) = gen_0;
226 : 100437 : t = a+2 + n0; l = lgpol(s1); s1 += 2;
227 [ + + ]: 1229731 : for (i=0; i < l; i++) gel(t,i) = gadd(gel(t,i), gel(s1,i));
228 : 5158523 : return gerepilecopy(av, normalizepol_lg(a, N+2));
229 : : }
230 : : /* spec function. nP = lgpol(P) */
231 : : static GEN
232 : 4857212 : cook_square(GEN P, long nP)
233 : : {
234 : : GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
235 : 4857212 : long n0, n3, i, j, n = nP - 1;
236 : : pari_sp av;
237 : :
238 [ + - ][ + - ]: 4857212 : if (n <= COOKSQUARE_LIMIT) return nP? karasquare(P, nP): pol_0(0);
239 : 0 : av = avma;
240 : :
241 : 0 : n0 = (n+1) >> 2; n3 = n+1 - 3*n0;
242 : 0 : p0 = P;
243 : 0 : p1 = p0+n0;
244 : 0 : p2 = p1+n0;
245 : 0 : p3 = p2+n0; /* lgpol(p0,p1,p2) = n0, lgpol(p3) = n3 */
246 : :
247 : 0 : q = cgetg(8,t_VEC) + 4;
248 : 0 : Q = cook_square(p0, n0);
249 : 0 : r = RgX_addspec(p0,n0, p2,n0);
250 : 0 : t = RgX_addspec(p1,n0, p3,n3);
251 : 0 : gel(q,-1) = RgX_sub(r,t);
252 : 0 : gel(q,1) = RgX_add(r,t);
253 : 0 : r = RgX_addspec(p0,n0, RgX_shiftspec(p2,n0, 2)+2,n0);
254 : 0 : t = gmul2n(RgX_addspec(p1,n0, RgX_shiftspec(p3,n3, 2)+2,n3), 1);
255 : 0 : gel(q,-2) = RgX_sub(r,t);
256 : 0 : gel(q,2) = RgX_add(r,t);
257 : 0 : r = RgX_addspec(p0,n0, RgX_s_mulspec(p2,n0, 9)+2,n0);
258 : 0 : t = gmulsg(3, RgX_addspec(p1,n0, RgX_s_mulspec(p3,n3, 9)+2,n3));
259 : 0 : gel(q,-3) = RgX_sub(r,t);
260 : 0 : gel(q,3) = RgX_add(r,t);
261 : :
262 : 0 : r = new_chunk(7);
263 : 0 : vp = cgetg(4,t_VEC);
264 : 0 : vm = cgetg(4,t_VEC);
265 [ # # ]: 0 : for (i=1; i<=3; i++)
266 : : {
267 : 0 : GEN a = gel(q,i), b = gel(q,-i);
268 : 0 : a = cook_square(a+2, lgpol(a));
269 : 0 : b = cook_square(b+2, lgpol(b));
270 : 0 : gel(vp,i) = RgX_add(b, a);
271 : 0 : gel(vm,i) = RgX_sub(b, a);
272 : : }
273 : 0 : gel(r,0) = Q;
274 : 0 : gel(r,1) = gdivgs(gsub(gsub(gmulgs(gel(vm,2),9),gel(vm,3)),
275 : 0 : gmulgs(gel(vm,1),45)),
276 : : 60);
277 : 0 : gel(r,2) = gdivgs(gadd(gadd(gmulgs(gel(vp,1),270),gmulgs(Q,-490)),
278 : 0 : gadd(gmulgs(gel(vp,2),-27),gmulgs(gel(vp,3),2))),
279 : : 360);
280 : 0 : gel(r,3) = gdivgs(gadd(gadd(gmulgs(gel(vm,1),13),gmulgs(gel(vm,2),-8)),
281 : 0 : gel(vm,3)),
282 : : 48);
283 : 0 : gel(r,4) = gdivgs(gadd(gadd(gmulgs(Q,56),gmulgs(gel(vp,1),-39)),
284 : 0 : gsub(gmulgs(gel(vp,2),12),gel(vp,3))),
285 : : 144);
286 : 0 : gel(r,5) = gdivgs(gsub(gadd(gmulgs(gel(vm,1),-5),gmulgs(gel(vm,2),4)),
287 : 0 : gel(vm,3)),
288 : : 240);
289 : 0 : gel(r,6) = gdivgs(gadd(gadd(gmulgs(Q,-20),gmulgs(gel(vp,1),15)),
290 : 0 : gadd(gmulgs(gel(vp,2),-6),gel(vp,3))),
291 : : 720);
292 : 0 : q = cgetg(2*n+3,t_POL); q[1] = evalsigne(1)|evalvarn(0);
293 : 0 : t = q+2;
294 [ # # ]: 0 : for (i=0; i<=2*n; i++) gel(t,i) = gen_0;
295 [ # # ]: 0 : for (i=0; i<=6; i++,t += n0)
296 : : {
297 : 0 : GEN h = gel(r,i);
298 : 0 : long d = lgpol(h);
299 : 0 : h += 2;
300 [ # # ]: 0 : for (j=0; j<d; j++) gel(t,j) = gadd(gel(t,j), gel(h,j));
301 : : }
302 : 4857212 : return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
303 : : }
304 : :
305 : : static GEN
306 : 2428606 : graeffe(GEN p)
307 : : {
308 : : GEN p0, p1, s0, s1;
309 : 2428606 : long n = degpol(p), n0, n1, i;
310 : :
311 [ - + ]: 2428606 : if (!n) return gcopy(p);
312 : 2428606 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
313 : 2428606 : p0 = new_chunk(n0);
314 : 2428606 : p1 = new_chunk(n1);
315 [ + + ]: 8192602 : for (i=0; i<n1; i++)
316 : : {
317 : 5763996 : p0[i] = p[2+(i<<1)];
318 : 5763996 : p1[i] = p[3+(i<<1)];
319 : : }
320 [ + + ]: 2428606 : if (n1 != n0)
321 : 1241331 : p0[i] = p[2+(i<<1)];
322 : 2428606 : s0 = cook_square(p0, n0);
323 : 2428606 : s1 = cook_square(p1, n1);
324 : 2428606 : return RgX_sub(s0, RgX_shift_shallow(s1,1));
325 : : }
326 : : GEN
327 : 4655 : ZX_graeffe(GEN p)
328 : : {
329 : 4655 : pari_sp av = avma;
330 : : GEN p0, p1, s0, s1;
331 : 4655 : long n = degpol(p);
332 : :
333 [ - + ]: 4655 : if (!n) return ZX_copy(p);
334 : 4655 : RgX_even_odd(p, &p0, &p1);
335 : : /* p = p0(x^2) + x p1(x^2) */
336 : 4655 : s0 = ZX_sqr(p0);
337 : 4655 : s1 = ZX_sqr(p1);
338 : 4655 : return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
339 : : }
340 : : GEN
341 : 14 : polgraeffe(GEN p)
342 : : {
343 : 14 : pari_sp av = avma;
344 : : GEN p0, p1, s0, s1;
345 : 14 : long n = degpol(p);
346 : :
347 [ - + ]: 14 : if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
348 : 14 : n = degpol(p);
349 [ - + ]: 14 : if (!n) return gcopy(p);
350 : 14 : RgX_even_odd(p, &p0, &p1);
351 : : /* p = p0(x^2) + x p1(x^2) */
352 : 14 : s0 = RgX_sqr(p0);
353 : 14 : s1 = RgX_sqr(p1);
354 : 14 : return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
355 : : }
356 : :
357 : : /********************************************************************/
358 : : /** **/
359 : : /** MODULUS OF ROOTS **/
360 : : /** **/
361 : : /********************************************************************/
362 : :
363 : : /* Quick approximation to log2(|x); first define y s.t. |y-x| < 2^-32 then
364 : : * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
365 : : * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
366 : : static double
367 : 13698183 : mydbllog2i(GEN x)
368 : : {
369 : : #ifdef LONG_IS_64BIT
370 : 11745592 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
371 : : #else
372 : 1952591 : const double W = 1/4294967296.; /*2^-32*/
373 : : #endif
374 : : GEN m;
375 : 13698183 : long lx = lgefint(x);
376 : : double l;
377 [ + + ]: 13698183 : if (lx == 2) return -pariINFINITY;
378 : 13600591 : m = int_MSW(x);
379 : 13600591 : l = (double)(ulong)*m;
380 [ + + ]: 13600591 : if (lx == 3) return log2(l);
381 : 6331544 : l += ((double)(ulong)*int_precW(m)) * W;
382 : : /* at least m = min(53,BIL) bits are correct in the mantissa, thus log2
383 : : * is correct with error < log(1 + 2^-m) ~ 2^-m. Adding the correct
384 : : * exponent BIL(lx-3) causes 1ulp further round-off error */
385 : 13698183 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
386 : : }
387 : :
388 : : /* return log(|x|) or -pariINFINITY */
389 : : static double
390 : 755872 : mydbllogr(GEN x) {
391 [ - + ]: 755872 : if (!signe(x)) return -pariINFINITY;
392 : 755872 : return LOG2*dbllog2r(x);
393 : : }
394 : :
395 : : /* return log2(|x|) or -pariINFINITY */
396 : : static double
397 : 7072970 : mydbllog2r(GEN x) {
398 [ + + ]: 7072970 : if (!signe(x)) return -pariINFINITY;
399 : 7072970 : return dbllog2r(x);
400 : : }
401 : : static double
402 [ + + ]: 3939708 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
403 : : double
404 : 18427826 : dbllog2(GEN z)
405 : : {
406 : : double x, y;
407 [ + + + ]: 18427826 : switch(typ(z))
408 : : {
409 : 10384893 : case t_INT: return mydbllog2i(z);
410 : 6073079 : case t_REAL: return mydbllog2r(z);
411 : : default: /*t_COMPLEX*/
412 : 1969854 : x = dbllog2mp(gel(z,1));
413 : 1969854 : y = dbllog2mp(gel(z,2));
414 [ + + ]: 1969854 : if (fabs(x-y) > 10) return maxdd(x,y);
415 : 18427826 : return x + 0.5*log2(1 + exp2(2*(y-x)));
416 : : }
417 : : }
418 : : static GEN /* beware overflow */
419 [ + + ]: 370748 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
420 : :
421 : : /* find s such that A_h <= 2^s <= 2 A_i for one h and all i < n = deg(p),
422 : : * with A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and p = sum p_i X^i */
423 : : static long
424 : 1873344 : findpower(GEN p)
425 : : {
426 : 1873344 : double x, L, mins = pariINFINITY;
427 : 1873344 : long n = degpol(p),i;
428 : :
429 : 1873344 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
430 [ + + ]: 9603974 : for (i=n-1; i>=0; i--)
431 : : {
432 : 7730630 : L += log2((double)(i+1) / (double)(n-i));
433 : 7730630 : x = dbllog2(gel(p,i+2));
434 [ + + ]: 7730630 : if (x != -pariINFINITY)
435 : : {
436 : 7671705 : double s = (L - x) / (double)(n-i);
437 [ + + ]: 7671705 : if (s < mins) mins = s;
438 : : }
439 : : }
440 : 1873344 : i = (long)ceil(mins);
441 [ + + ]: 1873344 : if (i - mins > 1 - 1e-12) i--;
442 : 1873344 : return i;
443 : : }
444 : :
445 : : /* returns the exponent for logmodulus(), from the Newton diagram */
446 : : static long
447 : 323074 : newton_polygon(GEN p, long k)
448 : : {
449 : 323074 : pari_sp av = avma;
450 : : double *logcoef, slope;
451 : 323074 : long n = degpol(p), i, j, h, l, *vertex;
452 : :
453 : 323074 : init_dalloc();
454 : 323074 : logcoef = (double*)stack_malloc((n+1)*sizeof(double));
455 : 323074 : vertex = (long*)new_chunk(n+1);
456 : :
457 : : /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
458 [ + + ]: 1782911 : for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
459 : 323074 : vertex[0] = 1; /* sentinel */
460 [ + + ]: 1369188 : for (i=0; i < n; i=h)
461 : : {
462 : 1046114 : slope = logcoef[i+1]-logcoef[i];
463 [ + + ]: 4989069 : for (j = h = i+1; j<=n; j++)
464 : : {
465 : 3942955 : double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
466 [ + + ]: 3942955 : if (slope < pij) { slope = pij; h = j; }
467 : : }
468 : 1046114 : vertex[h] = 1;
469 : : }
470 [ + + ]: 349496 : h = k; while (!vertex[h]) h++;
471 [ + + ]: 340836 : l = k-1; while (!vertex[l]) l--;
472 : 323074 : avma = av;
473 : 323074 : return (long)floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
474 : : }
475 : :
476 : : /* change z into z*2^e, where z is real or complex of real */
477 : : static void
478 : 2234047 : myshiftrc(GEN z, long e)
479 : : {
480 [ + + ]: 2234047 : if (typ(z)==t_COMPLEX)
481 : : {
482 [ + + ]: 464942 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
483 [ + + ]: 464942 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
484 : : }
485 : : else
486 [ + + ]: 1769105 : if (signe(z)) shiftr_inplace(z, e);
487 : 2234047 : }
488 : :
489 : : /* return z*2^e, where z is integer or complex of integer (destroy z) */
490 : : static GEN
491 : 8275741 : myshiftic(GEN z, long e)
492 : : {
493 [ + + ]: 8275741 : if (typ(z)==t_COMPLEX)
494 : : {
495 [ + + ]: 1299839 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
496 : 1299839 : gel(z,2) = mpshift(gel(z,2),e);
497 : 1299839 : return z;
498 : : }
499 [ + + ]: 8275741 : return signe(z)? mpshift(z,e): gen_0;
500 : : }
501 : :
502 : : /* as real_1 with precision in bits, not in words */
503 : : static GEN
504 : 381582 : myreal_1(long bit)
505 : : {
506 [ - + ]: 381582 : if (bit < 0) bit = 0;
507 : 381582 : return real_1(nbits2prec(bit));
508 : : }
509 : : static GEN
510 : 408081 : RgX_gtofp_bit(GEN q, long bit)
511 : : {
512 [ - + ]: 408081 : if (bit < 0) bit = 0;
513 : 408081 : return RgX_gtofp(q, nbits2prec(bit));
514 : : }
515 : :
516 : : static GEN
517 : 15211630 : mygprecrc(GEN x, long prec, long e)
518 : : {
519 : : GEN y;
520 [ + + + ]: 15211630 : switch(typ(x))
521 : : {
522 [ + + ]: 11732247 : case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
523 : : case t_COMPLEX:
524 : 2672597 : y = cgetg(3,t_COMPLEX);
525 : 2672597 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
526 : 2672597 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
527 : 2672597 : return y;
528 : 15211630 : default: return gcopy(x);
529 : : }
530 : : }
531 : :
532 : : /* gprec behaves badly with the zero for polynomials.
533 : : The second parameter in mygprec is the precision in base 2 */
534 : : static GEN
535 : 3959515 : mygprec(GEN x, long bit)
536 : : {
537 : : long lx, i, e, prec;
538 : : GEN y;
539 : :
540 [ - + ]: 3959515 : if (bit < 0) bit = 0; /* should rarely happen */
541 : 3959515 : e = gexpo(x) - bit;
542 : 3959515 : prec = nbits2prec(bit);
543 [ + + ]: 3959515 : switch(typ(x))
544 : : {
545 : : case t_POL:
546 : 2352115 : y = cgetg_copy(x, &lx); y[1] = x[1];
547 [ + + ]: 10536692 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
548 : 2352115 : break;
549 : :
550 : 1607400 : default: y = mygprecrc(x,prec,e);
551 : : }
552 : 3959515 : return y;
553 : : }
554 : :
555 : : /* normalize a polynomial p, that is change it with coefficients in Z[i],
556 : : after making product by 2^shift */
557 : : static GEN
558 : 965346 : pol_to_gaussint(GEN p, long shift)
559 : : {
560 : 965346 : long i, l = lg(p);
561 : 965346 : GEN q = cgetg(l, t_POL); q[1] = p[1];
562 [ + + ]: 7006862 : for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
563 : 965346 : return q;
564 : : }
565 : :
566 : : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
567 : : static GEN
568 : 765068 : eval_rel_pol(GEN p, long bit)
569 : : {
570 : : long i;
571 [ + + ]: 5469093 : for (i = 2; i < lg(p); i++)
572 [ + + ]: 4704025 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
573 : 765068 : return pol_to_gaussint(p, bit-gexpo(p)+1);
574 : : }
575 : :
576 : : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
577 : : static GEN
578 : 85145 : homothetie(GEN p, double lrho, long bit)
579 : : {
580 : : GEN q, r, t, iR;
581 : 85145 : long n = degpol(p), i;
582 : :
583 : 85145 : iR = mygprec(dblexp(-lrho),bit);
584 : 85145 : q = mygprec(p, bit);
585 : 85145 : r = cgetg(n+3,t_POL); r[1] = p[1];
586 : 85145 : t = iR; r[n+2] = q[n+2];
587 [ + + ]: 545371 : for (i=n-1; i>0; i--)
588 : : {
589 : 460226 : gel(r,i+2) = gmul(t, gel(q,i+2));
590 : 460226 : t = mulrr(t, iR);
591 : : }
592 : 85145 : gel(r,2) = gmul(t, gel(q,2)); return r;
593 : : }
594 : :
595 : : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) [ ~as above with R = 2^-e ]*/
596 : : static void
597 : 523352 : homothetie2n(GEN p, long e)
598 : : {
599 [ + + ]: 523352 : if (e)
600 : : {
601 : 396509 : long i,n = lg(p)-1;
602 [ + + ]: 2630556 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
603 : : }
604 : 523352 : }
605 : :
606 : : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
607 : : static void
608 : 1673066 : homothetie_gauss(GEN p, long e, long f)
609 : : {
610 [ + + ][ + + ]: 1673066 : if (e || f)
611 : : {
612 : 1532580 : long i, n = lg(p)-1;
613 [ + + ]: 9808321 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
614 : : }
615 : 1673066 : }
616 : :
617 : : /* Lower bound on the modulus of the largest root z_0
618 : : * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
619 : : static double
620 : 1873344 : lower_bound(GEN p, long *k, double eps)
621 : : {
622 : 1873344 : long n = degpol(p), i, j;
623 : 1873344 : pari_sp ltop = avma;
624 : : GEN a, s, S, ilc;
625 : : double r, R, rho;
626 : :
627 [ + + ]: 1873344 : if (n < 4) { *k = n; return 0.; }
628 : 877600 : S = cgetg(5,t_VEC);
629 : 877600 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
630 [ + + ]: 4388000 : for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
631 : : /* i = 1 split out from next loop for efficiency and initialization */
632 : 877600 : s = gel(a,1);
633 : 877600 : gel(S,1) = gneg(s); /* Newton sum S_i */
634 : 877600 : rho = r = gtodouble(gabs(s,3));
635 : 877600 : R = r / n;
636 [ + + ]: 3510400 : for (i=2; i<=4; i++)
637 : : {
638 : 2632800 : s = gmulsg(i,gel(a,i));
639 [ + + ]: 7898400 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
640 : 2632800 : gel(S,i) = gneg(s); /* Newton sum S_i */
641 : 2632800 : r = gtodouble(gabs(s,3));
642 [ + + ]: 2632800 : if (r > 0.)
643 : : {
644 : 2626849 : r = exp(log(r/n) / (double)i);
645 [ + + ]: 2626849 : if (r > R) R = r;
646 : : }
647 : : }
648 [ + + ][ + + ]: 877600 : if (R > 0. && eps < 1.2)
649 : 876830 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
650 : : else
651 : 770 : *k = n;
652 : 1873344 : avma = ltop; return R;
653 : : }
654 : :
655 : : /* log of modulus of the largest root of p with relative error tau. Assume
656 : : * P(0) != 0 and P non constant */
657 : : static double
658 : 200278 : logmax_modulus(GEN p, double tau)
659 : : {
660 : : GEN r,q,aux,gunr;
661 : 200278 : pari_sp av, ltop = avma;
662 : 200278 : long i,k,n=degpol(p),nn,bit,M,e;
663 [ - + ]: 200278 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
664 : :
665 : 200278 : r = cgeti(BIGDEFAULTPREC);
666 : 200278 : av = avma;
667 : :
668 : 200278 : eps = - 1/log(1.5*tau2); /* > 0 */
669 : 200278 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
670 : 200278 : gunr = myreal_1(bit+2*n);
671 : 200278 : aux = gdiv(gunr, gel(p,2+n));
672 : 200278 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
673 : 200278 : e = findpower(q);
674 : 200278 : homothetie2n(q,e);
675 : 200278 : affsi(e, r);
676 : 200278 : q = pol_to_gaussint(q, bit);
677 : 200278 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
678 : 200278 : nn = n;
679 : 200278 : for (i=0,e=0;;)
680 : : { /* nn = deg(q) */
681 : 1873344 : rho = lower_bound(q, &k, eps);
682 [ + + ]: 1873344 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
683 : 1873344 : affii(shifti(addis(r,e), 1), r);
684 [ + + ]: 1873344 : if (++i == M) break;
685 : :
686 : 5019198 : bit = (long) ((double)k * log2(1./tau2) +
687 : 3346132 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
688 : 1673066 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
689 : 1673066 : nn -= RgX_valrem(q, &q);
690 : 1673066 : set_karasquare_limit(gexpo(q));
691 : 1673066 : q = gerepileupto(av, graeffe(q));
692 [ + + ]: 1673066 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
693 : 1673066 : eps = -1/log(tau2); /* > 0 */
694 : 1673066 : e = findpower(q);
695 : 1673066 : }
696 [ + + ]: 200278 : if (!signe(r)) { avma = ltop; return 0.; }
697 : 187584 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
698 : 200278 : avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
699 : : }
700 : : /* assume P non constant */
701 : : GEN
702 : 1596 : logmax_modulus_bound(GEN P)
703 : : {
704 : 1596 : (void)RgX_valrem_inexact(P,&P);
705 [ + + ]: 1596 : return lg(P)==3? gen_0: dblexp(logmax_modulus(P, 0.01) + 0.01);
706 : : }
707 : :
708 : : /* log of modulus of the smallest root of p, with relative error tau */
709 : : static double
710 : 70621 : logmin_modulus(GEN p, double tau)
711 : : {
712 : 70621 : pari_sp av = avma;
713 : : double r;
714 : :
715 [ - + ]: 70621 : if (gequal0(gel(p,2))) return -pariINFINITY;
716 : 70621 : r = - logmax_modulus(RgX_recip_shallow(p),tau);
717 : 70621 : avma = av; return r;
718 : : }
719 : :
720 : : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
721 : : static double
722 : 35454 : logmodulus(GEN p, long k, double tau)
723 : : {
724 : : GEN q;
725 : 35454 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
726 : 35454 : pari_sp av, ltop=avma;
727 : 35454 : double r, tau2 = tau/6;
728 : :
729 : 35454 : bit = (long)(n * (2. + log2(3.*n/tau2)));
730 : 35454 : av = avma;
731 : 35454 : q = gprec_w(p, nbits2prec(bit));
732 : 35454 : q = RgX_gtofp_bit(q, bit);
733 : 35454 : e = newton_polygon(q,k);
734 : 35454 : r = (double)e;
735 : 35454 : homothetie2n(q,e);
736 : 35454 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
737 [ + + ]: 323074 : for (i=1; i<imax; i++)
738 : : {
739 : 287620 : q = eval_rel_pol(q,bit);
740 : 287620 : kk -= RgX_valrem(q, &q);
741 : 287620 : nn = degpol(q);
742 : :
743 : 287620 : set_karasquare_limit(bit);
744 : 287620 : q = gerepileupto(av, graeffe(q));
745 : 287620 : e = newton_polygon(q,kk);
746 : 287620 : r += e / exp2((double)i);
747 : 287620 : q = RgX_gtofp_bit(q, bit);
748 : 287620 : homothetie2n(q,e);
749 : :
750 [ - + ]: 287620 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
751 : 287620 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
752 : : }
753 : 35454 : avma = ltop; return -r * LOG2;
754 : : }
755 : :
756 : : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
757 : : * rmin < r_k < rmax. This information helps because we may reduce precision
758 : : * quicker */
759 : : static double
760 : 35454 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
761 : : {
762 : : GEN q;
763 : 35454 : long n = degpol(p), i, imax, imax2, bit;
764 : 35454 : pari_sp ltop = avma, av;
765 : 35454 : double lrho, aux, tau2 = tau/6.;
766 : :
767 : 35454 : aux = (lrmax - lrmin) / 2. + 4*tau2;
768 : 35454 : imax = (long) log2(log((double)n)/ aux);
769 [ + + ]: 35454 : if (imax <= 0) return logmodulus(p,k,tau);
770 : :
771 : 34587 : lrho = (lrmin + lrmax) / 2;
772 : 34587 : av = avma;
773 : 34587 : bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
774 : 34587 : q = homothetie(p, lrho, bit);
775 : 34587 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
776 [ - + ]: 34587 : if (imax > imax2) imax = imax2;
777 : :
778 [ + + ]: 108012 : for (i=0; i<imax; i++)
779 : : {
780 : 73425 : q = eval_rel_pol(q,bit);
781 : 73425 : set_karasquare_limit(bit);
782 : 73425 : q = gerepileupto(av, graeffe(q));
783 : 73425 : aux = 2*aux + 2*tau2;
784 : 73425 : tau2 *= 1.5;
785 : 73425 : bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
786 : 73425 : q = RgX_gtofp_bit(q, bit);
787 : : }
788 : 34587 : aux = exp2((double)imax);
789 : 34587 : aux = logmodulus(q,k, aux*tau/3.) / aux;
790 : 35454 : avma = ltop; return lrho + aux;
791 : : }
792 : :
793 : : static double
794 : 41030 : ind_maxlog2(GEN q)
795 : : {
796 : 41030 : long i, k = -1;
797 : 41030 : double L = - pariINFINITY;
798 [ + + ]: 114460 : for (i=0; i<=degpol(q); i++)
799 : : {
800 : 73430 : double d = dbllog2(gel(q,2+i));
801 [ + + ]: 73430 : if (d > L) { L = d; k = i; }
802 : : }
803 : 41030 : return k;
804 : : }
805 : :
806 : : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
807 : : * Assume that l <= k <= n-l */
808 : : static long
809 : 50558 : dual_modulus(GEN p, double lrho, double tau, long l)
810 : : {
811 : 50558 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
812 : 50558 : double tau2 = tau * 7./8.;
813 : 50558 : pari_sp av = avma;
814 : : GEN q;
815 : :
816 : 50558 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
817 : 50558 : q = homothetie(p, lrho, bit);
818 : 50558 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
819 : :
820 [ + + ]: 445053 : for (i=0; i<imax; i++)
821 : : {
822 : 404023 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
823 : 404023 : v = RgX_valrem(q, &q);
824 [ + + ]: 404023 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
825 : :
826 : 404023 : nn = degpol(q); delta_k += v;
827 [ + + ]: 404023 : if (!nn) return delta_k;
828 : :
829 : 394495 : set_karasquare_limit(bit);
830 : 394495 : q = gerepileupto(av, graeffe(q));
831 : 394495 : tau2 *= 7./4.;
832 : 394495 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
833 : : }
834 : 50558 : avma = av; return delta_k + (long)ind_maxlog2(q);
835 : : }
836 : :
837 : : /********************************************************************/
838 : : /** **/
839 : : /** FACTORS THROUGH CIRCLE INTEGRATION **/
840 : : /** **/
841 : : /********************************************************************/
842 : : /* l power of 2 */
843 : : static void
844 : 1559136 : fft(GEN Omega, GEN p, GEN f, long step, long l)
845 : : {
846 : : pari_sp ltop;
847 : : long i, l1, l2, l3, rapi, step4;
848 : : GEN f1, f2, f3, f02, f13, g02, g13, ff;
849 : :
850 [ + + ]: 1559136 : if (l == 2)
851 : : {
852 : 893508 : gel(f,0) = gadd(gel(p,0),gel(p,step));
853 : 893508 : gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
854 : : }
855 [ + + ]: 665628 : if (l == 4)
856 : : {
857 : 347356 : f1 = gadd(gel(p,0), gel(p,step<<1));
858 : 347356 : f2 = gsub(gel(p,0), gel(p,step<<1));
859 : 347356 : f3 = gadd(gel(p,step),gel(p,3*step));
860 : 347356 : f02= gsub(gel(p,step),gel(p,3*step));
861 : 347356 : f02 = mulcxI(f02);
862 : 347356 : gel(f,0) = gadd(f1, f3);
863 : 347356 : gel(f,1) = gadd(f2, f02);
864 : 347356 : gel(f,2) = gsub(f1, f3);
865 : 347356 : gel(f,3) = gsub(f2, f02); return;
866 : : }
867 : :
868 : 318272 : ltop = avma;
869 : 318272 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
870 : 318272 : fft(Omega,p, f, step4,l1);
871 : 318272 : fft(Omega,p+step, f+l1,step4,l1);
872 : 318272 : fft(Omega,p+(step<<1),f+l2,step4,l1);
873 : 318272 : fft(Omega,p+3*step, f+l3,step4,l1);
874 : :
875 : 318272 : ff = cgetg(l+1,t_VEC);
876 [ + + ]: 1277958 : for (i=0; i<l1; i++)
877 : : {
878 : 959686 : rapi = step*i;
879 : 959686 : f1 = gmul(gel(Omega,rapi), gel(f,i+l1));
880 : 959686 : f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
881 : 959686 : f3 = gmul(gel(Omega,3*rapi), gel(f,i+l3));
882 : :
883 : 959686 : f02 = gadd(gel(f,i),f2);
884 : 959686 : g02 = gsub(gel(f,i),f2);
885 : 959686 : f13 = gadd(f1,f3);
886 : 959686 : g13 = mulcxI(gsub(f1,f3));
887 : :
888 : 959686 : gel(ff,i+1) = gadd(f02, f13);
889 : 959686 : gel(ff,i+l1+1) = gadd(g02, g13);
890 : 959686 : gel(ff,i+l2+1) = gsub(f02, f13);
891 : 959686 : gel(ff,i+l3+1) = gsub(g02, g13);
892 : : }
893 : 318272 : ff = gerepilecopy(ltop,ff);
894 [ + + ]: 5397880 : for (i=0; i<l; i++) f[i] = ff[i+1];
895 : : }
896 : :
897 : : /* e(1/N) */
898 : : static GEN
899 : 101666 : RUgen(long N, long bit)
900 : : {
901 [ + + ]: 101666 : if (N == 2) return real_m1(nbits2prec(bit));
902 [ + + ]: 97911 : if (N == 4) return gen_I();
903 : 101666 : return expIr(divru(Pi2n(1, nbits2prec(bit)), N));
904 : : }
905 : :
906 : : /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
907 : : static GEN
908 : 48244 : initRU(long N, long bit)
909 : : {
910 : 48244 : GEN *RU, z = RUgen(N, bit);
911 : 48244 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
912 : :
913 : 48244 : RU = (GEN*)cgetg(N+1,t_VEC); RU++;
914 : :
915 : 48244 : RU[0] = myreal_1(bit);
916 : 48244 : RU[1] = z;
917 [ + + ]: 66346 : for (i=1; i<N8; i++)
918 : : {
919 : 18102 : GEN t = RU[i];
920 : 18102 : RU[i+1] = gmul(z, t);
921 : 18102 : RU[N4-i] = mkcomplex(gel(t,2), gel(t,1));
922 : : }
923 [ + + ]: 168712 : for (i=0; i<N4; i++) RU[i+N4] = mulcxI(RU[i]);
924 [ + + ]: 289180 : for (i=0; i<N2; i++) RU[i+N2] = gneg(RU[i]);
925 : 48244 : return (GEN)RU;
926 : : }
927 : :
928 : : /* as above, N arbitrary */
929 : : static GEN
930 : 5178 : initRUgen(long N, long bit)
931 : : {
932 : 5178 : GEN *RU = (GEN*)cgetg(N+1,t_VEC), z = RUgen(N,bit);
933 : 5178 : long i, k = (N+3)>>1;
934 : 5178 : RU[0] = gen_1;
935 : 5178 : RU[1] = z;
936 [ + + ]: 7857 : for (i=2; i<k; i++) RU[i] = gmul(z, RU[i-1]);
937 [ + + ]: 7233 : for ( ; i<N; i++) RU[i] = gconj(RU[N-i]);
938 : 5178 : return (GEN)RU;
939 : : }
940 : :
941 : : GEN
942 : 0 : FFTinit(long k, long prec)
943 : : {
944 [ # # ]: 0 : if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
945 : 0 : return initRU(1L << k, prec2nbits(prec)) - 1;
946 : : }
947 : :
948 : : GEN
949 : 0 : FFT(GEN x, GEN Omega)
950 : : {
951 : 0 : long i, l = lg(Omega), n = lg(x);
952 : : GEN y, z;
953 [ # # ]: 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
954 [ # # ]: 0 : if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
955 [ # # ]: 0 : if (n > l) pari_err_DIM("FFT");
956 : :
957 [ # # ]: 0 : if (n < l) {
958 : 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
959 [ # # ]: 0 : for (i = 1; i < n; i++) z[i] = x[i];
960 [ # # ]: 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
961 : : }
962 : 0 : else z = x;
963 : 0 : y = cgetg(l, t_VEC);
964 : 0 : fft(Omega+1, z+1, y+1, 1, l-1);
965 : 0 : return y;
966 : : }
967 : :
968 : : /* returns 1 if p has only real coefficients, 0 else */
969 : : static int
970 : 44451 : isreal(GEN p)
971 : : {
972 : : long i;
973 [ + + ]: 284024 : for (i = lg(p)-1; i > 1; i--)
974 [ + + ]: 250538 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
975 : 44451 : return 1;
976 : : }
977 : :
978 : : /* x non complex */
979 : : static GEN
980 : 33042 : abs_update_r(GEN x, double *mu) {
981 : 33042 : GEN y = gtofp(x, DEFAULTPREC);
982 [ + + ]: 33042 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
983 : 33042 : setabssign(y); return y;
984 : : }
985 : : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
986 : : static GEN
987 : 676504 : abs_update(GEN x, double *mu) {
988 : : GEN y, xr, yr;
989 : : double ly;
990 [ + + ]: 676504 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
991 : 644456 : xr = gel(x,1);
992 : 644456 : yr = gel(x,2);
993 [ + + ]: 644456 : if (gequal0(xr)) return abs_update_r(yr,mu);
994 [ + + ]: 644330 : if (gequal0(yr)) return abs_update_r(xr,mu);
995 : : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
996 : 643462 : xr = gtofp(xr, DEFAULTPREC);
997 : 643462 : yr = gtofp(yr, DEFAULTPREC);
998 : 643462 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
999 [ + + ]: 643462 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
1000 : 676504 : return y;
1001 : : }
1002 : :
1003 : : static void
1004 : 23247 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
1005 : : int polreal, double param, double param2)
1006 : : {
1007 : : GEN q, pc, Omega, A, RU, prim, g, ONE,TWO;
1008 : 23247 : long n = degpol(p), bit, NN, K, i, j, Lmax;
1009 : 23247 : pari_sp av2, av = avma;
1010 : :
1011 : 23247 : bit = gexpo(p) + (long)param2+8;
1012 [ + + ]: 46715 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
1013 [ + + ]: 23247 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
1014 [ + + ]: 23247 : K = NN/Lmax; if (K & 1) K++;
1015 : 23247 : NN = Lmax*K;
1016 [ + + ]: 23247 : if (polreal) K = K/2+1;
1017 : :
1018 : 23247 : Omega = initRU(Lmax,bit);
1019 : 23247 : prim = RUgen(NN, bit);
1020 : :
1021 : 23247 : q = mygprec(p,bit) + 2;
1022 : 23247 : A = cgetg(Lmax+1,t_VEC); A++;
1023 : 23247 : pc= cgetg(Lmax+1,t_VEC); pc++;
1024 [ + + ]: 185492 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
1025 [ + + ]: 83666 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
1026 : :
1027 : 23247 : *mu = pariINFINITY;
1028 : 23247 : g = real_0_bit(-bit);
1029 : 23247 : ONE = real_1(DEFAULTPREC);
1030 : 23247 : TWO = real2n(1, DEFAULTPREC);
1031 : 23247 : av2 = avma;
1032 : 23247 : RU = myreal_1(bit);
1033 [ + + ]: 88719 : for (i=0; i<K; i++)
1034 : : {
1035 [ + + ]: 65472 : if (i) {
1036 : 42225 : GEN z = RU;
1037 [ + + ]: 301118 : for (j=1; j<n; j++)
1038 : : {
1039 : 258893 : gel(pc,j) = gmul(gel(q,j),z);
1040 : 258893 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
1041 : : }
1042 : 42225 : gel(pc,n) = gmul(gel(q,n),z);
1043 : : }
1044 : :
1045 : 65472 : fft(Omega,pc,A,1,Lmax);
1046 [ + + ][ + + ]: 75308 : if (polreal && i>0 && i<K-1)
[ + + ]
1047 [ + + ]: 164708 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
1048 : : else
1049 [ + + ]: 577268 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(ONE, abs_update(gel(A,j),mu)));
1050 : 65472 : RU = gmul(RU, prim);
1051 [ - + ]: 65472 : if (gc_needed(av,1))
1052 : : {
1053 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
1054 : 0 : gerepileall(av2,2, &g,&RU);
1055 : : }
1056 : : }
1057 : 23247 : *gamma = mydbllog2r(divru(g,NN));
1058 : 23247 : *LMAX = Lmax; avma = av;
1059 : 23247 : }
1060 : :
1061 : : /* NN is a multiple of Lmax */
1062 : : static void
1063 : 24997 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
1064 : : {
1065 : : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
1066 : 24997 : long n = degpol(p), i, j, K;
1067 : : pari_sp ltop;
1068 : :
1069 : 24997 : Omega = initRU(Lmax,bit);
1070 : 24997 : prim = RUgen(NN, bit);
1071 : 24997 : RU = cgetg(n+2,t_VEC); RU++;
1072 : :
1073 [ + + ]: 24997 : K = NN/Lmax; if (polreal) K = K/2+1;
1074 : 24997 : q = mygprec(p,bit);
1075 : 24997 : qd = RgX_deriv(q);
1076 : :
1077 : 24997 : A = cgetg(Lmax+1,t_VEC); A++;
1078 : 24997 : B = cgetg(Lmax+1,t_VEC); B++;
1079 : 24997 : C = cgetg(Lmax+1,t_VEC); C++;
1080 : 24997 : pc = cgetg(Lmax+1,t_VEC); pc++;
1081 : 24997 : pd = cgetg(Lmax+1,t_VEC); pd++;
1082 [ + + ]: 95017 : pc[0] = q[2]; for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
1083 [ + + ]: 120014 : pd[0] = qd[2]; for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
1084 : :
1085 : 24997 : ltop = avma;
1086 : 24997 : W = cgetg(k+1,t_VEC);
1087 : 24997 : U = cgetg(k+1,t_VEC);
1088 [ + + ]: 75572 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
1089 : :
1090 : 24997 : gel(RU,0) = gen_1;
1091 : 24997 : prim2 = myreal_1(bit);
1092 [ + + ]: 80141 : for (i=0; i<K; i++)
1093 : : {
1094 : 55144 : gel(RU,1) = prim2;
1095 [ + + ]: 401624 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
1096 : : /* RU[j] = prim^(ij)= prim2^j */
1097 : :
1098 [ + + ]: 401624 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
1099 : 55144 : fft(Omega,pd,A,1,Lmax);
1100 [ + + ]: 456768 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
1101 : 55144 : fft(Omega,pc,B,1,Lmax);
1102 [ + + ]: 680128 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
1103 [ + + ]: 680128 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
1104 : 55144 : fft(Omega,B,A,1,Lmax);
1105 : 55144 : fft(Omega,C,B,1,Lmax);
1106 : :
1107 [ + + ]: 55144 : if (polreal) /* p has real coefficients */
1108 : : {
1109 [ + + ][ + + ]: 44252 : if (i>0 && i<K-1)
1110 : : {
1111 [ + + ]: 23696 : for (j=1; j<=k; j++)
1112 : : {
1113 : 19267 : gel(W,j) = gadd(gel(W,j), gshift(real_i(gmul(gel(A,j+1),gel(RU,j+1))),1));
1114 : 19267 : gel(U,j) = gadd(gel(U,j), gshift(real_i(gmul(gel(B,j),gel(RU,j))),1));
1115 : : }
1116 : : }
1117 : : else
1118 : : {
1119 [ + + ]: 116595 : for (j=1; j<=k; j++)
1120 : : {
1121 : 76772 : gel(W,j) = gadd(gel(W,j), real_i(gmul(gel(A,j+1),gel(RU,j+1))));
1122 : 76772 : gel(U,j) = gadd(gel(U,j), real_i(gmul(gel(B,j),gel(RU,j))));
1123 : : }
1124 : : }
1125 : : }
1126 : : else
1127 : : {
1128 [ + + ]: 43808 : for (j=1; j<=k; j++)
1129 : : {
1130 : 28487 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
1131 : 28487 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
1132 : : }
1133 : : }
1134 : 55144 : prim2 = gmul(prim2,prim);
1135 : 55144 : gerepileall(ltop,3, &W,&U,&prim2);
1136 : : }
1137 : :
1138 [ + + ]: 75572 : for (i=1; i<=k; i++)
1139 : : {
1140 : 50575 : aux=gel(W,i);
1141 [ + + ]: 117912 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
1142 : 50575 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
1143 : : }
1144 [ + + ]: 75572 : for (i=0; i<k; i++)
1145 : : {
1146 : 50575 : aux=gel(U,k-i);
1147 [ + + ]: 117912 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
1148 : 50575 : gel(H,i+2) = gdivgs(aux,NN);
1149 : : }
1150 : 24997 : }
1151 : :
1152 : : #define NEWTON_MAX 10
1153 : : static GEN
1154 : 126563 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
1155 : : {
1156 : 126563 : GEN H = HH, D, aux;
1157 : 126563 : pari_sp ltop = avma;
1158 : : long error, i, bit1, bit2;
1159 : :
1160 : 126563 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
1161 : 126563 : bit2 = bit + Sbit;
1162 [ + + ][ + + ]: 236557 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1163 : : {
1164 [ - + ]: 109994 : if (gc_needed(ltop,1))
1165 : : {
1166 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
1167 : 0 : gerepileall(ltop,2, &D,&H);
1168 : : }
1169 : 109994 : bit1 = -error + Sbit;
1170 : 109994 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
1171 : 109994 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
1172 : :
1173 [ + + ]: 109994 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
1174 : 109994 : H = RgX_add(mygprec(H,bit1), aux);
1175 : 109994 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
1176 [ + + ]: 109994 : error = gexpo(D); if (error < -bit1) error = -bit1;
1177 : : }
1178 [ + + ]: 126563 : if (error > -bit/2) return NULL; /* FAIL */
1179 : 126563 : return gerepilecopy(ltop,H);
1180 : : }
1181 : :
1182 : : /* return 0 if fails, 1 else */
1183 : : static long
1184 : 24997 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
1185 : : {
1186 : 24997 : GEN f0, FF, GG, r, HH = H;
1187 : 24997 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
1188 : 24997 : pari_sp av = avma;
1189 : :
1190 : 24997 : FF = *F; GG = RgX_divrem(p, FF, &r);
1191 [ - + ]: 24997 : error = gexpo(r); if (error <= -bit) error = 1-bit;
1192 : 24997 : normF = gexpo(FF);
1193 : 24997 : normG = gexpo(GG);
1194 [ + + ]: 24997 : enh = gexpo(H); if (enh < 0) enh = 0;
1195 : 24997 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
1196 : 24997 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
1197 : 24997 : bit2 = bit + Sbit;
1198 [ + + ][ + + ]: 151511 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1199 : : {
1200 [ + + ][ + + ]: 126563 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
1201 [ - + ]: 126563 : if (gc_needed(av,1))
1202 : : {
1203 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1204 : 0 : gerepileall(av,4, &FF,&GG,&r,&HH);
1205 : : }
1206 : :
1207 : 126563 : bit1 = -error + Sbit2;
1208 : 126563 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1209 : : 1-error, Sbit2);
1210 [ + + ]: 126563 : if (!HH) return 0; /* FAIL */
1211 : :
1212 : 126514 : bit1 = -error + Sbit;
1213 : 126514 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1214 : 126514 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1215 : :
1216 [ + + ]: 126514 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1217 : 126514 : FF = gadd(mygprec(FF,bit1),f0);
1218 : :
1219 [ + + ]: 126514 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1220 : 126514 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1221 [ + + ]: 126514 : error = gexpo(r); if (error < -bit1) error = -bit1;
1222 : : }
1223 [ + + ]: 24948 : if (error>-bit) return 0; /* FAIL */
1224 : 24997 : *F = FF; *G = GG; return 1;
1225 : : }
1226 : :
1227 : : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1228 : : where cd is the leading coefficient of p */
1229 : : static void
1230 : 23247 : split_fromU(GEN p, long k, double delta, long bit,
1231 : : GEN *F, GEN *G, double param, double param2)
1232 : : {
1233 : : GEN pp, FF, GG, H;
1234 : 23247 : long n = degpol(p), NN, bit2, Lmax;
1235 : 23247 : int polreal = isreal(p);
1236 : : pari_sp ltop;
1237 : : double mu, gamma;
1238 : :
1239 : 23247 : pp = gdiv(p, gel(p,2+n));
1240 : 23247 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1241 : :
1242 : 23247 : H = cgetg(k+2,t_POL); H[1] = p[1];
1243 : 23247 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1244 : 23247 : gel(FF,k+2) = gen_1;
1245 : :
1246 [ + + ]: 23247 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1247 : 23247 : NN *= Lmax; ltop = avma;
1248 : : for(;;)
1249 : : {
1250 : 24997 : bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
1251 : 24997 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1252 [ + + ]: 24997 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1253 : 1750 : NN <<= 1; avma = ltop;
1254 : 1750 : }
1255 : 23247 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1256 : 23247 : }
1257 : :
1258 : : static void
1259 : 23247 : optimize_split(GEN p, long k, double delta, long bit,
1260 : : GEN *F, GEN *G, double param, double param2)
1261 : : {
1262 : 23247 : long n = degpol(p);
1263 : : GEN FF, GG;
1264 : :
1265 [ + + ]: 23247 : if (k <= n/2)
1266 : 17867 : split_fromU(p,k,delta,bit,F,G,param,param2);
1267 : : else
1268 : : {
1269 : 5380 : split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1270 : 5380 : *F = RgX_recip_shallow(GG);
1271 : 5380 : *G = RgX_recip_shallow(FF);
1272 : : }
1273 : 23247 : }
1274 : :
1275 : : /********************************************************************/
1276 : : /** **/
1277 : : /** SEARCH FOR SEPARATING CIRCLE **/
1278 : : /** **/
1279 : : /********************************************************************/
1280 : :
1281 : : /* return p(2^e*x) *2^(-n*e) */
1282 : : static void
1283 : 0 : scalepol2n(GEN p, long e)
1284 : : {
1285 : 0 : long i,n=lg(p)-1;
1286 [ # # ]: 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1287 : 0 : }
1288 : :
1289 : : /* returns p(x/R)*R^n */
1290 : : static GEN
1291 : 196965 : scalepol(GEN p, GEN R, long bit)
1292 : : {
1293 : : GEN q,aux,gR;
1294 : : long i;
1295 : :
1296 : 196965 : aux = gR = mygprec(R,bit); q = mygprec(p,bit);
1297 [ + + ]: 944001 : for (i=lg(p)-2; i>=2; i--)
1298 : : {
1299 : 747036 : gel(q,i) = gmul(aux,gel(q,i));
1300 : 747036 : aux = gmul(aux,gR);
1301 : : }
1302 : 196965 : return q;
1303 : : }
1304 : :
1305 : : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1306 : : static GEN
1307 : 63612 : conformal_pol(GEN p, GEN a, long bit)
1308 : : {
1309 : 63612 : GEN z, r, ma = gneg(a), ca = gconj(a);
1310 : 63612 : long n = degpol(p), i;
1311 : 63612 : pari_sp av = avma;
1312 : :
1313 : 63612 : z = mkpoln(2, ca, negr(myreal_1(bit)));
1314 : 63612 : r = scalarpol(gel(p,2+n), 0);
1315 : 63612 : for (i=n-1; ; i--)
1316 : : {
1317 : 234520 : r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
1318 : 234520 : r = gadd(r, gmul(z, gel(p,2+i)));
1319 [ + + ]: 234520 : if (i == 0) return gerepileupto(av, r);
1320 : 170908 : z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1321 [ - + ]: 170908 : if (gc_needed(av,2))
1322 : : {
1323 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
1324 : 0 : gerepileall(av,2, &r,&z);
1325 : : }
1326 : 170908 : }
1327 : : }
1328 : :
1329 : : static const double UNDEF = -100000.;
1330 : :
1331 : : static double
1332 : 23247 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1333 : : {
1334 : 23247 : long i, n = degpol(p);
1335 : : double lrho, lrmin, lrmax;
1336 [ + + ]: 23247 : if (k > 1)
1337 : : {
1338 [ + - ][ + + ]: 24664 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1339 : 15532 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1340 : : }
1341 : : else /* k=1 */
1342 : 7715 : lrmin = logmin_modulus(p,aux);
1343 : 23247 : radii[k] = lrmin;
1344 : :
1345 [ + + ]: 23247 : if (k+1<n)
1346 : : {
1347 [ + - ][ + + ]: 45228 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1348 : 19922 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1349 : : }
1350 : : else /* k+1=n */
1351 : 3325 : lrmax = logmax_modulus(p,aux);
1352 : 23247 : radii[k+1] = lrmax;
1353 : :
1354 : 23247 : lrho = radii[k];
1355 [ + + ]: 55979 : for (i=k-1; i>=1; i--)
1356 : : {
1357 [ + + ][ + + ]: 32732 : if (radii[i] == UNDEF || radii[i] > lrho)
1358 : 22976 : radii[i] = lrho;
1359 : : else
1360 : 9756 : lrho = radii[i];
1361 : : }
1362 : 23247 : lrho = radii[k+1];
1363 [ + + ]: 106266 : for (i=k+1; i<=n; i++)
1364 : : {
1365 [ + + ][ + + ]: 83019 : if (radii[i] == UNDEF || radii[i] < lrho)
1366 : 51135 : radii[i] = lrho;
1367 : : else
1368 : 31884 : lrho = radii[i];
1369 : : }
1370 : 23247 : *delta = (lrmax - lrmin) / 2;
1371 [ + + ]: 23247 : if (*delta > 1.) *delta = 1.;
1372 : 23247 : return (lrmin + lrmax) / 2;
1373 : : }
1374 : :
1375 : : static void
1376 : 23247 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1377 : : {
1378 : 23247 : double t, param = 0., param2 = 0.;
1379 : : long i;
1380 [ + + ]: 162245 : for (i=1; i<=n; i++)
1381 : : {
1382 : 138998 : radii[i] -= lrho;
1383 : 138998 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1384 [ + + ]: 138998 : param += t; if (t > 1.) param2 += log2(t);
1385 : : }
1386 : 23247 : *par = param; *par2 = param2;
1387 : 23247 : }
1388 : :
1389 : : /* apply the conformal mapping then split from U */
1390 : : static void
1391 : 21204 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1392 : : double aux, GEN *F,GEN *G)
1393 : : {
1394 : 21204 : long bit2, n = degpol(p), i;
1395 : 21204 : pari_sp ltop = avma, av;
1396 : : GEN q, FF, GG, a, R;
1397 : : double lrho, delta, param, param2;
1398 : : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1399 : 21204 : bit2 = bit + (long)(n*3.4848775) + 1;
1400 : 21204 : a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
1401 : 21204 : a = divrs(a, -6);
1402 : 21204 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1403 : :
1404 : 21204 : av = avma;
1405 : 21204 : q = conformal_pol(mygprec(p,bit2), a, bit2);
1406 [ + + ]: 138464 : for (i=1; i<=n; i++)
1407 [ + + ]: 117260 : if (radii[i] != UNDEF) /* update array radii */
1408 : : {
1409 : 79368 : pari_sp av2 = avma;
1410 : 79368 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1411 : : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1412 : 79368 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1413 : 79368 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1414 : 79368 : avma = av2;
1415 : : }
1416 : 21204 : lrho = logradius(radii, q,k,aux/10., &delta);
1417 : 21204 : update_radius(n, radii, lrho, ¶m, ¶m2);
1418 : :
1419 : 21204 : bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
1420 : 21204 : R = mygprec(dblexp(-lrho), bit2);
1421 : 21204 : q = scalepol(q,R,bit2);
1422 : 21204 : gerepileall(av,2, &q,&R);
1423 : :
1424 : 21204 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1425 : 21204 : bit2 += n; R = invr(R);
1426 : 21204 : FF = scalepol(FF,R,bit2);
1427 : 21204 : GG = scalepol(GG,R,bit2);
1428 : :
1429 : 21204 : a = mygprec(a,bit2);
1430 : 21204 : FF = conformal_pol(FF,a,bit2);
1431 : 21204 : GG = conformal_pol(GG,a,bit2);
1432 : :
1433 : 21204 : a = invr(subsr(1, gnorm(a)));
1434 : 21204 : FF = RgX_Rg_mul(FF, powru(a,k));
1435 : 21204 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1436 : :
1437 : 21204 : *F = mygprec(FF,bit+n);
1438 : 21204 : *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1439 : 21204 : }
1440 : :
1441 : : /* split p, this time without scaling. returns in F and G two polynomials
1442 : : * such that |p-FG|< 2^(-bit)|p| */
1443 : : static void
1444 : 23247 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1445 : : {
1446 : : GEN q, FF, GG, R;
1447 : : double aux, delta, param, param2;
1448 : 23247 : long n = degpol(p), i, j, k, bit2;
1449 : : double lrmin, lrmax, lrho, *radii;
1450 : :
1451 : 23247 : init_dalloc();
1452 : 23247 : radii = (double*) stack_malloc((n+1) * sizeof(double));
1453 : :
1454 [ + + ]: 115751 : for (i=2; i<n; i++) radii[i] = UNDEF;
1455 : 23247 : aux = thickness/(double)(4 * n);
1456 : 23247 : lrmin = logmin_modulus(p, aux);
1457 : 23247 : lrmax = logmax_modulus(p, aux);
1458 : 23247 : radii[1] = lrmin;
1459 : 23247 : radii[n] = lrmax;
1460 : 23247 : i = 1; j = n;
1461 : 23247 : lrho = (lrmin + lrmax) / 2;
1462 : 23247 : k = dual_modulus(p, lrho, aux, 1);
1463 [ + + ][ + + ]: 23247 : if (5*k < n || (n < 2*k && 5*k < 4*n))
[ + + ]
1464 : 4988 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1465 : : else
1466 : 18259 : { lrmin = lrho; i=k; radii[i] = lrho; }
1467 [ + + ]: 50558 : while (j > i+1)
1468 : : {
1469 [ + + ]: 27311 : if (i+j == n+1)
1470 : 9796 : lrho = (lrmin + lrmax) / 2;
1471 : : else
1472 : : {
1473 : 17515 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1474 [ + + ]: 17515 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1475 : 12901 : else lrho = lrmin * kappa + lrmax;
1476 : 17515 : lrho /= 1+kappa;
1477 : : }
1478 : 27311 : aux = (lrmax - lrmin) / (4*(j-i));
1479 : 27311 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1480 [ + + ][ + + ]: 27311 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
[ + + ]
1481 : 19080 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1482 : : else
1483 : 8231 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1484 : : }
1485 : 23247 : aux = lrmax - lrmin;
1486 : :
1487 [ + + ]: 23247 : if (ctr)
1488 : : {
1489 : 21204 : lrho = (lrmax + lrmin) / 2;
1490 [ + + ]: 138464 : for (i=1; i<=n; i++)
1491 [ + + ]: 117260 : if (radii[i] != UNDEF) radii[i] -= lrho;
1492 : :
1493 : 21204 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1494 : 21204 : R = mygprec(dblexp(-lrho), bit2);
1495 : 21204 : q = scalepol(p,R,bit2);
1496 : 21204 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1497 : : }
1498 : : else
1499 : : {
1500 : 2043 : lrho = logradius(radii, p, k, aux/10., &delta);
1501 : 2043 : update_radius(n, radii, lrho, ¶m, ¶m2);
1502 : :
1503 : 2043 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1504 : 2043 : R = mygprec(dblexp(-lrho), bit2);
1505 : 2043 : q = scalepol(p,R,bit2);
1506 : 2043 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1507 : : }
1508 : 23247 : bit += n;
1509 : 23247 : bit2 += n; R = invr(mygprec(R,bit2));
1510 : 23247 : *F = mygprec(scalepol(FF,R,bit2), bit);
1511 : 23247 : *G = mygprec(scalepol(GG,R,bit2), bit);
1512 : 23247 : }
1513 : :
1514 : : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1515 : : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1516 : : * where the maximum modulus of the roots of p is <=1.
1517 : : * Assume sum of roots is 0. */
1518 : : static void
1519 : 21204 : split_1(GEN p, long bit, GEN *F, GEN *G)
1520 : : {
1521 : 21204 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1522 : : GEN TWO, ctr, q, qq, FF, GG, v, gr, r, newq;
1523 : : double lrmin, lrmax, lthick;
1524 : 21204 : const double LOG3 = 1.098613;
1525 : :
1526 : 21204 : lrmax = logmax_modulus(p, 0.01);
1527 : 21204 : gr = mygprec(dblexp(-lrmax), bit2);
1528 : 21204 : q = scalepol(p,gr,bit2);
1529 : :
1530 : 21204 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1531 : 21204 : TWO = myreal_1(bit2); setexpo(TWO,1);
1532 : 21204 : v = cgetg(5,t_VEC);
1533 : 21204 : gel(v,1) = TWO;
1534 : 21204 : gel(v,2) = negr(TWO);
1535 : 21204 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1536 : 21204 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1537 : 21204 : q = mygprec(q,bit2); lthick = 0;
1538 : 21204 : newq = ctr = NULL; /* -Wall */
1539 [ + + ]: 21204 : imax = polreal? 3: 4;
1540 [ + + ]: 40170 : for (i=1; i<=imax; i++)
1541 : : {
1542 : 39659 : qq = RgX_translate(q, gel(v,i));
1543 : 39659 : lrmin = logmin_modulus(qq,0.05);
1544 [ + + ]: 39659 : if (LOG3 > lrmin + lthick)
1545 : : {
1546 : 38656 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1547 [ + + ]: 38656 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1548 : : }
1549 [ + + ]: 39659 : if (lthick > LOG2) break;
1550 [ + + ][ + + ]: 22353 : if (polreal && i==2 && lthick > LOG3 - LOG2) break;
[ + + ]
1551 : : }
1552 : 21204 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
1553 : 21204 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1554 : 21204 : r = gneg(mygprec(ctr,bit2));
1555 : 21204 : FF = RgX_translate(FF,r);
1556 : 21204 : GG = RgX_translate(GG,r);
1557 : :
1558 : 21204 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1559 : 21204 : *F = scalepol(FF,gr,bit2);
1560 : 21204 : *G = scalepol(GG,gr,bit2);
1561 : 21204 : }
1562 : :
1563 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1564 : : where the maximum modulus of the roots of p is < 0.5 */
1565 : : static int
1566 : 21204 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1567 : : {
1568 : : GEN q, b, FF, GG;
1569 : 21204 : long n = degpol(p), k, bit2, eq;
1570 : 21204 : double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
1571 : :
1572 : : /* beware double overflow */
1573 [ + + ][ + - ]: 21204 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
[ - + ]
1574 : :
1575 [ + + ]: 21204 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1576 : 21204 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1577 : :
1578 : 21204 : q = mygprec(p,bit2);
1579 : 21204 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1580 : 21204 : q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
1581 : 21204 : k = 0;
1582 [ + - ][ - + ]: 21204 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1583 [ - + ]: 21204 : || gequal0(gel(q,k+2)))) k++;
1584 [ - + ]: 21204 : if (k > 0)
1585 : : {
1586 [ # # ]: 0 : if (k > n/2) k = n/2;
1587 : 0 : bit2 += k<<1;
1588 : 0 : FF = monomial(myreal_1(bit2), k, 0);
1589 : 0 : GG = RgX_shift_shallow(q, -k);
1590 : : }
1591 : : else
1592 : : {
1593 : 21204 : split_1(q,bit2,&FF,&GG);
1594 : 21204 : bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
1595 : 21204 : FF = mygprec(FF,bit2);
1596 : : }
1597 : 21204 : GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
1598 : 21204 : *F = RgX_translate(FF, b);
1599 : 21204 : *G = RgX_translate(GG, b); return 1;
1600 : : }
1601 : :
1602 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1603 : : * Assume max_modulus(p) < 2 */
1604 : : static void
1605 : 21204 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1606 : : {
1607 : : GEN FF, GG;
1608 : : long n, bit2, normp;
1609 : :
1610 [ - + ]: 42408 : if (split_0_2(p,bit,F,G)) return;
1611 : :
1612 : 0 : normp = gexpo(p);
1613 : 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1614 : 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1615 : 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1616 : 0 : scalepol2n(FF,-2);
1617 : 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1618 : 0 : *F = mygprec(FF,bit2);
1619 : 21204 : *G = mygprec(GG,bit2);
1620 : : }
1621 : :
1622 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1623 : : static void
1624 : 23247 : split_0(GEN p, long bit, GEN *F, GEN *G)
1625 : : {
1626 : 23247 : const double LOG1_9 = 0.6418539;
1627 : 23247 : long n = degpol(p), k = 0;
1628 : : GEN q;
1629 : :
1630 [ - + ][ # # ]: 23247 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1631 [ - + ]: 23247 : if (k > 0)
1632 : : {
1633 [ # # ]: 0 : if (k > n/2) k = n/2;
1634 : 0 : *F = monomial(myreal_1(bit), k, 0);
1635 : 0 : *G = RgX_shift_shallow(p, -k);
1636 : : }
1637 : : else
1638 : : {
1639 : 23247 : double lr = logmax_modulus(p, 0.05);
1640 [ + + ]: 23247 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1641 : : else
1642 : : {
1643 : 18396 : q = RgX_recip_shallow(p);
1644 : 18396 : lr = logmax_modulus(q,0.05);
1645 [ + + ]: 18396 : if (lr < LOG1_9)
1646 : : {
1647 : 16353 : split_0_1(q, bit, F, G);
1648 : 16353 : *F = RgX_recip_shallow(*F);
1649 : 16353 : *G = RgX_recip_shallow(*G);
1650 : : }
1651 : : else
1652 : 2043 : split_2(p,bit,NULL, 1.2837,F,G);
1653 : : }
1654 : : }
1655 : 23247 : }
1656 : :
1657 : : /********************************************************************/
1658 : : /** **/
1659 : : /** ERROR ESTIMATE FOR THE ROOTS **/
1660 : : /** **/
1661 : : /********************************************************************/
1662 : :
1663 : : static GEN
1664 : 75756 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1665 : : {
1666 : 75756 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1667 : : long i, j;
1668 : :
1669 : 75756 : d = cgetg(n+1,t_VEC);
1670 [ + + ]: 1164416 : for (i=1; i<=n; i++)
1671 : : {
1672 [ + + ]: 1088660 : if (i!=k)
1673 : : {
1674 : 1012904 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1675 : 1012904 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1676 : : }
1677 : : }
1678 : 75756 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1679 [ + + ]: 75756 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1680 : 75756 : eps = mulrr(rho, shatzle);
1681 : 75756 : aux = shiftr(powru(rho,n), err);
1682 : :
1683 [ + + ][ + - ]: 239779 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
[ + + ]
1684 : : {
1685 : 164023 : GEN prod = NULL; /* 1. */
1686 : 164023 : long m = n;
1687 : 164023 : epsbis = mulrr(eps, dbltor(1.25));
1688 [ + + ]: 2905977 : for (i=1; i<=n; i++)
1689 : : {
1690 [ + + ][ + + ]: 2741954 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1691 : : {
1692 : 2553259 : GEN dif = subrr(gel(d,i),eps);
1693 [ + + ]: 2553259 : prod = prod? mulrr(prod, dif): dif;
1694 : 2553259 : m--;
1695 : : }
1696 : : }
1697 [ + + ]: 164023 : eps2 = prod? divrr(aux, prod): aux;
1698 [ + + ]: 164023 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1699 : 164023 : rap = divrr(eps,eps2); eps = eps2;
1700 : : }
1701 : 75756 : return eps;
1702 : : }
1703 : :
1704 : : /* round a complex or real number x to an absolute value of 2^(-bit) */
1705 : : static GEN
1706 : 179058 : mygprec_absolute(GEN x, long bit)
1707 : : {
1708 : : long e;
1709 : : GEN y;
1710 : :
1711 [ + + + ]: 179058 : switch(typ(x))
1712 : : {
1713 : : case t_REAL:
1714 : 121000 : e = expo(x) + bit;
1715 [ + + ][ - + ]: 121000 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1716 : : case t_COMPLEX:
1717 [ + + ]: 52358 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1718 : 50944 : y = cgetg(3,t_COMPLEX);
1719 : 50944 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1720 : 50944 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1721 : 50944 : return y;
1722 : 179058 : default: return x;
1723 : : }
1724 : : }
1725 : :
1726 : : static long
1727 : 11582 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1728 : : {
1729 : 11582 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1730 : : GEN sigma, shatzle;
1731 : :
1732 : 11582 : err += (long)log2((double)n) + 1;
1733 [ - + ]: 11582 : if (err > -2) return 0;
1734 : 11582 : sigma = real2n(-err, LOWDEFAULTPREC);
1735 : : /* 2 / ((s - 1)^(1/n) - 1) */
1736 : 11582 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1737 [ + + ]: 87338 : for (i=1; i<=n; i++)
1738 : : {
1739 : 75756 : pari_sp av = avma;
1740 : 75756 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1741 : 75756 : long e = gexpo(x);
1742 [ + + ]: 75756 : avma = av; if (e > e_max) e_max = e;
1743 : 75756 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1744 : : }
1745 : 11582 : return e_max;
1746 : : }
1747 : :
1748 : : /********************************************************************/
1749 : : /** **/
1750 : : /** MAIN **/
1751 : : /** **/
1752 : : /********************************************************************/
1753 : : static GEN
1754 : 55683 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1755 : :
1756 : : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1757 : : * returns prod (x-roots_pol[i]) */
1758 : : static GEN
1759 : 58076 : split_complete(GEN p, long bit, GEN roots_pol)
1760 : : {
1761 : 58076 : long n = degpol(p);
1762 : : pari_sp ltop;
1763 : : GEN p1, F, G, a, b, m1, m2;
1764 : :
1765 [ + + ]: 58076 : if (n == 1)
1766 : : {
1767 : 13975 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1768 : 13975 : (void)append_clone(roots_pol,a); return p;
1769 : : }
1770 : 44101 : ltop = avma;
1771 [ + + ]: 44101 : if (n == 2)
1772 : : {
1773 : 20854 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1774 : 20854 : F = gsqrt(F, nbits2prec(bit));
1775 : 20854 : p1 = ginv(gmul2n(gel(p,4),1));
1776 : 20854 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1777 : 20854 : b = gmul(gsub(F,gel(p,3)), p1);
1778 : 20854 : a = append_clone(roots_pol,a);
1779 : 20854 : b = append_clone(roots_pol,b); avma = ltop;
1780 : 20854 : a = mygprec(a, 3*bit);
1781 : 20854 : b = mygprec(b, 3*bit);
1782 : 20854 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1783 : : }
1784 : 23247 : split_0(p,bit,&F,&G);
1785 : 23247 : m1 = split_complete(F,bit,roots_pol);
1786 : 23247 : m2 = split_complete(G,bit,roots_pol);
1787 : 58076 : return gerepileupto(ltop, gmul(m1,m2));
1788 : : }
1789 : :
1790 : : static GEN
1791 : 350226 : quicktofp(GEN x)
1792 : : {
1793 : 350226 : const long prec = DEFAULTPREC;
1794 [ + + - + : 350226 : switch(typ(x))
- ]
1795 : : {
1796 : 349512 : case t_INT: return itor(x, prec);
1797 : 693 : case t_REAL: return rtor(x, prec);
1798 : 0 : case t_FRAC: return fractor(x, prec);
1799 : : case t_COMPLEX: {
1800 : 21 : GEN a = gel(x,1), b = gel(x,2);
1801 : : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1802 [ - + ]: 21 : if (isintzero(a)) return cxcompotor(b, prec);
1803 [ - + ]: 21 : if (isintzero(b)) return cxcompotor(a, prec);
1804 : 21 : a = cxcompotor(a, prec);
1805 : 21 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1806 : : }
1807 : 0 : default: pari_err_TYPE("quicktofp",x);
1808 : 350226 : return NULL;/*not reached*/
1809 : : }
1810 : :
1811 : : }
1812 : :
1813 : : /* bound log_2 |largest root of p| (Fujiwara's bound) */
1814 : : double
1815 : 71671 : fujiwara_bound(GEN p)
1816 : : {
1817 : 71671 : pari_sp av = avma;
1818 : 71671 : long i, n = degpol(p);
1819 : : GEN cc;
1820 : : double loglc, Lmax;
1821 : :
1822 [ - + ]: 71671 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1823 : 71671 : loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
1824 : 71671 : cc = gel(p, 2);
1825 [ + + ]: 71671 : if (gequal0(cc))
1826 : 5682 : Lmax = -pariINFINITY-1;
1827 : : else
1828 : 65989 : Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
1829 [ + + ]: 334243 : for (i = 1; i < n; i++)
1830 : : {
1831 : 262572 : GEN y = gel(p,i+2);
1832 : : double L;
1833 [ + + ]: 262572 : if (gequal0(y)) continue;
1834 : 212566 : L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
1835 [ + + ]: 212566 : if (L > Lmax) Lmax = L;
1836 : : }
1837 : 71671 : avma = av; return Lmax + 1;
1838 : : }
1839 : :
1840 : : /* Fujiwara's bound, real roots. Based on the following remark: if
1841 : : * p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
1842 : : * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
1843 : : * of q is a bound for the positive roots of p. */
1844 : : double
1845 : 8111 : fujiwara_bound_real(GEN p, long sign)
1846 : : {
1847 : 8111 : pari_sp av = avma;
1848 : : GEN x;
1849 : 8111 : long n = degpol(p), i, signodd, signeven;
1850 : : double fb;
1851 [ - + ]: 8111 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1852 : 8111 : x = shallowcopy(p);
1853 [ + - ]: 8111 : if (gsigne(gel(x, n+2)) > 0)
1854 : : {
1855 : 8111 : signeven = 1;
1856 : 8111 : signodd = sign;
1857 : : }
1858 : : else
1859 : : {
1860 : 0 : signeven = -1;
1861 : 0 : signodd = -sign;
1862 : : }
1863 [ + + ]: 49032 : for (i = 0; i < n; i++)
1864 : : {
1865 [ + + ]: 40921 : if ((n - i) % 2)
1866 : : {
1867 [ + + ]: 21988 : if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0;
1868 : : }
1869 : : else
1870 : : {
1871 [ + + ]: 18933 : if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0;
1872 : : }
1873 : : }
1874 : 8111 : fb = fujiwara_bound(x);
1875 : 8111 : avma = av; return fb;
1876 : : }
1877 : :
1878 : : static GEN
1879 : 87380 : mygprecrc_special(GEN x, long prec, long e)
1880 : : {
1881 : : GEN y;
1882 [ + + + ]: 87380 : switch(typ(x))
1883 : : {
1884 : : case t_REAL:
1885 [ - + ]: 726 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1886 [ + - ]: 726 : return (prec > realprec(x))? rtor(x, prec): x;
1887 : : case t_COMPLEX:
1888 : 21 : y = cgetg(3,t_COMPLEX);
1889 : 21 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1890 : 21 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1891 : 21 : return y;
1892 : 87380 : default: return x;
1893 : : }
1894 : : }
1895 : :
1896 : : /* like mygprec but keep at least the same precision as before */
1897 : : static GEN
1898 : 11582 : mygprec_special(GEN x, long bit)
1899 : : {
1900 : : long lx, i, e, prec;
1901 : : GEN y;
1902 : :
1903 [ - + ]: 11582 : if (bit < 0) bit = 0; /* should not happen */
1904 : 11582 : e = gexpo(x) - bit;
1905 : 11582 : prec = nbits2prec(bit);
1906 [ + - ]: 11582 : switch(typ(x))
1907 : : {
1908 : : case t_POL:
1909 : 11582 : y = cgetg_copy(x, &lx); y[1] = x[1];
1910 [ + + ]: 98920 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1911 : 11582 : break;
1912 : :
1913 : 0 : default: y = mygprecrc_special(x,prec,e);
1914 : : }
1915 : 11582 : return y;
1916 : : }
1917 : :
1918 : : static GEN
1919 : 6404 : fix_roots1(GEN r)
1920 : : {
1921 : 6404 : long i, l = lg(r);
1922 : 6404 : GEN allr = cgetg(l, t_VEC);
1923 [ + + ]: 48637 : for (i=1; i<l; i++)
1924 : : {
1925 : 42233 : GEN t = gel(r,i);
1926 : 42233 : gel(allr,i) = gcopy(t); gunclone(t);
1927 : : }
1928 : 6404 : return allr;
1929 : : }
1930 : : static GEN
1931 : 11582 : fix_roots(GEN r, GEN *m, long h, long bit)
1932 : : {
1933 : : long i, j, k, l, prec;
1934 : : GEN allr, ro1;
1935 [ + + ]: 11582 : if (h == 1) return fix_roots1(r);
1936 : 5178 : ro1 = initRUgen(h, bit);
1937 : 5178 : prec = nbits2prec(bit);
1938 : 5178 : l = lg(r)-1;
1939 : 5178 : allr = cgetg(h*l+1, t_VEC);
1940 [ + + ]: 18628 : for (k=1,i=1; i<=l; i++)
1941 : : {
1942 : 13450 : GEN p2, p1 = gel(r,i);
1943 [ + + ]: 13450 : p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1944 [ + + ]: 46973 : for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1945 : 13450 : gunclone(p1);
1946 : : }
1947 : 5178 : *m = roots_to_pol(allr, 0);
1948 : 11582 : return allr;
1949 : : }
1950 : :
1951 : : static GEN
1952 : 11416 : all_roots(GEN p, long bit)
1953 : : {
1954 : : GEN lc, pd, q, roots_pol, m;
1955 : 11416 : long bit0, bit2, n = degpol(p), i, e, h;
1956 : : pari_sp av;
1957 : :
1958 : 11416 : pd = RgX_deflate_max(p, &h); lc = leading_term(pd);
1959 [ + + ]: 11416 : e = (long)(2 * fujiwara_bound(pd)); if (e < 0) e = 0;
1960 : 11416 : bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
1961 : 11416 : bit2 = bit0; e = 0;
1962 : 11416 : for (av=avma,i=1;; i++,avma=av)
1963 : : {
1964 : 11582 : roots_pol = vectrunc_init(n+1);
1965 : 11582 : bit2 += e + (n << i);
1966 : 11582 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1967 : 11582 : q[1] = evalsigne(1)|evalvarn(0);
1968 : 11582 : m = split_complete(q,bit2,roots_pol);
1969 : 11582 : roots_pol = fix_roots(roots_pol, &m, h, bit2);
1970 : 11582 : q = mygprec_special(p,bit2); lc = leading_term(q);
1971 : 11582 : q[1] = evalsigne(1)|evalvarn(0);
1972 [ + + ]: 11582 : if (h > 1) m = gmul(m,lc);
1973 : :
1974 : 11582 : e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
1975 [ - + ]: 11582 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
1976 [ + - ]: 11582 : if (e < 0)
1977 : : {
1978 : 11582 : e = bit + a_posteriori_errors(p,roots_pol,e);
1979 [ + + ]: 11582 : if (e < 0) return roots_pol;
1980 : : }
1981 [ - + ]: 166 : if (DEBUGLEVEL > 7)
1982 : 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1983 : 166 : }
1984 : : }
1985 : :
1986 : : INLINE int
1987 : 2174 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1988 : :
1989 : : static int
1990 : 715 : isexactpol(GEN p)
1991 : : {
1992 : 715 : long i,n = degpol(p);
1993 [ + + ]: 2651 : for (i=0; i<=n; i++)
1994 [ + + ]: 2174 : if (!isexactscalar(gel(p,i+2))) return 0;
1995 : 715 : return 1;
1996 : : }
1997 : :
1998 : : static long
1999 : 3035 : isvalidcoeff(GEN x)
2000 : : {
2001 [ + + + ]: 3035 : switch (typ(x))
2002 : : {
2003 : 3000 : case t_INT: case t_REAL: case t_FRAC: return 1;
2004 [ + - ][ + - ]: 21 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
2005 : : }
2006 : 3035 : return 0;
2007 : : }
2008 : :
2009 : : static void
2010 : 722 : checkvalidpol(GEN p)
2011 : : {
2012 : 722 : long i,n = lg(p);
2013 [ + + ]: 3694 : for (i=2; i<n; i++)
2014 [ + + ]: 2979 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
2015 : 715 : }
2016 : :
2017 : : static GEN
2018 : 477 : solve_exact_pol(GEN p, long bit)
2019 : : {
2020 : 477 : long i, j, k, m, n = degpol(p), iroots = 0;
2021 : 477 : GEN ex, factors, v = zerovec(n);
2022 : :
2023 : 477 : factors = ZX_squff(Q_primpart(p), &ex);
2024 [ + + ]: 954 : for (i=1; i<lg(factors); i++)
2025 : : {
2026 : 477 : GEN roots_fact = all_roots(gel(factors,i), bit);
2027 : 477 : n = degpol(gel(factors,i));
2028 : 477 : m = ex[i];
2029 [ + + ]: 1915 : for (j=1; j<=n; j++)
2030 [ + + ]: 2876 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
2031 : : }
2032 : 477 : return v;
2033 : : }
2034 : :
2035 : : /* return the roots of p with absolute error bit */
2036 : : static GEN
2037 : 715 : roots_com(GEN q, long bit)
2038 : : {
2039 : : GEN L, p;
2040 : 715 : long v = RgX_valrem_inexact(q, &p);
2041 [ - + ]: 715 : if (lg(p) == 3) L = cgetg(1,t_VEC); /* constant polynomial */
2042 [ + + ]: 715 : else L = isexactpol(p)? solve_exact_pol(p,bit): all_roots(p,bit);
2043 [ + + ]: 715 : if (v)
2044 : : {
2045 : 77 : GEN M, z, t = gel(q,2);
2046 : : long i, x, y, l, n;
2047 : :
2048 [ + + ]: 77 : if (isrationalzero(t)) x = -bit;
2049 : : else
2050 : : {
2051 : 7 : n = gexpo(t);
2052 : 7 : x = n / v; l = degpol(q);
2053 [ + + ]: 35 : for (i = v; i <= l; i++)
2054 : : {
2055 : 28 : t = gel(q,i+2);
2056 [ - + ]: 28 : if (isrationalzero(t)) continue;
2057 : 28 : y = (n - gexpo(t)) / i;
2058 [ + + ]: 28 : if (y < x) x = y;
2059 : : }
2060 : : }
2061 : 77 : z = real_0_bit(x); l = v + lg(L);
2062 : 77 : M = cgetg(l, t_VEC); L -= v;
2063 [ + + ]: 154 : for (i = 1; i <= v; i++) gel(M,i) = z;
2064 [ + + ]: 238 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
2065 : 77 : L = M;
2066 : : }
2067 : 715 : return L;
2068 : : }
2069 : :
2070 : : static GEN
2071 : 51398 : tocomplex(GEN x, long l, long bit)
2072 : : {
2073 : : GEN y;
2074 [ + + ]: 51398 : if (typ(x) == t_COMPLEX)
2075 : : {
2076 [ + + ]: 50324 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
2077 : 5344 : x = gel(x,2);
2078 : 5344 : y = cgetg(3,t_COMPLEX);
2079 : 5344 : gel(y,1) = real_0_bit(-bit);
2080 : 5344 : gel(y,2) = mygprecrc(x, l, -bit);
2081 : : }
2082 : : else
2083 : : {
2084 : 1074 : y = cgetg(3,t_COMPLEX);
2085 : 1074 : gel(y,1) = mygprecrc(x, l, -bit);
2086 : 1074 : gel(y,2) = real_0_bit(-bit);
2087 : : }
2088 : 51398 : return y;
2089 : : }
2090 : :
2091 : : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
2092 : : * up to 2^-e absolute error */
2093 : : static int
2094 : 151385 : cmp_complex_appr(void *E, GEN x, GEN y)
2095 : : {
2096 : 151385 : long e = (long)E;
2097 : : GEN z, xi, yi, xr, yr;
2098 : : long sxi, syi;
2099 [ + + ]: 151385 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
2100 : 58789 : else { xr = x; xi = NULL; sxi = 0; }
2101 [ + + ]: 151385 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
2102 : 52132 : else { yr = y; yi = NULL; syi = 0; }
2103 : : /* Compare absolute values of imaginary parts */
2104 [ + + ]: 151385 : if (!sxi)
2105 : : {
2106 [ + + ][ + - ]: 59996 : if (syi && expo(yi) >= e) return -1;
2107 : : /* |Im x| ~ |Im y| ~ 0 */
2108 : : }
2109 [ + + ]: 91389 : else if (!syi)
2110 : : {
2111 [ + - ][ + - ]: 2545 : if (sxi && expo(xi) >= e) return 1;
2112 : : /* |Im x| ~ |Im y| ~ 0 */
2113 : : }
2114 : : else
2115 : : {
2116 : : long sz;
2117 : 88844 : z = addrr_sign(xi, 1, yi, -1);
2118 : 88844 : sz = signe(z);
2119 [ + + ][ + - ]: 88844 : if (sz && expo(z) >= e) return (int)sz;
2120 : : }
2121 : : /* |Im x| ~ |Im y|, sort according to real parts */
2122 : 88052 : z = subrr(xr, yr);
2123 [ + + ]: 88052 : if (expo(z) >= e) return (int)signe(z);
2124 : : /* Re x ~ Re y. Place negative absolute value before positive */
2125 : 151385 : return (int) (sxi - syi);
2126 : : }
2127 : :
2128 : : static GEN
2129 : 11416 : clean_roots(GEN L, long l, long bit, long clean)
2130 : : {
2131 : 11416 : long i, n = lg(L), ex = 5 - bit;
2132 : 11416 : GEN res = cgetg(n,t_COL);
2133 [ + + ]: 85875 : for (i=1; i<n; i++)
2134 : : {
2135 : 74459 : GEN c = gel(L,i);
2136 [ + + ][ + + ]: 74459 : if (clean && isrealappr(c,ex))
2137 : : {
2138 [ - + ]: 23061 : if (typ(c) == t_COMPLEX) c = gel(c,1);
2139 : 23061 : c = mygprecrc(c, l, -bit);
2140 : : }
2141 : : else
2142 : 51398 : c = tocomplex(c, l, bit);
2143 : 74459 : gel(res,i) = c;
2144 : : }
2145 : 11416 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
2146 : 11416 : return res;
2147 : : }
2148 : :
2149 : : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
2150 : : static GEN
2151 : 743 : roots_aux(GEN p, long l, long clean)
2152 : : {
2153 : 743 : pari_sp av = avma;
2154 : : long bit;
2155 : : GEN L;
2156 : :
2157 [ + + ]: 743 : if (typ(p) != t_POL)
2158 : : {
2159 [ + + ]: 21 : if (gequal0(p)) pari_err_ROOTS0("roots");
2160 [ + + ]: 14 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
2161 : 7 : return cgetg(1,t_COL); /* constant polynomial */
2162 : : }
2163 [ - + ]: 722 : if (!signe(p)) pari_err_ROOTS0("roots");
2164 : 722 : checkvalidpol(p);
2165 [ - + ]: 715 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2166 [ - + ]: 715 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2167 : 715 : bit = prec2nbits(l);
2168 : 715 : L = roots_com(p, bit);
2169 : 722 : return gerepileupto(av, clean_roots(L, l, bit, clean));
2170 : : }
2171 : : GEN
2172 : 596 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
2173 : : /* clean up roots. If root is real replace it by its real part */
2174 : : GEN
2175 : 147 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
2176 : :
2177 : : /* private variant of conjvec. Allow non rational coefficients, shallow
2178 : : * function. */
2179 : : GEN
2180 : 77 : polmod_to_embed(GEN x, long prec)
2181 : : {
2182 : 77 : GEN v, T = gel(x,1), A = gel(x,2);
2183 : : long i, l;
2184 [ + - ][ - + ]: 77 : if (typ(A) != t_POL || varn(A) != varn(T))
2185 : : {
2186 : 0 : checkvalidpol(T);
2187 : 0 : return const_col(degpol(T), A);
2188 : : }
2189 : 77 : v = cleanroots(T,prec); l = lg(v);
2190 [ + + ]: 231 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
2191 : 77 : return v;
2192 : : }
2193 : :
2194 : : GEN
2195 : 10701 : QX_complex_roots(GEN p, long l)
2196 : : {
2197 : 10701 : pari_sp av = avma;
2198 : : long bit;
2199 : : GEN L;
2200 : :
2201 [ - + ]: 10701 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
2202 [ - + ]: 10701 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2203 [ - + ]: 10701 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2204 : 10701 : bit = prec2nbits(l);
2205 : 10701 : L = all_roots(Q_primpart(p), bit);
2206 : 10701 : return gerepileupto(av, clean_roots(L, l, bit, 1));
2207 : : }
2208 : :
2209 : : /********************************************************************/
2210 : : /** **/
2211 : : /** REAL ROOTS OF INTEGER POLYNOMIAL **/
2212 : : /** **/
2213 : : /********************************************************************/
2214 : :
2215 : : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
2216 : : * The inversion is implicit (we take coefficients backwards) */
2217 : : static long
2218 : 54245 : X2XP1(GEN P, long deg, GEN *Premapped)
2219 : : {
2220 : 54245 : pari_sp av = avma;
2221 : 54245 : GEN v = shallowcopy(P);
2222 : : long i, j, vlim, nb, s, s2;
2223 : : char flag;
2224 : :
2225 : 54245 : vlim = deg+2;
2226 : 54245 : nb = 0;
2227 : 54245 : i = 0;
2228 : : do
2229 : : {
2230 [ + + ]: 636145 : for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2231 : 54294 : s = -signe(gel(v, vlim));
2232 : 54294 : vlim--;
2233 : 54294 : i++;
2234 : : }
2235 [ + + ]: 54294 : while (!s);
2236 [ + + ][ + + ]: 54245 : if (vlim != deg + 1 && Premapped) setlg(v, vlim + 2);
2237 : :
2238 [ + + ]: 312501 : for (; i <= deg - 1; i++)
2239 : : {
2240 : 293127 : s2 = -signe(gel(v, 2));
2241 : 293127 : flag = (s2 == s);
2242 [ + + ]: 3376664 : for (j = 2; j < vlim; j++)
2243 : : {
2244 : 3083537 : gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2245 [ + + ]: 3083537 : if (flag) flag = (s2 != signe(gel(v, j+1)));
2246 : : }
2247 [ + + ]: 293127 : if (s == signe(gel(v, vlim)))
2248 : : {
2249 : 52515 : nb++;
2250 [ + + ]: 52515 : if (nb >= 2) { avma = av; return 2; }
2251 : 33793 : s = -s;
2252 : : }
2253 [ - + ]: 274405 : if (gc_needed(av, 3))
2254 : : {
2255 [ # # ]: 0 : if (!Premapped) setlg(v, vlim);
2256 : 0 : v = gerepileupto(av, v);
2257 [ # # ]: 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "X2XP1");
2258 : : }
2259 [ + + ][ + + ]: 274405 : if (flag && !Premapped) goto END;
2260 : 258256 : vlim--;
2261 : : }
2262 [ + + ]: 19374 : if (s == signe(gel(v, vlim))) nb++;
2263 : : END:
2264 [ + + ][ + + ]: 35523 : if (Premapped && nb == 1)
2265 : 5855 : *Premapped = gerepileupto(av, v);
2266 : : else
2267 : 29668 : avma = av;
2268 : 54245 : return nb;
2269 : : }
2270 : :
2271 : : static long
2272 : 0 : _intervalcmp(GEN x, GEN y)
2273 : : {
2274 [ # # ]: 0 : if (typ(x) == t_VEC) x = gel(x, 1);
2275 [ # # ]: 0 : if (typ(y) == t_VEC) y = gel(y, 1);
2276 : 0 : return gcmp(x, y);
2277 : : }
2278 : :
2279 : : static GEN
2280 : 509723 : _gen_nored(void *E, GEN x) { (void)E; return x; }
2281 : : static GEN
2282 : 2086255 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2283 : : static GEN
2284 : 461777 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2285 : : static GEN
2286 : 464334 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2287 : : static GEN
2288 : 595342 : _gen_one(void *E) { (void)E; return gen_1; }
2289 : : static GEN
2290 : 8326 : _gen_zero(void *E) { (void)E; return gen_0; }
2291 : :
2292 : : static struct bb_algebra mp_algebra = { _gen_nored,_mp_add,_mp_mul,_mp_sqr,_gen_one,_gen_zero };
2293 : :
2294 : : static GEN
2295 : 2478227 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
2296 : :
2297 : : /* Split the polynom P in two parts, each with unique sign.
2298 : : * Moreover compute the two parts of the derivative of P. */
2299 : : static long
2300 : 5855 : split_polynoms(GEN P, long deg, long s0, GEN *Pp, GEN *Pm, GEN *Pprimep, GEN *Pprimem)
2301 : : {
2302 : 5855 : long i, degneg, v = evalvarn(varn(P));
2303 [ + + ][ + - ]: 33821 : for(i=1; i <= deg; i++) if (signe(gel(P, i+2)) != s0) break;
2304 : 5855 : degneg = i;
2305 : 5855 : *Pm = cgetg(degneg + 2, t_POL);
2306 : 5855 : (*Pm)[1] = v;
2307 : 5855 : *Pprimem = cgetg(degneg + 1, t_POL);
2308 : 5855 : (*Pprimem)[1] = v;
2309 [ + + ]: 39676 : for(i=0; i < degneg; i++)
2310 : : {
2311 : 33821 : GEN elt = gel(P, i+2);
2312 : 33821 : gel(*Pm, i+2) = elt;
2313 [ + + ]: 33821 : if (i) gel(*Pprimem, i+1) = mului(i, elt);
2314 : : }
2315 : 5855 : *Pp = cgetg(deg - degneg + 3, t_POL);
2316 : 5855 : (*Pp)[1] = v;
2317 : 5855 : *Pprimep = cgetg(deg - degneg + 3, t_POL);
2318 : 5855 : (*Pprimep)[1] = v;
2319 [ + + ]: 37354 : for(i=degneg; i <= deg; i++)
2320 : : {
2321 : 31499 : GEN elt = gel(P, i+2);
2322 : 31499 : gel(*Pp, i+2-degneg) = elt;
2323 : 31499 : gel(*Pprimep, i+2-degneg) = mului(i, elt);
2324 : : }
2325 : 5855 : *Pm = normalizepol_lg(*Pm, degneg+2);
2326 : 5855 : *Pprimem = normalizepol_lg(*Pprimem, degneg+1);
2327 : 5855 : *Pp = normalizepol_lg(*Pp, deg-degneg+3);
2328 : 5855 : *Pprimep = normalizepol_lg(*Pprimep, deg-degneg+3);
2329 : 5855 : return degpol(*Pm) + 1;
2330 : : }
2331 : :
2332 : : static GEN
2333 : 200149 : bkeval_single_power(long d, GEN V)
2334 : : {
2335 : 200149 : long mp = lg(V) - 2;
2336 [ - + ]: 200149 : if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2337 : 200149 : return gel(V, d+1);
2338 : : }
2339 : :
2340 : : static GEN
2341 : 200149 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long deg, long degneg, long bitprec)
2342 : : {
2343 : 200149 : GEN vp = gen_bkeval_powers(Pp, deg-degneg, pows, NULL, &mp_algebra, _mp_cmul);
2344 : 200149 : GEN vm = gen_bkeval_powers(Pm, degneg-1, pows, NULL, &mp_algebra, _mp_cmul);
2345 : 200149 : GEN xa = bkeval_single_power(degneg, pows);
2346 : : GEN r;
2347 : 200149 : vp = gmul(vp, xa);
2348 : 200149 : r = gadd(vp, vm);
2349 [ + + ][ + + ]: 200149 : if (expo(vp) - (signe(r) ? expo(r) : 0) > prec2nbits(realprec(vp)) - bitprec)
2350 : 14572 : return NULL;
2351 : 200149 : return r;
2352 : : }
2353 : :
2354 : : /* Newton for polynom P. One solution between 0 and infinity, P' has also at
2355 : : * most one zero. P is almost certainly without zero coefficients. */
2356 : : static GEN
2357 : 5855 : polsolve(GEN P, long bitprec)
2358 : : {
2359 : 5855 : pari_sp av = avma;
2360 : : GEN Pp, Pm, Pprimep, Pprimem, pows;
2361 : : GEN Pprime, Pprime2;
2362 : 5855 : GEN ra, rb, rc, rcold = NULL, Pc, Ppc;
2363 : 5855 : long deg = degpol(P), degneg, rt;
2364 : 5855 : long s0, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
2365 : : long bitaddprec, addprec;
2366 : 5855 : long expoold = LONG_MAX, iter;
2367 [ - + ]: 5855 : if (deg == 1)
2368 : : {
2369 : 0 : ra = negr(divrr(itor(gel(P, 2), bitprec), itor(gel(P, 3), bitprec)));
2370 : 0 : return gerepileuptoleaf(av, ra);
2371 : : }
2372 : 5855 : Pprime = ZX_deriv(P); Pprime2 = ZX_deriv(Pprime);
2373 : 5855 : bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
2374 : 5855 : s0 = signe(gel(P, 2));
2375 : 5855 : degneg = split_polynoms(P, deg, s0, &Pp, &Pm, &Pprimep, &Pprimem);
2376 : 5855 : rt = maxss(degneg, brent_kung_optpow(maxss(deg-degneg, degneg-1), 2, 1));
2377 : : {
2378 : : /* Optimized Cauchy bound */
2379 : 5855 : GEN summin = gen_0, absmaj;
2380 : : long i;
2381 : :
2382 : 5855 : absmaj = gel(P, 2);
2383 [ + + ]: 33821 : for(i=1; i < degneg; i++)
2384 : : {
2385 : 27966 : GEN elt = gel(P, i+2);
2386 [ + + ]: 27966 : if (absi_cmp(absmaj, elt) < 0) absmaj = elt;
2387 : : }
2388 : 5855 : summin = gel(P, i+2);
2389 [ + + ]: 31499 : while (++i <= deg) summin = addii(summin, gel(P, i+2));
2390 : 5855 : rb = subsr(1, rdivii(absmaj, summin, cprec));
2391 : : do
2392 : : {
2393 : 5855 : pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2394 : 5855 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec);
2395 [ - + ]: 5855 : if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
2396 [ + - ]: 5855 : if (signe(Pc) != s0) break;
2397 : 0 : shiftr_inplace(rb,1);
2398 : : }
2399 : 0 : while (1);
2400 : : }
2401 : 5855 : ra = NULL;
2402 : 5855 : iter = 0;
2403 : : while (1)
2404 : : {
2405 : : GEN wdth;
2406 : 55576 : iter++;
2407 [ + + ]: 55576 : if (ra)
2408 : 36736 : rc = shiftr(addrr(ra, rb), -1);
2409 : : else
2410 : 18840 : rc = shiftr(rb, -1);
2411 : : do
2412 : : {
2413 : 55576 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2414 : 55576 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+2);
2415 [ + - ]: 55576 : if (Pc) break;
2416 : 0 : cprec++; rc = rtor(rc, cprec);
2417 : 0 : } while (1);
2418 [ + + ]: 55576 : if (signe(Pc) == s0)
2419 : 23773 : ra = rc;
2420 : : else
2421 : 31803 : rb = rc;
2422 [ + + ]: 55576 : if (!ra) continue;
2423 : 42591 : wdth = subrr(rb, ra);
2424 [ + + ]: 42591 : if (!(iter % 8))
2425 : : {
2426 : 6345 : GEN m1 = poleval(Pprime, ra), M2;
2427 [ + + ]: 6345 : if (signe(m1) == s0) continue;
2428 : 6156 : M2 = poleval(Pprime2, rb);
2429 [ + + ]: 6156 : if (absr_cmp(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2430 : 5855 : break;
2431 : : }
2432 [ - + ]: 36246 : else if (gexpo(wdth) <= -bitprec)
2433 : 0 : break;
2434 : 49721 : }
2435 : 5855 : rc = rb;
2436 : 5855 : iter = 0;
2437 : : while (1)
2438 : : {
2439 : : long exponew;
2440 : : GEN dist;
2441 : 69359 : iter++;
2442 : 69359 : rcold = rc;
2443 : 69359 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2444 : 69359 : Ppc = splitpoleval(Pprimep, Pprimem, pows, deg-1, degneg-1, bitaddprec+4);
2445 [ + - ]: 69359 : if (Ppc)
2446 : 69359 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+4);
2447 [ + - ][ + + ]: 69359 : if (!Ppc || !Pc)
2448 : : {
2449 [ + + ]: 14572 : if (cprec >= prec+addprec)
2450 : 1962 : cprec++;
2451 : : else
2452 : 12610 : cprec = minss(2*cprec, prec+addprec);
2453 : 14572 : rc = rtor(rc, cprec); /* Backtrack one step */
2454 : 14572 : continue;
2455 : : }
2456 : 54787 : dist = divrr(Pc, Ppc);
2457 : 54787 : rc = subrr(rc, dist);
2458 [ + - ][ - + ]: 54787 : if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2459 : : {
2460 [ # # ]: 0 : if (cprec >= prec+addprec) break;
2461 : 0 : cprec = minss(2*cprec, prec+addprec);
2462 : 0 : rc = rtor(rcold, cprec); /* Backtrack one step */
2463 : 0 : continue;
2464 : : }
2465 [ + + ]: 54787 : if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2466 : 46910 : exponew = expo(dist);
2467 [ + + ]: 46910 : if (exponew < -bitprec - 1)
2468 : : {
2469 [ + + ]: 6679 : if (cprec >= prec+addprec) break;
2470 : 824 : cprec = minss(2*cprec, prec+addprec);
2471 : 824 : rc = rtor(rc, cprec);
2472 : 824 : continue;
2473 : : }
2474 [ + + ]: 40231 : if (exponew > expoold - 2)
2475 : : {
2476 [ - + ]: 2022 : if (cprec >= prec+addprec) break;
2477 : 2022 : cprec = minss(2*cprec, prec+addprec);
2478 : 2022 : rc = rtor(rc, cprec);
2479 : 2022 : expoold = LONG_MAX;
2480 : 2022 : continue;
2481 : : }
2482 : 38209 : expoold = exponew;
2483 : 63504 : }
2484 : 5855 : return gerepileuptoleaf(av, rtor(rc, prec));
2485 : : }
2486 : :
2487 : : static GEN
2488 : 6781 : usp(GEN Q0, long deg, long *nb_donep, long flag, long bitprec)
2489 : : {
2490 : : pari_sp av;
2491 : : GEN Q, sol;
2492 : 6781 : long nb_todo, nbr = 0, ind, deg0, indf, i, k, nb, j;
2493 : 6781 : long listsize = 64, nb_done = 0;
2494 : : GEN c, Lc, Lk;
2495 : :
2496 : 6781 : av = avma;
2497 : :
2498 : 6781 : sol = const_col(deg, gen_0);
2499 : 6781 : deg0 = deg;
2500 : 6781 : Lc = const_vec(listsize, gen_0);
2501 : 6781 : Lk = cgetg(listsize+1, t_VECSMALL);
2502 : 6781 : c = gen_0;
2503 : 6781 : k = Lk[1] = 0;
2504 : 6781 : ind = 1; indf = 2;
2505 : 6781 : Q = gcopy(Q0);
2506 : :
2507 : 6781 : nb_todo = 1;
2508 [ + + ]: 61026 : while (nb_todo)
2509 : : {
2510 : 54245 : GEN nc = gel(Lc, ind), Qremapped;
2511 [ + + ]: 54245 : if (Lk[ind] == k + 1)
2512 : : {
2513 : 21324 : deg0 = deg;
2514 : 21324 : setlg(Q0, deg + 3);
2515 : 21324 : Q0 = ZX_rescale2n(Q0, 1);
2516 : 21324 : Q = Q_primpart(Q0);
2517 : 21324 : c = gen_0;
2518 : : }
2519 : :
2520 [ + + ]: 54245 : if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
2521 : :
2522 : 54245 : k = Lk[ind];
2523 : 54245 : c = nc;
2524 : 54245 : ind++;
2525 : 54245 : nb_todo--;
2526 : :
2527 [ + + ]: 54245 : if (equalii(gel(Q, 2), gen_0))
2528 : : {
2529 : 126 : GEN newsol = gmul2n(c, -k);
2530 [ + + ]: 168 : for (j = 1; j <= nbr; j++)
2531 [ + + ]: 77 : if (gequal(gel(sol, j), newsol)) break;
2532 [ + + ]: 126 : if (j > nbr) gel(sol, ++nbr) = newsol;
2533 : :
2534 : 126 : deg0--;
2535 [ + + ]: 504 : for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
2536 : 126 : setlg(Q, j);
2537 : : }
2538 : :
2539 [ + + ]: 54245 : nb = X2XP1(Q, deg0, flag == 1 ? &Qremapped : NULL);
2540 : 54245 : nb_done++;
2541 : :
2542 [ + + + ]: 54245 : switch (nb)
2543 : : {
2544 : : case 0:
2545 : 18708 : break;
2546 : : case 1:
2547 [ - + + ]: 11805 : switch(flag)
2548 : : {
2549 : : case 0:
2550 : : {
2551 : : GEN low, hi;
2552 : 0 : low = gmul2n(c, -k);
2553 : 0 : hi = gmul2n(addiu(c,1), -k);
2554 : 0 : gel(sol, ++nbr) = mkvec2(low, hi);
2555 : : }
2556 : 0 : break;
2557 : : case 1:
2558 : : { /* Caveat emptor: Qremapped is the reciprocal polynomial */
2559 : 5855 : GEN sr = polsolve(Qremapped, bitprec+1);
2560 : 5855 : sr = divrr(sr, addsr(1, sr));
2561 : 5855 : sr = gmul2n(addir(c, sr), -k);
2562 : 5855 : gel(sol, ++nbr) = rtor(sr, nbits2prec(bitprec));
2563 : : }
2564 : 5855 : break;
2565 : : default:
2566 : 5950 : gel(sol, ++nbr) = gen_0;
2567 : : }
2568 : 11805 : break;
2569 : :
2570 : : default:
2571 [ + + ]: 23732 : if (indf + 2 > listsize)
2572 : : {
2573 [ + - ]: 126 : if (ind>1)
2574 : : {
2575 [ + + ]: 280 : for (i = ind; i < indf; i++)
2576 : : {
2577 : 154 : gel(Lc, i-ind+1) = gel(Lc, i);
2578 : 154 : Lk[i-ind+1] = Lk[i];
2579 : : }
2580 : 126 : indf -= ind-1; ind = 1;
2581 : : }
2582 [ - + ]: 126 : if (indf + 2 > listsize)
2583 : : {
2584 : 0 : listsize *= 2;
2585 : 0 : Lc = vec_lengthen(Lc, listsize);
2586 : 0 : Lk = vecsmall_lengthen(Lk, listsize);
2587 : : }
2588 [ + + ]: 8036 : for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2589 : : }
2590 : 23732 : nc = shifti(c, 1);
2591 : 23732 : gel(Lc, indf) = nc;
2592 : 23732 : gel(Lc, indf + 1) = addis(nc, 1);
2593 : 23732 : Lk[indf] = Lk[indf + 1] = k + 1;
2594 : 23732 : indf += 2;
2595 : 23732 : nb_todo += 2;
2596 : : }
2597 : :
2598 [ - + ]: 54245 : if (gc_needed(av, 2))
2599 : : {
2600 : 0 : gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
2601 [ # # ]: 54245 : if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_uspensky", avma);
2602 : : }
2603 : : }
2604 : :
2605 : 6781 : setlg(sol, nbr+1);
2606 : 6781 : *nb_donep += nb_done;
2607 : 6781 : return gerepilecopy(av, sol);
2608 : : }
2609 : :
2610 : : static GEN
2611 : 1554 : ZX_uspensky_cst_pol(long nbz, long flag, long bitprec)
2612 : : {
2613 [ - - + ]: 1554 : switch(flag)
2614 : : {
2615 : 0 : case 0: return zerocol(nbz);
2616 [ # # ]: 0 : case 1: retconst_col(nbz, real_0_bit(bitprec));
2617 : 1554 : default: return utoi(nbz);
2618 : : }
2619 : : }
2620 : :
2621 : : GEN
2622 : 8450 : ZX_uspensky(GEN P, GEN ab, long flag, long bitprec)
2623 : : {
2624 : 8450 : pari_sp av = avma;
2625 : 8450 : GEN a, b, sol = NULL, Pcur;
2626 : : double fb;
2627 : 8450 : long nbz, deg, nb_done = 0;
2628 : :
2629 : 8450 : deg = degpol(P);
2630 [ + + ][ + - ]: 8450 : if (deg == 0) return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2631 [ + + ]: 7883 : if (ab)
2632 : : {
2633 [ + + ]: 3250 : if (typ(ab) == t_VEC)
2634 : : {
2635 [ - + ]: 119 : if (lg(ab) != 3) pari_err_DIM("ZX_uspensky");
2636 : 119 : a = gel(ab, 1);
2637 : 119 : b = gel(ab, 2);
2638 : : }
2639 : : else
2640 : : {
2641 : 3131 : a = ab;
2642 : 3131 : b = mkoo();
2643 : : }
2644 : : }
2645 : : else
2646 : : {
2647 : 4633 : a = mkmoo();
2648 : 4633 : b = mkoo();
2649 : : }
2650 [ + + + ]: 7883 : switch (gcmp(a, b))
2651 : : {
2652 [ - + ]: 7 : case 1: avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2653 : : case 0:
2654 [ + - ][ + + ]: 21 : if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
2655 [ + + ]: 14 : { avma = av; return flag <= 1 ? mkcolcopy(a): gen_1; }
2656 : : else
2657 [ + - ]: 7 : { avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
2658 : : }
2659 : 7855 : nbz = ZX_valrem(P, &Pcur);
2660 : 7855 : deg -= nbz;
2661 [ + + ]: 7855 : if (!nbz) Pcur = P;
2662 [ + + ][ + - ]: 7855 : if (nbz && (gsigne(a) > 0 || gsigne(b) < 0)) nbz = 0;
[ - + ]
2663 [ + + ]: 7855 : if (deg == 0) { avma = av; return ZX_uspensky_cst_pol(nbz, flag, bitprec); }
2664 [ + + ]: 7288 : if (deg == 1)
2665 : : {
2666 : 2739 : sol = gdiv(gneg(gel(Pcur, 2)), pollead(Pcur, -1));
2667 [ + + ][ - + ]: 2739 : if (gcmp(a, sol) > 0 || gcmp(sol, b) > 0)
2668 : : {
2669 : 987 : avma = av;
2670 : 987 : return ZX_uspensky_cst_pol(nbz, flag, bitprec);
2671 : : }
2672 [ + + ]: 1752 : if (flag >= 2) { avma = av; return utoi(nbz+1); }
2673 : 737 : sol = concat(zerocol(nbz), mkcol(sol));
2674 [ + - ]: 737 : if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
2675 : 737 : return gerepilecopy(av, sol);
2676 : : }
2677 [ - + + ]: 4549 : switch(flag)
2678 : : {
2679 : : case 0:
2680 : 0 : sol = zerocol(nbz);
2681 : 0 : break;
2682 : : case 1:
2683 : 1322 : sol = const_col(nbz, real_0_bit(bitprec));
2684 : 1322 : break;
2685 : : /* case 2: nothing */
2686 : : }
2687 : :
2688 [ + + ][ + + ]: 4549 : if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
[ + - ]
2689 : : {
2690 : 28 : fb = fujiwara_bound_real(Pcur, -1);
2691 [ + + ]: 28 : if (fb > -pariINFINITY)
2692 : 7 : a = negi(int2n((long)ceil(fb)));
2693 : : else
2694 : 21 : a = gen_0;
2695 : : }
2696 [ + + ][ + + ]: 4549 : if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
[ + + ]
2697 : : {
2698 : 14 : fb = fujiwara_bound_real(Pcur, 1);
2699 [ + - ]: 14 : if (fb > -pariINFINITY)
2700 : 14 : b = int2n((long)ceil(fb));
2701 : : else
2702 : 0 : b = gen_0;
2703 : : }
2704 : :
2705 [ + + ][ + + ]: 4549 : if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
2706 : : {
2707 : : pari_sp av1;
2708 : 91 : GEN den = lcmii(denom(a), denom(b)), diff, unscaledres, co, Pdiv;
2709 : : GEN ascaled;
2710 : : long i;
2711 [ + + ]: 91 : if (!is_pm1(den))
2712 : : {
2713 : 28 : Pcur = ZX_rescale(Pcur, den);
2714 : 28 : ascaled = gmul(a, den);
2715 : : }
2716 : : else
2717 : : {
2718 : 63 : den = NULL;
2719 : 63 : ascaled = a;
2720 : : }
2721 [ + + ]: 91 : diff = subii(den ? gmul(b,den) : b, ascaled);
2722 : 91 : Pcur = ZX_unscale(ZX_translate(Pcur, ascaled), diff);
2723 : 91 : av1 = avma;
2724 : 91 : Pdiv = cgetg(deg+2, t_POL);
2725 : 91 : Pdiv[1] = Pcur[1];
2726 : 91 : co = gel(Pcur, deg+2);
2727 [ + + ]: 364 : for (i = deg; --i >= 0; )
2728 : : {
2729 : 273 : gel(Pdiv, i+2) = co;
2730 : 273 : co = addii(co, gel(Pcur, i+2));
2731 : : }
2732 [ + + ]: 91 : if (!signe(co))
2733 : : {
2734 : 49 : Pcur = Pdiv;
2735 : 49 : deg--;
2736 [ + + ]: 49 : if (flag <= 1)
2737 : 14 : sol = concat(sol, b);
2738 : : else
2739 : 35 : nbz++;
2740 : : }
2741 : : else
2742 : 42 : avma = av1;
2743 : 91 : unscaledres = usp(Pcur, deg, &nb_done, flag, bitprec);
2744 [ + + ]: 91 : if (flag <= 1)
2745 : : {
2746 [ + + ]: 63 : for (i = 1; i < lg(unscaledres); i++)
2747 : : {
2748 : 35 : GEN z = gmul(diff, gel(unscaledres, i));
2749 [ - + ]: 35 : if (typ(z) == t_VEC)
2750 : : {
2751 : 0 : gel(z, 1) = gadd(ascaled, gel(z, 1));
2752 : 0 : gel(z, 2) = gadd(ascaled, gel(z, 2));
2753 : : }
2754 : : else
2755 : 35 : z = gadd(ascaled, z);
2756 [ + + ]: 35 : if (den) z = gdiv(z, den);
2757 : 35 : gel(unscaledres, i) = z;
2758 : : }
2759 : 28 : sol = concat(sol, unscaledres);
2760 : : }
2761 : : else
2762 : 63 : nbz += lg(unscaledres) - 1;
2763 : : }
2764 [ + + ][ + + ]: 4549 : if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(Pcur, 1)) > -pariINFINITY)
2765 : : {
2766 : : GEN Pcurp, unscaledres;
2767 : 3604 : long bp = (long)ceil(fb);
2768 [ + + ]: 3604 : if (bp < 0) bp = 0;
2769 : 3604 : Pcurp = ZX_unscale2n(Pcur, bp);
2770 : 3604 : unscaledres = usp(Pcurp, deg, &nb_done, flag, bitprec);
2771 [ + + ]: 3604 : if (flag <= 1)
2772 : 1287 : sol = concat(sol, gmul2n(unscaledres, bp));
2773 : : else
2774 : 2317 : nbz += lg(unscaledres) - 1;
2775 : : }
2776 [ + + ][ + + ]: 4549 : if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(Pcur,-1)) > -pariINFINITY)
2777 : : {
2778 : : GEN Pcurm, unscaledres;
2779 : 3086 : long i, bm = (long)ceil(fb);
2780 [ + + ]: 3086 : if (bm < 0) bm = 0;
2781 : 3086 : Pcurm = ZX_unscale(Pcur, gen_m1);
2782 : 3086 : Pcurm = ZX_unscale2n(Pcurm, bm);
2783 : 3086 : unscaledres = usp(Pcurm,deg,&nb_done,flag,bitprec);
2784 [ + + ]: 3086 : if (flag <= 1)
2785 : : {
2786 [ + + ]: 3168 : for (i = 1; i < lg(unscaledres); i++)
2787 : : {
2788 : 2133 : GEN z = gneg(gmul2n(gel(unscaledres, i), bm));
2789 [ - + ]: 2133 : if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
2790 : 2133 : gel(unscaledres, i) = z;
2791 : : }
2792 : 1035 : sol = concat(unscaledres, sol);
2793 : : }
2794 : : else
2795 : 2051 : nbz += lg(unscaledres) - 1;
2796 : : }
2797 : :
2798 [ - + ]: 4549 : if (DEBUGLEVEL > 4)
2799 : 0 : err_printf("ZX_uspensky: Number of visited nodes: %d\n", nb_done);
2800 : :
2801 [ + + ]: 4549 : if (flag >= 2) return utoi(nbz);
2802 [ + - ]: 1322 : if (flag)
2803 : 1322 : sol = sort(sol);
2804 : : else
2805 : 0 : sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
2806 : 8450 : return gerepileupto(av, sol);
2807 : : }
2808 : :
2809 : : /* x a scalar */
2810 : : static GEN
2811 : 35 : rootsdeg0(GEN x)
2812 : : {
2813 [ + + ]: 35 : if (!is_rational_t(typ(x))) pari_err_TYPE("realroots",x);
2814 [ + + ]: 28 : if (gequal0(x)) pari_err_ROOTS0("realroots");
2815 : 14 : return cgetg(1,t_COL); /* constant polynomial */
2816 : : }
2817 : : static void
2818 : 1190 : checkbound(GEN a)
2819 : : {
2820 [ + - ]: 1190 : switch(typ(a))
2821 : : {
2822 : 1190 : case t_INT: case t_FRAC: case t_INFINITY: break;
2823 : 0 : default: pari_err_TYPE("polrealroots", a);
2824 : : }
2825 : 1190 : }
2826 : : static GEN
2827 : 3214 : check_ab(GEN ab)
2828 : : {
2829 : : GEN a, b;
2830 [ + + ]: 3214 : if (!ab) return NULL;
2831 [ + - ][ - + ]: 595 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
2832 : 595 : a = gel(ab,1); checkbound(a);
2833 : 595 : b = gel(ab,2); checkbound(b);
2834 [ + + ][ + - ]: 595 : if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
[ + + ]
2835 [ + - ]: 476 : typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
2836 : 3214 : return ab;
2837 : : }
2838 : : GEN
2839 : 2668 : realroots(GEN P, GEN ab, long prec)
2840 : : {
2841 : 2668 : pari_sp av = avma;
2842 : 2668 : long nrr = 0;
2843 : 2668 : GEN sol = NULL, fa, ex;
2844 : : long i, j, k;
2845 : :
2846 : 2668 : ab = check_ab(ab);
2847 [ + + ]: 2668 : if (typ(P) != t_POL) return rootsdeg0(P);
2848 [ + + + ]: 2647 : switch(degpol(P))
2849 : : {
2850 : 7 : case -1: return rootsdeg0(gen_0);
2851 : 7 : case 0: return rootsdeg0(gel(P,2));
2852 : : }
2853 : 2633 : P = Q_primpart(P);
2854 [ - + ]: 2633 : if (!RgX_is_ZX(P)) pari_err_TYPE("realroots",P);
2855 : 2633 : fa = ZX_squff(P, &ex);
2856 [ + + ]: 5273 : for (i = 1; i < lg(fa); i++)
2857 : : {
2858 : 2640 : GEN Pi = gel(fa, i), soli, soli2 = NULL;
2859 : 2640 : long n, nrri = 0, h, nbz;
2860 [ + + ]: 2640 : if (ab)
2861 : 42 : h = 1;
2862 : : else
2863 : 2598 : Pi = RgX_deflate_max(Pi, &h);
2864 [ + + ]: 2640 : if (!signe(gel(Pi, 2)))
2865 : : {
2866 : 567 : Pi = RgX_shift_shallow(Pi, -1);
2867 : 567 : nbz = 1;
2868 : : }
2869 : : else
2870 : 2073 : nbz = 0;
2871 [ + + ]: 2640 : soli = ZX_uspensky(Pi, h%2 ? ab: gen_0, 1, prec2nbits(prec));
2872 : 2640 : n = lg(soli);
2873 [ + + ]: 2640 : if (!(h % 2)) soli2 = cgetg(n, t_COL);
2874 [ + + ]: 9295 : for (j = 1; j < n; j++)
2875 : : {
2876 : 6655 : GEN elt = gel(soli, j);
2877 [ + + ]: 6655 : if (typ(elt) != t_REAL)
2878 : : {
2879 : 63 : nrri++;
2880 : 63 : elt = gtofp(elt, prec);
2881 : 63 : gel(soli, j) = elt;
2882 : : }
2883 [ + + ]: 6655 : if (h > 1)
2884 : : { /* note: elt != 0 because we are square free */
2885 : : GEN r;
2886 [ + + ]: 1248 : if (h == 2)
2887 : 1234 : r = sqrtr(elt);
2888 : : else
2889 : : {
2890 [ + + ]: 14 : if (signe(elt) < 0)
2891 : 7 : r = negr(sqrtnr(negr(elt), h));
2892 : : else
2893 : 7 : r = sqrtnr(elt, h);
2894 : : }
2895 : 1248 : gel(soli, j) = r;
2896 [ + + ]: 1248 : if (!(h % 2)) gel(soli2, j) = negr(r);
2897 : : }
2898 : : }
2899 [ + + ]: 2640 : if (!(h % 2)) soli = shallowconcat(soli, soli2);
2900 [ + + ]: 2640 : if (nbz) soli = shallowconcat(soli, real_0(prec));
2901 [ + + ]: 5280 : for (k = 1; k <= ex[i]; k++)
2902 [ + + ]: 2640 : sol = sol ? shallowconcat(sol, soli) : soli;
2903 : 2640 : nrr += ex[i]*nrri;
2904 : : }
2905 : :
2906 [ - + ]: 2633 : if (DEBUGLEVEL > 4)
2907 : : {
2908 : 0 : err_printf("Number of real roots: %d\n", lg(sol)-1);
2909 : 0 : err_printf(" -- of which 2-integral: %ld\n", nrr);
2910 : : }
2911 : :
2912 : 2647 : return gerepileupto(av, sort(sol));
2913 : : }
2914 : :
2915 : : /* P non-constant, squarefree ZX */
2916 : : long
2917 : 5733 : ZX_sturm(GEN P)
2918 : : {
2919 : 5733 : pari_sp av = avma;
2920 : : long h, r;
2921 : 5733 : P = RgX_deflate_max(P, &h);
2922 [ + + ]: 5733 : if (odd(h))
2923 : 3353 : r = itos(ZX_uspensky(P, NULL, 2, 0));
2924 : : else
2925 : 2380 : r = 2*itos(ZX_uspensky(P, gen_0, 2, 0));
2926 : 5733 : avma = av; return r;
2927 : : }
2928 : : /* P non-constant, squarefree ZX */
2929 : : long
2930 : 546 : ZX_sturmpart(GEN P, GEN ab)
2931 : : {
2932 : 546 : pari_sp av = avma;
2933 : : long r;
2934 [ + + ]: 546 : if (!check_ab(ab)) return ZX_sturm(P);
2935 : 77 : r = itos(ZX_uspensky(P, ab, 2, 0));
2936 : 546 : avma = av; return r;
2937 : : }
|
