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 : 6441740 : addCC(GEN x, GEN y)
35 : : {
36 : : GEN z;
37 : :
38 [ + + ]: 6441740 : if (typ(x) == t_INT)
39 : : {
40 [ + + ]: 5778706 : if (typ(y) == t_INT) return addii(x,y);
41 : : /* ty == t_COMPLEX */
42 : 576 : z = cgetg(3,t_COMPLEX);
43 : 576 : gel(z,1) = addii(x, gel(y,1));
44 : 576 : gel(z,2) = icopy(gel(y,2)); return z;
45 : : }
46 : : /* tx == t_COMPLEX */
47 : 663034 : z = cgetg(3,t_COMPLEX);
48 [ + + ]: 663034 : if (typ(y) == t_INT)
49 : : {
50 : 9361 : gel(z,1) = addii(gel(x,1),y);
51 : 9361 : gel(z,2) = icopy(gel(x,2)); return z;
52 : : }
53 : : /* ty == t_COMPLEX */
54 : 653673 : gel(z,1) = addii(gel(x,1),gel(y,1));
55 : 6441740 : 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 : 12234516 : mulCC(GEN x, GEN y)
60 : : {
61 : : GEN z;
62 : :
63 [ + + ]: 12234516 : if (typ(x) == t_INT)
64 : : {
65 [ + + ]: 10614547 : if (typ(y) == t_INT) return mulii(x,y);
66 : : /* ty == t_COMPLEX */
67 : 2121 : z = cgetg(3,t_COMPLEX);
68 : 2121 : gel(z,1) = mulii(x, gel(y,1));
69 : 2121 : gel(z,2) = mulii(x, gel(y,2)); return z;
70 : : }
71 : : /* tx == t_COMPLEX */
72 : 1619969 : z = cgetg(3,t_COMPLEX);
73 [ + + ]: 1619969 : if (typ(y) == t_INT)
74 : : {
75 : 474298 : gel(z,1) = mulii(gel(x,1),y);
76 : 474298 : gel(z,2) = mulii(gel(x,2),y); return z;
77 : : }
78 : : /* ty == t_COMPLEX */
79 : : {
80 : 1145671 : pari_sp av = avma, tetpil;
81 : : GEN p1, p2;
82 : :
83 : 1145671 : p1 = mulii(gel(x,1),gel(y,1));
84 : 1145671 : p2 = mulii(gel(x,2),gel(y,2));
85 : 1145671 : y = mulii(addii(gel(x,1),gel(x,2)),
86 : 2291342 : addii(gel(y,1),gel(y,2)));
87 : 1145671 : x = addii(p1,p2); tetpil = avma;
88 : 1145671 : gel(z,1) = subii(p1,p2);
89 : 1145671 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
90 : 12234516 : return z;
91 : : }
92 : : }
93 : : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
94 : : static GEN
95 : 10227143 : sqrCC(GEN x)
96 : : {
97 : : GEN z;
98 : :
99 [ + + ]: 10227143 : if (typ(x) == t_INT) return sqri(x);
100 : : /* tx == t_COMPLEX */
101 : 1591951 : z = cgetg(3,t_COMPLEX);
102 : : {
103 : 1591951 : pari_sp av = avma, tetpil;
104 : : GEN y, p1, p2;
105 : :
106 : 1591951 : p1 = sqri(gel(x,1));
107 : 1591951 : p2 = sqri(gel(x,2));
108 : 1591951 : y = sqri(addii(gel(x,1),gel(x,2)));
109 : 1591951 : x = addii(p1,p2); tetpil = avma;
110 : 1591951 : gel(z,1) = subii(p1,p2);
111 : 1591951 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
112 : 10227143 : return z;
113 : : }
114 : : }
115 : :
116 : : static void
117 : 1726885 : set_karasquare_limit(long bit)
118 : : {
119 [ + + ]: 1726885 : if (bit<600) { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
120 [ + - ]: 351 : 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 : 1726885 : }
125 : :
126 : : /* assume lP > 0, lP = lgpol(P) */
127 : : static GEN
128 : 3571233 : CX_square_spec(GEN P, long lP)
129 : : {
130 : : GEN s, t;
131 : 3571233 : long i, j, l, nn, n = lP - 1;
132 : : pari_sp av;
133 : :
134 : 3571233 : nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
135 : 3571233 : gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
136 [ + + ]: 9364009 : for (i=1; i<=n; i++)
137 : : {
138 : 5792776 : av = avma; l = (i+1)>>1;
139 : 5792776 : t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
140 [ + + ]: 8004612 : for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
141 : 5792776 : t = gmul2n(t,1);
142 [ + + ]: 5792776 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
143 : 5792776 : gel(s,i+2) = gerepileupto(av, t);
144 : : }
145 : 3571233 : gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
146 [ + + ]: 6655910 : for ( ; i<nn; i++)
147 : : {
148 : 3084677 : av = avma; l = (i+1)>>1;
149 : 3084677 : t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
150 [ + + ]: 4229904 : for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
151 : 3084677 : t = gmul2n(t,1);
152 [ + + ]: 3084677 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
153 : 3084677 : gel(s,i+2) = gerepileupto(av, t);
154 : : }
155 : 3571233 : return normalizepol_lg(s, nn+3);
156 : : }
157 : : /* not stack clean */
158 : : static GEN
159 : 58904 : RgX_addspec(GEN x, long nx, GEN y, long ny)
160 : : {
161 : : GEN z, t;
162 : : long i;
163 [ + + ]: 58904 : if (nx == ny) {
164 : 33007 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
165 [ + + ]: 244321 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
166 : 33007 : return normalizepol_lg(z, nx+2);
167 : : }
168 [ + - ]: 25897 : if (ny < nx) {
169 : 25897 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
170 [ + + ]: 182407 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
171 [ + + ]: 51794 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
172 : 25897 : 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 : 58904 : 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 : 3630482 : karasquare(GEN P, long nP)
206 : : {
207 : : GEN Q, s0, s1, s2, a, t;
208 : 3630482 : long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
209 : : pari_sp av;
210 : :
211 [ + + ][ + + ]: 3630482 : if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
212 : 58904 : av = avma;
213 : 58904 : n0 = (n>>1) + 1; n1 = nP - n0;
214 : 58904 : s0 = karasquare(P, n0); Q = P + n0;
215 : 58904 : s2 = karasquare(Q, n1);
216 : 58904 : s1 = RgX_addspec(P, n0, Q, n1);
217 : 58904 : s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
218 : 58904 : N = (n<<1) + 1;
219 : 58904 : a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
220 : 58904 : t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
221 [ + + ]: 772732 : for (i=0; i < l; i++) gel(t,i) = gel(s0,i);
222 [ + + ]: 132518 : for ( ; i < N0; i++) gel(t,i) = gen_0;
223 : 58904 : t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
224 [ + + ]: 690203 : for (i=0; i < l; i++) gel(t,i) = gel(s2,i);
225 [ + + ]: 104349 : for ( ; i < N1; i++) gel(t,i) = gen_0;
226 : 58904 : t = a+2 + n0; l = lgpol(s1); s1 += 2;
227 [ + + ]: 728460 : for (i=0; i < l; i++) gel(t,i) = gadd(gel(t,i), gel(s1,i));
228 : 3630482 : return gerepilecopy(av, normalizepol_lg(a, N+2));
229 : : }
230 : : /* spec function. nP = lgpol(P) */
231 : : static GEN
232 : 3453770 : cook_square(GEN P, long nP)
233 : : {
234 : : GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
235 : 3453770 : long n0, n3, i, j, n = nP - 1;
236 : : pari_sp av;
237 : :
238 [ + - ][ + - ]: 3453770 : 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 : 3453770 : return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
303 : : }
304 : :
305 : : static GEN
306 : 1726885 : graeffe(GEN p)
307 : : {
308 : : GEN p0, p1, s0, s1;
309 : 1726885 : long n = degpol(p), n0, n1, i;
310 : :
311 [ - + ]: 1726885 : if (!n) return gcopy(p);
312 : 1726885 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
313 : 1726885 : p0 = new_chunk(n0);
314 : 1726885 : p1 = new_chunk(n1);
315 [ + + ]: 5790365 : for (i=0; i<n1; i++)
316 : : {
317 : 4063480 : p0[i] = p[2+(i<<1)];
318 : 4063480 : p1[i] = p[3+(i<<1)];
319 : : }
320 [ + + ]: 1726885 : if (n1 != n0)
321 : 847678 : p0[i] = p[2+(i<<1)];
322 : 1726885 : s0 = cook_square(p0, n0);
323 : 1726885 : s1 = cook_square(p1, n1);
324 : 1726885 : return RgX_sub(s0, RgX_shift_shallow(s1,1));
325 : : }
326 : : GEN
327 : 3085 : ZX_graeffe(GEN p)
328 : : {
329 : 3085 : pari_sp av = avma;
330 : : GEN p0, p1, s0, s1;
331 : 3085 : long n = degpol(p);
332 : :
333 [ - + ]: 3085 : if (!n) return ZX_copy(p);
334 : 3085 : RgX_even_odd(p, &p0, &p1);
335 : : /* p = p0(x^2) + x p1(x^2) */
336 : 3085 : s0 = ZX_sqr(p0);
337 : 3085 : s1 = ZX_sqr(p1);
338 : 3085 : return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
339 : : }
340 : : GEN
341 : 10 : polgraeffe(GEN p)
342 : : {
343 : 10 : pari_sp av = avma;
344 : : GEN p0, p1, s0, s1;
345 : 10 : long n = degpol(p);
346 : :
347 [ - + ]: 10 : if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
348 : 10 : n = degpol(p);
349 [ - + ]: 10 : if (!n) return gcopy(p);
350 : 10 : RgX_even_odd(p, &p0, &p1);
351 : : /* p = p0(x^2) + x p1(x^2) */
352 : 10 : s0 = RgX_sqr(p0);
353 : 10 : s1 = RgX_sqr(p1);
354 : 10 : 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 : 9664181 : mydbllog2i(GEN x)
368 : : {
369 : : #ifdef LONG_IS_64BIT
370 : 7751656 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
371 : : #else
372 : 1912525 : const double W = 1/4294967296.; /*2^-32*/
373 : : #endif
374 : : GEN m;
375 : 9664181 : long lx = lgefint(x);
376 : : double l;
377 [ + + ]: 9664181 : if (lx == 2) return -pariINFINITY;
378 : 9610093 : m = int_MSW(x);
379 : 9610093 : l = (double)(ulong)*m;
380 [ + + ]: 9610093 : if (lx == 3) return log2(l);
381 : 4543379 : 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 : 9664181 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
386 : : }
387 : :
388 : : /* return log(|x|) or -pariINFINITY */
389 : : static double
390 : 702958 : mydbllogr(GEN x) {
391 [ - + ]: 702958 : if (!signe(x)) return -pariINFINITY;
392 : 702958 : return LOG2*dbllog2r(x);
393 : : }
394 : :
395 : : /* return log2(|x|) or -pariINFINITY */
396 : : static double
397 : 4851332 : mydbllog2r(GEN x) {
398 [ + + ]: 4851332 : if (!signe(x)) return -pariINFINITY;
399 : 4851332 : return dbllog2r(x);
400 : : }
401 : : static double
402 [ + + ]: 2498526 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
403 : : double
404 : 13249344 : dbllog2(GEN z)
405 : : {
406 : : double x, y;
407 [ + + + ]: 13249344 : switch(typ(z))
408 : : {
409 : 7558453 : case t_INT: return mydbllog2i(z);
410 : 4441628 : case t_REAL: return mydbllog2r(z);
411 : : default: /*t_COMPLEX*/
412 : 1249263 : x = dbllog2mp(gel(z,1));
413 : 1249263 : y = dbllog2mp(gel(z,2));
414 [ + + ]: 1249263 : if (fabs(x-y) > 10) return maxdd(x,y);
415 : 13249344 : return x + 0.5*log2(1 + exp2(2*(y-x)));
416 : : }
417 : : }
418 : : static GEN /* beware overflow */
419 [ + + ]: 267212 : 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 : 1325203 : findpower(GEN p)
425 : : {
426 : 1325203 : double x, L, mins = pariINFINITY;
427 : 1325203 : long n = degpol(p),i;
428 : :
429 : 1325203 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
430 [ + + ]: 6692968 : for (i=n-1; i>=0; i--)
431 : : {
432 : 5367765 : L += log2((double)(i+1) / (double)(n-i));
433 : 5367765 : x = dbllog2(gel(p,i+2));
434 [ + + ]: 5367765 : if (x != -pariINFINITY)
435 : : {
436 : 5336147 : double s = (L - x) / (double)(n-i);
437 [ + + ]: 5336147 : if (s < mins) mins = s;
438 : : }
439 : : }
440 : 1325203 : i = (long)ceil(mins);
441 [ + + ]: 1325203 : if (i - mins > 1 - 1e-12) i--;
442 : 1325203 : return i;
443 : : }
444 : :
445 : : /* returns the exponent for logmodulus(), from the Newton diagram */
446 : : static long
447 : 231586 : newton_polygon(GEN p, long k)
448 : : {
449 : 231586 : pari_sp av = avma;
450 : : double *logcoef, slope;
451 : 231586 : long n = degpol(p), i, j, h, l, *vertex;
452 : :
453 : 231586 : init_dalloc();
454 : 231586 : logcoef = (double*)stack_malloc((n+1)*sizeof(double));
455 : 231586 : vertex = (long*)new_chunk(n+1);
456 : :
457 : : /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
458 [ + + ]: 1275037 : for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
459 : 231586 : vertex[0] = 1; /* sentinel */
460 [ + + ]: 985509 : for (i=0; i < n; i=h)
461 : : {
462 : 753923 : slope = logcoef[i+1]-logcoef[i];
463 [ + + ]: 3591192 : for (j = h = i+1; j<=n; j++)
464 : : {
465 : 2837269 : double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
466 [ + + ]: 2837269 : if (slope < pij) { slope = pij; h = j; }
467 : : }
468 : 753923 : vertex[h] = 1;
469 : : }
470 [ + + ]: 248723 : h = k; while (!vertex[h]) h++;
471 [ + + ]: 243037 : l = k-1; while (!vertex[l]) l--;
472 : 231586 : avma = av;
473 : 231586 : 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 : 1585646 : myshiftrc(GEN z, long e)
479 : : {
480 [ + + ]: 1585646 : if (typ(z)==t_COMPLEX)
481 : : {
482 [ + + ]: 291043 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
483 [ + + ]: 291043 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
484 : : }
485 : : else
486 [ + + ]: 1294603 : if (signe(z)) shiftr_inplace(z, e);
487 : 1585646 : }
488 : :
489 : : /* return z*2^e, where z is integer or complex of integer (destroy z) */
490 : : static GEN
491 : 5758558 : myshiftic(GEN z, long e)
492 : : {
493 [ + + ]: 5758558 : if (typ(z)==t_COMPLEX)
494 : : {
495 [ + + ]: 826171 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
496 : 826171 : gel(z,2) = mpshift(gel(z,2),e);
497 : 826171 : return z;
498 : : }
499 [ + + ]: 5758558 : 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 : 273238 : myreal_1(long bit)
505 : : {
506 [ - + ]: 273238 : if (bit < 0) bit = 0;
507 : 273238 : return real_1(nbits2prec(bit));
508 : : }
509 : : static GEN
510 : 292113 : RgX_gtofp_bit(GEN q, long bit)
511 : : {
512 [ - + ]: 292113 : if (bit < 0) bit = 0;
513 : 292113 : return RgX_gtofp(q, nbits2prec(bit));
514 : : }
515 : :
516 : : static GEN
517 : 10469473 : mygprecrc(GEN x, long prec, long e)
518 : : {
519 : : GEN y;
520 [ + + + ]: 10469473 : switch(typ(x))
521 : : {
522 [ + + ]: 8213200 : case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
523 : : case t_COMPLEX:
524 : 1676233 : y = cgetg(3,t_COMPLEX);
525 : 1676233 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
526 : 1676233 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
527 : 1676233 : return y;
528 : 10469473 : 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 : 2855268 : mygprec(GEN x, long bit)
536 : : {
537 : : long lx, i, e, prec;
538 : : GEN y;
539 : :
540 [ - + ]: 2855268 : if (bit < 0) bit = 0; /* should rarely happen */
541 : 2855268 : e = gexpo(x) - bit;
542 : 2855268 : prec = nbits2prec(bit);
543 [ + + ]: 2855268 : switch(typ(x))
544 : : {
545 : : case t_POL:
546 : 1706556 : y = cgetg_copy(x, &lx); y[1] = x[1];
547 [ + + ]: 7622997 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
548 : 1706556 : break;
549 : :
550 : 1148712 : default: y = mygprecrc(x,prec,e);
551 : : }
552 : 2855268 : 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 : 692505 : pol_to_gaussint(GEN p, long shift)
559 : : {
560 : 692505 : long i, l = lg(p);
561 : 692505 : GEN q = cgetg(l, t_POL); q[1] = p[1];
562 [ + + ]: 5005476 : for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
563 : 692505 : return q;
564 : : }
565 : :
566 : : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
567 : : static GEN
568 : 550623 : eval_rel_pol(GEN p, long bit)
569 : : {
570 : : long i;
571 [ + + ]: 3921683 : for (i = 2; i < lg(p); i++)
572 [ + + ]: 3371060 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
573 : 550623 : 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 : 61587 : homothetie(GEN p, double lrho, long bit)
579 : : {
580 : : GEN q, r, t, iR;
581 : 61587 : long n = degpol(p), i;
582 : :
583 : 61587 : iR = mygprec(dblexp(-lrho),bit);
584 : 61587 : q = mygprec(p, bit);
585 : 61587 : r = cgetg(n+3,t_POL); r[1] = p[1];
586 : 61587 : t = iR; r[n+2] = q[n+2];
587 [ + + ]: 393875 : for (i=n-1; i>0; i--)
588 : : {
589 : 332288 : gel(r,i+2) = gmul(t, gel(q,i+2));
590 : 332288 : t = mulrr(t, iR);
591 : : }
592 : 61587 : 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 : 373468 : homothetie2n(GEN p, long e)
598 : : {
599 [ + + ]: 373468 : if (e)
600 : : {
601 : 282029 : long i,n = lg(p)-1;
602 [ + + ]: 1867675 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
603 : : }
604 : 373468 : }
605 : :
606 : : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
607 : : static void
608 : 1183321 : homothetie_gauss(GEN p, long e, long f)
609 : : {
610 [ + + ][ + + ]: 1183321 : if (e || f)
611 : : {
612 : 1083651 : long i, n = lg(p)-1;
613 [ + + ]: 6842209 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
614 : : }
615 : 1183321 : }
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 : 1325203 : lower_bound(GEN p, long *k, double eps)
621 : : {
622 : 1325203 : long n = degpol(p), i, j;
623 : 1325203 : pari_sp ltop = avma;
624 : : GEN a, s, S, ilc;
625 : : double r, R, rho;
626 : :
627 [ + + ]: 1325203 : if (n < 4) { *k = n; return 0.; }
628 : 605425 : S = cgetg(5,t_VEC);
629 : 605425 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
630 [ + + ]: 3027125 : 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 : 605425 : s = gel(a,1);
633 : 605425 : gel(S,1) = gneg(s); /* Newton sum S_i */
634 : 605425 : rho = r = gtodouble(gabs(s,3));
635 : 605425 : R = r / n;
636 [ + + ]: 2421700 : for (i=2; i<=4; i++)
637 : : {
638 : 1816275 : s = gmulsg(i,gel(a,i));
639 [ + + ]: 5448825 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
640 : 1816275 : gel(S,i) = gneg(s); /* Newton sum S_i */
641 : 1816275 : r = gtodouble(gabs(s,3));
642 [ + + ]: 1816275 : if (r > 0.)
643 : : {
644 : 1813775 : r = exp(log(r/n) / (double)i);
645 [ + + ]: 1813775 : if (r > R) R = r;
646 : : }
647 : : }
648 [ + + ][ + + ]: 605425 : if (R > 0. && eps < 1.2)
649 : 605218 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
650 : : else
651 : 207 : *k = n;
652 : 1325203 : avma = ltop; return R;
653 : : }
654 : :
655 : : /* log of modulus of the largest root of p with relative error tau */
656 : : static double
657 : 141882 : logmax_modulus(GEN p, double tau)
658 : : {
659 : : GEN r,q,aux,gunr;
660 : 141882 : pari_sp av, ltop = avma;
661 : 141882 : long i,k,n=degpol(p),nn,bit,M,e;
662 [ - + ]: 141882 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
663 : :
664 : 141882 : r = cgeti(BIGDEFAULTPREC);
665 : 141882 : av = avma;
666 : :
667 : 141882 : eps = - 1/log(1.5*tau2); /* > 0 */
668 : 141882 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
669 : 141882 : gunr = myreal_1(bit+2*n);
670 : 141882 : aux = gdiv(gunr, gel(p,2+n));
671 : 141882 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
672 : 141882 : e = findpower(q);
673 : 141882 : homothetie2n(q,e);
674 : 141882 : affsi(e, r);
675 : 141882 : q = pol_to_gaussint(q, bit);
676 : 141882 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
677 : 141882 : nn = n;
678 : 141882 : for (i=0,e=0;;)
679 : : { /* nn = deg(q) */
680 : 1325203 : rho = lower_bound(q, &k, eps);
681 [ + + ]: 1325203 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
682 : 1325203 : affii(shifti(addis(r,e), 1), r);
683 [ + + ]: 1325203 : if (++i == M) break;
684 : :
685 : 3549963 : bit = (long) ((double)k * log2(1./tau2) +
686 : 2366642 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
687 : 1183321 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
688 : 1183321 : nn -= RgX_valrem(q, &q);
689 : 1183321 : set_karasquare_limit(gexpo(q));
690 : 1183321 : q = gerepileupto(av, graeffe(q));
691 [ + + ]: 1183321 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
692 : 1183321 : eps = -1/log(tau2); /* > 0 */
693 : 1183321 : e = findpower(q);
694 : 1183321 : }
695 [ + + ]: 141882 : if (!signe(r)) { avma = ltop; return 0.; }
696 : 132020 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
697 : 141882 : avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
698 : : }
699 : : GEN
700 : 300 : logmax_modulus_bound(GEN P)
701 : : {
702 : 300 : return dblexp(logmax_modulus(P, 0.01) + 0.01);
703 : : }
704 : :
705 : : /* log of modulus of the smallest root of p, with relative error tau */
706 : : static double
707 : 49696 : logmin_modulus(GEN p, double tau)
708 : : {
709 : 49696 : pari_sp av = avma;
710 : : double r;
711 : :
712 [ - + ]: 49696 : if (gequal0(gel(p,2))) return -pariINFINITY;
713 : 49696 : r = - logmax_modulus(RgX_recip_shallow(p),tau);
714 : 49696 : avma = av; return r;
715 : : }
716 : :
717 : : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
718 : : static double
719 : 25519 : logmodulus(GEN p, long k, double tau)
720 : : {
721 : : GEN q;
722 : 25519 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
723 : 25519 : pari_sp av, ltop=avma;
724 : 25519 : double r, tau2 = tau/6;
725 : :
726 : 25519 : bit = (long)(n * (2. + log2(3.*n/tau2)));
727 : 25519 : av = avma;
728 : 25519 : q = gprec_w(p, nbits2prec(bit));
729 : 25519 : q = RgX_gtofp_bit(q, bit);
730 : 25519 : e = newton_polygon(q,k);
731 : 25519 : r = (double)e;
732 : 25519 : homothetie2n(q,e);
733 : 25519 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
734 [ + + ]: 231586 : for (i=1; i<imax; i++)
735 : : {
736 : 206067 : q = eval_rel_pol(q,bit);
737 : 206067 : kk -= RgX_valrem(q, &q);
738 : 206067 : nn = degpol(q);
739 : :
740 : 206067 : set_karasquare_limit(bit);
741 : 206067 : q = gerepileupto(av, graeffe(q));
742 : 206067 : e = newton_polygon(q,kk);
743 : 206067 : r += e / exp2((double)i);
744 : 206067 : q = RgX_gtofp_bit(q, bit);
745 : 206067 : homothetie2n(q,e);
746 : :
747 [ - + ]: 206067 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
748 : 206067 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
749 : : }
750 : 25519 : avma = ltop; return -r * LOG2;
751 : : }
752 : :
753 : : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
754 : : * rmin < r_k < rmax. This information helps because we may reduce precision
755 : : * quicker */
756 : : static double
757 : 25519 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
758 : : {
759 : : GEN q;
760 : 25519 : long n = degpol(p), i, imax, imax2, bit;
761 : 25519 : pari_sp ltop = avma, av;
762 : 25519 : double lrho, aux, tau2 = tau/6.;
763 : :
764 : 25519 : aux = (lrmax - lrmin) / 2. + 4*tau2;
765 : 25519 : imax = (long) log2(log((double)n)/ aux);
766 [ + + ]: 25519 : if (imax <= 0) return logmodulus(p,k,tau);
767 : :
768 : 24788 : lrho = (lrmin + lrmax) / 2;
769 : 24788 : av = avma;
770 : 24788 : bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
771 : 24788 : q = homothetie(p, lrho, bit);
772 : 24788 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
773 [ - + ]: 24788 : if (imax > imax2) imax = imax2;
774 : :
775 [ + + ]: 77188 : for (i=0; i<imax; i++)
776 : : {
777 : 52400 : q = eval_rel_pol(q,bit);
778 : 52400 : set_karasquare_limit(bit);
779 : 52400 : q = gerepileupto(av, graeffe(q));
780 : 52400 : aux = 2*aux + 2*tau2;
781 : 52400 : tau2 *= 1.5;
782 : 52400 : bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
783 : 52400 : q = RgX_gtofp_bit(q, bit);
784 : : }
785 : 24788 : aux = exp2((double)imax);
786 : 24788 : aux = logmodulus(q,k, aux*tau/3.) / aux;
787 : 25519 : avma = ltop; return lrho + aux;
788 : : }
789 : :
790 : : static double
791 : 29740 : ind_maxlog2(GEN q)
792 : : {
793 : 29740 : long i, k = -1;
794 : 29740 : double L = - pariINFINITY;
795 [ + + ]: 83481 : for (i=0; i<=degpol(q); i++)
796 : : {
797 : 53741 : double d = dbllog2(gel(q,2+i));
798 [ + + ]: 53741 : if (d > L) { L = d; k = i; }
799 : : }
800 : 29740 : return k;
801 : : }
802 : :
803 : : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
804 : : * Assume that l <= k <= n-l */
805 : : static long
806 : 36799 : dual_modulus(GEN p, double lrho, double tau, long l)
807 : : {
808 : 36799 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
809 : 36799 : double tau2 = tau * 7./8.;
810 : 36799 : pari_sp av = avma;
811 : : GEN q;
812 : :
813 : 36799 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
814 : 36799 : q = homothetie(p, lrho, bit);
815 : 36799 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
816 : :
817 [ + + ]: 321896 : for (i=0; i<imax; i++)
818 : : {
819 : 292156 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
820 : 292156 : v = RgX_valrem(q, &q);
821 [ + + ]: 292156 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
822 : :
823 : 292156 : nn = degpol(q); delta_k += v;
824 [ + + ]: 292156 : if (!nn) return delta_k;
825 : :
826 : 285097 : set_karasquare_limit(bit);
827 : 285097 : q = gerepileupto(av, graeffe(q));
828 : 285097 : tau2 *= 7./4.;
829 : 285097 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
830 : : }
831 : 36799 : avma = av; return delta_k + (long)ind_maxlog2(q);
832 : : }
833 : :
834 : : /********************************************************************/
835 : : /** **/
836 : : /** FACTORS THROUGH CIRCLE INTEGRATION **/
837 : : /** **/
838 : : /********************************************************************/
839 : : /* l power of 2 */
840 : : static void
841 : 1145584 : fft(GEN Omega, GEN p, GEN f, long step, long l)
842 : : {
843 : : pari_sp ltop;
844 : : long i, l1, l2, l3, rapi, step4;
845 : : GEN f1, f2, f3, f02, f13, g02, g13, ff;
846 : :
847 [ + + ]: 1145584 : if (l == 2)
848 : : {
849 : 656188 : gel(f,0) = gadd(gel(p,0),gel(p,step));
850 : 656188 : gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
851 : : }
852 [ + + ]: 489396 : if (l == 4)
853 : : {
854 : 254980 : f1 = gadd(gel(p,0), gel(p,step<<1));
855 : 254980 : f2 = gsub(gel(p,0), gel(p,step<<1));
856 : 254980 : f3 = gadd(gel(p,step),gel(p,3*step));
857 : 254980 : f02= gsub(gel(p,step),gel(p,3*step));
858 : 254980 : f02 = mulcxI(f02);
859 : 254980 : gel(f,0) = gadd(f1, f3);
860 : 254980 : gel(f,1) = gadd(f2, f02);
861 : 254980 : gel(f,2) = gsub(f1, f3);
862 : 254980 : gel(f,3) = gsub(f2, f02); return;
863 : : }
864 : :
865 : 234416 : ltop = avma;
866 : 234416 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
867 : 234416 : fft(Omega,p, f, step4,l1);
868 : 234416 : fft(Omega,p+step, f+l1,step4,l1);
869 : 234416 : fft(Omega,p+(step<<1),f+l2,step4,l1);
870 : 234416 : fft(Omega,p+3*step, f+l3,step4,l1);
871 : :
872 : 234416 : ff = cgetg(l+1,t_VEC);
873 [ + + ]: 944286 : for (i=0; i<l1; i++)
874 : : {
875 : 709870 : rapi = step*i;
876 : 709870 : f1 = gmul(gel(Omega,rapi), gel(f,i+l1));
877 : 709870 : f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
878 : 709870 : f3 = gmul(gel(Omega,3*rapi), gel(f,i+l3));
879 : :
880 : 709870 : f02 = gadd(gel(f,i),f2);
881 : 709870 : g02 = gsub(gel(f,i),f2);
882 : 709870 : f13 = gadd(f1,f3);
883 : 709870 : g13 = mulcxI(gsub(f1,f3));
884 : :
885 : 709870 : gel(ff,i+1) = gadd(f02, f13);
886 : 709870 : gel(ff,i+l1+1) = gadd(g02, g13);
887 : 709870 : gel(ff,i+l2+1) = gsub(f02, f13);
888 : 709870 : gel(ff,i+l3+1) = gsub(g02, g13);
889 : : }
890 : 234416 : ff = gerepilecopy(ltop,ff);
891 [ + + ]: 3985064 : for (i=0; i<l; i++) f[i] = ff[i+1];
892 : : }
893 : :
894 : : /* e(1/N) */
895 : : static GEN
896 : 73559 : RUgen(long N, long bit)
897 : : {
898 [ + + ]: 73559 : if (N == 2) return real_m1(nbits2prec(bit));
899 [ + + ]: 71031 : if (N == 4) return gen_I();
900 : 73559 : return expIr(divru(Pi2n(1, nbits2prec(bit)), N));
901 : : }
902 : :
903 : : /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
904 : : static GEN
905 : 35136 : initRU(long N, long bit)
906 : : {
907 : 35136 : GEN *RU, z = RUgen(N, bit);
908 : 35136 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
909 : :
910 : 35136 : RU = (GEN*)cgetg(N+1,t_VEC); RU++;
911 : :
912 : 35136 : RU[0] = myreal_1(bit);
913 : 35136 : RU[1] = z;
914 [ + + ]: 48533 : for (i=1; i<N8; i++)
915 : : {
916 : 13397 : GEN t = RU[i];
917 : 13397 : RU[i+1] = gmul(z, t);
918 : 13397 : RU[N4-i] = mkcomplex(gel(t,2), gel(t,1));
919 : : }
920 [ + + ]: 122929 : for (i=0; i<N4; i++) RU[i+N4] = mulcxI(RU[i]);
921 [ + + ]: 210722 : for (i=0; i<N2; i++) RU[i+N2] = gneg(RU[i]);
922 : 35136 : return (GEN)RU;
923 : : }
924 : :
925 : : /* as above, N arbitrary */
926 : : static GEN
927 : 3287 : initRUgen(long N, long bit)
928 : : {
929 : 3287 : GEN *RU = (GEN*)cgetg(N+1,t_VEC), z = RUgen(N,bit);
930 : 3287 : long i, k = (N+3)>>1;
931 : 3287 : RU[0] = gen_1;
932 : 3287 : RU[1] = z;
933 [ + + ]: 4744 : for (i=2; i<k; i++) RU[i] = gmul(z, RU[i-1]);
934 [ + + ]: 4356 : for ( ; i<N; i++) RU[i] = gconj(RU[N-i]);
935 : 3287 : return (GEN)RU;
936 : : }
937 : :
938 : : GEN
939 : 0 : FFTinit(long k, long prec)
940 : : {
941 [ # # ]: 0 : if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
942 : 0 : return initRU(1L << k, prec2nbits(prec)) - 1;
943 : : }
944 : :
945 : : GEN
946 : 0 : FFT(GEN x, GEN Omega)
947 : : {
948 : 0 : long i, l = lg(Omega), n = lg(x);
949 : : GEN y, z;
950 [ # # ]: 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
951 [ # # ]: 0 : if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
952 [ # # ]: 0 : if (n > l) pari_err_DIM("FFT");
953 : :
954 [ # # ]: 0 : if (n < l) {
955 : 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
956 [ # # ]: 0 : for (i = 1; i < n; i++) z[i] = x[i];
957 [ # # ]: 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
958 : : }
959 : 0 : else z = x;
960 : 0 : y = cgetg(l, t_VEC);
961 : 0 : fft(Omega+1, z+1, y+1, 1, l-1);
962 : 0 : return y;
963 : : }
964 : :
965 : : /* returns 1 if p has only real coefficients, 0 else */
966 : : static int
967 : 32177 : isreal(GEN p)
968 : : {
969 : : long i;
970 [ + + ]: 210724 : for (i = lg(p)-1; i > 1; i--)
971 [ + + ]: 185541 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
972 : 32177 : return 1;
973 : : }
974 : :
975 : : /* x non complex */
976 : : static GEN
977 : 25116 : abs_update_r(GEN x, double *mu) {
978 : 25116 : GEN y = gtofp(x, DEFAULTPREC);
979 [ + + ]: 25116 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
980 : 25116 : setabssign(y); return y;
981 : : }
982 : : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
983 : : static GEN
984 : 478728 : abs_update(GEN x, double *mu) {
985 : : GEN y, xr, yr;
986 : : double ly;
987 [ + + ]: 478728 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
988 : 454176 : xr = gel(x,1);
989 : 454176 : yr = gel(x,2);
990 [ + + ]: 454176 : if (gequal0(xr)) return abs_update_r(yr,mu);
991 [ + + ]: 454106 : if (gequal0(yr)) return abs_update_r(xr,mu);
992 : : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
993 : 453612 : xr = gtofp(xr, DEFAULTPREC);
994 : 453612 : yr = gtofp(yr, DEFAULTPREC);
995 : 453612 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
996 [ + + ]: 453612 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
997 : 478728 : return y;
998 : : }
999 : :
1000 : : static void
1001 : 16906 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
1002 : : int polreal, double param, double param2)
1003 : : {
1004 : : GEN q, pc, Omega, A, RU, prim, g, ONE,TWO;
1005 : 16906 : long n = degpol(p), bit, NN, K, i, j, Lmax;
1006 : 16906 : pari_sp av2, av = avma, lim = stack_lim(av, 1);
1007 : :
1008 : 16906 : bit = gexpo(p) + (long)param2+8;
1009 [ + + ]: 33809 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
1010 [ + + ]: 16906 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
1011 [ + + ]: 16906 : K = NN/Lmax; if (K & 1) K++;
1012 : 16906 : NN = Lmax*K;
1013 [ + + ]: 16906 : if (polreal) K = K/2+1;
1014 : :
1015 : 16906 : Omega = initRU(Lmax,bit);
1016 : 16906 : prim = RUgen(NN, bit);
1017 : :
1018 : 16906 : q = mygprec(p,bit) + 2;
1019 : 16906 : A = cgetg(Lmax+1,t_VEC); A++;
1020 : 16906 : pc= cgetg(Lmax+1,t_VEC); pc++;
1021 [ + + ]: 134780 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
1022 [ + + ]: 60448 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
1023 : :
1024 : 16906 : *mu = pariINFINITY;
1025 : 16906 : g = real_0_bit(-bit);
1026 : 16906 : ONE = real_1(DEFAULTPREC);
1027 : 16906 : TWO = real2n(1, DEFAULTPREC);
1028 : 16906 : av2 = avma;
1029 : 16906 : RU = myreal_1(bit);
1030 [ + + ]: 63870 : for (i=0; i<K; i++)
1031 : : {
1032 [ + + ]: 46964 : if (i) {
1033 : 30058 : GEN z = RU;
1034 [ + + ]: 210425 : for (j=1; j<n; j++)
1035 : : {
1036 : 180367 : gel(pc,j) = gmul(gel(q,j),z);
1037 : 180367 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
1038 : : }
1039 : 30058 : gel(pc,n) = gmul(gel(q,n),z);
1040 : : }
1041 : :
1042 : 46964 : fft(Omega,pc,A,1,Lmax);
1043 [ + + ][ + + ]: 54218 : if (polreal && i>0 && i<K-1)
[ + + ]
1044 [ + + ]: 114110 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
1045 : : else
1046 [ + + ]: 411582 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(ONE, abs_update(gel(A,j),mu)));
1047 : 46964 : RU = gmul(RU, prim);
1048 [ - + ]: 46964 : if (low_stack(lim, stack_lim(av,1)))
1049 : : {
1050 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
1051 : 0 : gerepileall(av2,2, &g,&RU);
1052 : : }
1053 : : }
1054 : 16906 : *gamma = mydbllog2r(divru(g,NN));
1055 : 16906 : *LMAX = Lmax; avma = av;
1056 : 16906 : }
1057 : :
1058 : : /* NN is a multiple of Lmax */
1059 : : static void
1060 : 18230 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
1061 : : {
1062 : : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
1063 : 18230 : long n = degpol(p), i, j, K;
1064 : : pari_sp ltop;
1065 : :
1066 : 18230 : Omega = initRU(Lmax,bit);
1067 : 18230 : prim = RUgen(NN, bit);
1068 : 18230 : RU = cgetg(n+2,t_VEC); RU++;
1069 : :
1070 [ + + ]: 18230 : K = NN/Lmax; if (polreal) K = K/2+1;
1071 : 18230 : q = mygprec(p,bit);
1072 : 18230 : qd = RgX_deriv(q);
1073 : :
1074 : 18230 : A = cgetg(Lmax+1,t_VEC); A++;
1075 : 18230 : B = cgetg(Lmax+1,t_VEC); B++;
1076 : 18230 : C = cgetg(Lmax+1,t_VEC); C++;
1077 : 18230 : pc = cgetg(Lmax+1,t_VEC); pc++;
1078 : 18230 : pd = cgetg(Lmax+1,t_VEC); pd++;
1079 [ + + ]: 69534 : pc[0] = q[2]; for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
1080 [ + + ]: 87764 : pd[0] = qd[2]; for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
1081 : :
1082 : 18230 : ltop = avma;
1083 : 18230 : W = cgetg(k+1,t_VEC);
1084 : 18230 : U = cgetg(k+1,t_VEC);
1085 [ + + ]: 54535 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
1086 : :
1087 : 18230 : gel(RU,0) = gen_1;
1088 : 18230 : prim2 = myreal_1(bit);
1089 [ + + ]: 58469 : for (i=0; i<K; i++)
1090 : : {
1091 : 40239 : gel(RU,1) = prim2;
1092 [ + + ]: 296675 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
1093 : : /* RU[j] = prim^(ij)= prim2^j */
1094 : :
1095 [ + + ]: 296675 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
1096 : 40239 : fft(Omega,pd,A,1,Lmax);
1097 [ + + ]: 336914 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
1098 : 40239 : fft(Omega,pc,B,1,Lmax);
1099 [ + + ]: 503631 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
1100 [ + + ]: 503631 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
1101 : 40239 : fft(Omega,B,A,1,Lmax);
1102 : 40239 : fft(Omega,C,B,1,Lmax);
1103 : :
1104 [ + + ]: 40239 : if (polreal) /* p has real coefficients */
1105 : : {
1106 [ + + ][ + + ]: 33702 : if (i>0 && i<K-1)
1107 : : {
1108 [ + + ]: 16806 : for (j=1; j<=k; j++)
1109 : : {
1110 : 13503 : gel(W,j) = gadd(gel(W,j), gshift(real_i(gmul(gel(A,j+1),gel(RU,j+1))),1));
1111 : 13503 : gel(U,j) = gadd(gel(U,j), gshift(real_i(gmul(gel(B,j),gel(RU,j))),1));
1112 : : }
1113 : : }
1114 : : else
1115 : : {
1116 [ + + ]: 87571 : for (j=1; j<=k; j++)
1117 : : {
1118 : 57172 : gel(W,j) = gadd(gel(W,j), real_i(gmul(gel(A,j+1),gel(RU,j+1))));
1119 : 57172 : gel(U,j) = gadd(gel(U,j), real_i(gmul(gel(B,j),gel(RU,j))));
1120 : : }
1121 : : }
1122 : : }
1123 : : else
1124 : : {
1125 [ + + ]: 28121 : for (j=1; j<=k; j++)
1126 : : {
1127 : 18281 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
1128 : 18281 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
1129 : : }
1130 : : }
1131 : 40239 : prim2 = gmul(prim2,prim);
1132 : 40239 : gerepileall(ltop,3, &W,&U,&prim2);
1133 : : }
1134 : :
1135 [ + + ]: 54535 : for (i=1; i<=k; i++)
1136 : : {
1137 : 36305 : aux=gel(W,i);
1138 [ + + ]: 82664 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
1139 : 36305 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
1140 : : }
1141 [ + + ]: 54535 : for (i=0; i<k; i++)
1142 : : {
1143 : 36305 : aux=gel(U,k-i);
1144 [ + + ]: 82664 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
1145 : 36305 : gel(H,i+2) = gdivgs(aux,NN);
1146 : : }
1147 : 18230 : }
1148 : :
1149 : : #define NEWTON_MAX 10
1150 : : static GEN
1151 : 92129 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
1152 : : {
1153 : 92129 : GEN H = HH, D, aux;
1154 : 92129 : pari_sp ltop = avma, lim = stack_lim(ltop, 1);
1155 : : long error, i, bit1, bit2;
1156 : :
1157 : 92129 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
1158 : 92129 : bit2 = bit + Sbit;
1159 [ + + ][ + + ]: 171659 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1160 : : {
1161 [ - + ]: 79530 : if (low_stack(lim, stack_lim(ltop,1)))
1162 : : {
1163 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
1164 : 0 : gerepileall(ltop,2, &D,&H);
1165 : : }
1166 : 79530 : bit1 = -error + Sbit;
1167 : 79530 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
1168 : 79530 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
1169 : :
1170 [ + + ]: 79530 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
1171 : 79530 : H = RgX_add(mygprec(H,bit1), aux);
1172 : 79530 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
1173 [ + + ]: 79530 : error = gexpo(D); if (error < -bit1) error = -bit1;
1174 : : }
1175 [ + + ]: 92129 : if (error > -bit/2) return NULL; /* FAIL */
1176 : 92129 : return gerepilecopy(ltop,H);
1177 : : }
1178 : :
1179 : : /* return 0 if fails, 1 else */
1180 : : static long
1181 : 18230 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
1182 : : {
1183 : 18230 : GEN f0, FF, GG, r, HH = H;
1184 : 18230 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
1185 : 18230 : pari_sp av = avma, lim = stack_lim(av, 1);
1186 : :
1187 : 18230 : FF = *F; GG = RgX_divrem(p, FF, &r);
1188 [ - + ]: 18230 : error = gexpo(r); if (error <= -bit) error = 1-bit;
1189 : 18230 : normF = gexpo(FF);
1190 : 18230 : normG = gexpo(GG);
1191 [ + + ]: 18230 : enh = gexpo(H); if (enh < 0) enh = 0;
1192 : 18230 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
1193 : 18230 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
1194 : 18230 : bit2 = bit + Sbit;
1195 [ + + ][ + + ]: 110318 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1196 : : {
1197 [ + + ][ + + ]: 92129 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
1198 [ - + ]: 92129 : if (low_stack(lim, stack_lim(av,1)))
1199 : : {
1200 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1201 : 0 : gerepileall(av,4, &FF,&GG,&r,&HH);
1202 : : }
1203 : :
1204 : 92129 : bit1 = -error + Sbit2;
1205 : 92129 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1206 : : 1-error, Sbit2);
1207 [ + + ]: 92129 : if (!HH) return 0; /* FAIL */
1208 : :
1209 : 92088 : bit1 = -error + Sbit;
1210 : 92088 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1211 : 92088 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1212 : :
1213 [ + + ]: 92088 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1214 : 92088 : FF = gadd(mygprec(FF,bit1),f0);
1215 : :
1216 [ + + ]: 92088 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1217 : 92088 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1218 [ + + ]: 92088 : error = gexpo(r); if (error < -bit1) error = -bit1;
1219 : : }
1220 [ + + ]: 18189 : if (error>-bit) return 0; /* FAIL */
1221 : 18230 : *F = FF; *G = GG; return 1;
1222 : : }
1223 : :
1224 : : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1225 : : where cd is the leading coefficient of p */
1226 : : static void
1227 : 16906 : split_fromU(GEN p, long k, double delta, long bit,
1228 : : GEN *F, GEN *G, double param, double param2)
1229 : : {
1230 : : GEN pp, FF, GG, H;
1231 : 16906 : long n = degpol(p), NN, bit2, Lmax;
1232 : 16906 : int polreal = isreal(p);
1233 : : pari_sp ltop;
1234 : : double mu, gamma;
1235 : :
1236 : 16906 : pp = gdiv(p, gel(p,2+n));
1237 : 16906 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1238 : :
1239 : 16906 : H = cgetg(k+2,t_POL); H[1] = p[1];
1240 : 16906 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1241 : 16906 : gel(FF,k+2) = gen_1;
1242 : :
1243 [ + + ]: 16906 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1244 : 16906 : NN *= Lmax; ltop = avma;
1245 : : for(;;)
1246 : : {
1247 : 18230 : bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
1248 : 18230 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1249 [ + + ]: 18230 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1250 : 1324 : NN <<= 1; avma = ltop;
1251 : 1324 : }
1252 : 16906 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1253 : 16906 : }
1254 : :
1255 : : static void
1256 : 16906 : optimize_split(GEN p, long k, double delta, long bit,
1257 : : GEN *F, GEN *G, double param, double param2)
1258 : : {
1259 : 16906 : long n = degpol(p);
1260 : : GEN FF, GG;
1261 : :
1262 [ + + ]: 16906 : if (k <= n/2)
1263 : 12805 : split_fromU(p,k,delta,bit,F,G,param,param2);
1264 : : else
1265 : : {
1266 : 4101 : split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1267 : 4101 : *F = RgX_recip_shallow(GG);
1268 : 4101 : *G = RgX_recip_shallow(FF);
1269 : : }
1270 : 16906 : }
1271 : :
1272 : : /********************************************************************/
1273 : : /** **/
1274 : : /** SEARCH FOR SEPARATING CIRCLE **/
1275 : : /** **/
1276 : : /********************************************************************/
1277 : :
1278 : : /* return p(2^e*x) *2^(-n*e) */
1279 : : static void
1280 : 0 : scalepol2n(GEN p, long e)
1281 : : {
1282 : 0 : long i,n=lg(p)-1;
1283 [ # # ]: 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1284 : 0 : }
1285 : :
1286 : : /* returns p(x/R)*R^n */
1287 : : static GEN
1288 : 142344 : scalepol(GEN p, GEN R, long bit)
1289 : : {
1290 : : GEN q,aux,gR;
1291 : : long i;
1292 : :
1293 : 142344 : aux = gR = mygprec(R,bit); q = mygprec(p,bit);
1294 [ + + ]: 680296 : for (i=lg(p)-2; i>=2; i--)
1295 : : {
1296 : 537952 : gel(q,i) = gmul(aux,gel(q,i));
1297 : 537952 : aux = gmul(aux,gR);
1298 : : }
1299 : 142344 : return q;
1300 : : }
1301 : :
1302 : : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1303 : : static GEN
1304 : 45813 : conformal_pol(GEN p, GEN a, long bit)
1305 : : {
1306 : 45813 : GEN z, r, ma = gneg(a), ca = gconj(a);
1307 : 45813 : long n = degpol(p), i;
1308 : 45813 : pari_sp av = avma, lim = stack_lim(av,2);
1309 : :
1310 : 45813 : z = mkpoln(2, ca, negr(myreal_1(bit)));
1311 : 45813 : r = scalarpol(gel(p,2+n), 0);
1312 : 45813 : for (i=n-1; ; i--)
1313 : : {
1314 : 168008 : r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
1315 : 168008 : r = gadd(r, gmul(z, gel(p,2+i)));
1316 [ + + ]: 168008 : if (i == 0) return gerepileupto(av, r);
1317 : 122195 : z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1318 [ - + ]: 122195 : if (low_stack(lim, stack_lim(av,2)))
1319 : : {
1320 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
1321 : 0 : gerepileall(av,2, &r,&z);
1322 : : }
1323 : 122195 : }
1324 : : }
1325 : :
1326 : : static const double UNDEF = -100000.;
1327 : :
1328 : : static double
1329 : 16906 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1330 : : {
1331 : 16906 : long i, n = degpol(p);
1332 : : double lrho, lrmin, lrmax;
1333 [ + + ]: 16906 : if (k > 1)
1334 : : {
1335 [ + - ][ + + ]: 17865 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1336 : 11204 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1337 : : }
1338 : : else /* k=1 */
1339 : 5702 : lrmin = logmin_modulus(p,aux);
1340 : 16906 : radii[k] = lrmin;
1341 : :
1342 [ + + ]: 16906 : if (k+1<n)
1343 : : {
1344 [ + - ][ + + ]: 32521 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1345 : 14315 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1346 : : }
1347 : : else /* k+1=n */
1348 : 2591 : lrmax = logmax_modulus(p,aux);
1349 : 16906 : radii[k+1] = lrmax;
1350 : :
1351 : 16906 : lrho = radii[k];
1352 [ + + ]: 41036 : for (i=k-1; i>=1; i--)
1353 : : {
1354 [ + + ][ + + ]: 24130 : if (radii[i] == UNDEF || radii[i] > lrho)
1355 : 16936 : radii[i] = lrho;
1356 : : else
1357 : 7194 : lrho = radii[i];
1358 : : }
1359 : 16906 : lrho = radii[k+1];
1360 [ + + ]: 76838 : for (i=k+1; i<=n; i++)
1361 : : {
1362 [ + + ][ + + ]: 59932 : if (radii[i] == UNDEF || radii[i] < lrho)
1363 : 36676 : radii[i] = lrho;
1364 : : else
1365 : 23256 : lrho = radii[i];
1366 : : }
1367 : 16906 : *delta = (lrmax - lrmin) / 2;
1368 [ + + ]: 16906 : if (*delta > 1.) *delta = 1.;
1369 : 16906 : return (lrmin + lrmax) / 2;
1370 : : }
1371 : :
1372 : : static void
1373 : 16906 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1374 : : {
1375 : 16906 : double t, param = 0., param2 = 0.;
1376 : : long i;
1377 [ + + ]: 117874 : for (i=1; i<=n; i++)
1378 : : {
1379 : 100968 : radii[i] -= lrho;
1380 : 100968 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1381 [ + + ]: 100968 : param += t; if (t > 1.) param2 += log2(t);
1382 : : }
1383 : 16906 : *par = param; *par2 = param2;
1384 : 16906 : }
1385 : :
1386 : : /* apply the conformal mapping then split from U */
1387 : : static void
1388 : 15271 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1389 : : double aux, GEN *F,GEN *G)
1390 : : {
1391 : 15271 : long bit2, n = degpol(p), i;
1392 : 15271 : pari_sp ltop = avma, av;
1393 : : GEN q, FF, GG, a, R;
1394 : : double lrho, delta, param, param2;
1395 : : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1396 : 15271 : bit2 = bit + (long)(n*3.4848775) + 1;
1397 : 15271 : a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
1398 : 15271 : a = divrs(a, -6);
1399 : 15271 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1400 : :
1401 : 15271 : av = avma;
1402 : 15271 : q = conformal_pol(mygprec(p,bit2), a, bit2);
1403 [ + + ]: 99275 : for (i=1; i<=n; i++)
1404 [ + + ]: 84004 : if (radii[i] != UNDEF) /* update array radii */
1405 : : {
1406 : 56909 : pari_sp av2 = avma;
1407 : 56909 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1408 : : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1409 : 56909 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1410 : 56909 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1411 : 56909 : avma = av2;
1412 : : }
1413 : 15271 : lrho = logradius(radii, q,k,aux/10., &delta);
1414 : 15271 : update_radius(n, radii, lrho, ¶m, ¶m2);
1415 : :
1416 : 15271 : bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
1417 : 15271 : R = mygprec(dblexp(-lrho), bit2);
1418 : 15271 : q = scalepol(q,R,bit2);
1419 : 15271 : gerepileall(av,2, &q,&R);
1420 : :
1421 : 15271 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1422 : 15271 : bit2 += n; R = invr(R);
1423 : 15271 : FF = scalepol(FF,R,bit2);
1424 : 15271 : GG = scalepol(GG,R,bit2);
1425 : :
1426 : 15271 : a = mygprec(a,bit2);
1427 : 15271 : FF = conformal_pol(FF,a,bit2);
1428 : 15271 : GG = conformal_pol(GG,a,bit2);
1429 : :
1430 : 15271 : a = invr(subsr(1, gnorm(a)));
1431 : 15271 : FF = RgX_Rg_mul(FF, powru(a,k));
1432 : 15271 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1433 : :
1434 : 15271 : *F = mygprec(FF,bit+n);
1435 : 15271 : *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1436 : 15271 : }
1437 : :
1438 : : /* split p, this time without scaling. returns in F and G two polynomials
1439 : : * such that |p-FG|< 2^(-bit)|p| */
1440 : : static void
1441 : 16906 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1442 : : {
1443 : : GEN q, FF, GG, R;
1444 : : double aux, delta, param, param2;
1445 : 16906 : long n = degpol(p), i, j, k, bit2;
1446 : : double lrmin, lrmax, lrho, *radii;
1447 : :
1448 : 16906 : init_dalloc();
1449 : 16906 : radii = (double*) stack_malloc((n+1) * sizeof(double));
1450 : :
1451 [ + + ]: 84062 : for (i=2; i<n; i++) radii[i] = UNDEF;
1452 : 16906 : aux = thickness/(double)(4 * n);
1453 : 16906 : lrmin = logmin_modulus(p, aux);
1454 : 16906 : lrmax = logmax_modulus(p, aux);
1455 : 16906 : radii[1] = lrmin;
1456 : 16906 : radii[n] = lrmax;
1457 : 16906 : i = 1; j = n;
1458 : 16906 : lrho = (lrmin + lrmax) / 2;
1459 : 16906 : k = dual_modulus(p, lrho, aux, 1);
1460 [ + + ][ + + ]: 16906 : if (5*k < n || (n < 2*k && 5*k < 4*n))
[ + + ]
1461 : 3755 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1462 : : else
1463 : 13151 : { lrmin = lrho; i=k; radii[i] = lrho; }
1464 [ + + ]: 36799 : while (j > i+1)
1465 : : {
1466 [ + + ]: 19893 : if (i+j == n+1)
1467 : 7453 : lrho = (lrmin + lrmax) / 2;
1468 : : else
1469 : : {
1470 : 12440 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1471 [ + + ]: 12440 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1472 : 9093 : else lrho = lrmin * kappa + lrmax;
1473 : 12440 : lrho /= 1+kappa;
1474 : : }
1475 : 19893 : aux = (lrmax - lrmin) / (4*(j-i));
1476 : 19893 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1477 [ + + ][ + + ]: 19893 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
[ + + ]
1478 : 13698 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1479 : : else
1480 : 6195 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1481 : : }
1482 : 16906 : aux = lrmax - lrmin;
1483 : :
1484 [ + + ]: 16906 : if (ctr)
1485 : : {
1486 : 15271 : lrho = (lrmax + lrmin) / 2;
1487 [ + + ]: 99275 : for (i=1; i<=n; i++)
1488 [ + + ]: 84004 : if (radii[i] != UNDEF) radii[i] -= lrho;
1489 : :
1490 : 15271 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1491 : 15271 : R = mygprec(dblexp(-lrho), bit2);
1492 : 15271 : q = scalepol(p,R,bit2);
1493 : 15271 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1494 : : }
1495 : : else
1496 : : {
1497 : 1635 : lrho = logradius(radii, p, k, aux/10., &delta);
1498 : 1635 : update_radius(n, radii, lrho, ¶m, ¶m2);
1499 : :
1500 : 1635 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1501 : 1635 : R = mygprec(dblexp(-lrho), bit2);
1502 : 1635 : q = scalepol(p,R,bit2);
1503 : 1635 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1504 : : }
1505 : 16906 : bit += n;
1506 : 16906 : bit2 += n; R = invr(mygprec(R,bit2));
1507 : 16906 : *F = mygprec(scalepol(FF,R,bit2), bit);
1508 : 16906 : *G = mygprec(scalepol(GG,R,bit2), bit);
1509 : 16906 : }
1510 : :
1511 : : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1512 : : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1513 : : * where the maximum modulus of the roots of p is <=1.
1514 : : * Assume sum of roots is 0. */
1515 : : static void
1516 : 15271 : split_1(GEN p, long bit, GEN *F, GEN *G)
1517 : : {
1518 : 15271 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1519 : : GEN TWO, ctr, q, qq, FF, GG, v, gr, r, newq;
1520 : : double lrmin, lrmax, lthick;
1521 : 15271 : const double LOG3 = 1.098613;
1522 : :
1523 : 15271 : lrmax = logmax_modulus(p, 0.01);
1524 : 15271 : gr = mygprec(dblexp(-lrmax), bit2);
1525 : 15271 : q = scalepol(p,gr,bit2);
1526 : :
1527 : 15271 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1528 : 15271 : TWO = myreal_1(bit2); setexpo(TWO,1);
1529 : 15271 : v = cgetg(5,t_VEC);
1530 : 15271 : gel(v,1) = TWO;
1531 : 15271 : gel(v,2) = negr(TWO);
1532 : 15271 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1533 : 15271 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1534 : 15271 : q = mygprec(q,bit2); lthick = 0;
1535 : 15271 : newq = ctr = NULL; /* -Wall */
1536 [ + + ]: 15271 : imax = polreal? 3: 4;
1537 [ + + ]: 27439 : for (i=1; i<=imax; i++)
1538 : : {
1539 : 27088 : qq = RgX_translate(q, gel(v,i));
1540 : 27088 : lrmin = logmin_modulus(qq,0.05);
1541 [ + + ]: 27088 : if (LOG3 > lrmin + lthick)
1542 : : {
1543 : 26469 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1544 [ + + ]: 26469 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1545 : : }
1546 [ + + ]: 27088 : if (lthick > LOG2) break;
1547 [ + + ][ + + ]: 14296 : if (polreal && i==2 && lthick > LOG3 - LOG2) break;
[ + + ]
1548 : : }
1549 : 15271 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
1550 : 15271 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1551 : 15271 : r = gneg(mygprec(ctr,bit2));
1552 : 15271 : FF = RgX_translate(FF,r);
1553 : 15271 : GG = RgX_translate(GG,r);
1554 : :
1555 : 15271 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1556 : 15271 : *F = scalepol(FF,gr,bit2);
1557 : 15271 : *G = scalepol(GG,gr,bit2);
1558 : 15271 : }
1559 : :
1560 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1561 : : where the maximum modulus of the roots of p is < 0.5 */
1562 : : static int
1563 : 15271 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1564 : : {
1565 : : GEN q, b, FF, GG;
1566 : 15271 : long n = degpol(p), k, bit2, eq;
1567 : 15271 : double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
1568 : :
1569 : : /* beware double overflow */
1570 [ + + ][ + - ]: 15271 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
[ - + ]
1571 : :
1572 [ + + ]: 15271 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1573 : 15271 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1574 : :
1575 : 15271 : q = mygprec(p,bit2);
1576 : 15271 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1577 : 15271 : q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
1578 : 15271 : k = 0;
1579 [ + - ][ - + ]: 15271 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1580 [ - + ]: 15271 : || gequal0(gel(q,k+2)))) k++;
1581 [ - + ]: 15271 : if (k > 0)
1582 : : {
1583 [ # # ]: 0 : if (k > n/2) k = n/2;
1584 : 0 : bit2 += k<<1;
1585 : 0 : FF = monomial(myreal_1(bit2), k, 0);
1586 : 0 : GG = RgX_shift_shallow(q, -k);
1587 : : }
1588 : : else
1589 : : {
1590 : 15271 : split_1(q,bit2,&FF,&GG);
1591 : 15271 : bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
1592 : 15271 : FF = mygprec(FF,bit2);
1593 : : }
1594 : 15271 : GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
1595 : 15271 : *F = RgX_translate(FF, b);
1596 : 15271 : *G = RgX_translate(GG, b); return 1;
1597 : : }
1598 : :
1599 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1600 : : * Assume max_modulus(p) < 2 */
1601 : : static void
1602 : 15271 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1603 : : {
1604 : : GEN FF, GG;
1605 : : long n, bit2, normp;
1606 : :
1607 [ - + ]: 30542 : if (split_0_2(p,bit,F,G)) return;
1608 : :
1609 : 0 : normp = gexpo(p);
1610 : 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1611 : 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1612 : 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1613 : 0 : scalepol2n(FF,-2);
1614 : 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1615 : 0 : *F = mygprec(FF,bit2);
1616 : 15271 : *G = mygprec(GG,bit2);
1617 : : }
1618 : :
1619 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1620 : : static void
1621 : 16906 : split_0(GEN p, long bit, GEN *F, GEN *G)
1622 : : {
1623 : 16906 : const double LOG1_9 = 0.6418539;
1624 : 16906 : long n = degpol(p), k = 0;
1625 : : GEN q;
1626 : :
1627 [ - + ][ # # ]: 16906 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1628 [ - + ]: 16906 : if (k > 0)
1629 : : {
1630 [ # # ]: 0 : if (k > n/2) k = n/2;
1631 : 0 : *F = monomial(myreal_1(bit), k, 0);
1632 : 0 : *G = RgX_shift_shallow(p, -k);
1633 : : }
1634 : : else
1635 : : {
1636 : 16906 : double lr = logmax_modulus(p, 0.05);
1637 [ + + ]: 16906 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1638 : : else
1639 : : {
1640 : 13743 : q = RgX_recip_shallow(p);
1641 : 13743 : lr = logmax_modulus(q,0.05);
1642 [ + + ]: 13743 : if (lr < LOG1_9)
1643 : : {
1644 : 12108 : split_0_1(q, bit, F, G);
1645 : 12108 : *F = RgX_recip_shallow(*F);
1646 : 12108 : *G = RgX_recip_shallow(*G);
1647 : : }
1648 : : else
1649 : 1635 : split_2(p,bit,NULL, 1.2837,F,G);
1650 : : }
1651 : : }
1652 : 16906 : }
1653 : :
1654 : : /********************************************************************/
1655 : : /** **/
1656 : : /** ERROR ESTIMATE FOR THE ROOTS **/
1657 : : /** **/
1658 : : /********************************************************************/
1659 : :
1660 : : static GEN
1661 : 52754 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1662 : : {
1663 : 52754 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1664 : : long i, j;
1665 : :
1666 : 52754 : d = cgetg(n+1,t_VEC);
1667 [ + + ]: 827302 : for (i=1; i<=n; i++)
1668 : : {
1669 [ + + ]: 774548 : if (i!=k)
1670 : : {
1671 : 721794 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1672 : 721794 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1673 : : }
1674 : : }
1675 : 52754 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1676 [ + + ]: 52754 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1677 : 52754 : eps = mulrr(rho, shatzle);
1678 : 52754 : aux = shiftr(powru(rho,n), err);
1679 : :
1680 [ + + ][ + - ]: 167879 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
[ + + ]
1681 : : {
1682 : 115125 : GEN prod = NULL; /* 1. */
1683 : 115125 : long m = n;
1684 : 115125 : epsbis = mulrr(eps, dbltor(1.25));
1685 [ + + ]: 2090005 : for (i=1; i<=n; i++)
1686 : : {
1687 [ + + ][ + + ]: 1974880 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1688 : : {
1689 : 1841663 : GEN dif = subrr(gel(d,i),eps);
1690 [ + + ]: 1841663 : prod = prod? mulrr(prod, dif): dif;
1691 : 1841663 : m--;
1692 : : }
1693 : : }
1694 [ + + ]: 115125 : eps2 = prod? divrr(aux, prod): aux;
1695 [ + + ]: 115125 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1696 : 115125 : rap = divrr(eps,eps2); eps = eps2;
1697 : : }
1698 : 52754 : return eps;
1699 : : }
1700 : :
1701 : : /* round a complex or real number x to an absolute value of 2^(-bit) */
1702 : : static GEN
1703 : 114830 : mygprec_absolute(GEN x, long bit)
1704 : : {
1705 : : long e;
1706 : : GEN y;
1707 : :
1708 [ + + + ]: 114830 : switch(typ(x))
1709 : : {
1710 : : case t_REAL:
1711 : 80300 : e = expo(x) + bit;
1712 [ + + ][ - + ]: 80300 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1713 : : case t_COMPLEX:
1714 [ + + ]: 31464 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1715 : 30612 : y = cgetg(3,t_COMPLEX);
1716 : 30612 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1717 : 30612 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1718 : 30612 : return y;
1719 : 114830 : default: return x;
1720 : : }
1721 : : }
1722 : :
1723 : : static long
1724 : 8127 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1725 : : {
1726 : 8127 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1727 : : GEN sigma, shatzle;
1728 : :
1729 : 8127 : err += (long)log2((double)n) + 1;
1730 [ - + ]: 8127 : if (err > -2) return 0;
1731 : 8127 : sigma = real2n(-err, LOWDEFAULTPREC);
1732 : : /* 2 / ((s - 1)^(1/n) - 1) */
1733 : 8127 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1734 [ + + ]: 60881 : for (i=1; i<=n; i++)
1735 : : {
1736 : 52754 : pari_sp av = avma;
1737 : 52754 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1738 : 52754 : long e = gexpo(x);
1739 [ + + ]: 52754 : avma = av; if (e > e_max) e_max = e;
1740 : 52754 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1741 : : }
1742 : 8127 : return e_max;
1743 : : }
1744 : :
1745 : : /********************************************************************/
1746 : : /** **/
1747 : : /** MAIN **/
1748 : : /** **/
1749 : : /********************************************************************/
1750 : : static GEN
1751 : 39795 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1752 : :
1753 : : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1754 : : * returns prod (x-roots_pol[i]) */
1755 : : static GEN
1756 : 41939 : split_complete(GEN p, long bit, GEN roots_pol)
1757 : : {
1758 : 41939 : long n = degpol(p);
1759 : : pari_sp ltop;
1760 : : GEN p1, F, G, a, b, m1, m2;
1761 : :
1762 [ + + ]: 41939 : if (n == 1)
1763 : : {
1764 : 10271 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1765 : 10271 : (void)append_clone(roots_pol,a); return p;
1766 : : }
1767 : 31668 : ltop = avma;
1768 [ + + ]: 31668 : if (n == 2)
1769 : : {
1770 : 14762 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1771 : 14762 : F = gsqrt(F, nbits2prec(bit));
1772 : 14762 : p1 = ginv(gmul2n(gel(p,4),1));
1773 : 14762 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1774 : 14762 : b = gmul(gsub(F,gel(p,3)), p1);
1775 : 14762 : a = append_clone(roots_pol,a);
1776 : 14762 : b = append_clone(roots_pol,b); avma = ltop;
1777 : 14762 : a = mygprec(a, 3*bit);
1778 : 14762 : b = mygprec(b, 3*bit);
1779 : 14762 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1780 : : }
1781 : 16906 : split_0(p,bit,&F,&G);
1782 : 16906 : m1 = split_complete(F,bit,roots_pol);
1783 : 16906 : m2 = split_complete(G,bit,roots_pol);
1784 : 41939 : return gerepileupto(ltop, gmul(m1,m2));
1785 : : }
1786 : :
1787 : : static GEN
1788 : 210538 : quickabs(GEN x)
1789 : : {
1790 : 210538 : const long prec = DEFAULTPREC;
1791 : : GEN y;
1792 [ + + - + : 210538 : switch(typ(x))
- ]
1793 : : {
1794 : 210028 : case t_INT: y = itor(x, prec); setabssign(y); return y;
1795 : 495 : case t_REAL: y = rtor(x, prec); setabssign(y); return y;
1796 : 0 : case t_FRAC: y = fractor(x, prec); setabssign(y); return y;
1797 : : case t_COMPLEX: {
1798 : 15 : GEN a = gel(x,1), b = gel(x,2);
1799 : : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1800 [ - + ]: 15 : if (isintzero(a)) return cxcompotor(b, prec);
1801 [ - + ]: 15 : if (isintzero(b)) return cxcompotor(a, prec);
1802 : 15 : a = cxcompotor(a, prec);
1803 : 15 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1804 : : }
1805 : 0 : default: pari_err_TYPE("quickabs",x);
1806 : 210538 : return NULL;/*not reached*/
1807 : : }
1808 : :
1809 : : }
1810 : :
1811 : : /* bound ln |largest root of p| */
1812 : : double
1813 : 43217 : cauchy_bound(GEN p)
1814 : : {
1815 : 43217 : pari_sp av = avma;
1816 : 43217 : long i, n = degpol(p);
1817 : : GEN invlc;
1818 : 43217 : double Lmax = -pariINFINITY;
1819 : :
1820 [ - + ]: 43217 : if (n <= 0) pari_err_CONSTPOL("cauchy_bound");
1821 : 43217 : invlc = invr( quickabs(gel(p,n+2)) ); /* 1 / |lc(p)| */
1822 [ + + ]: 227728 : for (i = 0; i < n; i++)
1823 : : {
1824 : 184511 : GEN y = gel(p,i+2);
1825 : : double L;
1826 [ + + ]: 184511 : if (gequal0(y)) continue;
1827 : 167321 : L = mydbllogr(mulrr(quickabs(y), invlc)) / (n-i);
1828 [ + + ]: 167321 : if (L > Lmax) Lmax = L;
1829 : : }
1830 : 43217 : avma = av; return Lmax + LOG2;
1831 : : }
1832 : :
1833 : : static GEN
1834 : 60911 : mygprecrc_special(GEN x, long prec, long e)
1835 : : {
1836 : : GEN y;
1837 [ + + + ]: 60911 : switch(typ(x))
1838 : : {
1839 : : case t_REAL:
1840 [ - + ]: 522 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1841 [ + - ]: 522 : return (prec > realprec(x))? rtor(x, prec): x;
1842 : : case t_COMPLEX:
1843 : 15 : y = cgetg(3,t_COMPLEX);
1844 : 15 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1845 : 15 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1846 : 15 : return y;
1847 : 60911 : default: return x;
1848 : : }
1849 : : }
1850 : :
1851 : : /* like mygprec but keep at least the same precision as before */
1852 : : static GEN
1853 : 8127 : mygprec_special(GEN x, long bit)
1854 : : {
1855 : : long lx, i, e, prec;
1856 : : GEN y;
1857 : :
1858 [ - + ]: 8127 : if (bit < 0) bit = 0; /* should not happen */
1859 : 8127 : e = gexpo(x) - bit;
1860 : 8127 : prec = nbits2prec(bit);
1861 [ + - ]: 8127 : switch(typ(x))
1862 : : {
1863 : : case t_POL:
1864 : 8127 : y = cgetg_copy(x, &lx); y[1] = x[1];
1865 [ + + ]: 69008 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1866 : 8127 : break;
1867 : :
1868 : 0 : default: y = mygprecrc_special(x,prec,e);
1869 : : }
1870 : 8127 : return y;
1871 : : }
1872 : :
1873 : : static GEN
1874 : 4840 : fix_roots1(GEN r)
1875 : : {
1876 : 4840 : long i, l = lg(r);
1877 : 4840 : GEN allr = cgetg(l, t_VEC);
1878 [ + + ]: 35505 : for (i=1; i<l; i++)
1879 : : {
1880 : 30665 : GEN t = gel(r,i);
1881 : 30665 : gel(allr,i) = gcopy(t); gunclone(t);
1882 : : }
1883 : 4840 : return allr;
1884 : : }
1885 : : static GEN
1886 : 8127 : fix_roots(GEN r, GEN *m, long h, long bit)
1887 : : {
1888 : : long i, j, k, l, prec;
1889 : : GEN allr, ro1;
1890 [ + + ]: 8127 : if (h == 1) return fix_roots1(r);
1891 : 3287 : ro1 = initRUgen(h, bit);
1892 : 3287 : prec = nbits2prec(bit);
1893 : 3287 : l = lg(r)-1;
1894 : 3287 : allr = cgetg(h*l+1, t_VEC);
1895 [ + + ]: 12417 : for (k=1,i=1; i<=l; i++)
1896 : : {
1897 : 9130 : GEN p2, p1 = gel(r,i);
1898 [ + + ]: 9130 : p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1899 [ + + ]: 31219 : for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1900 : 9130 : gunclone(p1);
1901 : : }
1902 : 3287 : *m = roots_to_pol(allr, 0);
1903 : 8127 : return allr;
1904 : : }
1905 : :
1906 : : static GEN
1907 : 8009 : all_roots(GEN p, long bit)
1908 : : {
1909 : : GEN lc, pd, q, roots_pol, m;
1910 : 8009 : long bit0, bit2, n = degpol(p), i, e, h;
1911 : : pari_sp av;
1912 : :
1913 : 8009 : pd = RgX_deflate_max(p, &h); lc = leading_term(pd);
1914 [ + + ]: 8009 : e = (long)((2/LOG2) * cauchy_bound(pd)); if (e < 0) e = 0;
1915 : 8009 : bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
1916 : 8009 : bit2 = bit0; e = 0;
1917 : 8009 : for (av=avma,i=1;; i++,avma=av)
1918 : : {
1919 : 8127 : roots_pol = vectrunc_init(n+1);
1920 : 8127 : bit2 += e + (n << i);
1921 : 8127 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1922 : 8127 : q[1] = evalsigne(1)|evalvarn(0);
1923 : 8127 : m = split_complete(q,bit2,roots_pol);
1924 : 8127 : roots_pol = fix_roots(roots_pol, &m, h, bit2);
1925 : 8127 : q = mygprec_special(p,bit2); lc = leading_term(q);
1926 : 8127 : q[1] = evalsigne(1)|evalvarn(0);
1927 [ + + ]: 8127 : if (h > 1) m = gmul(m,lc);
1928 : :
1929 : 8127 : e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
1930 [ - + ]: 8127 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
1931 [ + - ]: 8127 : if (e < 0)
1932 : : {
1933 : 8127 : e = bit + a_posteriori_errors(p,roots_pol,e);
1934 [ + + ]: 8127 : if (e < 0) return roots_pol;
1935 : : }
1936 [ - + ]: 118 : if (DEBUGLEVEL > 7)
1937 : 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1938 : 118 : }
1939 : : }
1940 : :
1941 : : INLINE int
1942 : 1513 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1943 : :
1944 : : static int
1945 : 502 : isexactpol(GEN p)
1946 : : {
1947 : 502 : long i,n = degpol(p);
1948 [ + + ]: 1845 : for (i=0; i<=n; i++)
1949 [ + + ]: 1513 : if (!isexactscalar(gel(p,i+2))) return 0;
1950 : 502 : return 1;
1951 : : }
1952 : :
1953 : : static long
1954 : 2148 : isvalidcoeff(GEN x)
1955 : : {
1956 [ + + + ]: 2148 : switch (typ(x))
1957 : : {
1958 : 2123 : case t_INT: case t_REAL: case t_FRAC: return 1;
1959 [ + - ][ + - ]: 15 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
1960 : : }
1961 : 2148 : return 0;
1962 : : }
1963 : :
1964 : : static void
1965 : 507 : checkvalidpol(GEN p)
1966 : : {
1967 : 507 : long i,n = lg(p);
1968 [ + + ]: 2605 : for (i=2; i<n; i++)
1969 [ + + ]: 2103 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
1970 : 502 : }
1971 : :
1972 : : static GEN
1973 : 332 : solve_exact_pol(GEN p, long bit)
1974 : : {
1975 : 332 : long i, j, k, m, n = degpol(p), iroots = 0;
1976 : 332 : GEN ex, factors, v = zerovec(n);
1977 : :
1978 : 332 : factors = ZX_squff(Q_primpart(p), &ex);
1979 [ + + ]: 664 : for (i=1; i<lg(factors); i++)
1980 : : {
1981 : 332 : GEN roots_fact = all_roots(gel(factors,i), bit);
1982 : 332 : n = degpol(gel(factors,i));
1983 : 332 : m = ex[i];
1984 [ + + ]: 1328 : for (j=1; j<=n; j++)
1985 [ + + ]: 1992 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
1986 : : }
1987 : 332 : return v;
1988 : : }
1989 : :
1990 : : /* return the roots of p with absolute error bit */
1991 : : static GEN
1992 : 502 : roots_com(GEN q, long bit)
1993 : : {
1994 : : GEN L, p;
1995 : 502 : long v = RgX_valrem_inexact(q, &p);
1996 [ - + ]: 502 : if (lg(p) == 3) L = cgetg(1,t_VEC); /* constant polynomial */
1997 [ + + ]: 502 : else L = isexactpol(p)? solve_exact_pol(p,bit): all_roots(p,bit);
1998 [ + + ]: 502 : if (v)
1999 : : {
2000 : 70 : GEN M, z, t = gel(q,2);
2001 : : long i, x, y, l, n;
2002 : :
2003 [ + - ]: 70 : if (isrationalzero(t)) x = -bit;
2004 : : else
2005 : : {
2006 : 0 : n = gexpo(t);
2007 : 0 : x = n / v; l = degpol(q);
2008 [ # # ]: 0 : for (i = v; i <= l; i++)
2009 : : {
2010 : 0 : t = gel(q,i+2);
2011 [ # # ]: 0 : if (isrationalzero(t)) continue;
2012 : 0 : y = (n - gexpo(t)) / i;
2013 [ # # ]: 0 : if (y < x) x = y;
2014 : : }
2015 : : }
2016 : 70 : z = real_0_bit(x); l = v + lg(L);
2017 : 70 : M = cgetg(l, t_VEC); L -= v;
2018 [ + + ]: 140 : for (i = 1; i <= v; i++) gel(M,i) = z;
2019 [ + + ]: 210 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
2020 : 70 : L = M;
2021 : : }
2022 : 502 : return L;
2023 : : }
2024 : :
2025 : : static GEN
2026 : 30952 : tocomplex(GEN x, long l, long bit)
2027 : : {
2028 : : GEN y;
2029 [ + + ]: 30952 : if (typ(x) == t_COMPLEX)
2030 : : {
2031 [ + + ]: 30176 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
2032 : 2810 : x = gel(x,2);
2033 : 2810 : y = cgetg(3,t_COMPLEX);
2034 : 2810 : gel(y,1) = real_0_bit(-bit);
2035 : 2810 : gel(y,2) = mygprecrc(x, l, -bit);
2036 : : }
2037 : : else
2038 : : {
2039 : 776 : y = cgetg(3,t_COMPLEX);
2040 : 776 : gel(y,1) = mygprecrc(x, l, -bit);
2041 : 776 : gel(y,2) = real_0_bit(-bit);
2042 : : }
2043 : 30952 : return y;
2044 : : }
2045 : :
2046 : : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
2047 : : * up to 2^-e absolute error */
2048 : : static int
2049 : 105992 : cmp_complex_appr(void *E, GEN x, GEN y)
2050 : : {
2051 : 105992 : long e = (long)E;
2052 : : GEN z, xi, yi, xr, yr;
2053 : : long sxi, syi;
2054 [ + + ]: 105992 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
2055 : 47712 : else { xr = x; xi = NULL; sxi = 0; }
2056 [ + + ]: 105992 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
2057 : 43987 : else { yr = y; yi = NULL; syi = 0; }
2058 : : /* Compare absolute values of imaginary parts */
2059 [ + + ]: 105992 : if (!sxi)
2060 : : {
2061 [ + + ][ + - ]: 48573 : if (syi && expo(yi) >= e) return -1;
2062 : : /* |Im x| ~ |Im y| ~ 0 */
2063 : : }
2064 [ + + ]: 57419 : else if (!syi)
2065 : : {
2066 [ + - ][ + - ]: 1902 : if (sxi && expo(xi) >= e) return 1;
2067 : : /* |Im x| ~ |Im y| ~ 0 */
2068 : : }
2069 : : else
2070 : : {
2071 : : long sz;
2072 : 55517 : z = addrr_sign(xi, 1, yi, -1);
2073 : 55517 : sz = signe(z);
2074 [ + + ][ + - ]: 55517 : if (sz && expo(z) >= e) return (int)sz;
2075 : : }
2076 : : /* |Im x| ~ |Im y|, sort according to real parts */
2077 : 64724 : z = subrr(xr, yr);
2078 [ + + ]: 64724 : if (expo(z) >= e) return (int)signe(z);
2079 : : /* Re x ~ Re y. Place negative absolute value before positive */
2080 : 105992 : return (int) (sxi - syi);
2081 : : }
2082 : :
2083 : : static GEN
2084 : 8009 : clean_roots(GEN L, long l, long bit, long clean)
2085 : : {
2086 : 8009 : long i, n = lg(L), ex = 5 - bit;
2087 : 8009 : GEN res = cgetg(n,t_COL);
2088 [ + + ]: 59863 : for (i=1; i<n; i++)
2089 : : {
2090 : 51854 : GEN c = gel(L,i);
2091 [ + + ][ + + ]: 51854 : if (clean && isrealappr(c,ex))
2092 : : {
2093 [ - + ]: 20902 : if (typ(c) == t_COMPLEX) c = gel(c,1);
2094 : 20902 : c = mygprecrc(c, l, -bit);
2095 : : }
2096 : : else
2097 : 30952 : c = tocomplex(c, l, bit);
2098 : 51854 : gel(res,i) = c;
2099 : : }
2100 : 8009 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
2101 : 8009 : return res;
2102 : : }
2103 : :
2104 : : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
2105 : : static GEN
2106 : 507 : roots_aux(GEN p, long l, long clean)
2107 : : {
2108 : 507 : pari_sp av = avma;
2109 : : long bit;
2110 : : GEN L;
2111 : :
2112 [ - + ]: 507 : if (typ(p) != t_POL)
2113 : : {
2114 [ # # ]: 0 : if (gequal0(p)) pari_err_ROOTS0("roots");
2115 [ # # ]: 0 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
2116 : 0 : return cgetg(1,t_COL); /* constant polynomial */
2117 : : }
2118 [ - + ]: 507 : if (!signe(p)) pari_err_ROOTS0("roots");
2119 : 507 : checkvalidpol(p);
2120 [ - + ]: 502 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2121 [ - + ]: 502 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2122 : 502 : bit = prec2nbits(l);
2123 : 502 : L = roots_com(p, bit);
2124 : 502 : return gerepileupto(av, clean_roots(L, l, bit, clean));
2125 : : }
2126 : : GEN
2127 : 392 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
2128 : : /* clean up roots. If root is real replace it by its real part */
2129 : : GEN
2130 : 115 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
2131 : :
2132 : : /* private variant of conjvec. Allow non rational coefficients, shallow
2133 : : * function. */
2134 : : GEN
2135 : 55 : polmod_to_embed(GEN x, long prec)
2136 : : {
2137 : 55 : GEN v, T = gel(x,1), A = gel(x,2);
2138 : : long i, l;
2139 [ + - ][ - + ]: 55 : if (typ(A) != t_POL || varn(A) != varn(T))
2140 : : {
2141 : 0 : checkvalidpol(T);
2142 : 0 : return const_col(degpol(T), A);
2143 : : }
2144 : 55 : v = cleanroots(T,prec); l = lg(v);
2145 [ + + ]: 165 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
2146 : 55 : return v;
2147 : : }
2148 : :
2149 : : GEN
2150 : 7507 : QX_complex_roots(GEN p, long l)
2151 : : {
2152 : 7507 : pari_sp av = avma;
2153 : : long bit;
2154 : : GEN L;
2155 : :
2156 [ - + ]: 7507 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
2157 [ - + ]: 7507 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2158 [ - + ]: 7507 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2159 : 7507 : bit = prec2nbits(l);
2160 : 7507 : L = all_roots(Q_primpart(p), bit);
2161 : 7507 : return gerepileupto(av, clean_roots(L, l, bit, 1));
2162 : : }
2163 : :
2164 : : /********************************************************************/
2165 : : /** **/
2166 : : /** REAL ROOTS OF INTEGER POLYNOMIAL **/
2167 : : /** **/
2168 : : /********************************************************************/
2169 : :
2170 : : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
2171 : : * The inversion is implicit (we take coefficients backwards) */
2172 : : static long
2173 : 190 : X2XP1(GEN P, unsigned long deg, GEN *Premapped)
2174 : : {
2175 : 190 : pari_sp av = avma, lim = stack_lim(av, 3);
2176 : 190 : GEN v = shallowcopy(P);
2177 : : long i, j, vlim, nb, s, s2;
2178 : : char flag;
2179 : :
2180 : 190 : vlim = deg+2;
2181 : 190 : nb = 0;
2182 : 190 : i = 0;
2183 : : do
2184 : : {
2185 [ + + ]: 1055 : for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2186 : 210 : s = -signe(gel(v, vlim));
2187 : 210 : vlim--;
2188 : 210 : i++;
2189 : : }
2190 [ + + ]: 210 : while (!s);
2191 [ + + ][ + - ]: 190 : if (vlim != deg + 1 && Premapped) setlg(v, vlim + 2);
2192 : :
2193 [ + + ]: 625 : for (; i <= deg - 1; i++)
2194 : : {
2195 : 500 : s2 = -signe(gel(v, 2));
2196 : 500 : flag = (s2 == s);
2197 [ + + ]: 1620 : for (j = 2; j < vlim; j++)
2198 : : {
2199 : 1120 : gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2200 [ + + ]: 1120 : if (flag) flag = (s2 != signe(gel(v, j+1)));
2201 : : }
2202 : :
2203 [ - + ]: 500 : if (low_stack(lim, stack_lim(av, 3)))
2204 : : {
2205 [ # # ]: 0 : if (!Premapped) setlg(v, vlim);
2206 : 0 : v = gerepileupto(av, v);
2207 [ # # ]: 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "X2XP1");
2208 : : }
2209 : :
2210 [ + + ]: 500 : if (s == signe(gel(v, vlim)))
2211 : : {
2212 : 185 : nb++;
2213 [ + + ]: 185 : if (nb >= 2) { avma = av; return 2; }
2214 : 120 : s = -s;
2215 : : }
2216 : :
2217 [ + + ][ - + ]: 435 : if (flag && !Premapped) goto END;
2218 : 435 : vlim--;
2219 : : }
2220 [ + + ]: 125 : if (s == signe(gel(v, vlim))) nb++;
2221 : : END:
2222 [ + - ][ + + ]: 125 : if (Premapped && nb == 1)
2223 : 45 : *Premapped = gerepileupto(av, v);
2224 : : else
2225 : 80 : avma = av;
2226 : 190 : return nb;
2227 : : }
2228 : :
2229 : : static long
2230 : 0 : _intervalcmp(GEN x, GEN y)
2231 : : {
2232 [ # # ]: 0 : if (typ(x) == t_VEC) x = gel(x, 1);
2233 [ # # ]: 0 : if (typ(y) == t_VEC) y = gel(y, 1);
2234 : 0 : return gcmp(x, y);
2235 : : }
2236 : :
2237 : : static GEN
2238 : 2434 : _gen_nored(void *E, GEN x) { (void)E; return x; }
2239 : : static GEN
2240 : 3032 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2241 : : static GEN
2242 : 656 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2243 : : static GEN
2244 : 867 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2245 : : static GEN
2246 : 3060 : _gen_one(void *E) { (void)E; return gen_1; }
2247 : : static GEN
2248 : 70 : _gen_zero(void *E) { (void)E; return gen_0; }
2249 : :
2250 : : static struct bb_algebra mp_algebra = { _gen_nored,_mp_add,_mp_mul,_mp_sqr,_gen_one,_gen_zero };
2251 : :
2252 : : static GEN
2253 : 5196 : _mp_cmul(void *E, GEN P, long a, GEN x) { return mpmul(gel(P, a+2), x); }
2254 : :
2255 : : /* Split the polynom P in two parts, each with unique sign.
2256 : : * Moreover compute the two parts of the derivative of P. */
2257 : : static long
2258 : 45 : split_polynoms(GEN P, long deg, long s0, GEN *Pp, GEN *Pm, GEN *Pprimep, GEN *Pprimem)
2259 : : {
2260 : 45 : long i, degneg, v = evalvarn(varn(P));
2261 [ + + ][ + - ]: 115 : for(i=1; i <= deg; i++) if (signe(gel(P, i+2)) != s0) break;
2262 : 45 : degneg = i;
2263 : 45 : *Pm = cgetg(degneg + 2, t_POL);
2264 : 45 : (*Pm)[1] = v;
2265 : 45 : *Pprimem = cgetg(degneg + 1, t_POL);
2266 : 45 : (*Pprimem)[1] = v;
2267 [ + + ]: 160 : for(i=0; i < degneg; i++)
2268 : : {
2269 : 115 : GEN elt = gel(P, i+2);
2270 : 115 : gel(*Pm, i+2) = elt;
2271 [ + + ]: 115 : if (i) gel(*Pprimem, i+1) = mului(i, elt);
2272 : : }
2273 : 45 : *Pp = cgetg(deg - degneg + 3, t_POL);
2274 : 45 : (*Pp)[1] = v;
2275 : 45 : *Pprimep = cgetg(deg - degneg + 3, t_POL);
2276 : 45 : (*Pprimep)[1] = v;
2277 [ + + ]: 150 : for(i=degneg; i <= deg; i++)
2278 : : {
2279 : 105 : GEN elt = gel(P, i+2);
2280 : 105 : gel(*Pp, i+2-degneg) = elt;
2281 : 105 : gel(*Pprimep, i+2-degneg) = mului(i, elt);
2282 : : }
2283 : 45 : *Pm = normalizepol_lg(*Pm, degneg+2);
2284 : 45 : *Pprimem = normalizepol_lg(*Pprimem, degneg+1);
2285 : 45 : *Pp = normalizepol_lg(*Pp, deg-degneg+3);
2286 : 45 : *Pprimep = normalizepol_lg(*Pprimep, deg-degneg+3);
2287 : 45 : return degpol(*Pm) + 1;
2288 : : }
2289 : :
2290 : : static GEN
2291 : 1117 : bkeval_single_power(long d, GEN V)
2292 : : {
2293 : 1117 : long mp = lg(V) - 2;
2294 [ - + ]: 1117 : if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2295 : 1117 : return gel(V, d+1);
2296 : : }
2297 : :
2298 : : static GEN
2299 : 1117 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long deg, long degneg, long bitprec)
2300 : : {
2301 : 1117 : GEN vp = gen_bkeval_powers(Pp, deg-degneg, pows, NULL, &mp_algebra, _mp_cmul);
2302 : 1117 : GEN vm = gen_bkeval_powers(Pm, degneg-1, pows, NULL, &mp_algebra, _mp_cmul);
2303 : 1117 : GEN xa = bkeval_single_power(degneg, pows);
2304 : : GEN r;
2305 : 1117 : vp = gmul(vp, xa);
2306 : 1117 : r = gadd(vp, vm);
2307 [ + + ][ + + ]: 1117 : if (expo(vp) - (signe(r) ? expo(r) : 0) > prec2nbits(realprec(vp)) - bitprec)
2308 : 44 : return NULL;
2309 : 1117 : return r;
2310 : : }
2311 : :
2312 : : /* Newton for polynom P. One solution between 0 and infinity, P' has also at
2313 : : * most one zero. P is almost certainly without zero coefficients. */
2314 : : static GEN
2315 : 45 : polsolve(GEN P, long bitprec)
2316 : : {
2317 : 45 : pari_sp av = avma;
2318 : : GEN Pp, Pm, Pprimep, Pprimem, pows;
2319 : : GEN Pprime, Pprime2;
2320 : 45 : GEN ra, rb, rc, rcold = NULL, Pc, Ppc;
2321 : 45 : long deg = degpol(P), degneg, rt;
2322 : 45 : long s0, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
2323 : : long bitaddprec, addprec;
2324 : 45 : long expoold = LONG_MAX, iter;
2325 [ - + ]: 45 : if (deg == 1)
2326 : : {
2327 : 0 : ra = negr(divrr(itor(gel(P, 2), bitprec), itor(gel(P, 3), bitprec)));
2328 : 0 : return gerepileuptoleaf(av, ra);
2329 : : }
2330 : 45 : Pprime = ZX_deriv(P); Pprime2 = ZX_deriv(Pprime);
2331 : 45 : bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
2332 : 45 : s0 = signe(gel(P, 2));
2333 : 45 : degneg = split_polynoms(P, deg, s0, &Pp, &Pm, &Pprimep, &Pprimem);
2334 : 45 : rt = maxss(degneg, brent_kung_optpow(maxss(deg-degneg, degneg-1), 2, 1));
2335 : : {
2336 : : /* Optimized Cauchy bound */
2337 : 45 : GEN summin = gen_0, absmaj;
2338 : : long i;
2339 : :
2340 : 45 : absmaj = gel(P, 2);
2341 [ + + ]: 115 : for(i=1; i < degneg; i++)
2342 : : {
2343 : 70 : GEN elt = gel(P, i+2);
2344 [ + + ]: 70 : if (absi_cmp(absmaj, elt) < 0) absmaj = elt;
2345 : : }
2346 : 45 : summin = gel(P, i+2);
2347 [ + + ]: 105 : while (++i <= deg) summin = addii(summin, gel(P, i+2));
2348 : 45 : rb = subsr(1, rdivii(absmaj, summin, cprec));
2349 : : do
2350 : : {
2351 : 45 : pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2352 : 45 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec);
2353 [ - + ]: 45 : if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
2354 [ + - ]: 45 : if (signe(Pc) != s0) break;
2355 : 0 : shiftr_inplace(rb,1);
2356 : : }
2357 : 0 : while (1);
2358 : : }
2359 : 45 : ra = NULL;
2360 : 45 : iter = 0;
2361 : : while (1)
2362 : : {
2363 : : GEN wdth;
2364 : 360 : iter++;
2365 [ + + ]: 360 : if (ra)
2366 : 290 : rc = shiftr(addrr(ra, rb), -1);
2367 : : else
2368 : 70 : rc = shiftr(rb, -1);
2369 : : do
2370 : : {
2371 : 360 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2372 : 360 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+2);
2373 [ + - ]: 360 : if (Pc) break;
2374 : 0 : cprec++; rc = rtor(rc, cprec);
2375 : 0 : } while (1);
2376 [ + + ]: 360 : if (signe(Pc) == s0)
2377 : 195 : ra = rc;
2378 : : else
2379 : 165 : rb = rc;
2380 [ + + ]: 360 : if (!ra) continue;
2381 : 335 : wdth = subrr(rb, ra);
2382 [ + + ]: 335 : if (!(iter % 8))
2383 : : {
2384 : 45 : GEN m1 = poleval(Pprime, ra), M2;
2385 [ - + ]: 45 : if (signe(m1) == s0) continue;
2386 : 45 : M2 = poleval(Pprime2, rb);
2387 [ - + ]: 45 : if (absr_cmp(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2388 : 45 : break;
2389 : : }
2390 [ - + ]: 290 : else if (gexpo(wdth) <= -bitprec)
2391 : 0 : break;
2392 : 315 : }
2393 : 45 : rc = rb;
2394 : 45 : iter = 0;
2395 : : while (1)
2396 : : {
2397 : : long exponew;
2398 : : GEN dist;
2399 : 356 : iter++;
2400 : 356 : rcold = rc;
2401 : 356 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2402 : 356 : Ppc = splitpoleval(Pprimep, Pprimem, pows, deg-1, degneg-1, bitaddprec+4);
2403 [ + - ]: 356 : if (Ppc)
2404 : 356 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+4);
2405 [ + - ][ + + ]: 356 : if (!Ppc || !Pc)
2406 : : {
2407 [ - + ]: 44 : if (cprec >= prec+addprec)
2408 : 0 : cprec++;
2409 : : else
2410 : 44 : cprec = minss(2*cprec, prec+addprec);
2411 : 44 : rc = rtor(rc, cprec); /* Backtrack one step */
2412 : 44 : continue;
2413 : : }
2414 : 312 : dist = divrr(Pc, Ppc);
2415 : 312 : rc = subrr(rc, dist);
2416 [ + - ][ - + ]: 312 : if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2417 : : {
2418 [ # # ]: 0 : if (cprec >= prec+addprec) break;
2419 : 0 : cprec = minss(2*cprec, prec+addprec);
2420 : 0 : rc = rtor(rcold, cprec); /* Backtrack one step */
2421 : 0 : continue;
2422 : : }
2423 [ + + ]: 312 : if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2424 : 251 : exponew = expo(dist);
2425 [ + + ]: 251 : if (exponew < -bitprec - 1)
2426 : : {
2427 [ + + ]: 75 : if (cprec >= prec+addprec) break;
2428 : 30 : cprec = minss(2*cprec, prec+addprec);
2429 : 30 : rc = rtor(rc, cprec);
2430 : 30 : continue;
2431 : : }
2432 [ + + ]: 176 : if (exponew > expoold - 2)
2433 : : {
2434 [ - + ]: 16 : if (cprec >= prec+addprec) break;
2435 : 16 : cprec = minss(2*cprec, prec+addprec);
2436 : 16 : rc = rtor(rc, cprec);
2437 : 16 : expoold = LONG_MAX;
2438 : 16 : continue;
2439 : : }
2440 : 160 : expoold = exponew;
2441 : 311 : }
2442 : 45 : return gerepileuptoleaf(av, rtor(rc, prec));
2443 : : }
2444 : :
2445 : : static GEN
2446 : 20 : usp(GEN Q0, long deg, long *nb_donep, long flag, long bitprec)
2447 : : {
2448 : : pari_sp av, lim;
2449 : : GEN Q, sol;
2450 : 20 : long nb_todo, nbr = 0, ind, deg0, indf, i, k, nb, j;
2451 : 20 : long listsize = 64, nb_done = 0;
2452 : : GEN c, Lc, Lk;
2453 : :
2454 : 20 : av = avma;
2455 : :
2456 : 20 : sol = const_col(deg, gen_0);
2457 : 20 : deg0 = deg;
2458 : 20 : Lc = const_vec(listsize, gen_0);
2459 : 20 : Lk = cgetg(listsize+1, t_VECSMALL);
2460 : 20 : c = gen_0;
2461 : 20 : k = Lk[1] = 0;
2462 : 20 : ind = 1; indf = 2;
2463 : 20 : Q = gcopy(Q0);
2464 : :
2465 : 20 : nb_todo = 1;
2466 : 20 : lim = stack_lim(av, 2);
2467 [ + + ]: 210 : while (nb_todo)
2468 : : {
2469 : 190 : GEN nc = gel(Lc, ind), Qremapped;
2470 [ - + ]: 190 : if (DEBUGLEVEL > 1)
2471 : 0 : err_printf("Checking interval [%Ps/2^%ld, %Ps/2^%ld)\n",
2472 : 0 : nc, Lk[ind], addis(nc, 1), Lk[ind]);
2473 : :
2474 [ + + ]: 190 : if (Lk[ind] == k + 1)
2475 : : {
2476 [ - + ]: 70 : if (DEBUGLEVEL > 2) err_printf("depth %d\n", Lk[ind]);
2477 : 70 : deg0 = deg;
2478 : 70 : setlg(Q0, deg + 3);
2479 : 70 : Q0 = ZX_rescale2n(Q0, 1);
2480 : 70 : Q = Q_primpart(Q0);
2481 : 70 : c = gen_0;
2482 : : }
2483 : :
2484 [ + + ]: 190 : if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
2485 : :
2486 [ - + ][ # # ]: 190 : if (DEBUGLEVEL > 2 && Lk[ind] == 19) err_printf("Q=%Ps\n", Q);
2487 : :
2488 : 190 : k = Lk[ind];
2489 : 190 : c = nc;
2490 : 190 : ind++;
2491 : 190 : nb_todo--;
2492 : :
2493 [ + + ]: 190 : if (equalii(gel(Q, 2), gen_0))
2494 : : {
2495 : 20 : GEN newsol = gmul2n(c, -k);
2496 [ + + ]: 35 : for (j = 1; j <= nbr; j++)
2497 [ + + ]: 20 : if (gequal(gel(sol, j), newsol)) break;
2498 [ + + ]: 20 : if (j > nbr) gel(sol, ++nbr) = newsol;
2499 : :
2500 : 20 : deg0--;
2501 [ + + ]: 105 : for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
2502 : 20 : setlg(Q, j);
2503 [ - + ]: 20 : if (DEBUGLEVEL > 1)
2504 : 0 : err_printf("2-root %Ps, Q=%Ps\n", newsol, Q);
2505 : : }
2506 [ - + ]: 170 : else if (DEBUGLEVEL > 2) err_printf("Q=%Ps\n", Q);
2507 : :
2508 [ + - ]: 190 : nb = X2XP1(Q, deg0, flag ? &Qremapped : NULL);
2509 : 190 : nb_done++;
2510 : :
2511 [ - + ]: 190 : if (DEBUGLEVEL > 1) err_printf("nb = %ld\n", nb);
2512 : :
2513 [ + + + ]: 190 : switch (nb)
2514 : : {
2515 : : case 0:
2516 [ - + ]: 60 : if (DEBUGLEVEL > 2) err_printf("No root in this interval\n");
2517 : 60 : break;
2518 : : case 1:
2519 [ - + ]: 45 : if (DEBUGLEVEL > 2) err_printf("A simple root in this interval\n");
2520 [ - + ]: 45 : if (!flag)
2521 : : {
2522 : : GEN low, hi;
2523 : 0 : low = gmul2n(c, -k);
2524 : 0 : hi = gmul2n(addiu(c,1), -k);
2525 : 0 : gel(sol, ++nbr) = mkvec2(low, hi);
2526 : : }
2527 : : else
2528 : : { /* Caveat emptor: Qremapped is the reciprocal polynomial */
2529 : 45 : GEN sr = polsolve(Qremapped, bitprec);
2530 : 45 : gel(sol, ++nbr) = rtor(gmul2n(addir(c, divrr(sr, addsr(1, sr))), -k), nbits2prec(bitprec));
2531 : : }
2532 : 45 : break;
2533 : :
2534 : : default:
2535 [ - + ]: 85 : if (DEBUGLEVEL > 2)
2536 : 0 : err_printf("Maybe at least %ld roots, retry\n", nb);
2537 [ - + ]: 85 : if (indf + 2 > listsize)
2538 : : {
2539 [ # # ]: 0 : if (ind>1)
2540 : : {
2541 [ # # ]: 0 : for (i = ind; i < indf; i++)
2542 : : {
2543 : 0 : gel(Lc, i-ind+1) = gel(Lc, i);
2544 : 0 : Lk[i-ind+1] = Lk[i];
2545 : : }
2546 : 0 : indf -= ind-1; ind = 1;
2547 : : }
2548 [ # # ]: 0 : if (indf + 2 > listsize)
2549 : : {
2550 : 0 : listsize *= 2;
2551 : 0 : Lc = vec_lengthen(Lc, listsize);
2552 : 0 : Lk = vecsmall_lengthen(Lk, listsize);
2553 : : }
2554 [ # # ]: 0 : for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2555 : : }
2556 : 85 : nc = shifti(c, 1);
2557 : 85 : gel(Lc, indf) = nc;
2558 : 85 : gel(Lc, indf + 1) = addis(nc, 1);
2559 : 85 : Lk[indf] = Lk[indf + 1] = k + 1;
2560 : 85 : indf += 2;
2561 : 85 : nb_todo += 2;
2562 : : }
2563 : :
2564 [ - + ]: 190 : if (low_stack(lim, stack_lim(av, 2)))
2565 : : {
2566 : 0 : gerepileall(av, 5, &Q0, &Q, &sol, &Lc, &Lk);
2567 [ # # ]: 190 : if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_uspensky", avma);
2568 : : }
2569 : : }
2570 : :
2571 : 20 : setlg(sol, nbr+1);
2572 : 20 : *nb_donep += nb_done;
2573 : 20 : return gerepilecopy(av, sol);
2574 : : }
2575 : :
2576 : : GEN
2577 : 20 : ZX_uspensky(GEN P, long onlypos, long flag, long bitprec)
2578 : : {
2579 : 20 : pari_sp av = avma;
2580 : : GEN sol, Pcur;
2581 : 20 : long i, nbr, deg, b, nb_done = 0;
2582 : :
2583 : 20 : deg = degpol(P);
2584 [ - + ]: 20 : if (deg == 0) return cgetg(1, t_COL);
2585 [ + + ]: 20 : if (deg == 1)
2586 : : {
2587 : 5 : GEN a = pollead(P, -1), b = gel(P,2);
2588 [ - + ][ # # ]: 5 : if (onlypos && signe(a) == signe(b)) return cgetg(1, t_COL);
2589 : 5 : retmkcol(gdiv(gneg(b), a));
2590 : : }
2591 : 15 : nbr = ZX_valrem(P, &Pcur);
2592 : 15 : deg -= nbr;
2593 [ + - ]: 15 : if (!nbr) Pcur = P;
2594 : :
2595 : 15 : b = (long)ceil(cauchy_bound(P)/LOG2);
2596 [ - + ]: 15 : if (DEBUGLEVEL > 1) err_printf("b = %d\n", b);
2597 : 15 : Pcur = ZX_unscale2n(Pcur, b); /* Maps roots to [-1;1] */
2598 : 15 : sol = usp(Pcur, deg, &nb_done, flag, bitprec);
2599 [ + + ]: 15 : if (!onlypos)
2600 : 5 : sol = concat(sol, gneg(usp(ZX_unscale(Pcur, gen_m1), deg, &nb_done, flag, bitprec)));
2601 [ + + ]: 75 : for (i = 1; i < lg(sol); i++)
2602 : : {
2603 : 60 : GEN soli = gel(sol, i);
2604 [ - + ]: 60 : if (typ(soli) == t_VEC)
2605 : : {
2606 : 0 : gel(soli, 1) = gmul2n(gel(soli, 1), b);
2607 : 0 : gel(soli, 2) = gmul2n(gel(soli, 2), b);
2608 : : }
2609 : : else
2610 : 60 : gel(sol, i) = gmul2n(soli, b);
2611 : : }
2612 [ - + ]: 15 : if (nbr)
2613 : 0 : sol = concat(zerovec(nbr), sol);
2614 : :
2615 [ - + ]: 15 : if (DEBUGLEVEL > 2)
2616 : 0 : err_printf("Number of visited nodes: %d\n", nb_done);
2617 : :
2618 [ + - ]: 15 : if (flag)
2619 : : {
2620 : 15 : long prec = nbits2prec(bitprec);
2621 [ - + ]: 15 : for (i = 1; i <= nbr; i++)
2622 [ # # ]: 0 : if (typ(gel(sol, i)) != t_REAL)
2623 : 0 : gel(sol, i) = gtofp(gel(sol, i), prec);
2624 : 15 : return gerepileupto(av, sort(sol));
2625 : : }
2626 : 20 : return gerepileupto(av, gen_sort(sol, (void *)_intervalcmp, cmp_nodata));
2627 : : }
2628 : :
2629 : : /* x a scalar */
2630 : : static GEN
2631 : 20 : rootsdeg0(GEN x)
2632 : : {
2633 [ + + ]: 20 : if (gequal0(x)) pari_err_ROOTS0("realroots");
2634 [ + + ]: 15 : if (!isvalidcoeff(x)) pari_err_TYPE("realroots",x);
2635 : 10 : return cgetg(1,t_COL); /* constant polynomial */
2636 : : }
2637 : : GEN
2638 : 35 : realroots(GEN P, long prec)
2639 : : {
2640 : 35 : pari_sp av = avma;
2641 : 35 : long nrr = 0, v = varn(P);
2642 : 35 : GEN sol = NULL, x, fa, ex;
2643 : : long i, j, k;
2644 : :
2645 [ + + ]: 35 : if (typ(P) != t_POL) return rootsdeg0(P);
2646 [ - + + ]: 20 : switch(degpol(P))
2647 : : {
2648 : 0 : case -1: return rootsdeg0(gen_0);
2649 : 5 : case 0: return rootsdeg0(gel(P,2));
2650 : : }
2651 : 15 : x = pol_x(v);
2652 : 15 : P = Q_primpart(P);
2653 [ - + ]: 15 : if (!RgX_is_ZX(P)) pari_err_TYPE("realroots",P);
2654 : 15 : fa = ZX_squff(P, &ex);
2655 [ + + ]: 35 : for (i = 1; i < lg(fa); i++)
2656 : : {
2657 : 20 : GEN Pi = gel(fa, i), soli, soli2 = NULL;
2658 : 20 : long n, nrri = 0, h;
2659 : 20 : Pi = RgX_deflate_max(Pi, &h);
2660 : 20 : soli = ZX_uspensky(Pi, !(h%2), 1, prec2nbits(prec));
2661 : 20 : n = lg(soli);
2662 [ + + ]: 20 : if (!(h % 2)) soli2 = cgetg(n, t_COL);
2663 [ + + ]: 85 : for (j = 1; j < n; j++)
2664 : : {
2665 : 65 : GEN elt = gel(soli, j);
2666 [ + - ]: 65 : if (typ(elt) != t_VEC)
2667 : : {
2668 : 65 : nrri++;
2669 : 65 : Pi = gdiv(Pi, gsub(gmul(denom(elt), x), numer(elt)));
2670 : 65 : gel(soli, j) = gtofp(elt, prec);
2671 : : }
2672 [ + + ]: 65 : if (h > 1)
2673 : : {
2674 : : /* note: elt != 0 because we are square free */
2675 : 45 : GEN sqrtroot, elt = gel(soli, j);
2676 [ + - ]: 45 : if (h == 2)
2677 : 45 : sqrtroot = gsqrt(elt, prec);
2678 : : else
2679 : : {
2680 [ # # ]: 0 : if (signe(elt) < 0)
2681 : 0 : sqrtroot = mpneg(gsqrtn(mpneg(elt), utoipos(h), NULL, prec));
2682 : : else
2683 : 0 : sqrtroot = gsqrtn(elt, utoipos(h), NULL, prec);
2684 : : }
2685 : 45 : gel(soli, j) = sqrtroot;
2686 [ + - ]: 45 : if (!(h % 2)) gel(soli2, j) = mpneg(sqrtroot);
2687 : : }
2688 : : }
2689 [ + + ]: 85 : for (j = 1; j < n; j++)
2690 : : {
2691 : 65 : GEN interval = gel(soli, j);
2692 [ + - ]: 65 : if (typ(interval) != t_VEC) continue;
2693 : 0 : gel(soli, j) = polsolve(NULL,prec2nbits(prec));
2694 : : }
2695 [ + + ]: 20 : if (!(h % 2)) soli = shallowconcat(soli, soli2);
2696 [ + + ]: 40 : for (k = 1; k <= ex[i]; k++)
2697 [ + + ]: 20 : sol = sol ? shallowconcat(sol, soli) : soli;
2698 : 20 : nrr += ex[i]*nrri;
2699 : : }
2700 : :
2701 [ - + ]: 15 : if (DEBUGLEVEL > 2)
2702 : : {
2703 : 0 : err_printf("Number of real roots: %d\n", lg(sol)-1);
2704 : 0 : err_printf(" -- of which 2-integral: %ld\n", nrr);
2705 : : }
2706 : :
2707 : 25 : return gerepilecopy(av, sort(sol));
2708 : : }
|
