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 : 6053982 : addCC(GEN x, GEN y)
35 : : {
36 : : GEN z;
37 : :
38 [ + + ]: 6053982 : if (typ(x) == t_INT)
39 : : {
40 [ + + ]: 5391654 : 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 : 662328 : z = cgetg(3,t_COMPLEX);
48 [ + + ]: 662328 : if (typ(y) == t_INT)
49 : : {
50 : 9270 : gel(z,1) = addii(gel(x,1),y);
51 : 9270 : gel(z,2) = icopy(gel(x,2)); return z;
52 : : }
53 : : /* ty == t_COMPLEX */
54 : 653058 : gel(z,1) = addii(gel(x,1),gel(y,1));
55 : 6053982 : 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 : 11492894 : mulCC(GEN x, GEN y)
60 : : {
61 : : GEN z;
62 : :
63 [ + + ]: 11492894 : if (typ(x) == t_INT)
64 : : {
65 [ + + ]: 9876067 : 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 : 1616827 : z = cgetg(3,t_COMPLEX);
73 [ + + ]: 1616827 : if (typ(y) == t_INT)
74 : : {
75 : 472861 : gel(z,1) = mulii(gel(x,1),y);
76 : 472861 : gel(z,2) = mulii(gel(x,2),y); return z;
77 : : }
78 : : /* ty == t_COMPLEX */
79 : : {
80 : 1143966 : pari_sp av = avma, tetpil;
81 : : GEN p1, p2;
82 : :
83 : 1143966 : p1 = mulii(gel(x,1),gel(y,1));
84 : 1143966 : p2 = mulii(gel(x,2),gel(y,2));
85 : 1143966 : y = mulii(addii(gel(x,1),gel(x,2)),
86 : 2287932 : addii(gel(y,1),gel(y,2)));
87 : 1143966 : x = addii(p1,p2); tetpil = avma;
88 : 1143966 : gel(z,1) = subii(p1,p2);
89 : 1143966 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
90 : 11492894 : return z;
91 : : }
92 : : }
93 : : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
94 : : static GEN
95 : 9536203 : sqrCC(GEN x)
96 : : {
97 : : GEN z;
98 : :
99 [ + + ]: 9536203 : if (typ(x) == t_INT) return sqri(x);
100 : : /* tx == t_COMPLEX */
101 : 1587125 : z = cgetg(3,t_COMPLEX);
102 : : {
103 : 1587125 : pari_sp av = avma, tetpil;
104 : : GEN y, p1, p2;
105 : :
106 : 1587125 : p1 = sqri(gel(x,1));
107 : 1587125 : p2 = sqri(gel(x,2));
108 : 1587125 : y = sqri(addii(gel(x,1),gel(x,2)));
109 : 1587125 : x = addii(p1,p2); tetpil = avma;
110 : 1587125 : gel(z,1) = subii(p1,p2);
111 : 1587125 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
112 : 9536203 : return z;
113 : : }
114 : : }
115 : :
116 : : static void
117 : 1602115 : set_karasquare_limit(long bit)
118 : : {
119 [ + + ]: 1602115 : 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 : 1602115 : }
125 : :
126 : : /* assume lP > 0, lP = lgpol(P) */
127 : : static GEN
128 : 3316843 : CX_square_spec(GEN P, long lP)
129 : : {
130 : : GEN s, t;
131 : 3316843 : long i, j, l, nn, n = lP - 1;
132 : : pari_sp av;
133 : :
134 : 3316843 : nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
135 : 3316843 : gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
136 [ + + ]: 8755755 : for (i=1; i<=n; i++)
137 : : {
138 : 5438912 : av = avma; l = (i+1)>>1;
139 : 5438912 : t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
140 [ + + ]: 7516390 : for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
141 : 5438912 : t = gmul2n(t,1);
142 [ + + ]: 5438912 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
143 : 5438912 : gel(s,i+2) = gerepileupto(av, t);
144 : : }
145 : 3316843 : gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
146 [ + + ]: 6219360 : for ( ; i<nn; i++)
147 : : {
148 : 2902517 : av = avma; l = (i+1)>>1;
149 : 2902517 : t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
150 [ + + ]: 3976504 : for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
151 : 2902517 : t = gmul2n(t,1);
152 [ + + ]: 2902517 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
153 : 2902517 : gel(s,i+2) = gerepileupto(av, t);
154 : : }
155 : 3316843 : return normalizepol_lg(s, nn+3);
156 : : }
157 : : /* not stack clean */
158 : : static GEN
159 : 56479 : RgX_addspec(GEN x, long nx, GEN y, long ny)
160 : : {
161 : : GEN z, t;
162 : : long i;
163 [ + + ]: 56479 : if (nx == ny) {
164 : 30942 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
165 [ + + ]: 229646 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
166 : 30942 : return normalizepol_lg(z, nx+2);
167 : : }
168 [ + - ]: 25537 : if (ny < nx) {
169 : 25537 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
170 [ + + ]: 179432 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
171 [ + + ]: 51074 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
172 : 25537 : 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 : 56479 : 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 : 3373667 : karasquare(GEN P, long nP)
206 : : {
207 : : GEN Q, s0, s1, s2, a, t;
208 : 3373667 : long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
209 : : pari_sp av;
210 : :
211 [ + + ][ + + ]: 3373667 : if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
212 : 56479 : av = avma;
213 : 56479 : n0 = (n>>1) + 1; n1 = nP - n0;
214 : 56479 : s0 = karasquare(P, n0); Q = P + n0;
215 : 56479 : s2 = karasquare(Q, n1);
216 : 56479 : s1 = RgX_addspec(P, n0, Q, n1);
217 : 56479 : s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
218 : 56479 : N = (n<<1) + 1;
219 : 56479 : a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
220 : 56479 : t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
221 [ + + ]: 741562 : for (i=0; i < l; i++) gel(t,i) = gel(s0,i);
222 [ + + ]: 127668 : for ( ; i < N0; i++) gel(t,i) = gen_0;
223 : 56479 : t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
224 [ + + ]: 661293 : for (i=0; i < l; i++) gel(t,i) = gel(s2,i);
225 [ + + ]: 100384 : for ( ; i < N1; i++) gel(t,i) = gen_0;
226 : 56479 : t = a+2 + n0; l = lgpol(s1); s1 += 2;
227 [ + + ]: 698420 : for (i=0; i < l; i++) gel(t,i) = gadd(gel(t,i), gel(s1,i));
228 : 3373667 : return gerepilecopy(av, normalizepol_lg(a, N+2));
229 : : }
230 : : /* spec function. nP = lgpol(P) */
231 : : static GEN
232 : 3204230 : cook_square(GEN P, long nP)
233 : : {
234 : : GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
235 : 3204230 : long n0, n3, i, j, n = nP - 1;
236 : : pari_sp av;
237 : :
238 [ + - ][ + - ]: 3204230 : 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 : 3204230 : return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
303 : : }
304 : :
305 : : static GEN
306 : 1602115 : graeffe(GEN p)
307 : : {
308 : : GEN p0, p1, s0, s1;
309 : 1602115 : long n = degpol(p), n0, n1, i;
310 : :
311 [ - + ]: 1602115 : if (!n) return gcopy(p);
312 : 1602115 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
313 : 1602115 : p0 = new_chunk(n0);
314 : 1602115 : p1 = new_chunk(n1);
315 [ + + ]: 5392127 : for (i=0; i<n1; i++)
316 : : {
317 : 3790012 : p0[i] = p[2+(i<<1)];
318 : 3790012 : p1[i] = p[3+(i<<1)];
319 : : }
320 [ + + ]: 1602115 : if (n1 != n0)
321 : 801945 : p0[i] = p[2+(i<<1)];
322 : 1602115 : s0 = cook_square(p0, n0);
323 : 1602115 : s1 = cook_square(p1, n1);
324 : 1602115 : return RgX_sub(s0, RgX_shift_shallow(s1,1));
325 : : }
326 : : GEN
327 : 3200 : ZX_graeffe(GEN p)
328 : : {
329 : 3200 : pari_sp av = avma;
330 : : GEN p0, p1, s0, s1;
331 : 3200 : long n = degpol(p);
332 : :
333 [ - + ]: 3200 : if (!n) return ZX_copy(p);
334 : 3200 : RgX_even_odd(p, &p0, &p1);
335 : : /* p = p0(x^2) + x p1(x^2) */
336 : 3200 : s0 = ZX_sqr(p0);
337 : 3200 : s1 = ZX_sqr(p1);
338 : 3200 : 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 : 9103677 : mydbllog2i(GEN x)
368 : : {
369 : : #ifdef LONG_IS_64BIT
370 : 7295948 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
371 : : #else
372 : 1807729 : const double W = 1/4294967296.; /*2^-32*/
373 : : #endif
374 : : GEN m;
375 : 9103677 : long lx = lgefint(x);
376 : : double l;
377 [ + + ]: 9103677 : if (lx == 2) return -pariINFINITY;
378 : 9050899 : m = int_MSW(x);
379 : 9050899 : l = (double)(ulong)*m;
380 [ + + ]: 9050899 : if (lx == 3) return log2(l);
381 : 4292404 : 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 : 9103677 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
386 : : }
387 : :
388 : : /* return log(|x|) or -pariINFINITY */
389 : : static double
390 : 496525 : mydbllogr(GEN x) {
391 [ - + ]: 496525 : if (!signe(x)) return -pariINFINITY;
392 : 496525 : return LOG2*dbllog2r(x);
393 : : }
394 : :
395 : : /* return log2(|x|) or -pariINFINITY */
396 : : static double
397 : 5002056 : mydbllog2r(GEN x) {
398 [ + + ]: 5002056 : if (!signe(x)) return -pariINFINITY;
399 : 5002056 : return dbllog2r(x);
400 : : }
401 : : static double
402 [ + + ]: 2491346 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
403 : : double
404 : 12623131 : dbllog2(GEN z)
405 : : {
406 : : double x, y;
407 [ + + + ]: 12623131 : switch(typ(z))
408 : : {
409 : 7004175 : case t_INT: return mydbllog2i(z);
410 : 4373283 : case t_REAL: return mydbllog2r(z);
411 : : default: /*t_COMPLEX*/
412 : 1245673 : x = dbllog2mp(gel(z,1));
413 : 1245673 : y = dbllog2mp(gel(z,2));
414 [ + + ]: 1245673 : if (fabs(x-y) > 10) return maxdd(x,y);
415 : 12623131 : return x + 0.5*log2(1 + exp2(2*(y-x)));
416 : : }
417 : : }
418 : : static GEN /* beware overflow */
419 [ + + ]: 245839 : 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 : 1232797 : findpower(GEN p)
425 : : {
426 : 1232797 : double x, L, mins = pariINFINITY;
427 : 1232797 : long n = degpol(p),i;
428 : :
429 : 1232797 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
430 [ + + ]: 6274125 : for (i=n-1; i>=0; i--)
431 : : {
432 : 5041328 : L += log2((double)(i+1) / (double)(n-i));
433 : 5041328 : x = dbllog2(gel(p,i+2));
434 [ + + ]: 5041328 : if (x != -pariINFINITY)
435 : : {
436 : 5011347 : double s = (L - x) / (double)(n-i);
437 [ + + ]: 5011347 : if (s < mins) mins = s;
438 : : }
439 : : }
440 : 1232797 : i = (long)ceil(mins);
441 [ + + ]: 1232797 : if (i - mins > 1 - 1e-12) i--;
442 : 1232797 : return i;
443 : : }
444 : :
445 : : /* returns the exponent for logmodulus(), from the Newton diagram */
446 : : static long
447 : 213675 : newton_polygon(GEN p, long k)
448 : : {
449 : 213675 : pari_sp av = avma;
450 : : double *logcoef, slope;
451 : 213675 : long n = degpol(p), i, j, h, l, *vertex;
452 : :
453 : 213675 : init_dalloc();
454 : 213675 : logcoef = (double*)stack_malloc((n+1)*sizeof(double));
455 : 213675 : vertex = (long*)new_chunk(n+1);
456 : :
457 : : /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
458 [ + + ]: 1178130 : for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
459 : 213675 : vertex[0] = 1; /* sentinel */
460 [ + + ]: 906996 : for (i=0; i < n; i=h)
461 : : {
462 : 693321 : slope = logcoef[i+1]-logcoef[i];
463 [ + + ]: 3295675 : for (j = h = i+1; j<=n; j++)
464 : : {
465 : 2602354 : double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
466 [ + + ]: 2602354 : if (slope < pij) { slope = pij; h = j; }
467 : : }
468 : 693321 : vertex[h] = 1;
469 : : }
470 [ + + ]: 230898 : h = k; while (!vertex[h]) h++;
471 [ + + ]: 225090 : l = k-1; while (!vertex[l]) l--;
472 : 213675 : avma = av;
473 : 213675 : 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 : 1475417 : myshiftrc(GEN z, long e)
479 : : {
480 [ + + ]: 1475417 : if (typ(z)==t_COMPLEX)
481 : : {
482 [ + + ]: 289831 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
483 [ + + ]: 289831 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
484 : : }
485 : : else
486 [ + + ]: 1185586 : if (signe(z)) shiftr_inplace(z, e);
487 : 1475417 : }
488 : :
489 : : /* return z*2^e, where z is integer or complex of integer (destroy z) */
490 : : static GEN
491 : 5401111 : myshiftic(GEN z, long e)
492 : : {
493 [ + + ]: 5401111 : if (typ(z)==t_COMPLEX)
494 : : {
495 [ + + ]: 824526 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
496 : 824526 : gel(z,2) = mpshift(gel(z,2),e);
497 : 824526 : return z;
498 : : }
499 [ + + ]: 5401111 : 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 : 252427 : myreal_1(long bit)
505 : : {
506 [ - + ]: 252427 : if (bit < 0) bit = 0;
507 : 252427 : return real_1(nbits2prec(bit));
508 : : }
509 : : static GEN
510 : 269287 : RgX_gtofp_bit(GEN q, long bit)
511 : : {
512 [ - + ]: 269287 : if (bit < 0) bit = 0;
513 : 269287 : return RgX_gtofp(q, nbits2prec(bit));
514 : : }
515 : :
516 : : static GEN
517 : 9830130 : mygprecrc(GEN x, long prec, long e)
518 : : {
519 : : GEN y;
520 [ + + + ]: 9830130 : switch(typ(x))
521 : : {
522 [ + + ]: 7633907 : case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
523 : : case t_COMPLEX:
524 : 1672066 : y = cgetg(3,t_COMPLEX);
525 : 1672066 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
526 : 1672066 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
527 : 1672066 : return y;
528 : 9830130 : 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 : 2611274 : mygprec(GEN x, long bit)
536 : : {
537 : : long lx, i, e, prec;
538 : : GEN y;
539 : :
540 [ - + ]: 2611274 : if (bit < 0) bit = 0; /* should rarely happen */
541 : 2611274 : e = gexpo(x) - bit;
542 : 2611274 : prec = nbits2prec(bit);
543 [ + + ]: 2611274 : switch(typ(x))
544 : : {
545 : : case t_POL:
546 : 1538339 : y = cgetg_copy(x, &lx); y[1] = x[1];
547 [ + + ]: 6904246 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
548 : 1538339 : break;
549 : :
550 : 1072935 : default: y = mygprecrc(x,prec,e);
551 : : }
552 : 2611274 : 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 : 639658 : pol_to_gaussint(GEN p, long shift)
559 : : {
560 : 639658 : long i, l = lg(p);
561 : 639658 : GEN q = cgetg(l, t_POL); q[1] = p[1];
562 [ + + ]: 4636289 : for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
563 : 639658 : return q;
564 : : }
565 : :
566 : : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
567 : : static GEN
568 : 507725 : eval_rel_pol(GEN p, long bit)
569 : : {
570 : : long i;
571 [ + + ]: 3626671 : for (i = 2; i < lg(p); i++)
572 [ + + ]: 3118946 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
573 : 507725 : 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 : 56543 : homothetie(GEN p, double lrho, long bit)
579 : : {
580 : : GEN q, r, t, iR;
581 : 56543 : long n = degpol(p), i;
582 : :
583 : 56543 : iR = mygprec(dblexp(-lrho),bit);
584 : 56543 : q = mygprec(p, bit);
585 : 56543 : r = cgetg(n+3,t_POL); r[1] = p[1];
586 : 56543 : t = iR; r[n+2] = q[n+2];
587 [ + + ]: 362650 : for (i=n-1; i>0; i--)
588 : : {
589 : 306107 : gel(r,i+2) = gmul(t, gel(q,i+2));
590 : 306107 : t = mulrr(t, iR);
591 : : }
592 : 56543 : 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 : 345608 : homothetie2n(GEN p, long e)
598 : : {
599 [ + + ]: 345608 : if (e)
600 : : {
601 : 262178 : long i,n = lg(p)-1;
602 [ + + ]: 1737595 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
603 : : }
604 : 345608 : }
605 : :
606 : : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
607 : : static void
608 : 1100864 : homothetie_gauss(GEN p, long e, long f)
609 : : {
610 [ + + ][ + + ]: 1100864 : if (e || f)
611 : : {
612 : 1008414 : long i, n = lg(p)-1;
613 [ + + ]: 6409525 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
614 : : }
615 : 1100864 : }
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 : 1232797 : lower_bound(GEN p, long *k, double eps)
621 : : {
622 : 1232797 : long n = degpol(p), i, j;
623 : 1232797 : pari_sp ltop = avma;
624 : : GEN a, s, S, ilc;
625 : : double r, R, rho;
626 : :
627 [ + + ]: 1232797 : if (n < 4) { *k = n; return 0.; }
628 : 573525 : S = cgetg(5,t_VEC);
629 : 573525 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
630 [ + + ]: 2867625 : 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 : 573525 : s = gel(a,1);
633 : 573525 : gel(S,1) = gneg(s); /* Newton sum S_i */
634 : 573525 : rho = r = gtodouble(gabs(s,3));
635 : 573525 : R = r / n;
636 [ + + ]: 2294100 : for (i=2; i<=4; i++)
637 : : {
638 : 1720575 : s = gmulsg(i,gel(a,i));
639 [ + + ]: 5161725 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
640 : 1720575 : gel(S,i) = gneg(s); /* Newton sum S_i */
641 : 1720575 : r = gtodouble(gabs(s,3));
642 [ + + ]: 1720575 : if (r > 0.)
643 : : {
644 : 1718095 : r = exp(log(r/n) / (double)i);
645 [ + + ]: 1718095 : if (r > R) R = r;
646 : : }
647 : : }
648 [ + + ][ + + ]: 573525 : if (R > 0. && eps < 1.2)
649 : 573318 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
650 : : else
651 : 207 : *k = n;
652 : 1232797 : avma = ltop; return R;
653 : : }
654 : :
655 : : /* log of modulus of the largest root of p with relative error tau. Assume
656 : : * P(0) != 0 and P non constant */
657 : : static double
658 : 131933 : logmax_modulus(GEN p, double tau)
659 : : {
660 : : GEN r,q,aux,gunr;
661 : 131933 : pari_sp av, ltop = avma;
662 : 131933 : long i,k,n=degpol(p),nn,bit,M,e;
663 [ - + ]: 131933 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
664 : :
665 : 131933 : r = cgeti(BIGDEFAULTPREC);
666 : 131933 : av = avma;
667 : :
668 : 131933 : eps = - 1/log(1.5*tau2); /* > 0 */
669 : 131933 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
670 : 131933 : gunr = myreal_1(bit+2*n);
671 : 131933 : aux = gdiv(gunr, gel(p,2+n));
672 : 131933 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
673 : 131933 : e = findpower(q);
674 : 131933 : homothetie2n(q,e);
675 : 131933 : affsi(e, r);
676 : 131933 : q = pol_to_gaussint(q, bit);
677 : 131933 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
678 : 131933 : nn = n;
679 : 131933 : for (i=0,e=0;;)
680 : : { /* nn = deg(q) */
681 : 1232797 : rho = lower_bound(q, &k, eps);
682 [ + + ]: 1232797 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
683 : 1232797 : affii(shifti(addis(r,e), 1), r);
684 [ + + ]: 1232797 : if (++i == M) break;
685 : :
686 : 3302592 : bit = (long) ((double)k * log2(1./tau2) +
687 : 2201728 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
688 : 1100864 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
689 : 1100864 : nn -= RgX_valrem(q, &q);
690 : 1100864 : set_karasquare_limit(gexpo(q));
691 : 1100864 : q = gerepileupto(av, graeffe(q));
692 [ + + ]: 1100864 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
693 : 1100864 : eps = -1/log(tau2); /* > 0 */
694 : 1100864 : e = findpower(q);
695 : 1100864 : }
696 [ + + ]: 131933 : if (!signe(r)) { avma = ltop; return 0.; }
697 : 123431 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
698 : 131933 : avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
699 : : }
700 : : /* assume P non constant */
701 : : GEN
702 : 320 : logmax_modulus_bound(GEN P)
703 : : {
704 : 320 : (void)RgX_valrem_inexact(P,&P);
705 [ + + ]: 320 : return lg(P)==3? gen_0: dblexp(logmax_modulus(P, 0.01) + 0.01);
706 : : }
707 : :
708 : : /* log of modulus of the smallest root of p, with relative error tau */
709 : : static double
710 : 46593 : logmin_modulus(GEN p, double tau)
711 : : {
712 : 46593 : pari_sp av = avma;
713 : : double r;
714 : :
715 [ - + ]: 46593 : if (gequal0(gel(p,2))) return -pariINFINITY;
716 : 46593 : r = - logmax_modulus(RgX_recip_shallow(p),tau);
717 : 46593 : avma = av; return r;
718 : : }
719 : :
720 : : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
721 : : static double
722 : 23485 : logmodulus(GEN p, long k, double tau)
723 : : {
724 : : GEN q;
725 : 23485 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
726 : 23485 : pari_sp av, ltop=avma;
727 : 23485 : double r, tau2 = tau/6;
728 : :
729 : 23485 : bit = (long)(n * (2. + log2(3.*n/tau2)));
730 : 23485 : av = avma;
731 : 23485 : q = gprec_w(p, nbits2prec(bit));
732 : 23485 : q = RgX_gtofp_bit(q, bit);
733 : 23485 : e = newton_polygon(q,k);
734 : 23485 : r = (double)e;
735 : 23485 : homothetie2n(q,e);
736 : 23485 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
737 [ + + ]: 213675 : for (i=1; i<imax; i++)
738 : : {
739 : 190190 : q = eval_rel_pol(q,bit);
740 : 190190 : kk -= RgX_valrem(q, &q);
741 : 190190 : nn = degpol(q);
742 : :
743 : 190190 : set_karasquare_limit(bit);
744 : 190190 : q = gerepileupto(av, graeffe(q));
745 : 190190 : e = newton_polygon(q,kk);
746 : 190190 : r += e / exp2((double)i);
747 : 190190 : q = RgX_gtofp_bit(q, bit);
748 : 190190 : homothetie2n(q,e);
749 : :
750 [ - + ]: 190190 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
751 : 190190 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
752 : : }
753 : 23485 : avma = ltop; return -r * LOG2;
754 : : }
755 : :
756 : : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
757 : : * rmin < r_k < rmax. This information helps because we may reduce precision
758 : : * quicker */
759 : : static double
760 : 23485 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
761 : : {
762 : : GEN q;
763 : 23485 : long n = degpol(p), i, imax, imax2, bit;
764 : 23485 : pari_sp ltop = avma, av;
765 : 23485 : double lrho, aux, tau2 = tau/6.;
766 : :
767 : 23485 : aux = (lrmax - lrmin) / 2. + 4*tau2;
768 : 23485 : imax = (long) log2(log((double)n)/ aux);
769 [ + + ]: 23485 : if (imax <= 0) return logmodulus(p,k,tau);
770 : :
771 : 22899 : lrho = (lrmin + lrmax) / 2;
772 : 22899 : av = avma;
773 : 22899 : bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
774 : 22899 : q = homothetie(p, lrho, bit);
775 : 22899 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
776 [ - + ]: 22899 : if (imax > imax2) imax = imax2;
777 : :
778 [ + + ]: 71672 : for (i=0; i<imax; i++)
779 : : {
780 : 48773 : q = eval_rel_pol(q,bit);
781 : 48773 : set_karasquare_limit(bit);
782 : 48773 : q = gerepileupto(av, graeffe(q));
783 : 48773 : aux = 2*aux + 2*tau2;
784 : 48773 : tau2 *= 1.5;
785 : 48773 : bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
786 : 48773 : q = RgX_gtofp_bit(q, bit);
787 : : }
788 : 22899 : aux = exp2((double)imax);
789 : 22899 : aux = logmodulus(q,k, aux*tau/3.) / aux;
790 : 23485 : avma = ltop; return lrho + aux;
791 : : }
792 : :
793 : : static double
794 : 27170 : ind_maxlog2(GEN q)
795 : : {
796 : 27170 : long i, k = -1;
797 : 27170 : double L = - pariINFINITY;
798 [ + + ]: 75959 : for (i=0; i<=degpol(q); i++)
799 : : {
800 : 48789 : double d = dbllog2(gel(q,2+i));
801 [ + + ]: 48789 : if (d > L) { L = d; k = i; }
802 : : }
803 : 27170 : return k;
804 : : }
805 : :
806 : : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
807 : : * Assume that l <= k <= n-l */
808 : : static long
809 : 33644 : dual_modulus(GEN p, double lrho, double tau, long l)
810 : : {
811 : 33644 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
812 : 33644 : double tau2 = tau * 7./8.;
813 : 33644 : pari_sp av = avma;
814 : : GEN q;
815 : :
816 : 33644 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
817 : 33644 : q = homothetie(p, lrho, bit);
818 : 33644 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
819 : :
820 [ + + ]: 295932 : for (i=0; i<imax; i++)
821 : : {
822 : 268762 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
823 : 268762 : v = RgX_valrem(q, &q);
824 [ + + ]: 268762 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
825 : :
826 : 268762 : nn = degpol(q); delta_k += v;
827 [ + + ]: 268762 : if (!nn) return delta_k;
828 : :
829 : 262288 : set_karasquare_limit(bit);
830 : 262288 : q = gerepileupto(av, graeffe(q));
831 : 262288 : tau2 *= 7./4.;
832 : 262288 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
833 : : }
834 : 33644 : avma = av; return delta_k + (long)ind_maxlog2(q);
835 : : }
836 : :
837 : : /********************************************************************/
838 : : /** **/
839 : : /** FACTORS THROUGH CIRCLE INTEGRATION **/
840 : : /** **/
841 : : /********************************************************************/
842 : : /* l power of 2 */
843 : : static void
844 : 1035892 : fft(GEN Omega, GEN p, GEN f, long step, long l)
845 : : {
846 : : pari_sp ltop;
847 : : long i, l1, l2, l3, rapi, step4;
848 : : GEN f1, f2, f3, f02, f13, g02, g13, ff;
849 : :
850 [ + + ]: 1035892 : if (l == 2)
851 : : {
852 : 597728 : gel(f,0) = gadd(gel(p,0),gel(p,step));
853 : 597728 : gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
854 : : }
855 [ + + ]: 438164 : if (l == 4)
856 : : {
857 : 226380 : f1 = gadd(gel(p,0), gel(p,step<<1));
858 : 226380 : f2 = gsub(gel(p,0), gel(p,step<<1));
859 : 226380 : f3 = gadd(gel(p,step),gel(p,3*step));
860 : 226380 : f02= gsub(gel(p,step),gel(p,3*step));
861 : 226380 : f02 = mulcxI(f02);
862 : 226380 : gel(f,0) = gadd(f1, f3);
863 : 226380 : gel(f,1) = gadd(f2, f02);
864 : 226380 : gel(f,2) = gsub(f1, f3);
865 : 226380 : gel(f,3) = gsub(f2, f02); return;
866 : : }
867 : :
868 : 211784 : ltop = avma;
869 : 211784 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
870 : 211784 : fft(Omega,p, f, step4,l1);
871 : 211784 : fft(Omega,p+step, f+l1,step4,l1);
872 : 211784 : fft(Omega,p+(step<<1),f+l2,step4,l1);
873 : 211784 : fft(Omega,p+3*step, f+l3,step4,l1);
874 : :
875 : 211784 : ff = cgetg(l+1,t_VEC);
876 [ + + ]: 848356 : for (i=0; i<l1; i++)
877 : : {
878 : 636572 : rapi = step*i;
879 : 636572 : f1 = gmul(gel(Omega,rapi), gel(f,i+l1));
880 : 636572 : f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
881 : 636572 : f3 = gmul(gel(Omega,3*rapi), gel(f,i+l3));
882 : :
883 : 636572 : f02 = gadd(gel(f,i),f2);
884 : 636572 : g02 = gsub(gel(f,i),f2);
885 : 636572 : f13 = gadd(f1,f3);
886 : 636572 : g13 = mulcxI(gsub(f1,f3));
887 : :
888 : 636572 : gel(ff,i+1) = gadd(f02, f13);
889 : 636572 : gel(ff,i+l1+1) = gadd(g02, g13);
890 : 636572 : gel(ff,i+l2+1) = gsub(f02, f13);
891 : 636572 : gel(ff,i+l3+1) = gsub(g02, g13);
892 : : }
893 : 211784 : ff = gerepilecopy(ltop,ff);
894 [ + + ]: 3582180 : for (i=0; i<l; i++) f[i] = ff[i+1];
895 : : }
896 : :
897 : : /* e(1/N) */
898 : : static GEN
899 : 66985 : RUgen(long N, long bit)
900 : : {
901 [ + + ]: 66985 : if (N == 2) return real_m1(nbits2prec(bit));
902 [ + + ]: 64905 : if (N == 4) return gen_I();
903 : 66985 : return expIr(divru(Pi2n(1, nbits2prec(bit)), N));
904 : : }
905 : :
906 : : /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
907 : : static GEN
908 : 32033 : initRU(long N, long bit)
909 : : {
910 : 32033 : GEN *RU, z = RUgen(N, bit);
911 : 32033 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
912 : :
913 : 32033 : RU = (GEN*)cgetg(N+1,t_VEC); RU++;
914 : :
915 : 32033 : RU[0] = myreal_1(bit);
916 : 32033 : RU[1] = z;
917 [ + + ]: 44161 : for (i=1; i<N8; i++)
918 : : {
919 : 12128 : GEN t = RU[i];
920 : 12128 : RU[i+1] = gmul(z, t);
921 : 12128 : RU[N4-i] = mkcomplex(gel(t,2), gel(t,1));
922 : : }
923 [ + + ]: 112247 : for (i=0; i<N4; i++) RU[i+N4] = mulcxI(RU[i]);
924 [ + + ]: 192461 : for (i=0; i<N2; i++) RU[i+N2] = gneg(RU[i]);
925 : 32033 : return (GEN)RU;
926 : : }
927 : :
928 : : /* as above, N arbitrary */
929 : : static GEN
930 : 2919 : initRUgen(long N, long bit)
931 : : {
932 : 2919 : GEN *RU = (GEN*)cgetg(N+1,t_VEC), z = RUgen(N,bit);
933 : 2919 : long i, k = (N+3)>>1;
934 : 2919 : RU[0] = gen_1;
935 : 2919 : RU[1] = z;
936 [ + + ]: 4456 : for (i=2; i<k; i++) RU[i] = gmul(z, RU[i-1]);
937 [ + + ]: 4013 : for ( ; i<N; i++) RU[i] = gconj(RU[N-i]);
938 : 2919 : return (GEN)RU;
939 : : }
940 : :
941 : : GEN
942 : 0 : FFTinit(long k, long prec)
943 : : {
944 [ # # ]: 0 : if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
945 : 0 : return initRU(1L << k, prec2nbits(prec)) - 1;
946 : : }
947 : :
948 : : GEN
949 : 0 : FFT(GEN x, GEN Omega)
950 : : {
951 : 0 : long i, l = lg(Omega), n = lg(x);
952 : : GEN y, z;
953 [ # # ]: 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
954 [ # # ]: 0 : if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
955 [ # # ]: 0 : if (n > l) pari_err_DIM("FFT");
956 : :
957 [ # # ]: 0 : if (n < l) {
958 : 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
959 [ # # ]: 0 : for (i = 1; i < n; i++) z[i] = x[i];
960 [ # # ]: 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
961 : : }
962 : 0 : else z = x;
963 : 0 : y = cgetg(l, t_VEC);
964 : 0 : fft(Omega+1, z+1, y+1, 1, l-1);
965 : 0 : return y;
966 : : }
967 : :
968 : : /* returns 1 if p has only real coefficients, 0 else */
969 : : static int
970 : 29569 : isreal(GEN p)
971 : : {
972 : : long i;
973 [ + + ]: 191364 : for (i = lg(p)-1; i > 1; i--)
974 [ + + ]: 168768 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
975 : 29569 : return 1;
976 : : }
977 : :
978 : : /* x non complex */
979 : : static GEN
980 : 22242 : abs_update_r(GEN x, double *mu) {
981 : 22242 : GEN y = gtofp(x, DEFAULTPREC);
982 [ + + ]: 22242 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
983 : 22242 : setabssign(y); return y;
984 : : }
985 : : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
986 : : static GEN
987 : 443808 : abs_update(GEN x, double *mu) {
988 : : GEN y, xr, yr;
989 : : double ly;
990 [ + + ]: 443808 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
991 : 422100 : xr = gel(x,1);
992 : 422100 : yr = gel(x,2);
993 [ + + ]: 422100 : if (gequal0(xr)) return abs_update_r(yr,mu);
994 [ + + ]: 422050 : if (gequal0(yr)) return abs_update_r(xr,mu);
995 : : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
996 : 421566 : xr = gtofp(xr, DEFAULTPREC);
997 : 421566 : yr = gtofp(yr, DEFAULTPREC);
998 : 421566 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
999 [ + + ]: 421566 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
1000 : 443808 : return y;
1001 : : }
1002 : :
1003 : : static void
1004 : 15462 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
1005 : : int polreal, double param, double param2)
1006 : : {
1007 : : GEN q, pc, Omega, A, RU, prim, g, ONE,TWO;
1008 : 15462 : long n = degpol(p), bit, NN, K, i, j, Lmax;
1009 : 15462 : pari_sp av2, av = avma, lim = stack_lim(av, 1);
1010 : :
1011 : 15462 : bit = gexpo(p) + (long)param2+8;
1012 [ + + ]: 31087 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
1013 [ + + ]: 15462 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
1014 [ + + ]: 15462 : K = NN/Lmax; if (K & 1) K++;
1015 : 15462 : NN = Lmax*K;
1016 [ + + ]: 15462 : if (polreal) K = K/2+1;
1017 : :
1018 : 15462 : Omega = initRU(Lmax,bit);
1019 : 15462 : prim = RUgen(NN, bit);
1020 : :
1021 : 15462 : q = mygprec(p,bit) + 2;
1022 : 15462 : A = cgetg(Lmax+1,t_VEC); A++;
1023 : 15462 : pc= cgetg(Lmax+1,t_VEC); pc++;
1024 [ + + ]: 123512 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
1025 [ + + ]: 55792 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
1026 : :
1027 : 15462 : *mu = pariINFINITY;
1028 : 15462 : g = real_0_bit(-bit);
1029 : 15462 : ONE = real_1(DEFAULTPREC);
1030 : 15462 : TWO = real2n(1, DEFAULTPREC);
1031 : 15462 : av2 = avma;
1032 : 15462 : RU = myreal_1(bit);
1033 [ + + ]: 58770 : for (i=0; i<K; i++)
1034 : : {
1035 [ + + ]: 43308 : if (i) {
1036 : 27846 : GEN z = RU;
1037 [ + + ]: 195788 : for (j=1; j<n; j++)
1038 : : {
1039 : 167942 : gel(pc,j) = gmul(gel(q,j),z);
1040 : 167942 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
1041 : : }
1042 : 27846 : gel(pc,n) = gmul(gel(q,n),z);
1043 : : }
1044 : :
1045 : 43308 : fft(Omega,pc,A,1,Lmax);
1046 [ + + ][ + + ]: 49790 : if (polreal && i>0 && i<K-1)
[ + + ]
1047 [ + + ]: 104466 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
1048 : : else
1049 [ + + ]: 382650 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(ONE, abs_update(gel(A,j),mu)));
1050 : 43308 : RU = gmul(RU, prim);
1051 [ - + ]: 43308 : if (low_stack(lim, stack_lim(av,1)))
1052 : : {
1053 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
1054 : 0 : gerepileall(av2,2, &g,&RU);
1055 : : }
1056 : : }
1057 : 15462 : *gamma = mydbllog2r(divru(g,NN));
1058 : 15462 : *LMAX = Lmax; avma = av;
1059 : 15462 : }
1060 : :
1061 : : /* NN is a multiple of Lmax */
1062 : : static void
1063 : 16571 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
1064 : : {
1065 : : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
1066 : 16571 : long n = degpol(p), i, j, K;
1067 : : pari_sp ltop;
1068 : :
1069 : 16571 : Omega = initRU(Lmax,bit);
1070 : 16571 : prim = RUgen(NN, bit);
1071 : 16571 : RU = cgetg(n+2,t_VEC); RU++;
1072 : :
1073 [ + + ]: 16571 : K = NN/Lmax; if (polreal) K = K/2+1;
1074 : 16571 : q = mygprec(p,bit);
1075 : 16571 : qd = RgX_deriv(q);
1076 : :
1077 : 16571 : A = cgetg(Lmax+1,t_VEC); A++;
1078 : 16571 : B = cgetg(Lmax+1,t_VEC); B++;
1079 : 16571 : C = cgetg(Lmax+1,t_VEC); C++;
1080 : 16571 : pc = cgetg(Lmax+1,t_VEC); pc++;
1081 : 16571 : pd = cgetg(Lmax+1,t_VEC); pd++;
1082 [ + + ]: 63362 : pc[0] = q[2]; for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
1083 [ + + ]: 79933 : pd[0] = qd[2]; for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
1084 : :
1085 : 16571 : ltop = avma;
1086 : 16571 : W = cgetg(k+1,t_VEC);
1087 : 16571 : U = cgetg(k+1,t_VEC);
1088 [ + + ]: 49767 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
1089 : :
1090 : 16571 : gel(RU,0) = gen_1;
1091 : 16571 : prim2 = myreal_1(bit);
1092 [ + + ]: 52933 : for (i=0; i<K; i++)
1093 : : {
1094 : 36362 : gel(RU,1) = prim2;
1095 [ + + ]: 265574 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
1096 : : /* RU[j] = prim^(ij)= prim2^j */
1097 : :
1098 [ + + ]: 265574 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
1099 : 36362 : fft(Omega,pd,A,1,Lmax);
1100 [ + + ]: 301936 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
1101 : 36362 : fft(Omega,pc,B,1,Lmax);
1102 [ + + ]: 450654 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
1103 [ + + ]: 450654 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
1104 : 36362 : fft(Omega,B,A,1,Lmax);
1105 : 36362 : fft(Omega,C,B,1,Lmax);
1106 : :
1107 [ + + ]: 36362 : if (polreal) /* p has real coefficients */
1108 : : {
1109 [ + + ][ + + ]: 29310 : if (i>0 && i<K-1)
1110 : : {
1111 [ + + ]: 14602 : for (j=1; j<=k; j++)
1112 : : {
1113 : 11858 : gel(W,j) = gadd(gel(W,j), gshift(real_i(gmul(gel(A,j+1),gel(RU,j+1))),1));
1114 : 11858 : gel(U,j) = gadd(gel(U,j), gshift(real_i(gmul(gel(B,j),gel(RU,j))),1));
1115 : : }
1116 : : }
1117 : : else
1118 : : {
1119 [ + + ]: 77610 : for (j=1; j<=k; j++)
1120 : : {
1121 : 51044 : gel(W,j) = gadd(gel(W,j), real_i(gmul(gel(A,j+1),gel(RU,j+1))));
1122 : 51044 : gel(U,j) = gadd(gel(U,j), real_i(gmul(gel(B,j),gel(RU,j))));
1123 : : }
1124 : : }
1125 : : }
1126 : : else
1127 : : {
1128 [ + + ]: 27987 : for (j=1; j<=k; j++)
1129 : : {
1130 : 18191 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
1131 : 18191 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
1132 : : }
1133 : : }
1134 : 36362 : prim2 = gmul(prim2,prim);
1135 : 36362 : gerepileall(ltop,3, &W,&U,&prim2);
1136 : : }
1137 : :
1138 [ + + ]: 49767 : for (i=1; i<=k; i++)
1139 : : {
1140 : 33196 : aux=gel(W,i);
1141 [ + + ]: 75989 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
1142 : 33196 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
1143 : : }
1144 [ + + ]: 49767 : for (i=0; i<k; i++)
1145 : : {
1146 : 33196 : aux=gel(U,k-i);
1147 [ + + ]: 75989 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
1148 : 33196 : gel(H,i+2) = gdivgs(aux,NN);
1149 : : }
1150 : 16571 : }
1151 : :
1152 : : #define NEWTON_MAX 10
1153 : : static GEN
1154 : 82634 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
1155 : : {
1156 : 82634 : GEN H = HH, D, aux;
1157 : 82634 : pari_sp ltop = avma, lim = stack_lim(ltop, 1);
1158 : : long error, i, bit1, bit2;
1159 : :
1160 : 82634 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
1161 : 82634 : bit2 = bit + Sbit;
1162 [ + + ][ + + ]: 153861 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1163 : : {
1164 [ - + ]: 71227 : if (low_stack(lim, stack_lim(ltop,1)))
1165 : : {
1166 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
1167 : 0 : gerepileall(ltop,2, &D,&H);
1168 : : }
1169 : 71227 : bit1 = -error + Sbit;
1170 : 71227 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
1171 : 71227 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
1172 : :
1173 [ + + ]: 71227 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
1174 : 71227 : H = RgX_add(mygprec(H,bit1), aux);
1175 : 71227 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
1176 [ + + ]: 71227 : error = gexpo(D); if (error < -bit1) error = -bit1;
1177 : : }
1178 [ + + ]: 82634 : if (error > -bit/2) return NULL; /* FAIL */
1179 : 82634 : return gerepilecopy(ltop,H);
1180 : : }
1181 : :
1182 : : /* return 0 if fails, 1 else */
1183 : : static long
1184 : 16571 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
1185 : : {
1186 : 16571 : GEN f0, FF, GG, r, HH = H;
1187 : 16571 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
1188 : 16571 : pari_sp av = avma, lim = stack_lim(av, 1);
1189 : :
1190 : 16571 : FF = *F; GG = RgX_divrem(p, FF, &r);
1191 [ - + ]: 16571 : error = gexpo(r); if (error <= -bit) error = 1-bit;
1192 : 16571 : normF = gexpo(FF);
1193 : 16571 : normG = gexpo(GG);
1194 [ + + ]: 16571 : enh = gexpo(H); if (enh < 0) enh = 0;
1195 : 16571 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
1196 : 16571 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
1197 : 16571 : bit2 = bit + Sbit;
1198 [ + + ][ + - ]: 99169 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1199 : : {
1200 [ + + ][ + + ]: 82634 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
1201 [ - + ]: 82634 : if (low_stack(lim, stack_lim(av,1)))
1202 : : {
1203 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1204 : 0 : gerepileall(av,4, &FF,&GG,&r,&HH);
1205 : : }
1206 : :
1207 : 82634 : bit1 = -error + Sbit2;
1208 : 82634 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1209 : : 1-error, Sbit2);
1210 [ + + ]: 82634 : if (!HH) return 0; /* FAIL */
1211 : :
1212 : 82598 : bit1 = -error + Sbit;
1213 : 82598 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1214 : 82598 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1215 : :
1216 [ + + ]: 82598 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1217 : 82598 : FF = gadd(mygprec(FF,bit1),f0);
1218 : :
1219 [ + + ]: 82598 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1220 : 82598 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1221 [ + + ]: 82598 : error = gexpo(r); if (error < -bit1) error = -bit1;
1222 : : }
1223 [ + + ]: 16535 : if (error>-bit) return 0; /* FAIL */
1224 : 16571 : *F = FF; *G = GG; return 1;
1225 : : }
1226 : :
1227 : : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1228 : : where cd is the leading coefficient of p */
1229 : : static void
1230 : 15462 : split_fromU(GEN p, long k, double delta, long bit,
1231 : : GEN *F, GEN *G, double param, double param2)
1232 : : {
1233 : : GEN pp, FF, GG, H;
1234 : 15462 : long n = degpol(p), NN, bit2, Lmax;
1235 : 15462 : int polreal = isreal(p);
1236 : : pari_sp ltop;
1237 : : double mu, gamma;
1238 : :
1239 : 15462 : pp = gdiv(p, gel(p,2+n));
1240 : 15462 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1241 : :
1242 : 15462 : H = cgetg(k+2,t_POL); H[1] = p[1];
1243 : 15462 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1244 : 15462 : gel(FF,k+2) = gen_1;
1245 : :
1246 [ + + ]: 15462 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1247 : 15462 : NN *= Lmax; ltop = avma;
1248 : : for(;;)
1249 : : {
1250 : 16571 : bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
1251 : 16571 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1252 [ + + ]: 16571 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1253 : 1109 : NN <<= 1; avma = ltop;
1254 : 1109 : }
1255 : 15462 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1256 : 15462 : }
1257 : :
1258 : : static void
1259 : 15462 : optimize_split(GEN p, long k, double delta, long bit,
1260 : : GEN *F, GEN *G, double param, double param2)
1261 : : {
1262 : 15462 : long n = degpol(p);
1263 : : GEN FF, GG;
1264 : :
1265 [ + + ]: 15462 : if (k <= n/2)
1266 : 11866 : split_fromU(p,k,delta,bit,F,G,param,param2);
1267 : : else
1268 : : {
1269 : 3596 : split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1270 : 3596 : *F = RgX_recip_shallow(GG);
1271 : 3596 : *G = RgX_recip_shallow(FF);
1272 : : }
1273 : 15462 : }
1274 : :
1275 : : /********************************************************************/
1276 : : /** **/
1277 : : /** SEARCH FOR SEPARATING CIRCLE **/
1278 : : /** **/
1279 : : /********************************************************************/
1280 : :
1281 : : /* return p(2^e*x) *2^(-n*e) */
1282 : : static void
1283 : 0 : scalepol2n(GEN p, long e)
1284 : : {
1285 : 0 : long i,n=lg(p)-1;
1286 [ # # ]: 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1287 : 0 : }
1288 : :
1289 : : /* returns p(x/R)*R^n */
1290 : : static GEN
1291 : 131028 : scalepol(GEN p, GEN R, long bit)
1292 : : {
1293 : : GEN q,aux,gR;
1294 : : long i;
1295 : :
1296 : 131028 : aux = gR = mygprec(R,bit); q = mygprec(p,bit);
1297 [ + + ]: 628740 : for (i=lg(p)-2; i>=2; i--)
1298 : : {
1299 : 497712 : gel(q,i) = gmul(aux,gel(q,i));
1300 : 497712 : aux = gmul(aux,gR);
1301 : : }
1302 : 131028 : return q;
1303 : : }
1304 : :
1305 : : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1306 : : static GEN
1307 : 42321 : conformal_pol(GEN p, GEN a, long bit)
1308 : : {
1309 : 42321 : GEN z, r, ma = gneg(a), ca = gconj(a);
1310 : 42321 : long n = degpol(p), i;
1311 : 42321 : pari_sp av = avma, lim = stack_lim(av,2);
1312 : :
1313 : 42321 : z = mkpoln(2, ca, negr(myreal_1(bit)));
1314 : 42321 : r = scalarpol(gel(p,2+n), 0);
1315 : 42321 : for (i=n-1; ; i--)
1316 : : {
1317 : 156268 : r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
1318 : 156268 : r = gadd(r, gmul(z, gel(p,2+i)));
1319 [ + + ]: 156268 : if (i == 0) return gerepileupto(av, r);
1320 : 113947 : z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1321 [ - + ]: 113947 : if (low_stack(lim, stack_lim(av,2)))
1322 : : {
1323 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
1324 : 0 : gerepileall(av,2, &r,&z);
1325 : : }
1326 : 113947 : }
1327 : : }
1328 : :
1329 : : static const double UNDEF = -100000.;
1330 : :
1331 : : static double
1332 : 15462 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1333 : : {
1334 : 15462 : long i, n = degpol(p);
1335 : : double lrho, lrmin, lrmax;
1336 [ + + ]: 15462 : if (k > 1)
1337 : : {
1338 [ + - ][ + + ]: 16169 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1339 : 10229 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1340 : : }
1341 : : else /* k=1 */
1342 : 5233 : lrmin = logmin_modulus(p,aux);
1343 : 15462 : radii[k] = lrmin;
1344 : :
1345 [ + + ]: 15462 : if (k+1<n)
1346 : : {
1347 [ + - ][ + + ]: 30231 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1348 : 13256 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1349 : : }
1350 : : else /* k+1=n */
1351 : 2206 : lrmax = logmax_modulus(p,aux);
1352 : 15462 : radii[k+1] = lrmax;
1353 : :
1354 : 15462 : lrho = radii[k];
1355 [ + + ]: 37013 : for (i=k-1; i>=1; i--)
1356 : : {
1357 [ + + ][ + + ]: 21551 : if (radii[i] == UNDEF || radii[i] > lrho)
1358 : 15044 : radii[i] = lrho;
1359 : : else
1360 : 6507 : lrho = radii[i];
1361 : : }
1362 : 15462 : lrho = radii[k+1];
1363 [ + + ]: 71037 : for (i=k+1; i<=n; i++)
1364 : : {
1365 [ + + ][ + + ]: 55575 : if (radii[i] == UNDEF || radii[i] < lrho)
1366 : 34311 : radii[i] = lrho;
1367 : : else
1368 : 21264 : lrho = radii[i];
1369 : : }
1370 : 15462 : *delta = (lrmax - lrmin) / 2;
1371 [ + + ]: 15462 : if (*delta > 1.) *delta = 1.;
1372 : 15462 : return (lrmin + lrmax) / 2;
1373 : : }
1374 : :
1375 : : static void
1376 : 15462 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1377 : : {
1378 : 15462 : double t, param = 0., param2 = 0.;
1379 : : long i;
1380 [ + + ]: 108050 : for (i=1; i<=n; i++)
1381 : : {
1382 : 92588 : radii[i] -= lrho;
1383 : 92588 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1384 [ + + ]: 92588 : param += t; if (t > 1.) param2 += log2(t);
1385 : : }
1386 : 15462 : *par = param; *par2 = param2;
1387 : 15462 : }
1388 : :
1389 : : /* apply the conformal mapping then split from U */
1390 : : static void
1391 : 14107 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1392 : : double aux, GEN *F,GEN *G)
1393 : : {
1394 : 14107 : long bit2, n = degpol(p), i;
1395 : 14107 : pari_sp ltop = avma, av;
1396 : : GEN q, FF, GG, a, R;
1397 : : double lrho, delta, param, param2;
1398 : : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1399 : 14107 : bit2 = bit + (long)(n*3.4848775) + 1;
1400 : 14107 : a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
1401 : 14107 : a = divrs(a, -6);
1402 : 14107 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1403 : :
1404 : 14107 : av = avma;
1405 : 14107 : q = conformal_pol(mygprec(p,bit2), a, bit2);
1406 [ + + ]: 92241 : for (i=1; i<=n; i++)
1407 [ + + ]: 78134 : if (radii[i] != UNDEF) /* update array radii */
1408 : : {
1409 : 52717 : pari_sp av2 = avma;
1410 : 52717 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1411 : : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1412 : 52717 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1413 : 52717 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1414 : 52717 : avma = av2;
1415 : : }
1416 : 14107 : lrho = logradius(radii, q,k,aux/10., &delta);
1417 : 14107 : update_radius(n, radii, lrho, ¶m, ¶m2);
1418 : :
1419 : 14107 : bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
1420 : 14107 : R = mygprec(dblexp(-lrho), bit2);
1421 : 14107 : q = scalepol(q,R,bit2);
1422 : 14107 : gerepileall(av,2, &q,&R);
1423 : :
1424 : 14107 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1425 : 14107 : bit2 += n; R = invr(R);
1426 : 14107 : FF = scalepol(FF,R,bit2);
1427 : 14107 : GG = scalepol(GG,R,bit2);
1428 : :
1429 : 14107 : a = mygprec(a,bit2);
1430 : 14107 : FF = conformal_pol(FF,a,bit2);
1431 : 14107 : GG = conformal_pol(GG,a,bit2);
1432 : :
1433 : 14107 : a = invr(subsr(1, gnorm(a)));
1434 : 14107 : FF = RgX_Rg_mul(FF, powru(a,k));
1435 : 14107 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1436 : :
1437 : 14107 : *F = mygprec(FF,bit+n);
1438 : 14107 : *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1439 : 14107 : }
1440 : :
1441 : : /* split p, this time without scaling. returns in F and G two polynomials
1442 : : * such that |p-FG|< 2^(-bit)|p| */
1443 : : static void
1444 : 15462 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1445 : : {
1446 : : GEN q, FF, GG, R;
1447 : : double aux, delta, param, param2;
1448 : 15462 : long n = degpol(p), i, j, k, bit2;
1449 : : double lrmin, lrmax, lrho, *radii;
1450 : :
1451 : 15462 : init_dalloc();
1452 : 15462 : radii = (double*) stack_malloc((n+1) * sizeof(double));
1453 : :
1454 [ + + ]: 77126 : for (i=2; i<n; i++) radii[i] = UNDEF;
1455 : 15462 : aux = thickness/(double)(4 * n);
1456 : 15462 : lrmin = logmin_modulus(p, aux);
1457 : 15462 : lrmax = logmax_modulus(p, aux);
1458 : 15462 : radii[1] = lrmin;
1459 : 15462 : radii[n] = lrmax;
1460 : 15462 : i = 1; j = n;
1461 : 15462 : lrho = (lrmin + lrmax) / 2;
1462 : 15462 : k = dual_modulus(p, lrho, aux, 1);
1463 [ + + ][ + + ]: 15462 : if (5*k < n || (n < 2*k && 5*k < 4*n))
[ + + ]
1464 : 3416 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1465 : : else
1466 : 12046 : { lrmin = lrho; i=k; radii[i] = lrho; }
1467 [ + + ]: 33644 : while (j > i+1)
1468 : : {
1469 [ + + ]: 18182 : if (i+j == n+1)
1470 : 6638 : lrho = (lrmin + lrmax) / 2;
1471 : : else
1472 : : {
1473 : 11544 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1474 [ + + ]: 11544 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1475 : 8419 : else lrho = lrmin * kappa + lrmax;
1476 : 11544 : lrho /= 1+kappa;
1477 : : }
1478 : 18182 : aux = (lrmax - lrmin) / (4*(j-i));
1479 : 18182 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1480 [ + + ][ + + ]: 18182 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
[ + + ]
1481 : 12668 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1482 : : else
1483 : 5514 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1484 : : }
1485 : 15462 : aux = lrmax - lrmin;
1486 : :
1487 [ + + ]: 15462 : if (ctr)
1488 : : {
1489 : 14107 : lrho = (lrmax + lrmin) / 2;
1490 [ + + ]: 92241 : for (i=1; i<=n; i++)
1491 [ + + ]: 78134 : if (radii[i] != UNDEF) radii[i] -= lrho;
1492 : :
1493 : 14107 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1494 : 14107 : R = mygprec(dblexp(-lrho), bit2);
1495 : 14107 : q = scalepol(p,R,bit2);
1496 : 14107 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1497 : : }
1498 : : else
1499 : : {
1500 : 1355 : lrho = logradius(radii, p, k, aux/10., &delta);
1501 : 1355 : update_radius(n, radii, lrho, ¶m, ¶m2);
1502 : :
1503 : 1355 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1504 : 1355 : R = mygprec(dblexp(-lrho), bit2);
1505 : 1355 : q = scalepol(p,R,bit2);
1506 : 1355 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1507 : : }
1508 : 15462 : bit += n;
1509 : 15462 : bit2 += n; R = invr(mygprec(R,bit2));
1510 : 15462 : *F = mygprec(scalepol(FF,R,bit2), bit);
1511 : 15462 : *G = mygprec(scalepol(GG,R,bit2), bit);
1512 : 15462 : }
1513 : :
1514 : : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1515 : : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1516 : : * where the maximum modulus of the roots of p is <=1.
1517 : : * Assume sum of roots is 0. */
1518 : : static void
1519 : 14107 : split_1(GEN p, long bit, GEN *F, GEN *G)
1520 : : {
1521 : 14107 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1522 : : GEN TWO, ctr, q, qq, FF, GG, v, gr, r, newq;
1523 : : double lrmin, lrmax, lthick;
1524 : 14107 : const double LOG3 = 1.098613;
1525 : :
1526 : 14107 : lrmax = logmax_modulus(p, 0.01);
1527 : 14107 : gr = mygprec(dblexp(-lrmax), bit2);
1528 : 14107 : q = scalepol(p,gr,bit2);
1529 : :
1530 : 14107 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1531 : 14107 : TWO = myreal_1(bit2); setexpo(TWO,1);
1532 : 14107 : v = cgetg(5,t_VEC);
1533 : 14107 : gel(v,1) = TWO;
1534 : 14107 : gel(v,2) = negr(TWO);
1535 : 14107 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1536 : 14107 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1537 : 14107 : q = mygprec(q,bit2); lthick = 0;
1538 : 14107 : newq = ctr = NULL; /* -Wall */
1539 [ + + ]: 14107 : imax = polreal? 3: 4;
1540 [ + + ]: 26249 : for (i=1; i<=imax; i++)
1541 : : {
1542 : 25898 : qq = RgX_translate(q, gel(v,i));
1543 : 25898 : lrmin = logmin_modulus(qq,0.05);
1544 [ + + ]: 25898 : if (LOG3 > lrmin + lthick)
1545 : : {
1546 : 25279 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1547 [ + + ]: 25279 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1548 : : }
1549 [ + + ]: 25898 : if (lthick > LOG2) break;
1550 [ + + ][ + + ]: 14273 : if (polreal && i==2 && lthick > LOG3 - LOG2) break;
[ + + ]
1551 : : }
1552 : 14107 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
1553 : 14107 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1554 : 14107 : r = gneg(mygprec(ctr,bit2));
1555 : 14107 : FF = RgX_translate(FF,r);
1556 : 14107 : GG = RgX_translate(GG,r);
1557 : :
1558 : 14107 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1559 : 14107 : *F = scalepol(FF,gr,bit2);
1560 : 14107 : *G = scalepol(GG,gr,bit2);
1561 : 14107 : }
1562 : :
1563 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1564 : : where the maximum modulus of the roots of p is < 0.5 */
1565 : : static int
1566 : 14107 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1567 : : {
1568 : : GEN q, b, FF, GG;
1569 : 14107 : long n = degpol(p), k, bit2, eq;
1570 : 14107 : double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
1571 : :
1572 : : /* beware double overflow */
1573 [ + + ][ + - ]: 14107 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
[ - + ]
1574 : :
1575 [ + + ]: 14107 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1576 : 14107 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1577 : :
1578 : 14107 : q = mygprec(p,bit2);
1579 : 14107 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1580 : 14107 : q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
1581 : 14107 : k = 0;
1582 [ + - ][ - + ]: 14107 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1583 [ - + ]: 14107 : || gequal0(gel(q,k+2)))) k++;
1584 [ - + ]: 14107 : if (k > 0)
1585 : : {
1586 [ # # ]: 0 : if (k > n/2) k = n/2;
1587 : 0 : bit2 += k<<1;
1588 : 0 : FF = monomial(myreal_1(bit2), k, 0);
1589 : 0 : GG = RgX_shift_shallow(q, -k);
1590 : : }
1591 : : else
1592 : : {
1593 : 14107 : split_1(q,bit2,&FF,&GG);
1594 : 14107 : bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
1595 : 14107 : FF = mygprec(FF,bit2);
1596 : : }
1597 : 14107 : GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
1598 : 14107 : *F = RgX_translate(FF, b);
1599 : 14107 : *G = RgX_translate(GG, b); return 1;
1600 : : }
1601 : :
1602 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1603 : : * Assume max_modulus(p) < 2 */
1604 : : static void
1605 : 14107 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1606 : : {
1607 : : GEN FF, GG;
1608 : : long n, bit2, normp;
1609 : :
1610 [ - + ]: 28214 : if (split_0_2(p,bit,F,G)) return;
1611 : :
1612 : 0 : normp = gexpo(p);
1613 : 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1614 : 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1615 : 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1616 : 0 : scalepol2n(FF,-2);
1617 : 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1618 : 0 : *F = mygprec(FF,bit2);
1619 : 14107 : *G = mygprec(GG,bit2);
1620 : : }
1621 : :
1622 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1623 : : static void
1624 : 15462 : split_0(GEN p, long bit, GEN *F, GEN *G)
1625 : : {
1626 : 15462 : const double LOG1_9 = 0.6418539;
1627 : 15462 : long n = degpol(p), k = 0;
1628 : : GEN q;
1629 : :
1630 [ - + ][ # # ]: 15462 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1631 [ - + ]: 15462 : if (k > 0)
1632 : : {
1633 [ # # ]: 0 : if (k > n/2) k = n/2;
1634 : 0 : *F = monomial(myreal_1(bit), k, 0);
1635 : 0 : *G = RgX_shift_shallow(p, -k);
1636 : : }
1637 : : else
1638 : : {
1639 : 15462 : double lr = logmax_modulus(p, 0.05);
1640 [ + + ]: 15462 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1641 : : else
1642 : : {
1643 : 12509 : q = RgX_recip_shallow(p);
1644 : 12509 : lr = logmax_modulus(q,0.05);
1645 [ + + ]: 12509 : if (lr < LOG1_9)
1646 : : {
1647 : 11154 : split_0_1(q, bit, F, G);
1648 : 11154 : *F = RgX_recip_shallow(*F);
1649 : 11154 : *G = RgX_recip_shallow(*G);
1650 : : }
1651 : : else
1652 : 1355 : split_2(p,bit,NULL, 1.2837,F,G);
1653 : : }
1654 : : }
1655 : 15462 : }
1656 : :
1657 : : /********************************************************************/
1658 : : /** **/
1659 : : /** ERROR ESTIMATE FOR THE ROOTS **/
1660 : : /** **/
1661 : : /********************************************************************/
1662 : :
1663 : : static GEN
1664 : 48051 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1665 : : {
1666 : 48051 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1667 : : long i, j;
1668 : :
1669 : 48051 : d = cgetg(n+1,t_VEC);
1670 [ + + ]: 777094 : for (i=1; i<=n; i++)
1671 : : {
1672 [ + + ]: 729043 : if (i!=k)
1673 : : {
1674 : 680992 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1675 : 680992 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1676 : : }
1677 : : }
1678 : 48051 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1679 [ + + ]: 48051 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1680 : 48051 : eps = mulrr(rho, shatzle);
1681 : 48051 : aux = shiftr(powru(rho,n), err);
1682 : :
1683 [ + + ][ + - ]: 153770 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
[ + + ]
1684 : : {
1685 : 105719 : GEN prod = NULL; /* 1. */
1686 : 105719 : long m = n;
1687 : 105719 : epsbis = mulrr(eps, dbltor(1.25));
1688 [ + + ]: 1989589 : for (i=1; i<=n; i++)
1689 : : {
1690 [ + + ][ + + ]: 1883870 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1691 : : {
1692 : 1760057 : GEN dif = subrr(gel(d,i),eps);
1693 [ + + ]: 1760057 : prod = prod? mulrr(prod, dif): dif;
1694 : 1760057 : m--;
1695 : : }
1696 : : }
1697 [ + + ]: 105719 : eps2 = prod? divrr(aux, prod): aux;
1698 [ + + ]: 105719 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1699 : 105719 : rap = divrr(eps,eps2); eps = eps2;
1700 : : }
1701 : 48051 : return eps;
1702 : : }
1703 : :
1704 : : /* round a complex or real number x to an absolute value of 2^(-bit) */
1705 : : static GEN
1706 : 111076 : mygprec_absolute(GEN x, long bit)
1707 : : {
1708 : : long e;
1709 : : GEN y;
1710 : :
1711 [ + + + ]: 111076 : switch(typ(x))
1712 : : {
1713 : : case t_REAL:
1714 : 75963 : e = expo(x) + bit;
1715 [ + + ][ - + ]: 75963 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1716 : : case t_COMPLEX:
1717 [ + + ]: 31947 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1718 : 31078 : y = cgetg(3,t_COMPLEX);
1719 : 31078 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1720 : 31078 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1721 : 31078 : return y;
1722 : 111076 : default: return x;
1723 : : }
1724 : : }
1725 : :
1726 : : static long
1727 : 6839 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1728 : : {
1729 : 6839 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1730 : : GEN sigma, shatzle;
1731 : :
1732 : 6839 : err += (long)log2((double)n) + 1;
1733 [ - + ]: 6839 : if (err > -2) return 0;
1734 : 6839 : sigma = real2n(-err, LOWDEFAULTPREC);
1735 : : /* 2 / ((s - 1)^(1/n) - 1) */
1736 : 6839 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1737 [ + + ]: 54890 : for (i=1; i<=n; i++)
1738 : : {
1739 : 48051 : pari_sp av = avma;
1740 : 48051 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1741 : 48051 : long e = gexpo(x);
1742 [ + + ]: 48051 : avma = av; if (e > e_max) e_max = e;
1743 : 48051 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1744 : : }
1745 : 6839 : return e_max;
1746 : : }
1747 : :
1748 : : /********************************************************************/
1749 : : /** **/
1750 : : /** MAIN **/
1751 : : /** **/
1752 : : /********************************************************************/
1753 : : static GEN
1754 : 35625 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1755 : :
1756 : : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1757 : : * returns prod (x-roots_pol[i]) */
1758 : : static GEN
1759 : 37763 : split_complete(GEN p, long bit, GEN roots_pol)
1760 : : {
1761 : 37763 : long n = degpol(p);
1762 : : pari_sp ltop;
1763 : : GEN p1, F, G, a, b, m1, m2;
1764 : :
1765 [ + + ]: 37763 : if (n == 1)
1766 : : {
1767 : 8977 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1768 : 8977 : (void)append_clone(roots_pol,a); return p;
1769 : : }
1770 : 28786 : ltop = avma;
1771 [ + + ]: 28786 : if (n == 2)
1772 : : {
1773 : 13324 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1774 : 13324 : F = gsqrt(F, nbits2prec(bit));
1775 : 13324 : p1 = ginv(gmul2n(gel(p,4),1));
1776 : 13324 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1777 : 13324 : b = gmul(gsub(F,gel(p,3)), p1);
1778 : 13324 : a = append_clone(roots_pol,a);
1779 : 13324 : b = append_clone(roots_pol,b); avma = ltop;
1780 : 13324 : a = mygprec(a, 3*bit);
1781 : 13324 : b = mygprec(b, 3*bit);
1782 : 13324 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1783 : : }
1784 : 15462 : split_0(p,bit,&F,&G);
1785 : 15462 : m1 = split_complete(F,bit,roots_pol);
1786 : 15462 : m2 = split_complete(G,bit,roots_pol);
1787 : 37763 : return gerepileupto(ltop, gmul(m1,m2));
1788 : : }
1789 : :
1790 : : static GEN
1791 : 221467 : quicktofp(GEN x)
1792 : : {
1793 : 221467 : const long prec = DEFAULTPREC;
1794 [ + + - + : 221467 : switch(typ(x))
- ]
1795 : : {
1796 : 220989 : case t_INT: return itor(x, prec);
1797 : 463 : case t_REAL: return rtor(x, prec);
1798 : 0 : case t_FRAC: return fractor(x, prec);
1799 : : case t_COMPLEX: {
1800 : 15 : GEN a = gel(x,1), b = gel(x,2);
1801 : : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1802 [ - + ]: 15 : if (isintzero(a)) return cxcompotor(b, prec);
1803 [ - + ]: 15 : if (isintzero(b)) return cxcompotor(a, prec);
1804 : 15 : a = cxcompotor(a, prec);
1805 : 15 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1806 : : }
1807 : 0 : default: pari_err_TYPE("quicktofp",x);
1808 : 221467 : return NULL;/*not reached*/
1809 : : }
1810 : :
1811 : : }
1812 : :
1813 : : /* bound log_2 |largest root of p| (Fujiwara's bound) */
1814 : : double
1815 : 46977 : fujiwara_bound(GEN p)
1816 : : {
1817 : 46977 : pari_sp av = avma;
1818 : 46977 : long i, n = degpol(p);
1819 : : GEN cc;
1820 : : double loglc, Lmax;
1821 : :
1822 [ - + ]: 46977 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1823 : 46977 : loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
1824 : 46977 : cc = gel(p, 2);
1825 [ + + ]: 46977 : if (gequal0(cc))
1826 : 2719 : Lmax = -pariINFINITY-1;
1827 : : else
1828 : 44258 : Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
1829 [ + + ]: 204975 : for (i = 1; i < n; i++)
1830 : : {
1831 : 157998 : GEN y = gel(p,i+2);
1832 : : double L;
1833 [ + + ]: 157998 : if (gequal0(y)) continue;
1834 : 130232 : L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
1835 [ + + ]: 130232 : if (L > Lmax) Lmax = L;
1836 : : }
1837 : 46977 : avma = av; return Lmax + 1;
1838 : : }
1839 : :
1840 : : /* Fujiwara's bound, real roots. Based on the following remark: if
1841 : : * p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
1842 : : * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
1843 : : * of q is a bound for the positive roots of p. */
1844 : : double
1845 : 4800 : fujiwara_bound_real(GEN p, long sign)
1846 : : {
1847 : 4800 : pari_sp av = avma;
1848 : : GEN x;
1849 : 4800 : long n = degpol(p), i, signodd, signeven;
1850 : : double fb;
1851 [ - + ]: 4800 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1852 : 4800 : x = shallowcopy(p);
1853 [ + - ]: 4800 : if (gsigne(gel(x, n+2)) > 0)
1854 : : {
1855 : 4800 : signeven = 1;
1856 : 4800 : signodd = sign;
1857 : : }
1858 : : else
1859 : : {
1860 : 0 : signeven = -1;
1861 : 0 : signodd = -sign;
1862 : : }
1863 [ + + ]: 28498 : for (i = 0; i < n; i++)
1864 : : {
1865 [ + + ]: 23698 : if ((n - i) % 2)
1866 : : {
1867 [ + + ]: 12768 : if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0;
1868 : : }
1869 : : else
1870 : : {
1871 [ + + ]: 10930 : if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0;
1872 : : }
1873 : : }
1874 : 4800 : fb = fujiwara_bound(x);
1875 : 4800 : avma = av; return fb;
1876 : : }
1877 : :
1878 : : static GEN
1879 : 54920 : mygprecrc_special(GEN x, long prec, long e)
1880 : : {
1881 : : GEN y;
1882 [ + + + ]: 54920 : switch(typ(x))
1883 : : {
1884 : : case t_REAL:
1885 [ - + ]: 490 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1886 [ + - ]: 490 : return (prec > realprec(x))? rtor(x, prec): x;
1887 : : case t_COMPLEX:
1888 : 15 : y = cgetg(3,t_COMPLEX);
1889 : 15 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1890 : 15 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1891 : 15 : return y;
1892 : 54920 : default: return x;
1893 : : }
1894 : : }
1895 : :
1896 : : /* like mygprec but keep at least the same precision as before */
1897 : : static GEN
1898 : 6839 : mygprec_special(GEN x, long bit)
1899 : : {
1900 : : long lx, i, e, prec;
1901 : : GEN y;
1902 : :
1903 [ - + ]: 6839 : if (bit < 0) bit = 0; /* should not happen */
1904 : 6839 : e = gexpo(x) - bit;
1905 : 6839 : prec = nbits2prec(bit);
1906 [ + - ]: 6839 : switch(typ(x))
1907 : : {
1908 : : case t_POL:
1909 : 6839 : y = cgetg_copy(x, &lx); y[1] = x[1];
1910 [ + + ]: 61729 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1911 : 6839 : break;
1912 : :
1913 : 0 : default: y = mygprecrc_special(x,prec,e);
1914 : : }
1915 : 6839 : return y;
1916 : : }
1917 : :
1918 : : static GEN
1919 : 3920 : fix_roots1(GEN r)
1920 : : {
1921 : 3920 : long i, l = lg(r);
1922 : 3920 : GEN allr = cgetg(l, t_VEC);
1923 [ + + ]: 31053 : for (i=1; i<l; i++)
1924 : : {
1925 : 27133 : GEN t = gel(r,i);
1926 : 27133 : gel(allr,i) = gcopy(t); gunclone(t);
1927 : : }
1928 : 3920 : return allr;
1929 : : }
1930 : : static GEN
1931 : 6839 : fix_roots(GEN r, GEN *m, long h, long bit)
1932 : : {
1933 : : long i, j, k, l, prec;
1934 : : GEN allr, ro1;
1935 [ + + ]: 6839 : if (h == 1) return fix_roots1(r);
1936 : 2919 : ro1 = initRUgen(h, bit);
1937 : 2919 : prec = nbits2prec(bit);
1938 : 2919 : l = lg(r)-1;
1939 : 2919 : allr = cgetg(h*l+1, t_VEC);
1940 [ + + ]: 11411 : for (k=1,i=1; i<=l; i++)
1941 : : {
1942 : 8492 : GEN p2, p1 = gel(r,i);
1943 [ + + ]: 8492 : p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1944 [ + + ]: 29410 : for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1945 : 8492 : gunclone(p1);
1946 : : }
1947 : 2919 : *m = roots_to_pol(allr, 0);
1948 : 6839 : return allr;
1949 : : }
1950 : :
1951 : : static GEN
1952 : 6721 : all_roots(GEN p, long bit)
1953 : : {
1954 : : GEN lc, pd, q, roots_pol, m;
1955 : 6721 : long bit0, bit2, n = degpol(p), i, e, h;
1956 : : pari_sp av;
1957 : :
1958 : 6721 : pd = RgX_deflate_max(p, &h); lc = leading_term(pd);
1959 [ + + ]: 6721 : e = (long)(2 * fujiwara_bound(pd)); if (e < 0) e = 0;
1960 : 6721 : bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
1961 : 6721 : bit2 = bit0; e = 0;
1962 : 6721 : for (av=avma,i=1;; i++,avma=av)
1963 : : {
1964 : 6839 : roots_pol = vectrunc_init(n+1);
1965 : 6839 : bit2 += e + (n << i);
1966 : 6839 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1967 : 6839 : q[1] = evalsigne(1)|evalvarn(0);
1968 : 6839 : m = split_complete(q,bit2,roots_pol);
1969 : 6839 : roots_pol = fix_roots(roots_pol, &m, h, bit2);
1970 : 6839 : q = mygprec_special(p,bit2); lc = leading_term(q);
1971 : 6839 : q[1] = evalsigne(1)|evalvarn(0);
1972 [ + + ]: 6839 : if (h > 1) m = gmul(m,lc);
1973 : :
1974 : 6839 : e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
1975 [ - + ]: 6839 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
1976 [ + - ]: 6839 : if (e < 0)
1977 : : {
1978 : 6839 : e = bit + a_posteriori_errors(p,roots_pol,e);
1979 [ + + ]: 6839 : if (e < 0) return roots_pol;
1980 : : }
1981 [ - + ]: 118 : if (DEBUGLEVEL > 7)
1982 : 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1983 : 118 : }
1984 : : }
1985 : :
1986 : : INLINE int
1987 : 1691 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1988 : :
1989 : : static int
1990 : 525 : isexactpol(GEN p)
1991 : : {
1992 : 525 : long i,n = degpol(p);
1993 [ + + ]: 2078 : for (i=0; i<=n; i++)
1994 [ + + ]: 1691 : if (!isexactscalar(gel(p,i+2))) return 0;
1995 : 525 : return 1;
1996 : : }
1997 : :
1998 : : static long
1999 : 2294 : isvalidcoeff(GEN x)
2000 : : {
2001 [ + + + ]: 2294 : switch (typ(x))
2002 : : {
2003 : 2269 : case t_INT: case t_REAL: case t_FRAC: return 1;
2004 [ + - ][ + - ]: 15 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
2005 : : }
2006 : 2294 : return 0;
2007 : : }
2008 : :
2009 : : static void
2010 : 530 : checkvalidpol(GEN p)
2011 : : {
2012 : 530 : long i,n = lg(p);
2013 [ + + ]: 2779 : for (i=2; i<n; i++)
2014 [ + + ]: 2254 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
2015 : 525 : }
2016 : :
2017 : : static GEN
2018 : 387 : solve_exact_pol(GEN p, long bit)
2019 : : {
2020 : 387 : long i, j, k, m, n = degpol(p), iroots = 0;
2021 : 387 : GEN ex, factors, v = zerovec(n);
2022 : :
2023 : 387 : factors = ZX_squff(Q_primpart(p), &ex);
2024 [ + + ]: 774 : for (i=1; i<lg(factors); i++)
2025 : : {
2026 : 387 : GEN roots_fact = all_roots(gel(factors,i), bit);
2027 : 387 : n = degpol(gel(factors,i));
2028 : 387 : m = ex[i];
2029 [ + + ]: 1538 : for (j=1; j<=n; j++)
2030 [ + + ]: 2302 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
2031 : : }
2032 : 387 : return v;
2033 : : }
2034 : :
2035 : : /* return the roots of p with absolute error bit */
2036 : : static GEN
2037 : 525 : roots_com(GEN q, long bit)
2038 : : {
2039 : : GEN L, p;
2040 : 525 : long v = RgX_valrem_inexact(q, &p);
2041 [ - + ]: 525 : if (lg(p) == 3) L = cgetg(1,t_VEC); /* constant polynomial */
2042 [ + + ]: 525 : else L = isexactpol(p)? solve_exact_pol(p,bit): all_roots(p,bit);
2043 [ + + ]: 525 : if (v)
2044 : : {
2045 : 75 : GEN M, z, t = gel(q,2);
2046 : : long i, x, y, l, n;
2047 : :
2048 [ + + ]: 75 : if (isrationalzero(t)) x = -bit;
2049 : : else
2050 : : {
2051 : 5 : n = gexpo(t);
2052 : 5 : x = n / v; l = degpol(q);
2053 [ + + ]: 25 : for (i = v; i <= l; i++)
2054 : : {
2055 : 20 : t = gel(q,i+2);
2056 [ - + ]: 20 : if (isrationalzero(t)) continue;
2057 : 20 : y = (n - gexpo(t)) / i;
2058 [ + + ]: 20 : if (y < x) x = y;
2059 : : }
2060 : : }
2061 : 75 : z = real_0_bit(x); l = v + lg(L);
2062 : 75 : M = cgetg(l, t_VEC); L -= v;
2063 [ + + ]: 150 : for (i = 1; i <= v; i++) gel(M,i) = z;
2064 [ + + ]: 230 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
2065 : 75 : L = M;
2066 : : }
2067 : 525 : return L;
2068 : : }
2069 : :
2070 : : static GEN
2071 : 31476 : tocomplex(GEN x, long l, long bit)
2072 : : {
2073 : : GEN y;
2074 [ + + ]: 31476 : if (typ(x) == t_COMPLEX)
2075 : : {
2076 [ + + ]: 30642 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
2077 : 2910 : x = gel(x,2);
2078 : 2910 : y = cgetg(3,t_COMPLEX);
2079 : 2910 : gel(y,1) = real_0_bit(-bit);
2080 : 2910 : gel(y,2) = mygprecrc(x, l, -bit);
2081 : : }
2082 : : else
2083 : : {
2084 : 834 : y = cgetg(3,t_COMPLEX);
2085 : 834 : gel(y,1) = mygprecrc(x, l, -bit);
2086 : 834 : gel(y,2) = real_0_bit(-bit);
2087 : : }
2088 : 31476 : return y;
2089 : : }
2090 : :
2091 : : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
2092 : : * up to 2^-e absolute error */
2093 : : static int
2094 : 99064 : cmp_complex_appr(void *E, GEN x, GEN y)
2095 : : {
2096 : 99064 : long e = (long)E;
2097 : : GEN z, xi, yi, xr, yr;
2098 : : long sxi, syi;
2099 [ + + ]: 99064 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
2100 : 40497 : else { xr = x; xi = NULL; sxi = 0; }
2101 [ + + ]: 99064 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
2102 : 36646 : else { yr = y; yi = NULL; syi = 0; }
2103 : : /* Compare absolute values of imaginary parts */
2104 [ + + ]: 99064 : if (!sxi)
2105 : : {
2106 [ + + ][ + - ]: 41510 : if (syi && expo(yi) >= e) return -1;
2107 : : /* |Im x| ~ |Im y| ~ 0 */
2108 : : }
2109 [ + + ]: 57554 : else if (!syi)
2110 : : {
2111 [ + - ][ + - ]: 1914 : if (sxi && expo(xi) >= e) return 1;
2112 : : /* |Im x| ~ |Im y| ~ 0 */
2113 : : }
2114 : : else
2115 : : {
2116 : : long sz;
2117 : 55640 : z = addrr_sign(xi, 1, yi, -1);
2118 : 55640 : sz = signe(z);
2119 [ + + ][ + - ]: 55640 : if (sz && expo(z) >= e) return (int)sz;
2120 : : }
2121 : : /* |Im x| ~ |Im y|, sort according to real parts */
2122 : 57758 : z = subrr(xr, yr);
2123 [ + + ]: 57758 : if (expo(z) >= e) return (int)signe(z);
2124 : : /* Re x ~ Re y. Place negative absolute value before positive */
2125 : 99064 : return (int) (sxi - syi);
2126 : : }
2127 : :
2128 : : static GEN
2129 : 6721 : clean_roots(GEN L, long l, long bit, long clean)
2130 : : {
2131 : 6721 : long i, n = lg(L), ex = 5 - bit;
2132 : 6721 : GEN res = cgetg(n,t_COL);
2133 [ + + ]: 53877 : for (i=1; i<n; i++)
2134 : : {
2135 : 47156 : GEN c = gel(L,i);
2136 [ + + ][ + + ]: 47156 : if (clean && isrealappr(c,ex))
2137 : : {
2138 [ - + ]: 15680 : if (typ(c) == t_COMPLEX) c = gel(c,1);
2139 : 15680 : c = mygprecrc(c, l, -bit);
2140 : : }
2141 : : else
2142 : 31476 : c = tocomplex(c, l, bit);
2143 : 47156 : gel(res,i) = c;
2144 : : }
2145 : 6721 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
2146 : 6721 : return res;
2147 : : }
2148 : :
2149 : : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
2150 : : static GEN
2151 : 545 : roots_aux(GEN p, long l, long clean)
2152 : : {
2153 : 545 : pari_sp av = avma;
2154 : : long bit;
2155 : : GEN L;
2156 : :
2157 [ + + ]: 545 : if (typ(p) != t_POL)
2158 : : {
2159 [ + + ]: 15 : if (gequal0(p)) pari_err_ROOTS0("roots");
2160 [ + + ]: 10 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
2161 : 5 : return cgetg(1,t_COL); /* constant polynomial */
2162 : : }
2163 [ - + ]: 530 : if (!signe(p)) pari_err_ROOTS0("roots");
2164 : 530 : checkvalidpol(p);
2165 [ - + ]: 525 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2166 [ - + ]: 525 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2167 : 525 : bit = prec2nbits(l);
2168 : 525 : L = roots_com(p, bit);
2169 : 530 : return gerepileupto(av, clean_roots(L, l, bit, clean));
2170 : : }
2171 : : GEN
2172 : 430 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
2173 : : /* clean up roots. If root is real replace it by its real part */
2174 : : GEN
2175 : 115 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
2176 : :
2177 : : /* private variant of conjvec. Allow non rational coefficients, shallow
2178 : : * function. */
2179 : : GEN
2180 : 55 : polmod_to_embed(GEN x, long prec)
2181 : : {
2182 : 55 : GEN v, T = gel(x,1), A = gel(x,2);
2183 : : long i, l;
2184 [ + - ][ - + ]: 55 : if (typ(A) != t_POL || varn(A) != varn(T))
2185 : : {
2186 : 0 : checkvalidpol(T);
2187 : 0 : return const_col(degpol(T), A);
2188 : : }
2189 : 55 : v = cleanroots(T,prec); l = lg(v);
2190 [ + + ]: 165 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
2191 : 55 : return v;
2192 : : }
2193 : :
2194 : : GEN
2195 : 6196 : QX_complex_roots(GEN p, long l)
2196 : : {
2197 : 6196 : pari_sp av = avma;
2198 : : long bit;
2199 : : GEN L;
2200 : :
2201 [ - + ]: 6196 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
2202 [ - + ]: 6196 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2203 [ - + ]: 6196 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2204 : 6196 : bit = prec2nbits(l);
2205 : 6196 : L = all_roots(Q_primpart(p), bit);
2206 : 6196 : return gerepileupto(av, clean_roots(L, l, bit, 1));
2207 : : }
2208 : :
2209 : : /********************************************************************/
2210 : : /** **/
2211 : : /** REAL ROOTS OF INTEGER POLYNOMIAL **/
2212 : : /** **/
2213 : : /********************************************************************/
2214 : :
2215 : : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
2216 : : * The inversion is implicit (we take coefficients backwards) */
2217 : : static long
2218 : 34824 : X2XP1(GEN P, long deg, GEN *Premapped)
2219 : : {
2220 : 34824 : pari_sp av = avma, lim = stack_lim(av, 3);
2221 : 34824 : GEN v = shallowcopy(P);
2222 : : long i, j, vlim, nb, s, s2;
2223 : : char flag;
2224 : :
2225 : 34824 : vlim = deg+2;
2226 : 34824 : nb = 0;
2227 : 34824 : i = 0;
2228 : : do
2229 : : {
2230 [ + + ]: 392864 : for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2231 : 34859 : s = -signe(gel(v, vlim));
2232 : 34859 : vlim--;
2233 : 34859 : i++;
2234 : : }
2235 [ + + ]: 34859 : while (!s);
2236 [ + + ][ + + ]: 34824 : if (vlim != deg + 1 && Premapped) setlg(v, vlim + 2);
2237 : :
2238 [ + + ]: 192628 : for (; i <= deg - 1; i++)
2239 : : {
2240 : 179607 : s2 = -signe(gel(v, 2));
2241 : 179607 : flag = (s2 == s);
2242 [ + + ]: 1879915 : for (j = 2; j < vlim; j++)
2243 : : {
2244 : 1700308 : gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2245 [ + + ]: 1700308 : if (flag) flag = (s2 != signe(gel(v, j+1)));
2246 : : }
2247 [ + + ]: 179607 : if (s == signe(gel(v, vlim)))
2248 : : {
2249 : 33902 : nb++;
2250 [ + + ]: 33902 : if (nb >= 2) { avma = av; return 2; }
2251 : 21806 : s = -s;
2252 : : }
2253 [ - + ]: 167511 : if (low_stack(lim, stack_lim(av, 3)))
2254 : : {
2255 [ # # ]: 0 : if (!Premapped) setlg(v, vlim);
2256 : 0 : v = gerepileupto(av, v);
2257 [ # # ]: 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "X2XP1");
2258 : : }
2259 [ + + ][ + + ]: 167511 : if (flag && !Premapped) goto END;
2260 : 157804 : vlim--;
2261 : : }
2262 [ + + ]: 13021 : if (s == signe(gel(v, vlim))) nb++;
2263 : : END:
2264 [ + + ][ + + ]: 22728 : if (Premapped && nb == 1)
2265 : 4079 : *Premapped = gerepileupto(av, v);
2266 : : else
2267 : 18649 : avma = av;
2268 : 34824 : return nb;
2269 : : }
2270 : :
2271 : : static long
2272 : 0 : _intervalcmp(GEN x, GEN y)
2273 : : {
2274 [ # # ]: 0 : if (typ(x) == t_VEC) x = gel(x, 1);
2275 [ # # ]: 0 : if (typ(y) == t_VEC) y = gel(y, 1);
2276 : 0 : return gcmp(x, y);
2277 : : }
2278 : :
2279 : : static GEN
2280 : 324817 : _gen_nored(void *E, GEN x) { (void)E; return x; }
2281 : : static GEN
2282 : 1159653 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2283 : : static GEN
2284 : 254973 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2285 : : static GEN
2286 : 263137 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2287 : : static GEN
2288 : 383133 : _gen_one(void *E) { (void)E; return gen_1; }
2289 : : static GEN
2290 : 5161 : _gen_zero(void *E) { (void)E; return gen_0; }
2291 : :
2292 : : static struct bb_algebra mp_algebra = { _gen_nored,_mp_add,_mp_mul,_mp_sqr,_gen_one,_gen_zero };
2293 : :
2294 : : static GEN
2295 : 1414450 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
2296 : :
2297 : : /* Split the polynom P in two parts, each with unique sign.
2298 : : * Moreover compute the two parts of the derivative of P. */
2299 : : static long
2300 : 4079 : split_polynoms(GEN P, long deg, long s0, GEN *Pp, GEN *Pm, GEN *Pprimep, GEN *Pprimem)
2301 : : {
2302 : 4079 : long i, degneg, v = evalvarn(varn(P));
2303 [ + + ][ + - ]: 20549 : for(i=1; i <= deg; i++) if (signe(gel(P, i+2)) != s0) break;
2304 : 4079 : degneg = i;
2305 : 4079 : *Pm = cgetg(degneg + 2, t_POL);
2306 : 4079 : (*Pm)[1] = v;
2307 : 4079 : *Pprimem = cgetg(degneg + 1, t_POL);
2308 : 4079 : (*Pprimem)[1] = v;
2309 [ + + ]: 24628 : for(i=0; i < degneg; i++)
2310 : : {
2311 : 20549 : GEN elt = gel(P, i+2);
2312 : 20549 : gel(*Pm, i+2) = elt;
2313 [ + + ]: 20549 : if (i) gel(*Pprimem, i+1) = mului(i, elt);
2314 : : }
2315 : 4079 : *Pp = cgetg(deg - degneg + 3, t_POL);
2316 : 4079 : (*Pp)[1] = v;
2317 : 4079 : *Pprimep = cgetg(deg - degneg + 3, t_POL);
2318 : 4079 : (*Pprimep)[1] = v;
2319 [ + + ]: 22944 : for(i=degneg; i <= deg; i++)
2320 : : {
2321 : 18865 : GEN elt = gel(P, i+2);
2322 : 18865 : gel(*Pp, i+2-degneg) = elt;
2323 : 18865 : gel(*Pprimep, i+2-degneg) = mului(i, elt);
2324 : : }
2325 : 4079 : *Pm = normalizepol_lg(*Pm, degneg+2);
2326 : 4079 : *Pprimem = normalizepol_lg(*Pprimem, degneg+1);
2327 : 4079 : *Pp = normalizepol_lg(*Pp, deg-degneg+3);
2328 : 4079 : *Pprimep = normalizepol_lg(*Pprimep, deg-degneg+3);
2329 : 4079 : return degpol(*Pm) + 1;
2330 : : }
2331 : :
2332 : : static GEN
2333 : 129979 : bkeval_single_power(long d, GEN V)
2334 : : {
2335 : 129979 : long mp = lg(V) - 2;
2336 [ - + ]: 129979 : if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2337 : 129979 : return gel(V, d+1);
2338 : : }
2339 : :
2340 : : static GEN
2341 : 129979 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long deg, long degneg, long bitprec)
2342 : : {
2343 : 129979 : GEN vp = gen_bkeval_powers(Pp, deg-degneg, pows, NULL, &mp_algebra, _mp_cmul);
2344 : 129979 : GEN vm = gen_bkeval_powers(Pm, degneg-1, pows, NULL, &mp_algebra, _mp_cmul);
2345 : 129979 : GEN xa = bkeval_single_power(degneg, pows);
2346 : : GEN r;
2347 : 129979 : vp = gmul(vp, xa);
2348 : 129979 : r = gadd(vp, vm);
2349 [ + + ][ + + ]: 129979 : if (expo(vp) - (signe(r) ? expo(r) : 0) > prec2nbits(realprec(vp)) - bitprec)
2350 : 8745 : return NULL;
2351 : 129979 : return r;
2352 : : }
2353 : :
2354 : : /* Newton for polynom P. One solution between 0 and infinity, P' has also at
2355 : : * most one zero. P is almost certainly without zero coefficients. */
2356 : : static GEN
2357 : 4079 : polsolve(GEN P, long bitprec)
2358 : : {
2359 : 4079 : pari_sp av = avma;
2360 : : GEN Pp, Pm, Pprimep, Pprimem, pows;
2361 : : GEN Pprime, Pprime2;
2362 : 4079 : GEN ra, rb, rc, rcold = NULL, Pc, Ppc;
2363 : 4079 : long deg = degpol(P), degneg, rt;
2364 : 4079 : long s0, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
2365 : : long bitaddprec, addprec;
2366 : 4079 : long expoold = LONG_MAX, iter;
2367 [ - + ]: 4079 : if (deg == 1)
2368 : : {
2369 : 0 : ra = negr(divrr(itor(gel(P, 2), bitprec), itor(gel(P, 3), bitprec)));
2370 : 0 : return gerepileuptoleaf(av, ra);
2371 : : }
2372 : 4079 : Pprime = ZX_deriv(P); Pprime2 = ZX_deriv(Pprime);
2373 : 4079 : bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
2374 : 4079 : s0 = signe(gel(P, 2));
2375 : 4079 : degneg = split_polynoms(P, deg, s0, &Pp, &Pm, &Pprimep, &Pprimem);
2376 : 4079 : rt = maxss(degneg, brent_kung_optpow(maxss(deg-degneg, degneg-1), 2, 1));
2377 : : {
2378 : : /* Optimized Cauchy bound */
2379 : 4079 : GEN summin = gen_0, absmaj;
2380 : : long i;
2381 : :
2382 : 4079 : absmaj = gel(P, 2);
2383 [ + + ]: 20549 : for(i=1; i < degneg; i++)
2384 : : {
2385 : 16470 : GEN elt = gel(P, i+2);
2386 [ + + ]: 16470 : if (absi_cmp(absmaj, elt) < 0) absmaj = elt;
2387 : : }
2388 : 4079 : summin = gel(P, i+2);
2389 [ + + ]: 18865 : while (++i <= deg) summin = addii(summin, gel(P, i+2));
2390 : 4079 : rb = subsr(1, rdivii(absmaj, summin, cprec));
2391 : : do
2392 : : {
2393 : 4079 : pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2394 : 4079 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec);
2395 [ - + ]: 4079 : if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
2396 [ + - ]: 4079 : if (signe(Pc) != s0) break;
2397 : 0 : shiftr_inplace(rb,1);
2398 : : }
2399 : 0 : while (1);
2400 : : }
2401 : 4079 : ra = NULL;
2402 : 4079 : iter = 0;
2403 : : while (1)
2404 : : {
2405 : : GEN wdth;
2406 : 37552 : iter++;
2407 [ + + ]: 37552 : if (ra)
2408 : 25648 : rc = shiftr(addrr(ra, rb), -1);
2409 : : else
2410 : 11904 : rc = shiftr(rb, -1);
2411 : : do
2412 : : {
2413 : 37552 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2414 : 37552 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+2);
2415 [ + - ]: 37552 : if (Pc) break;
2416 : 0 : cprec++; rc = rtor(rc, cprec);
2417 : 0 : } while (1);
2418 [ + + ]: 37552 : if (signe(Pc) == s0)
2419 : 16526 : ra = rc;
2420 : : else
2421 : 21026 : rb = rc;
2422 [ + + ]: 37552 : if (!ra) continue;
2423 : 29727 : wdth = subrr(rb, ra);
2424 [ + + ]: 29727 : if (!(iter % 8))
2425 : : {
2426 : 4364 : GEN m1 = poleval(Pprime, ra), M2;
2427 [ + + ]: 4364 : if (signe(m1) == s0) continue;
2428 : 4264 : M2 = poleval(Pprime2, rb);
2429 [ + + ]: 4264 : if (absr_cmp(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2430 : 4079 : break;
2431 : : }
2432 [ - + ]: 25363 : else if (gexpo(wdth) <= -bitprec)
2433 : 0 : break;
2434 : 33473 : }
2435 : 4079 : rc = rb;
2436 : 4079 : iter = 0;
2437 : : while (1)
2438 : : {
2439 : : long exponew;
2440 : : GEN dist;
2441 : 44174 : iter++;
2442 : 44174 : rcold = rc;
2443 : 44174 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2444 : 44174 : Ppc = splitpoleval(Pprimep, Pprimem, pows, deg-1, degneg-1, bitaddprec+4);
2445 [ + - ]: 44174 : if (Ppc)
2446 : 44174 : Pc = splitpoleval(Pp, Pm, pows, deg, degneg, bitaddprec+4);
2447 [ + - ][ + + ]: 44174 : if (!Ppc || !Pc)
2448 : : {
2449 [ + + ]: 8745 : if (cprec >= prec+addprec)
2450 : 903 : cprec++;
2451 : : else
2452 : 7842 : cprec = minss(2*cprec, prec+addprec);
2453 : 8745 : rc = rtor(rc, cprec); /* Backtrack one step */
2454 : 8745 : continue;
2455 : : }
2456 : 35429 : dist = divrr(Pc, Ppc);
2457 : 35429 : rc = subrr(rc, dist);
2458 [ + - ][ - + ]: 35429 : if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2459 : : {
2460 [ # # ]: 0 : if (cprec >= prec+addprec) break;
2461 : 0 : cprec = minss(2*cprec, prec+addprec);
2462 : 0 : rc = rtor(rcold, cprec); /* Backtrack one step */
2463 : 0 : continue;
2464 : : }
2465 [ + + ]: 35429 : if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2466 : 30230 : exponew = expo(dist);
2467 [ + + ]: 30230 : if (exponew < -bitprec - 1)
2468 : : {
2469 [ + + ]: 5411 : if (cprec >= prec+addprec) break;
2470 : 1332 : cprec = minss(2*cprec, prec+addprec);
2471 : 1332 : rc = rtor(rc, cprec);
2472 : 1332 : continue;
2473 : : }
2474 [ + + ]: 24819 : if (exponew > expoold - 2)
2475 : : {
2476 [ - + ]: 1120 : if (cprec >= prec+addprec) break;
2477 : 1120 : cprec = minss(2*cprec, prec+addprec);
2478 : 1120 : rc = rtor(rc, cprec);
2479 : 1120 : expoold = LONG_MAX;
2480 : 1120 : continue;
2481 : : }
2482 : 23699 : expoold = exponew;
2483 : 40095 : }
2484 : 4079 : return gerepileuptoleaf(av, rtor(rc, prec));
2485 : : }
2486 : :
2487 : : static GEN
2488 : 4164 : usp(GEN Q0, long deg, long *nb_donep, long flag, long bitprec)
2489 : : {
2490 : : pari_sp av, lim;
2491 : : GEN Q, sol;
2492 : 4164 : long nb_todo, nbr = 0, ind, deg0, indf, i, k, nb, j;
2493 : 4164 : long listsize = 64, nb_done = 0;
2494 : : GEN c, Lc, Lk;
2495 : :
2496 : 4164 : av = avma;
2497 : :
2498 : 4164 : sol = const_col(deg, gen_0);
2499 : 4164 : deg0 = deg;
2500 : 4164 : Lc = const_vec(listsize, gen_0);
2501 : 4164 : Lk = cgetg(listsize+1, t_VECSMALL);
2502 : 4164 : c = gen_0;
2503 : 4164 : k = Lk[1] = 0;
2504 : 4164 : ind = 1; indf = 2;
2505 : 4164 : Q = gcopy(Q0);
2506 : :
2507 : 4164 : nb_todo = 1;
2508 : 4164 : lim = stack_lim(av, 2);
2509 [ + + ]: 38988 : while (nb_todo)
2510 : : {
2511 : 34824 : GEN nc = gel(Lc, ind), Qremapped;
2512 [ + + ]: 34824 : if (Lk[ind] == k + 1)
2513 : : {
2514 : 13990 : deg0 = deg;
2515 : 13990 : setlg(Q0, deg + 3);
2516 : 13990 : Q0 = ZX_rescale2n(Q0, 1);
2517 : 13990 : Q = Q_primpart(Q0);
2518 : 13990 : c = gen_0;
2519 : : }
2520 : :
2521 [ + + ]: 34824 : if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
2522 : :
2523 : 34824 : k = Lk[ind];
2524 : 34824 : c = nc;
2525 : 34824 : ind++;
2526 : 34824 : nb_todo--;
2527 : :
2528 [ + + ]: 34824 : if (equalii(gel(Q, 2), gen_0))
2529 : : {
2530 : 90 : GEN newsol = gmul2n(c, -k);
2531 [ + + ]: 120 : for (j = 1; j <= nbr; j++)
2532 [ + + ]: 55 : if (gequal(gel(sol, j), newsol)) break;
2533 [ + + ]: 90 : if (j > nbr) gel(sol, ++nbr) = newsol;
2534 : :
2535 : 90 : deg0--;
2536 [ + + ]: 360 : for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
2537 : 90 : setlg(Q, j);
2538 : : }
2539 : :
2540 [ + + ]: 34824 : nb = X2XP1(Q, deg0, flag == 1 ? &Qremapped : NULL);
2541 : 34824 : nb_done++;
2542 : :
2543 [ + + + ]: 34824 : switch (nb)
2544 : : {
2545 : : case 0:
2546 : 11920 : break;
2547 : : case 1:
2548 [ - + + ]: 7574 : switch(flag)
2549 : : {
2550 : : case 0:
2551 : : {
2552 : : GEN low, hi;
2553 : 0 : low = gmul2n(c, -k);
2554 : 0 : hi = gmul2n(addiu(c,1), -k);
2555 : 0 : gel(sol, ++nbr) = mkvec2(low, hi);
2556 : : }
2557 : 0 : break;
2558 : : case 1:
2559 : : { /* Caveat emptor: Qremapped is the reciprocal polynomial */
2560 : 4079 : GEN sr = polsolve(Qremapped, bitprec+1);
2561 : 4079 : sr = divrr(sr, addsr(1, sr));
2562 : 4079 : sr = gmul2n(addir(c, sr), -k);
2563 : 4079 : gel(sol, ++nbr) = rtor(sr, nbits2prec(bitprec));
2564 : : }
2565 : 4079 : break;
2566 : : default:
2567 : 3495 : gel(sol, ++nbr) = gen_0;
2568 : : }
2569 : 7574 : break;
2570 : :
2571 : : default:
2572 [ + + ]: 15330 : if (indf + 2 > listsize)
2573 : : {
2574 [ + - ]: 90 : if (ind>1)
2575 : : {
2576 [ + + ]: 200 : for (i = ind; i < indf; i++)
2577 : : {
2578 : 110 : gel(Lc, i-ind+1) = gel(Lc, i);
2579 : 110 : Lk[i-ind+1] = Lk[i];
2580 : : }
2581 : 90 : indf -= ind-1; ind = 1;
2582 : : }
2583 [ - + ]: 90 : if (indf + 2 > listsize)
2584 : : {
2585 : 0 : listsize *= 2;
2586 : 0 : Lc = vec_lengthen(Lc, listsize);
2587 : 0 : Lk = vecsmall_lengthen(Lk, listsize);
2588 : : }
2589 [ + + ]: 5740 : for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2590 : : }
2591 : 15330 : nc = shifti(c, 1);
2592 : 15330 : gel(Lc, indf) = nc;
2593 : 15330 : gel(Lc, indf + 1) = addis(nc, 1);
2594 : 15330 : Lk[indf] = Lk[indf + 1] = k + 1;
2595 : 15330 : indf += 2;
2596 : 15330 : nb_todo += 2;
2597 : : }
2598 : :
2599 [ - + ]: 34824 : if (low_stack(lim, stack_lim(av, 2)))
2600 : : {
2601 : 0 : gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
2602 [ # # ]: 34824 : if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_uspensky", avma);
2603 : : }
2604 : : }
2605 : :
2606 : 4164 : setlg(sol, nbr+1);
2607 : 4164 : *nb_donep += nb_done;
2608 : 4164 : return gerepilecopy(av, sol);
2609 : : }
2610 : :
2611 : : GEN
2612 : 4224 : ZX_uspensky(GEN P, GEN ab, long flag, long bitprec)
2613 : : {
2614 : 4224 : pari_sp av = avma;
2615 : 4224 : GEN a, b, sol = NULL, Pcur;
2616 : : double fb;
2617 : 4224 : long nbz, deg, nb_done = 0;
2618 : :
2619 : 4224 : deg = degpol(P);
2620 [ + + ][ + - ]: 4224 : if (deg == 0) return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2621 [ + + ]: 4129 : if (ab)
2622 : : {
2623 [ + + ]: 1654 : if (typ(ab) == t_VEC)
2624 : : {
2625 [ - + ]: 85 : if (lg(ab) != 3) pari_err_DIM("ZX_uspensky");
2626 : 85 : a = gel(ab, 1);
2627 : 85 : b = gel(ab, 2);
2628 : : }
2629 : : else
2630 : : {
2631 : 1569 : a = ab;
2632 : 1569 : b = mkoo();
2633 : : }
2634 : : }
2635 : : else
2636 : : {
2637 : 2475 : a = mkmoo();
2638 : 2475 : b = mkoo();
2639 : : }
2640 [ + + + ]: 4129 : switch (gcmp(a, b))
2641 : : {
2642 [ - + ]: 5 : case 1: avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2643 : : case 0:
2644 [ + - ][ + + ]: 15 : if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
2645 [ + + ]: 10 : { avma = av; return flag <= 1 ? mkcolcopy(a): gen_1; }
2646 : : else
2647 [ + - ]: 5 : { avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
2648 : : }
2649 [ + + ]: 4109 : if (deg == 1)
2650 : : {
2651 : 1429 : sol = gdiv(gneg(gel(P, 2)), pollead(P, -1));
2652 [ + + ][ - + ]: 1429 : if (gcmp(a, sol) > 0 || gcmp(sol, b) > 0)
2653 [ - + ]: 380 : { avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0; }
2654 [ + + ]: 1049 : if (flag >= 2) { avma = av; return gen_1; }
2655 : 479 : return gerepilecopy(av, mkcol(sol));
2656 : : }
2657 : 2680 : nbz = ZX_valrem(P, &Pcur);
2658 : 2680 : deg -= nbz;
2659 [ + - ]: 2680 : if (!nbz) Pcur = P;
2660 [ - + ][ # # ]: 2680 : if (nbz && (gsigne(a) > 0 || gsigne(b) < 0)) nbz = 0;
[ # # ]
2661 [ - + + ]: 2680 : switch(flag)
2662 : : {
2663 : : case 0:
2664 : 0 : sol = zerocol(nbz);
2665 : 0 : break;
2666 : : case 1:
2667 : 1020 : sol = const_col(nbz, real_0(bitprec));
2668 : 1020 : break;
2669 : : /* case 2: nothing */
2670 : : }
2671 : :
2672 [ + + ][ + + ]: 2680 : if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
[ + - ]
2673 : : {
2674 : 20 : fb = fujiwara_bound_real(Pcur, -1);
2675 [ + + ]: 20 : if (fb > -pariINFINITY)
2676 : 5 : a = negi(int2n((long)ceil(fb)));
2677 : : else
2678 : 15 : a = gen_0;
2679 : : }
2680 [ + + ][ + + ]: 2680 : if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
[ + + ]
2681 : : {
2682 : 10 : fb = fujiwara_bound_real(Pcur, 1);
2683 [ + - ]: 10 : if (fb > -pariINFINITY)
2684 : 10 : b = int2n((long)ceil(fb));
2685 : : else
2686 : 0 : b = gen_0;
2687 : : }
2688 : :
2689 [ + + ][ + + ]: 2680 : if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
2690 : : {
2691 : : pari_sp av1;
2692 : 65 : GEN den = lcmii(denom(a), denom(b)), diff, unscaledres, co, Pdiv;
2693 : : GEN ascaled;
2694 : : long i;
2695 [ + + ]: 65 : if (!is_pm1(den))
2696 : : {
2697 : 20 : Pcur = ZX_rescale(Pcur, den);
2698 : 20 : ascaled = gmul(a, den);
2699 : : }
2700 : : else
2701 : : {
2702 : 45 : den = NULL;
2703 : 45 : ascaled = a;
2704 : : }
2705 [ + + ]: 65 : diff = subii(den ? gmul(b,den) : b, ascaled);
2706 : 65 : Pcur = ZX_unscale(ZX_translate(Pcur, ascaled), diff);
2707 : 65 : av1 = avma;
2708 : 65 : Pdiv = cgetg(deg+2, t_POL);
2709 : 65 : Pdiv[1] = Pcur[1];
2710 : 65 : co = gel(Pcur, deg+2);
2711 [ + + ]: 260 : for (i = deg; --i >= 0; )
2712 : : {
2713 : 195 : gel(Pdiv, i+2) = co;
2714 : 195 : co = addii(co, gel(Pcur, i+2));
2715 : : }
2716 [ + + ]: 65 : if (!signe(co))
2717 : : {
2718 : 35 : Pcur = Pdiv;
2719 : 35 : deg--;
2720 [ + + ]: 35 : if (flag <= 1)
2721 : 10 : sol = concat(sol, b);
2722 : : else
2723 : 25 : nbz++;
2724 : : }
2725 : : else
2726 : 30 : avma = av1;
2727 : 65 : unscaledres = usp(Pcur, deg, &nb_done, flag, bitprec);
2728 [ + + ]: 65 : if (flag <= 1)
2729 : : {
2730 [ + + ]: 45 : for (i = 1; i < lg(unscaledres); i++)
2731 : : {
2732 : 25 : GEN z = gmul(diff, gel(unscaledres, i));
2733 [ - + ]: 25 : if (typ(z) == t_VEC)
2734 : : {
2735 : 0 : gel(z, 1) = gadd(ascaled, gel(z, 1));
2736 : 0 : gel(z, 2) = gadd(ascaled, gel(z, 2));
2737 : : }
2738 : : else
2739 : 25 : z = gadd(ascaled, z);
2740 [ + + ]: 25 : if (den) z = gdiv(z, den);
2741 : 25 : gel(unscaledres, i) = z;
2742 : : }
2743 : 20 : sol = concat(sol, unscaledres);
2744 : : }
2745 : : else
2746 : 45 : nbz += lg(unscaledres) - 1;
2747 : : }
2748 [ + + ][ + + ]: 2680 : if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(Pcur, 1)) > -pariINFINITY)
2749 : : {
2750 : : GEN Pcurp, unscaledres;
2751 : 2255 : long bp = (long)ceil(fb);
2752 : 2255 : Pcurp = ZX_unscale2n(Pcur, bp);
2753 : 2255 : unscaledres = usp(Pcurp, deg, &nb_done, flag, bitprec);
2754 [ + + ]: 2255 : if (flag <= 1)
2755 : 990 : sol = concat(sol, gmul2n(unscaledres, bp));
2756 : : else
2757 : 1265 : nbz += lg(unscaledres) - 1;
2758 : : }
2759 [ + + ][ + + ]: 2680 : if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(Pcur,-1)) > -pariINFINITY)
2760 : : {
2761 : : GEN Pcurm, unscaledres;
2762 : 1844 : long i, bm = (long)ceil(fb);
2763 : 1844 : Pcurm = ZX_unscale(Pcur, gen_m1);
2764 : 1844 : Pcurm = ZX_unscale2n(Pcurm, bm);
2765 : 1844 : unscaledres = usp(Pcurm,deg,&nb_done,flag,bitprec);
2766 [ + + ]: 1844 : if (flag <= 1)
2767 : : {
2768 [ + + ]: 2324 : for (i = 1; i < lg(unscaledres); i++)
2769 : : {
2770 : 1499 : GEN z = gneg(gmul2n(gel(unscaledres, i), bm));
2771 [ - + ]: 1499 : if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
2772 : 1499 : gel(unscaledres, i) = z;
2773 : : }
2774 : 825 : sol = concat(unscaledres, sol);
2775 : : }
2776 : : else
2777 : 1019 : nbz += lg(unscaledres) - 1;
2778 : : }
2779 : :
2780 [ - + ]: 2680 : if (DEBUGLEVEL > 4)
2781 : 0 : err_printf("ZX_uspensky: Number of visited nodes: %d\n", nb_done);
2782 : :
2783 [ + + ]: 2680 : if (flag >= 2) return utoi(nbz);
2784 [ + - ]: 1020 : if (flag)
2785 : 1020 : sol = sort(sol);
2786 : : else
2787 : 0 : sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
2788 : 4224 : return gerepileupto(av, sol);
2789 : : }
2790 : :
2791 : : /* x a scalar */
2792 : : static GEN
2793 : 25 : rootsdeg0(GEN x)
2794 : : {
2795 [ + + ]: 25 : if (!is_rational_t(typ(x))) pari_err_TYPE("realroots",x);
2796 [ + + ]: 20 : if (gequal0(x)) pari_err_ROOTS0("realroots");
2797 : 10 : return cgetg(1,t_COL); /* constant polynomial */
2798 : : }
2799 : : static void
2800 : 200 : checkbound(GEN a)
2801 : : {
2802 [ + - ]: 200 : switch(typ(a))
2803 : : {
2804 : 200 : case t_INT: case t_FRAC: case t_INFINITY: break;
2805 : 0 : default: pari_err_TYPE("polrealroots", a);
2806 : : }
2807 : 200 : }
2808 : : static GEN
2809 : 1689 : check_ab(GEN ab)
2810 : : {
2811 : : GEN a, b;
2812 [ + + ]: 1689 : if (!ab) return NULL;
2813 [ + - ][ - + ]: 100 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
2814 : 100 : a = gel(ab,1); checkbound(a);
2815 : 100 : b = gel(ab,2); checkbound(b);
2816 [ + + ][ + - ]: 100 : if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
[ + + ]
2817 [ + - ]: 15 : typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
2818 : 1689 : return ab;
2819 : : }
2820 : : GEN
2821 : 1624 : realroots(GEN P, GEN ab, long prec)
2822 : : {
2823 : 1624 : pari_sp av = avma;
2824 : 1624 : long nrr = 0;
2825 : 1624 : GEN sol = NULL, fa, ex;
2826 : : long i, j, k;
2827 : :
2828 : 1624 : ab = check_ab(ab);
2829 [ + + ]: 1624 : if (typ(P) != t_POL) return rootsdeg0(P);
2830 [ + + + ]: 1609 : switch(degpol(P))
2831 : : {
2832 : 5 : case -1: return rootsdeg0(gen_0);
2833 : 5 : case 0: return rootsdeg0(gel(P,2));
2834 : : }
2835 : 1599 : P = Q_primpart(P);
2836 [ - + ]: 1599 : if (!RgX_is_ZX(P)) pari_err_TYPE("realroots",P);
2837 : 1599 : fa = ZX_squff(P, &ex);
2838 [ + + ]: 3203 : for (i = 1; i < lg(fa); i++)
2839 : : {
2840 : 1604 : GEN Pi = gel(fa, i), soli, soli2 = NULL;
2841 : 1604 : long n, nrri = 0, h, nbz;
2842 [ + + ]: 1604 : if (ab)
2843 : 30 : h = 1;
2844 : : else
2845 : 1574 : Pi = RgX_deflate_max(Pi, &h);
2846 [ + + ]: 1604 : if (!signe(gel(Pi, 2)))
2847 : : {
2848 : 95 : Pi = RgX_shift_shallow(Pi, -1);
2849 : 95 : nbz = 1;
2850 : : }
2851 : : else
2852 : 1509 : nbz = 0;
2853 [ + + ]: 1604 : soli = ZX_uspensky(Pi, h%2 ? ab: gen_0, 1, prec2nbits(prec));
2854 : 1604 : n = lg(soli);
2855 [ + + ]: 1604 : if (!(h % 2)) soli2 = cgetg(n, t_COL);
2856 [ + + ]: 6207 : for (j = 1; j < n; j++)
2857 : : {
2858 : 4603 : GEN elt = gel(soli, j);
2859 [ + + ]: 4603 : if (typ(elt) != t_REAL)
2860 : : {
2861 : 524 : nrri++;
2862 : 524 : elt = gtofp(elt, prec);
2863 : 524 : gel(soli, j) = elt;
2864 : : }
2865 [ + + ]: 4603 : if (h > 1)
2866 : : { /* note: elt != 0 because we are square free */
2867 : : GEN r;
2868 [ + + ]: 834 : if (h == 2)
2869 : 824 : r = sqrtr(elt);
2870 : : else
2871 : : {
2872 [ + + ]: 10 : if (signe(elt) < 0)
2873 : 5 : r = negr(sqrtnr(negr(elt), h));
2874 : : else
2875 : 5 : r = sqrtnr(elt, h);
2876 : : }
2877 : 834 : gel(soli, j) = r;
2878 [ + + ]: 834 : if (!(h % 2)) gel(soli2, j) = negr(r);
2879 : : }
2880 : : }
2881 [ + + ]: 1604 : if (!(h % 2)) soli = shallowconcat(soli, soli2);
2882 [ + + ]: 1604 : if (nbz) soli = shallowconcat(soli, real_0(prec));
2883 [ + + ]: 3208 : for (k = 1; k <= ex[i]; k++)
2884 [ + + ]: 1604 : sol = sol ? shallowconcat(sol, soli) : soli;
2885 : 1604 : nrr += ex[i]*nrri;
2886 : : }
2887 : :
2888 [ - + ]: 1599 : if (DEBUGLEVEL > 4)
2889 : : {
2890 : 0 : err_printf("Number of real roots: %d\n", lg(sol)-1);
2891 : 0 : err_printf(" -- of which 2-integral: %ld\n", nrr);
2892 : : }
2893 : :
2894 : 1609 : return gerepileupto(av, sort(sol));
2895 : : }
2896 : :
2897 : : /* P non-constant, squarefree ZX */
2898 : : long
2899 : 2565 : ZX_sturm(GEN P)
2900 : : {
2901 : 2565 : pari_sp av = avma;
2902 : : long h, r;
2903 : 2565 : P = RgX_deflate_max(P, &h);
2904 [ + + ]: 2565 : if (odd(h))
2905 : 1515 : r = itos(ZX_uspensky(P, NULL, 2, 0));
2906 : : else
2907 : 1050 : r = 2*itos(ZX_uspensky(P, gen_0, 2, 0));
2908 : 2565 : avma = av; return r;
2909 : : }
2910 : : /* P non-constant, squarefree ZX */
2911 : : long
2912 : 65 : ZX_sturmpart(GEN P, GEN ab)
2913 : : {
2914 : 65 : pari_sp av = avma;
2915 : : long r;
2916 [ + + ]: 65 : if (!check_ab(ab)) return ZX_sturm(P);
2917 : 55 : r = itos(ZX_uspensky(P, ab, 2, 0));
2918 : 65 : avma = av; return r;
2919 : : }
|
