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 : 5152492 : addCC(GEN x, GEN y)
35 : : {
36 : : GEN z;
37 : :
38 [ + + ]: 5152492 : if (typ(x) == t_INT)
39 : : {
40 [ + + ]: 4623700 : if (typ(y) == t_INT) return addii(x,y);
41 : : /* ty == t_COMPLEX */
42 : 447 : z = cgetg(3,t_COMPLEX);
43 : 447 : gel(z,1) = addii(x, gel(y,1));
44 : 447 : gel(z,2) = icopy(gel(y,2)); return z;
45 : : }
46 : : /* tx == t_COMPLEX */
47 : 528792 : z = cgetg(3,t_COMPLEX);
48 [ + + ]: 528792 : if (typ(y) == t_INT)
49 : : {
50 : 7438 : gel(z,1) = addii(gel(x,1),y);
51 : 7438 : gel(z,2) = icopy(gel(x,2)); return z;
52 : : }
53 : : /* ty == t_COMPLEX */
54 : 521354 : gel(z,1) = addii(gel(x,1),gel(y,1));
55 : 5152492 : 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 : 9785469 : mulCC(GEN x, GEN y)
60 : : {
61 : : GEN z;
62 : :
63 [ + + ]: 9785469 : if (typ(x) == t_INT)
64 : : {
65 [ + + ]: 8492675 : if (typ(y) == t_INT) return mulii(x,y);
66 : : /* ty == t_COMPLEX */
67 : 1686 : z = cgetg(3,t_COMPLEX);
68 : 1686 : gel(z,1) = mulii(x, gel(y,1));
69 : 1686 : gel(z,2) = mulii(x, gel(y,2)); return z;
70 : : }
71 : : /* tx == t_COMPLEX */
72 : 1292794 : z = cgetg(3,t_COMPLEX);
73 [ + + ]: 1292794 : if (typ(y) == t_INT)
74 : : {
75 : 378580 : gel(z,1) = mulii(gel(x,1),y);
76 : 378580 : gel(z,2) = mulii(gel(x,2),y); return z;
77 : : }
78 : : /* ty == t_COMPLEX */
79 : : {
80 : 914214 : pari_sp av = avma, tetpil;
81 : : GEN p1, p2;
82 : :
83 : 914214 : p1 = mulii(gel(x,1),gel(y,1));
84 : 914214 : p2 = mulii(gel(x,2),gel(y,2));
85 : 914214 : y = mulii(addii(gel(x,1),gel(x,2)),
86 : 1828428 : addii(gel(y,1),gel(y,2)));
87 : 914214 : x = addii(p1,p2); tetpil = avma;
88 : 914214 : gel(z,1) = subii(p1,p2);
89 : 914214 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
90 : 9785469 : return z;
91 : : }
92 : : }
93 : : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
94 : : static GEN
95 : 8179467 : sqrCC(GEN x)
96 : : {
97 : : GEN z;
98 : :
99 [ + + ]: 8179467 : if (typ(x) == t_INT) return sqri(x);
100 : : /* tx == t_COMPLEX */
101 : 1271078 : z = cgetg(3,t_COMPLEX);
102 : : {
103 : 1271078 : pari_sp av = avma, tetpil;
104 : : GEN y, p1, p2;
105 : :
106 : 1271078 : p1 = sqri(gel(x,1));
107 : 1271078 : p2 = sqri(gel(x,2));
108 : 1271078 : y = sqri(addii(gel(x,1),gel(x,2)));
109 : 1271078 : x = addii(p1,p2); tetpil = avma;
110 : 1271078 : gel(z,1) = subii(p1,p2);
111 : 1271078 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
112 : 8179467 : return z;
113 : : }
114 : : }
115 : :
116 : : static void
117 : 1381043 : set_karasquare_limit(long bit)
118 : : {
119 [ + + ]: 1381043 : if (bit<600) { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
120 [ + - ]: 281 : 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 : 1381043 : }
125 : :
126 : : /* assume lP > 0, lP = lgpol(P) */
127 : : static GEN
128 : 2856084 : CX_square_spec(GEN P, long lP)
129 : : {
130 : : GEN s, t;
131 : 2856084 : long i, j, l, nn, n = lP - 1;
132 : : pari_sp av;
133 : :
134 : 2856084 : nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
135 : 2856084 : gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
136 [ + + ]: 7489061 : for (i=1; i<=n; i++)
137 : : {
138 : 4632977 : av = avma; l = (i+1)>>1;
139 : 4632977 : t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
140 [ + + ]: 6402125 : for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
141 : 4632977 : t = gmul2n(t,1);
142 [ + + ]: 4632977 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
143 : 4632977 : gel(s,i+2) = gerepileupto(av, t);
144 : : }
145 : 2856084 : gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
146 [ + + ]: 5323383 : for ( ; i<nn; i++)
147 : : {
148 : 2467299 : av = avma; l = (i+1)>>1;
149 : 2467299 : t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
150 [ + + ]: 3383344 : for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
151 : 2467299 : t = gmul2n(t,1);
152 [ + + ]: 2467299 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
153 : 2467299 : gel(s,i+2) = gerepileupto(av, t);
154 : : }
155 : 2856084 : return normalizepol_lg(s, nn+3);
156 : : }
157 : : /* not stack clean */
158 : : static GEN
159 : 47137 : RgX_addspec(GEN x, long nx, GEN y, long ny)
160 : : {
161 : : GEN z, t;
162 : : long i;
163 [ + + ]: 47137 : if (nx == ny) {
164 : 26411 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
165 [ + + ]: 195503 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
166 : 26411 : return normalizepol_lg(z, nx+2);
167 : : }
168 [ + - ]: 20726 : if (ny < nx) {
169 : 20726 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
170 [ + + ]: 145960 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
171 [ + + ]: 41452 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
172 : 20726 : 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 : 47137 : 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 : 2903497 : karasquare(GEN P, long nP)
206 : : {
207 : : GEN Q, s0, s1, s2, a, t;
208 : 2903497 : long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
209 : : pari_sp av;
210 : :
211 [ + + ][ + + ]: 2903497 : if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
212 : 47137 : av = avma;
213 : 47137 : n0 = (n>>1) + 1; n1 = nP - n0;
214 : 47137 : s0 = karasquare(P, n0); Q = P + n0;
215 : 47137 : s2 = karasquare(Q, n1);
216 : 47137 : s1 = RgX_addspec(P, n0, Q, n1);
217 : 47137 : s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
218 : 47137 : N = (n<<1) + 1;
219 : 47137 : a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
220 : 47137 : t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
221 [ + + ]: 618336 : for (i=0; i < l; i++) gel(t,i) = gel(s0,i);
222 [ + + ]: 106042 : for ( ; i < N0; i++) gel(t,i) = gen_0;
223 : 47137 : t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
224 [ + + ]: 552296 : for (i=0; i < l; i++) gel(t,i) = gel(s2,i);
225 [ + + ]: 83493 : for ( ; i < N1; i++) gel(t,i) = gen_0;
226 : 47137 : t = a+2 + n0; l = lgpol(s1); s1 += 2;
227 [ + + ]: 582910 : for (i=0; i < l; i++) gel(t,i) = gadd(gel(t,i), gel(s1,i));
228 : 2903497 : return gerepilecopy(av, normalizepol_lg(a, N+2));
229 : : }
230 : : /* spec function. nP = lgpol(P) */
231 : : static GEN
232 : 2762086 : cook_square(GEN P, long nP)
233 : : {
234 : : GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
235 : 2762086 : long n0, n3, i, j, n = nP - 1;
236 : : pari_sp av;
237 : :
238 [ + - ][ + - ]: 2762086 : 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 : 2762086 : return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
303 : : }
304 : :
305 : : static GEN
306 : 1381043 : graeffe(GEN p)
307 : : {
308 : : GEN p0, p1, s0, s1;
309 : 1381043 : long n = degpol(p), n0, n1, i;
310 : :
311 [ - + ]: 1381043 : if (!n) return gcopy(p);
312 : 1381043 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
313 : 1381043 : p0 = new_chunk(n0);
314 : 1381043 : p1 = new_chunk(n1);
315 [ + + ]: 4631162 : for (i=0; i<n1; i++)
316 : : {
317 : 3250119 : p0[i] = p[2+(i<<1)];
318 : 3250119 : p1[i] = p[3+(i<<1)];
319 : : }
320 [ + + ]: 1381043 : if (n1 != n0)
321 : 677251 : p0[i] = p[2+(i<<1)];
322 : 1381043 : s0 = cook_square(p0, n0);
323 : 1381043 : s1 = cook_square(p1, n1);
324 : 1381043 : return RgX_sub(s0, RgX_shift_shallow(s1,1));
325 : : }
326 : : GEN
327 : 2468 : ZX_graeffe(GEN p)
328 : : {
329 : 2468 : pari_sp av = avma;
330 : : GEN p0, p1, s0, s1;
331 : 2468 : long n = degpol(p);
332 : :
333 [ - + ]: 2468 : if (!n) return ZX_copy(p);
334 : 2468 : RgX_even_odd(p, &p0, &p1);
335 : : /* p = p0(x^2) + x p1(x^2) */
336 : 2468 : s0 = ZX_sqr(p0);
337 : 2468 : s1 = ZX_sqr(p1);
338 : 2468 : return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
339 : : }
340 : : GEN
341 : 8 : polgraeffe(GEN p)
342 : : {
343 : 8 : pari_sp av = avma;
344 : : GEN p0, p1, s0, s1;
345 : 8 : long n = degpol(p);
346 : :
347 [ - + ]: 8 : if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
348 : 8 : n = degpol(p);
349 [ - + ]: 8 : if (!n) return gcopy(p);
350 : 8 : RgX_even_odd(p, &p0, &p1);
351 : : /* p = p0(x^2) + x p1(x^2) */
352 : 8 : s0 = RgX_sqr(p0);
353 : 8 : s1 = RgX_sqr(p1);
354 : 8 : 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 : 7716513 : mydbllog2i(GEN x)
368 : : {
369 : : #ifdef LONG_IS_64BIT
370 : 5806245 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
371 : : #else
372 : 1910268 : const double W = 1/4294967296.; /*2^-32*/
373 : : #endif
374 : : GEN m;
375 : 7716513 : long lx = lgefint(x);
376 : : double l;
377 [ + + ]: 7716513 : if (lx == 2) return -pariINFINITY;
378 : 7673264 : m = int_MSW(x);
379 : 7673264 : l = (double)(ulong)*m;
380 [ + + ]: 7673264 : if (lx == 3) return log2(l);
381 : 3768286 : 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 : 7716513 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
386 : : }
387 : :
388 : : /* return log(|x|) or -pariINFINITY */
389 : : static double
390 : 562286 : mydbllogr(GEN x) {
391 [ - + ]: 562286 : if (!signe(x)) return -pariINFINITY;
392 : 562286 : return LOG2*dbllog2r(x);
393 : : }
394 : :
395 : : /* return log2(|x|) or -pariINFINITY */
396 : : static double
397 : 3872311 : mydbllog2r(GEN x) {
398 [ + + ]: 3872311 : if (!signe(x)) return -pariINFINITY;
399 : 3872311 : return dbllog2r(x);
400 : : }
401 : : static double
402 [ + + ]: 1994472 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
403 : : double
404 : 10578066 : dbllog2(GEN z)
405 : : {
406 : : double x, y;
407 [ + + + ]: 10578066 : switch(typ(z))
408 : : {
409 : 6035405 : case t_INT: return mydbllog2i(z);
410 : 3545425 : case t_REAL: return mydbllog2r(z);
411 : : default: /*t_COMPLEX*/
412 : 997236 : x = dbllog2mp(gel(z,1));
413 : 997236 : y = dbllog2mp(gel(z,2));
414 [ + + ]: 997236 : if (fabs(x-y) > 10) return maxdd(x,y);
415 : 10578066 : return x + 0.5*log2(1 + exp2(2*(y-x)));
416 : : }
417 : : }
418 : : static GEN /* beware overflow */
419 [ + + ]: 213717 : 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 : 1059785 : findpower(GEN p)
425 : : {
426 : 1059785 : double x, L, mins = pariINFINITY;
427 : 1059785 : long n = degpol(p),i;
428 : :
429 : 1059785 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
430 [ + + ]: 5352203 : for (i=n-1; i>=0; i--)
431 : : {
432 : 4292418 : L += log2((double)(i+1) / (double)(n-i));
433 : 4292418 : x = dbllog2(gel(p,i+2));
434 [ + + ]: 4292418 : if (x != -pariINFINITY)
435 : : {
436 : 4267147 : double s = (L - x) / (double)(n-i);
437 [ + + ]: 4267147 : if (s < mins) mins = s;
438 : : }
439 : : }
440 : 1059785 : i = (long)ceil(mins);
441 [ + + ]: 1059785 : if (i - mins > 1 - 1e-12) i--;
442 : 1059785 : return i;
443 : : }
444 : :
445 : : /* returns the exponent for logmodulus(), from the Newton diagram */
446 : : static long
447 : 185171 : newton_polygon(GEN p, long k)
448 : : {
449 : 185171 : pari_sp av = avma;
450 : : double *logcoef, slope;
451 : 185171 : long n = degpol(p), i, j, h, l, *vertex;
452 : :
453 : 185171 : init_dalloc();
454 : 185171 : logcoef = (double*)stack_malloc((n+1)*sizeof(double));
455 : 185171 : vertex = (long*)new_chunk(n+1);
456 : :
457 : : /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
458 [ + + ]: 1019742 : for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
459 : 185171 : vertex[0] = 1; /* sentinel */
460 [ + + ]: 788159 : for (i=0; i < n; i=h)
461 : : {
462 : 602988 : slope = logcoef[i+1]-logcoef[i];
463 [ + + ]: 2872609 : for (j = h = i+1; j<=n; j++)
464 : : {
465 : 2269621 : double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
466 [ + + ]: 2269621 : if (slope < pij) { slope = pij; h = j; }
467 : : }
468 : 602988 : vertex[h] = 1;
469 : : }
470 [ + + ]: 198890 : h = k; while (!vertex[h]) h++;
471 [ + + ]: 194346 : l = k-1; while (!vertex[l]) l--;
472 : 185171 : avma = av;
473 : 185171 : 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 : 1267534 : myshiftrc(GEN z, long e)
479 : : {
480 [ + + ]: 1267534 : if (typ(z)==t_COMPLEX)
481 : : {
482 [ + + ]: 231682 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
483 [ + + ]: 231682 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
484 : : }
485 : : else
486 [ + + ]: 1035852 : if (signe(z)) shiftr_inplace(z, e);
487 : 1267534 : }
488 : :
489 : : /* return z*2^e, where z is integer or complex of integer (destroy z) */
490 : : static GEN
491 : 4605445 : myshiftic(GEN z, long e)
492 : : {
493 [ + + ]: 4605445 : if (typ(z)==t_COMPLEX)
494 : : {
495 [ + + ]: 660043 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
496 : 660043 : gel(z,2) = mpshift(gel(z,2),e);
497 : 660043 : return z;
498 : : }
499 [ + + ]: 4605445 : 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 : 218530 : myreal_1(long bit)
505 : : {
506 [ - + ]: 218530 : if (bit < 0) bit = 0;
507 : 218530 : return real_1(nbits2prec(bit));
508 : : }
509 : : static GEN
510 : 233557 : RgX_gtofp_bit(GEN q, long bit)
511 : : {
512 [ - + ]: 233557 : if (bit < 0) bit = 0;
513 : 233557 : return RgX_gtofp(q, nbits2prec(bit));
514 : : }
515 : :
516 : : static GEN
517 : 8380093 : mygprecrc(GEN x, long prec, long e)
518 : : {
519 : : GEN y;
520 [ + + + ]: 8380093 : switch(typ(x))
521 : : {
522 [ + + ]: 6574735 : case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
523 : : case t_COMPLEX:
524 : 1340569 : y = cgetg(3,t_COMPLEX);
525 : 1340569 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
526 : 1340569 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
527 : 1340569 : return y;
528 : 8380093 : 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 : 2286219 : mygprec(GEN x, long bit)
536 : : {
537 : : long lx, i, e, prec;
538 : : GEN y;
539 : :
540 [ - + ]: 2286219 : if (bit < 0) bit = 0; /* should rarely happen */
541 : 2286219 : e = gexpo(x) - bit;
542 : 2286219 : prec = nbits2prec(bit);
543 [ + + ]: 2286219 : switch(typ(x))
544 : : {
545 : : case t_POL:
546 : 1367810 : y = cgetg_copy(x, &lx); y[1] = x[1];
547 [ + + ]: 6106972 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
548 : 1367810 : break;
549 : :
550 : 918409 : default: y = mygprecrc(x,prec,e);
551 : : }
552 : 2286219 : 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 : 553833 : pol_to_gaussint(GEN p, long shift)
559 : : {
560 : 553833 : long i, l = lg(p);
561 : 553833 : GEN q = cgetg(l, t_POL); q[1] = p[1];
562 [ + + ]: 4003637 : for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
563 : 553833 : return q;
564 : : }
565 : :
566 : : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
567 : : static GEN
568 : 440371 : eval_rel_pol(GEN p, long bit)
569 : : {
570 : : long i;
571 [ + + ]: 3136886 : for (i = 2; i < lg(p); i++)
572 [ + + ]: 2696515 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
573 : 440371 : 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 : 49251 : homothetie(GEN p, double lrho, long bit)
579 : : {
580 : : GEN q, r, t, iR;
581 : 49251 : long n = degpol(p), i;
582 : :
583 : 49251 : iR = mygprec(dblexp(-lrho),bit);
584 : 49251 : q = mygprec(p, bit);
585 : 49251 : r = cgetg(n+3,t_POL); r[1] = p[1];
586 : 49251 : t = iR; r[n+2] = q[n+2];
587 [ + + ]: 315041 : for (i=n-1; i>0; i--)
588 : : {
589 : 265790 : gel(r,i+2) = gmul(t, gel(q,i+2));
590 : 265790 : t = mulrr(t, iR);
591 : : }
592 : 49251 : 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 : 298633 : homothetie2n(GEN p, long e)
598 : : {
599 [ + + ]: 298633 : if (e)
600 : : {
601 : 225415 : long i,n = lg(p)-1;
602 [ + + ]: 1492949 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
603 : : }
604 : 298633 : }
605 : :
606 : : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
607 : : static void
608 : 946323 : homothetie_gauss(GEN p, long e, long f)
609 : : {
610 [ + + ][ + + ]: 946323 : if (e || f)
611 : : {
612 : 866697 : long i, n = lg(p)-1;
613 [ + + ]: 5472142 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
614 : : }
615 : 946323 : }
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 : 1059785 : lower_bound(GEN p, long *k, double eps)
621 : : {
622 : 1059785 : long n = degpol(p), i, j;
623 : 1059785 : pari_sp ltop = avma;
624 : : GEN a, s, S, ilc;
625 : : double r, R, rho;
626 : :
627 [ + + ]: 1059785 : if (n < 4) { *k = n; return 0.; }
628 : 483842 : S = cgetg(5,t_VEC);
629 : 483842 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
630 [ + + ]: 2419210 : 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 : 483842 : s = gel(a,1);
633 : 483842 : gel(S,1) = gneg(s); /* Newton sum S_i */
634 : 483842 : rho = r = gtodouble(gabs(s,3));
635 : 483842 : R = r / n;
636 [ + + ]: 1935368 : for (i=2; i<=4; i++)
637 : : {
638 : 1451526 : s = gmulsg(i,gel(a,i));
639 [ + + ]: 4354578 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
640 : 1451526 : gel(S,i) = gneg(s); /* Newton sum S_i */
641 : 1451526 : r = gtodouble(gabs(s,3));
642 [ + + ]: 1451526 : if (r > 0.)
643 : : {
644 : 1449529 : r = exp(log(r/n) / (double)i);
645 [ + + ]: 1449529 : if (r > R) R = r;
646 : : }
647 : : }
648 [ + + ][ + + ]: 483842 : if (R > 0. && eps < 1.2)
649 : 483677 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
650 : : else
651 : 165 : *k = n;
652 : 1059785 : avma = ltop; return R;
653 : : }
654 : :
655 : : /* log of modulus of the largest root of p with relative error tau */
656 : : static double
657 : 113462 : logmax_modulus(GEN p, double tau)
658 : : {
659 : : GEN r,q,aux,gunr;
660 : 113462 : pari_sp av, ltop = avma;
661 : 113462 : long i,k,n=degpol(p),nn,bit,M,e;
662 [ - + ]: 113462 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
663 : :
664 : 113462 : r = cgeti(BIGDEFAULTPREC);
665 : 113462 : av = avma;
666 : :
667 : 113462 : eps = - 1/log(1.5*tau2); /* > 0 */
668 : 113462 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
669 : 113462 : gunr = myreal_1(bit+2*n);
670 : 113462 : aux = gdiv(gunr, gel(p,2+n));
671 : 113462 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
672 : 113462 : e = findpower(q);
673 : 113462 : homothetie2n(q,e);
674 : 113462 : affsi(e, r);
675 : 113462 : q = pol_to_gaussint(q, bit);
676 : 113462 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
677 : 113462 : nn = n;
678 : 113462 : for (i=0,e=0;;)
679 : : { /* nn = deg(q) */
680 : 1059785 : rho = lower_bound(q, &k, eps);
681 [ + + ]: 1059785 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
682 : 1059785 : affii(shifti(addis(r,e), 1), r);
683 [ + + ]: 1059785 : if (++i == M) break;
684 : :
685 : 2838969 : bit = (long) ((double)k * log2(1./tau2) +
686 : 1892646 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
687 : 946323 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
688 : 946323 : nn -= RgX_valrem(q, &q);
689 : 946323 : set_karasquare_limit(gexpo(q));
690 : 946323 : q = gerepileupto(av, graeffe(q));
691 [ + + ]: 946323 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
692 : 946323 : eps = -1/log(tau2); /* > 0 */
693 : 946323 : e = findpower(q);
694 : 946323 : }
695 [ + + ]: 113462 : if (!signe(r)) { avma = ltop; return 0.; }
696 : 105560 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
697 : 113462 : avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
698 : : }
699 : : GEN
700 : 240 : logmax_modulus_bound(GEN P)
701 : : {
702 : 240 : return dblexp(logmax_modulus(P, 0.01) + 0.01);
703 : : }
704 : :
705 : : /* log of modulus of the smallest root of p, with relative error tau */
706 : : static double
707 : 39733 : logmin_modulus(GEN p, double tau)
708 : : {
709 : 39733 : pari_sp av = avma;
710 : : double r;
711 : :
712 [ - + ]: 39733 : if (gequal0(gel(p,2))) return -pariINFINITY;
713 : 39733 : r = - logmax_modulus(RgX_recip_shallow(p),tau);
714 : 39733 : avma = av; return r;
715 : : }
716 : :
717 : : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
718 : : static double
719 : 20405 : logmodulus(GEN p, long k, double tau)
720 : : {
721 : : GEN q;
722 : 20405 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
723 : 20405 : pari_sp av, ltop=avma;
724 : 20405 : double r, tau2 = tau/6;
725 : :
726 : 20405 : bit = (long)(n * (2. + log2(3.*n/tau2)));
727 : 20405 : av = avma;
728 : 20405 : q = gprec_w(p, nbits2prec(bit));
729 : 20405 : q = RgX_gtofp_bit(q, bit);
730 : 20405 : e = newton_polygon(q,k);
731 : 20405 : r = (double)e;
732 : 20405 : homothetie2n(q,e);
733 : 20405 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
734 [ + + ]: 185171 : for (i=1; i<imax; i++)
735 : : {
736 : 164766 : q = eval_rel_pol(q,bit);
737 : 164766 : kk -= RgX_valrem(q, &q);
738 : 164766 : nn = degpol(q);
739 : :
740 : 164766 : set_karasquare_limit(bit);
741 : 164766 : q = gerepileupto(av, graeffe(q));
742 : 164766 : e = newton_polygon(q,kk);
743 : 164766 : r += e / exp2((double)i);
744 : 164766 : q = RgX_gtofp_bit(q, bit);
745 : 164766 : homothetie2n(q,e);
746 : :
747 [ - + ]: 164766 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
748 : 164766 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
749 : : }
750 : 20405 : avma = ltop; return -r * LOG2;
751 : : }
752 : :
753 : : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
754 : : * rmin < r_k < rmax. This information helps because we may reduce precision
755 : : * quicker */
756 : : static double
757 : 20405 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
758 : : {
759 : : GEN q;
760 : 20405 : long n = degpol(p), i, imax, imax2, bit;
761 : 20405 : pari_sp ltop = avma, av;
762 : 20405 : double lrho, aux, tau2 = tau/6.;
763 : :
764 : 20405 : aux = (lrmax - lrmin) / 2. + 4*tau2;
765 : 20405 : imax = (long) log2(log((double)n)/ aux);
766 [ + + ]: 20405 : if (imax <= 0) return logmodulus(p,k,tau);
767 : :
768 : 19819 : lrho = (lrmin + lrmax) / 2;
769 : 19819 : av = avma;
770 : 19819 : bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
771 : 19819 : q = homothetie(p, lrho, bit);
772 : 19819 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
773 [ - + ]: 19819 : if (imax > imax2) imax = imax2;
774 : :
775 [ + + ]: 61733 : for (i=0; i<imax; i++)
776 : : {
777 : 41914 : q = eval_rel_pol(q,bit);
778 : 41914 : set_karasquare_limit(bit);
779 : 41914 : q = gerepileupto(av, graeffe(q));
780 : 41914 : aux = 2*aux + 2*tau2;
781 : 41914 : tau2 *= 1.5;
782 : 41914 : bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
783 : 41914 : q = RgX_gtofp_bit(q, bit);
784 : : }
785 : 19819 : aux = exp2((double)imax);
786 : 19819 : aux = logmodulus(q,k, aux*tau/3.) / aux;
787 : 20405 : avma = ltop; return lrho + aux;
788 : : }
789 : :
790 : : static double
791 : 23781 : ind_maxlog2(GEN q)
792 : : {
793 : 23781 : long i, k = -1;
794 : 23781 : double L = - pariINFINITY;
795 [ + + ]: 66783 : for (i=0; i<=degpol(q); i++)
796 : : {
797 : 43002 : double d = dbllog2(gel(q,2+i));
798 [ + + ]: 43002 : if (d > L) { L = d; k = i; }
799 : : }
800 : 23781 : return k;
801 : : }
802 : :
803 : : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
804 : : * Assume that l <= k <= n-l */
805 : : static long
806 : 29432 : dual_modulus(GEN p, double lrho, double tau, long l)
807 : : {
808 : 29432 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
809 : 29432 : double tau2 = tau * 7./8.;
810 : 29432 : pari_sp av = avma;
811 : : GEN q;
812 : :
813 : 29432 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
814 : 29432 : q = homothetie(p, lrho, bit);
815 : 29432 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
816 : :
817 [ + + ]: 257472 : for (i=0; i<imax; i++)
818 : : {
819 : 233691 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
820 : 233691 : v = RgX_valrem(q, &q);
821 [ + + ]: 233691 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
822 : :
823 : 233691 : nn = degpol(q); delta_k += v;
824 [ + + ]: 233691 : if (!nn) return delta_k;
825 : :
826 : 228040 : set_karasquare_limit(bit);
827 : 228040 : q = gerepileupto(av, graeffe(q));
828 : 228040 : tau2 *= 7./4.;
829 : 228040 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
830 : : }
831 : 29432 : avma = av; return delta_k + (long)ind_maxlog2(q);
832 : : }
833 : :
834 : : /********************************************************************/
835 : : /** **/
836 : : /** FACTORS THROUGH CIRCLE INTEGRATION **/
837 : : /** **/
838 : : /********************************************************************/
839 : : /* l power of 2 */
840 : : static void
841 : 916404 : fft(GEN Omega, GEN p, GEN f, long step, long l)
842 : : {
843 : : pari_sp ltop;
844 : : long i, l1, l2, l3, rapi, step4;
845 : : GEN f1, f2, f3, f02, f13, g02, g13, ff;
846 : :
847 [ + + ]: 916404 : if (l == 2)
848 : : {
849 : 524856 : gel(f,0) = gadd(gel(p,0),gel(p,step));
850 : 524856 : gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
851 : : }
852 [ + + ]: 391548 : if (l == 4)
853 : : {
854 : 204043 : f1 = gadd(gel(p,0), gel(p,step<<1));
855 : 204043 : f2 = gsub(gel(p,0), gel(p,step<<1));
856 : 204043 : f3 = gadd(gel(p,step),gel(p,3*step));
857 : 204043 : f02= gsub(gel(p,step),gel(p,3*step));
858 : 204043 : f02 = mulcxI(f02);
859 : 204043 : gel(f,0) = gadd(f1, f3);
860 : 204043 : gel(f,1) = gadd(f2, f02);
861 : 204043 : gel(f,2) = gsub(f1, f3);
862 : 204043 : gel(f,3) = gsub(f2, f02); return;
863 : : }
864 : :
865 : 187505 : ltop = avma;
866 : 187505 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
867 : 187505 : fft(Omega,p, f, step4,l1);
868 : 187505 : fft(Omega,p+step, f+l1,step4,l1);
869 : 187505 : fft(Omega,p+(step<<1),f+l2,step4,l1);
870 : 187505 : fft(Omega,p+3*step, f+l3,step4,l1);
871 : :
872 : 187505 : ff = cgetg(l+1,t_VEC);
873 [ + + ]: 755337 : for (i=0; i<l1; i++)
874 : : {
875 : 567832 : rapi = step*i;
876 : 567832 : f1 = gmul(gel(Omega,rapi), gel(f,i+l1));
877 : 567832 : f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
878 : 567832 : f3 = gmul(gel(Omega,3*rapi), gel(f,i+l3));
879 : :
880 : 567832 : f02 = gadd(gel(f,i),f2);
881 : 567832 : g02 = gsub(gel(f,i),f2);
882 : 567832 : f13 = gadd(f1,f3);
883 : 567832 : g13 = mulcxI(gsub(f1,f3));
884 : :
885 : 567832 : gel(ff,i+1) = gadd(f02, f13);
886 : 567832 : gel(ff,i+l1+1) = gadd(g02, g13);
887 : 567832 : gel(ff,i+l2+1) = gsub(f02, f13);
888 : 567832 : gel(ff,i+l3+1) = gsub(g02, g13);
889 : : }
890 : 187505 : ff = gerepilecopy(ltop,ff);
891 [ + + ]: 3187732 : for (i=0; i<l; i++) f[i] = ff[i+1];
892 : : }
893 : :
894 : : /* e(1/N) */
895 : : static GEN
896 : 58817 : RUgen(long N, long bit)
897 : : {
898 [ + + ]: 58817 : if (N == 2) return real_m1(nbits2prec(bit));
899 [ + + ]: 56819 : if (N == 4) return gen_I();
900 : 58817 : return expIr(divru(Pi2n(1, nbits2prec(bit)), N));
901 : : }
902 : :
903 : : /* N=2^k. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
904 : : static GEN
905 : 28106 : initRU(long N, long bit)
906 : : {
907 : 28106 : GEN *RU, z = RUgen(N, bit);
908 : 28106 : long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
909 : :
910 : 28106 : RU = (GEN*)cgetg(N+1,t_VEC); RU++;
911 : :
912 : 28106 : RU[0] = myreal_1(bit);
913 : 28106 : RU[1] = z;
914 [ + + ]: 38822 : for (i=1; i<N8; i++)
915 : : {
916 : 10716 : GEN t = RU[i];
917 : 10716 : RU[i+1] = gmul(z, t);
918 : 10716 : RU[N4-i] = mkcomplex(gel(t,2), gel(t,1));
919 : : }
920 [ + + ]: 98324 : for (i=0; i<N4; i++) RU[i+N4] = mulcxI(RU[i]);
921 [ + + ]: 168542 : for (i=0; i<N2; i++) RU[i+N2] = gneg(RU[i]);
922 : 28106 : return (GEN)RU;
923 : : }
924 : :
925 : : /* as above, N arbitrary */
926 : : static GEN
927 : 2605 : initRUgen(long N, long bit)
928 : : {
929 : 2605 : GEN *RU = (GEN*)cgetg(N+1,t_VEC), z = RUgen(N,bit);
930 : 2605 : long i, k = (N+3)>>1;
931 : 2605 : RU[0] = gen_1;
932 : 2605 : RU[1] = z;
933 [ + + ]: 3770 : for (i=2; i<k; i++) RU[i] = gmul(z, RU[i-1]);
934 [ + + ]: 3460 : for ( ; i<N; i++) RU[i] = gconj(RU[N-i]);
935 : 2605 : return (GEN)RU;
936 : : }
937 : :
938 : : GEN
939 : 0 : FFTinit(long k, long prec)
940 : : {
941 [ # # ]: 0 : if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
942 : 0 : return initRU(1L << k, prec2nbits(prec)) - 1;
943 : : }
944 : :
945 : : GEN
946 : 0 : FFT(GEN x, GEN Omega)
947 : : {
948 : 0 : long i, l = lg(Omega), n = lg(x);
949 : : GEN y, z;
950 [ # # ]: 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
951 [ # # ]: 0 : if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
952 [ # # ]: 0 : if (n > l) pari_err_DIM("FFT");
953 : :
954 [ # # ]: 0 : if (n < l) {
955 : 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
956 [ # # ]: 0 : for (i = 1; i < n; i++) z[i] = x[i];
957 [ # # ]: 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
958 : : }
959 : 0 : else z = x;
960 : 0 : y = cgetg(l, t_VEC);
961 : 0 : fft(Omega+1, z+1, y+1, 1, l-1);
962 : 0 : return y;
963 : : }
964 : :
965 : : /* returns 1 if p has only real coefficients, 0 else */
966 : : static int
967 : 25736 : isreal(GEN p)
968 : : {
969 : : long i;
970 [ + + ]: 168642 : for (i = lg(p)-1; i > 1; i--)
971 [ + + ]: 148489 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
972 : 25736 : return 1;
973 : : }
974 : :
975 : : /* x non complex */
976 : : static GEN
977 : 20106 : abs_update_r(GEN x, double *mu) {
978 : 20106 : GEN y = gtofp(x, DEFAULTPREC);
979 [ + + ]: 20106 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
980 : 20106 : setabssign(y); return y;
981 : : }
982 : : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
983 : : static GEN
984 : 383036 : abs_update(GEN x, double *mu) {
985 : : GEN y, xr, yr;
986 : : double ly;
987 [ + + ]: 383036 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
988 : 363380 : xr = gel(x,1);
989 : 363380 : yr = gel(x,2);
990 [ + + ]: 363380 : if (gequal0(xr)) return abs_update_r(yr,mu);
991 [ + + ]: 363324 : if (gequal0(yr)) return abs_update_r(xr,mu);
992 : : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
993 : 362930 : xr = gtofp(xr, DEFAULTPREC);
994 : 362930 : yr = gtofp(yr, DEFAULTPREC);
995 : 362930 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
996 [ + + ]: 362930 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
997 : 383036 : return y;
998 : : }
999 : :
1000 : : static void
1001 : 13522 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
1002 : : int polreal, double param, double param2)
1003 : : {
1004 : : GEN q, pc, Omega, A, RU, prim, g, ONE,TWO;
1005 : 13522 : long n = degpol(p), bit, NN, K, i, j, Lmax;
1006 : 13522 : pari_sp av2, av = avma, lim = stack_lim(av, 1);
1007 : :
1008 : 13522 : bit = gexpo(p) + (long)param2+8;
1009 [ + + ]: 27038 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
1010 [ + + ]: 13522 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
1011 [ + + ]: 13522 : K = NN/Lmax; if (K & 1) K++;
1012 : 13522 : NN = Lmax*K;
1013 [ + + ]: 13522 : if (polreal) K = K/2+1;
1014 : :
1015 : 13522 : Omega = initRU(Lmax,bit);
1016 : 13522 : prim = RUgen(NN, bit);
1017 : :
1018 : 13522 : q = mygprec(p,bit) + 2;
1019 : 13522 : A = cgetg(Lmax+1,t_VEC); A++;
1020 : 13522 : pc= cgetg(Lmax+1,t_VEC); pc++;
1021 [ + + ]: 107810 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
1022 [ + + ]: 48326 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
1023 : :
1024 : 13522 : *mu = pariINFINITY;
1025 : 13522 : g = real_0_bit(-bit);
1026 : 13522 : ONE = real_1(DEFAULTPREC);
1027 : 13522 : TWO = real2n(1, DEFAULTPREC);
1028 : 13522 : av2 = avma;
1029 : 13522 : RU = myreal_1(bit);
1030 [ + + ]: 51114 : for (i=0; i<K; i++)
1031 : : {
1032 [ + + ]: 37592 : if (i) {
1033 : 24070 : GEN z = RU;
1034 [ + + ]: 168430 : for (j=1; j<n; j++)
1035 : : {
1036 : 144360 : gel(pc,j) = gmul(gel(q,j),z);
1037 : 144360 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
1038 : : }
1039 : 24070 : gel(pc,n) = gmul(gel(q,n),z);
1040 : : }
1041 : :
1042 : 37592 : fft(Omega,pc,A,1,Lmax);
1043 [ + + ][ + + ]: 43396 : if (polreal && i>0 && i<K-1)
[ + + ]
1044 [ + + ]: 91280 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
1045 : : else
1046 [ + + ]: 329348 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(ONE, abs_update(gel(A,j),mu)));
1047 : 37592 : RU = gmul(RU, prim);
1048 [ - + ]: 37592 : if (low_stack(lim, stack_lim(av,1)))
1049 : : {
1050 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
1051 : 0 : gerepileall(av2,2, &g,&RU);
1052 : : }
1053 : : }
1054 : 13522 : *gamma = mydbllog2r(divru(g,NN));
1055 : 13522 : *LMAX = Lmax; avma = av;
1056 : 13522 : }
1057 : :
1058 : : /* NN is a multiple of Lmax */
1059 : : static void
1060 : 14584 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
1061 : : {
1062 : : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
1063 : 14584 : long n = degpol(p), i, j, K;
1064 : : pari_sp ltop;
1065 : :
1066 : 14584 : Omega = initRU(Lmax,bit);
1067 : 14584 : prim = RUgen(NN, bit);
1068 : 14584 : RU = cgetg(n+2,t_VEC); RU++;
1069 : :
1070 [ + + ]: 14584 : K = NN/Lmax; if (polreal) K = K/2+1;
1071 : 14584 : q = mygprec(p,bit);
1072 : 14584 : qd = RgX_deriv(q);
1073 : :
1074 : 14584 : A = cgetg(Lmax+1,t_VEC); A++;
1075 : 14584 : B = cgetg(Lmax+1,t_VEC); B++;
1076 : 14584 : C = cgetg(Lmax+1,t_VEC); C++;
1077 : 14584 : pc = cgetg(Lmax+1,t_VEC); pc++;
1078 : 14584 : pd = cgetg(Lmax+1,t_VEC); pd++;
1079 [ + + ]: 55599 : pc[0] = q[2]; for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
1080 [ + + ]: 70183 : pd[0] = qd[2]; for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
1081 : :
1082 : 14584 : ltop = avma;
1083 : 14584 : W = cgetg(k+1,t_VEC);
1084 : 14584 : U = cgetg(k+1,t_VEC);
1085 [ + + ]: 43621 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
1086 : :
1087 : 14584 : gel(RU,0) = gen_1;
1088 : 14584 : prim2 = myreal_1(bit);
1089 [ + + ]: 46782 : for (i=0; i<K; i++)
1090 : : {
1091 : 32198 : gel(RU,1) = prim2;
1092 [ + + ]: 237381 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
1093 : : /* RU[j] = prim^(ij)= prim2^j */
1094 : :
1095 [ + + ]: 237381 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
1096 : 32198 : fft(Omega,pd,A,1,Lmax);
1097 [ + + ]: 269579 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
1098 : 32198 : fft(Omega,pc,B,1,Lmax);
1099 [ + + ]: 402910 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
1100 [ + + ]: 402910 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
1101 : 32198 : fft(Omega,B,A,1,Lmax);
1102 : 32198 : fft(Omega,C,B,1,Lmax);
1103 : :
1104 [ + + ]: 32198 : if (polreal) /* p has real coefficients */
1105 : : {
1106 [ + + ][ + + ]: 26990 : if (i>0 && i<K-1)
1107 : : {
1108 [ + + ]: 13459 : for (j=1; j<=k; j++)
1109 : : {
1110 : 10812 : gel(W,j) = gadd(gel(W,j), gshift(real_i(gmul(gel(A,j+1),gel(RU,j+1))),1));
1111 : 10812 : gel(U,j) = gadd(gel(U,j), gshift(real_i(gmul(gel(B,j),gel(RU,j))),1));
1112 : : }
1113 : : }
1114 : : else
1115 : : {
1116 [ + + ]: 70105 : for (j=1; j<=k; j++)
1117 : : {
1118 : 45762 : gel(W,j) = gadd(gel(W,j), real_i(gmul(gel(A,j+1),gel(RU,j+1))));
1119 : 45762 : gel(U,j) = gadd(gel(U,j), real_i(gmul(gel(B,j),gel(RU,j))));
1120 : : }
1121 : : }
1122 : : }
1123 : : else
1124 : : {
1125 [ + + ]: 22444 : for (j=1; j<=k; j++)
1126 : : {
1127 : 14589 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
1128 : 14589 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
1129 : : }
1130 : : }
1131 : 32198 : prim2 = gmul(prim2,prim);
1132 : 32198 : gerepileall(ltop,3, &W,&U,&prim2);
1133 : : }
1134 : :
1135 [ + + ]: 43621 : for (i=1; i<=k; i++)
1136 : : {
1137 : 29037 : aux=gel(W,i);
1138 [ + + ]: 66119 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
1139 : 29037 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
1140 : : }
1141 [ + + ]: 43621 : for (i=0; i<k; i++)
1142 : : {
1143 : 29037 : aux=gel(U,k-i);
1144 [ + + ]: 66119 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
1145 : 29037 : gel(H,i+2) = gdivgs(aux,NN);
1146 : : }
1147 : 14584 : }
1148 : :
1149 : : #define NEWTON_MAX 10
1150 : : static GEN
1151 : 73919 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
1152 : : {
1153 : 73919 : GEN H = HH, D, aux;
1154 : 73919 : pari_sp ltop = avma, lim = stack_lim(ltop, 1);
1155 : : long error, i, bit1, bit2;
1156 : :
1157 : 73919 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
1158 : 73919 : bit2 = bit + Sbit;
1159 [ + + ][ + + ]: 137645 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1160 : : {
1161 [ - + ]: 63726 : if (low_stack(lim, stack_lim(ltop,1)))
1162 : : {
1163 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
1164 : 0 : gerepileall(ltop,2, &D,&H);
1165 : : }
1166 : 63726 : bit1 = -error + Sbit;
1167 : 63726 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
1168 : 63726 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
1169 : :
1170 [ + + ]: 63726 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
1171 : 63726 : H = RgX_add(mygprec(H,bit1), aux);
1172 : 63726 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
1173 [ + + ]: 63726 : error = gexpo(D); if (error < -bit1) error = -bit1;
1174 : : }
1175 [ + + ]: 73919 : if (error > -bit/2) return NULL; /* FAIL */
1176 : 73919 : return gerepilecopy(ltop,H);
1177 : : }
1178 : :
1179 : : /* return 0 if fails, 1 else */
1180 : : static long
1181 : 14584 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
1182 : : {
1183 : 14584 : GEN f0, FF, GG, r, HH = H;
1184 : 14584 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
1185 : 14584 : pari_sp av = avma, lim = stack_lim(av, 1);
1186 : :
1187 : 14584 : FF = *F; GG = RgX_divrem(p, FF, &r);
1188 [ - + ]: 14584 : error = gexpo(r); if (error <= -bit) error = 1-bit;
1189 : 14584 : normF = gexpo(FF);
1190 : 14584 : normG = gexpo(GG);
1191 [ + + ]: 14584 : enh = gexpo(H); if (enh < 0) enh = 0;
1192 : 14584 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
1193 : 14584 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
1194 : 14584 : bit2 = bit + Sbit;
1195 [ + + ][ + + ]: 88470 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
[ + + ]
1196 : : {
1197 [ + + ][ + + ]: 73919 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
1198 [ - + ]: 73919 : if (low_stack(lim, stack_lim(av,1)))
1199 : : {
1200 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1201 : 0 : gerepileall(av,4, &FF,&GG,&r,&HH);
1202 : : }
1203 : :
1204 : 73919 : bit1 = -error + Sbit2;
1205 : 73919 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1206 : : 1-error, Sbit2);
1207 [ + + ]: 73919 : if (!HH) return 0; /* FAIL */
1208 : :
1209 : 73886 : bit1 = -error + Sbit;
1210 : 73886 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1211 : 73886 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1212 : :
1213 [ + + ]: 73886 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1214 : 73886 : FF = gadd(mygprec(FF,bit1),f0);
1215 : :
1216 [ + + ]: 73886 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1217 : 73886 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1218 [ + + ]: 73886 : error = gexpo(r); if (error < -bit1) error = -bit1;
1219 : : }
1220 [ + + ]: 14551 : if (error>-bit) return 0; /* FAIL */
1221 : 14584 : *F = FF; *G = GG; return 1;
1222 : : }
1223 : :
1224 : : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1225 : : where cd is the leading coefficient of p */
1226 : : static void
1227 : 13522 : split_fromU(GEN p, long k, double delta, long bit,
1228 : : GEN *F, GEN *G, double param, double param2)
1229 : : {
1230 : : GEN pp, FF, GG, H;
1231 : 13522 : long n = degpol(p), NN, bit2, Lmax;
1232 : 13522 : int polreal = isreal(p);
1233 : : pari_sp ltop;
1234 : : double mu, gamma;
1235 : :
1236 : 13522 : pp = gdiv(p, gel(p,2+n));
1237 : 13522 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1238 : :
1239 : 13522 : H = cgetg(k+2,t_POL); H[1] = p[1];
1240 : 13522 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1241 : 13522 : gel(FF,k+2) = gen_1;
1242 : :
1243 [ + + ]: 13522 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1244 : 13522 : NN *= Lmax; ltop = avma;
1245 : : for(;;)
1246 : : {
1247 : 14584 : bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
1248 : 14584 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1249 [ + + ]: 14584 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1250 : 1062 : NN <<= 1; avma = ltop;
1251 : 1062 : }
1252 : 13522 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1253 : 13522 : }
1254 : :
1255 : : static void
1256 : 13522 : optimize_split(GEN p, long k, double delta, long bit,
1257 : : GEN *F, GEN *G, double param, double param2)
1258 : : {
1259 : 13522 : long n = degpol(p);
1260 : : GEN FF, GG;
1261 : :
1262 [ + + ]: 13522 : if (k <= n/2)
1263 : 10235 : split_fromU(p,k,delta,bit,F,G,param,param2);
1264 : : else
1265 : : {
1266 : 3287 : split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1267 : 3287 : *F = RgX_recip_shallow(GG);
1268 : 3287 : *G = RgX_recip_shallow(FF);
1269 : : }
1270 : 13522 : }
1271 : :
1272 : : /********************************************************************/
1273 : : /** **/
1274 : : /** SEARCH FOR SEPARATING CIRCLE **/
1275 : : /** **/
1276 : : /********************************************************************/
1277 : :
1278 : : /* return p(2^e*x) *2^(-n*e) */
1279 : : static void
1280 : 0 : scalepol2n(GEN p, long e)
1281 : : {
1282 : 0 : long i,n=lg(p)-1;
1283 [ # # ]: 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1284 : 0 : }
1285 : :
1286 : : /* returns p(x/R)*R^n */
1287 : : static GEN
1288 : 113850 : scalepol(GEN p, GEN R, long bit)
1289 : : {
1290 : : GEN q,aux,gR;
1291 : : long i;
1292 : :
1293 : 113850 : aux = gR = mygprec(R,bit); q = mygprec(p,bit);
1294 [ + + ]: 544190 : for (i=lg(p)-2; i>=2; i--)
1295 : : {
1296 : 430340 : gel(q,i) = gmul(aux,gel(q,i));
1297 : 430340 : aux = gmul(aux,gR);
1298 : : }
1299 : 113850 : return q;
1300 : : }
1301 : :
1302 : : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1303 : : static GEN
1304 : 36642 : conformal_pol(GEN p, GEN a, long bit)
1305 : : {
1306 : 36642 : GEN z, r, ma = gneg(a), ca = gconj(a);
1307 : 36642 : long n = degpol(p), i;
1308 : 36642 : pari_sp av = avma, lim = stack_lim(av,2);
1309 : :
1310 : 36642 : z = mkpoln(2, ca, negr(myreal_1(bit)));
1311 : 36642 : r = scalarpol(gel(p,2+n), 0);
1312 : 36642 : for (i=n-1; ; i--)
1313 : : {
1314 : 134404 : r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
1315 : 134404 : r = gadd(r, gmul(z, gel(p,2+i)));
1316 [ + + ]: 134404 : if (i == 0) return gerepileupto(av, r);
1317 : 97762 : z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1318 [ - + ]: 97762 : if (low_stack(lim, stack_lim(av,2)))
1319 : : {
1320 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
1321 : 0 : gerepileall(av,2, &r,&z);
1322 : : }
1323 : 97762 : }
1324 : : }
1325 : :
1326 : : static const double UNDEF = -100000.;
1327 : :
1328 : : static double
1329 : 13522 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1330 : : {
1331 : 13522 : long i, n = degpol(p);
1332 : : double lrho, lrmin, lrmax;
1333 [ + + ]: 13522 : if (k > 1)
1334 : : {
1335 [ + - ][ + + ]: 14292 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1336 : 8962 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1337 : : }
1338 : : else /* k=1 */
1339 : 4560 : lrmin = logmin_modulus(p,aux);
1340 : 13522 : radii[k] = lrmin;
1341 : :
1342 [ + + ]: 13522 : if (k+1<n)
1343 : : {
1344 [ + - ][ + + ]: 26020 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1345 : 11443 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1346 : : }
1347 : : else /* k+1=n */
1348 : 2079 : lrmax = logmax_modulus(p,aux);
1349 : 13522 : radii[k+1] = lrmax;
1350 : :
1351 : 13522 : lrho = radii[k];
1352 [ + + ]: 32827 : for (i=k-1; i>=1; i--)
1353 : : {
1354 [ + + ][ + + ]: 19305 : if (radii[i] == UNDEF || radii[i] > lrho)
1355 : 13541 : radii[i] = lrho;
1356 : : else
1357 : 5764 : lrho = radii[i];
1358 : : }
1359 : 13522 : lrho = radii[k+1];
1360 [ + + ]: 61461 : for (i=k+1; i<=n; i++)
1361 : : {
1362 [ + + ][ + + ]: 47939 : if (radii[i] == UNDEF || radii[i] < lrho)
1363 : 29333 : radii[i] = lrho;
1364 : : else
1365 : 18606 : lrho = radii[i];
1366 : : }
1367 : 13522 : *delta = (lrmax - lrmin) / 2;
1368 [ + + ]: 13522 : if (*delta > 1.) *delta = 1.;
1369 : 13522 : return (lrmin + lrmax) / 2;
1370 : : }
1371 : :
1372 : : static void
1373 : 13522 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1374 : : {
1375 : 13522 : double t, param = 0., param2 = 0.;
1376 : : long i;
1377 [ + + ]: 94288 : for (i=1; i<=n; i++)
1378 : : {
1379 : 80766 : radii[i] -= lrho;
1380 : 80766 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1381 [ + + ]: 80766 : param += t; if (t > 1.) param2 += log2(t);
1382 : : }
1383 : 13522 : *par = param; *par2 = param2;
1384 : 13522 : }
1385 : :
1386 : : /* apply the conformal mapping then split from U */
1387 : : static void
1388 : 12214 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1389 : : double aux, GEN *F,GEN *G)
1390 : : {
1391 : 12214 : long bit2, n = degpol(p), i;
1392 : 12214 : pari_sp ltop = avma, av;
1393 : : GEN q, FF, GG, a, R;
1394 : : double lrho, delta, param, param2;
1395 : : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1396 : 12214 : bit2 = bit + (long)(n*3.4848775) + 1;
1397 : 12214 : a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
1398 : 12214 : a = divrs(a, -6);
1399 : 12214 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1400 : :
1401 : 12214 : av = avma;
1402 : 12214 : q = conformal_pol(mygprec(p,bit2), a, bit2);
1403 [ + + ]: 79416 : for (i=1; i<=n; i++)
1404 [ + + ]: 67202 : if (radii[i] != UNDEF) /* update array radii */
1405 : : {
1406 : 45510 : pari_sp av2 = avma;
1407 : 45510 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1408 : : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1409 : 45510 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1410 : 45510 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1411 : 45510 : avma = av2;
1412 : : }
1413 : 12214 : lrho = logradius(radii, q,k,aux/10., &delta);
1414 : 12214 : update_radius(n, radii, lrho, ¶m, ¶m2);
1415 : :
1416 : 12214 : bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
1417 : 12214 : R = mygprec(dblexp(-lrho), bit2);
1418 : 12214 : q = scalepol(q,R,bit2);
1419 : 12214 : gerepileall(av,2, &q,&R);
1420 : :
1421 : 12214 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1422 : 12214 : bit2 += n; R = invr(R);
1423 : 12214 : FF = scalepol(FF,R,bit2);
1424 : 12214 : GG = scalepol(GG,R,bit2);
1425 : :
1426 : 12214 : a = mygprec(a,bit2);
1427 : 12214 : FF = conformal_pol(FF,a,bit2);
1428 : 12214 : GG = conformal_pol(GG,a,bit2);
1429 : :
1430 : 12214 : a = invr(subsr(1, gnorm(a)));
1431 : 12214 : FF = RgX_Rg_mul(FF, powru(a,k));
1432 : 12214 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1433 : :
1434 : 12214 : *F = mygprec(FF,bit+n);
1435 : 12214 : *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1436 : 12214 : }
1437 : :
1438 : : /* split p, this time without scaling. returns in F and G two polynomials
1439 : : * such that |p-FG|< 2^(-bit)|p| */
1440 : : static void
1441 : 13522 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1442 : : {
1443 : : GEN q, FF, GG, R;
1444 : : double aux, delta, param, param2;
1445 : 13522 : long n = degpol(p), i, j, k, bit2;
1446 : : double lrmin, lrmax, lrho, *radii;
1447 : :
1448 : 13522 : init_dalloc();
1449 : 13522 : radii = (double*) stack_malloc((n+1) * sizeof(double));
1450 : :
1451 [ + + ]: 67244 : for (i=2; i<n; i++) radii[i] = UNDEF;
1452 : 13522 : aux = thickness/(double)(4 * n);
1453 : 13522 : lrmin = logmin_modulus(p, aux);
1454 : 13522 : lrmax = logmax_modulus(p, aux);
1455 : 13522 : radii[1] = lrmin;
1456 : 13522 : radii[n] = lrmax;
1457 : 13522 : i = 1; j = n;
1458 : 13522 : lrho = (lrmin + lrmax) / 2;
1459 : 13522 : k = dual_modulus(p, lrho, aux, 1);
1460 [ + + ][ + + ]: 13522 : if (5*k < n || (n < 2*k && 5*k < 4*n))
[ + + ]
1461 : 3005 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1462 : : else
1463 : 10517 : { lrmin = lrho; i=k; radii[i] = lrho; }
1464 [ + + ]: 29432 : while (j > i+1)
1465 : : {
1466 [ + + ]: 15910 : if (i+j == n+1)
1467 : 5967 : lrho = (lrmin + lrmax) / 2;
1468 : : else
1469 : : {
1470 : 9943 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1471 [ + + ]: 9943 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1472 : 7267 : else lrho = lrmin * kappa + lrmax;
1473 : 9943 : lrho /= 1+kappa;
1474 : : }
1475 : 15910 : aux = (lrmax - lrmin) / (4*(j-i));
1476 : 15910 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1477 [ + + ][ + + ]: 15910 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
[ + + ]
1478 : 10946 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1479 : : else
1480 : 4964 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1481 : : }
1482 : 13522 : aux = lrmax - lrmin;
1483 : :
1484 [ + + ]: 13522 : if (ctr)
1485 : : {
1486 : 12214 : lrho = (lrmax + lrmin) / 2;
1487 [ + + ]: 79416 : for (i=1; i<=n; i++)
1488 [ + + ]: 67202 : if (radii[i] != UNDEF) radii[i] -= lrho;
1489 : :
1490 : 12214 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1491 : 12214 : R = mygprec(dblexp(-lrho), bit2);
1492 : 12214 : q = scalepol(p,R,bit2);
1493 : 12214 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1494 : : }
1495 : : else
1496 : : {
1497 : 1308 : lrho = logradius(radii, p, k, aux/10., &delta);
1498 : 1308 : update_radius(n, radii, lrho, ¶m, ¶m2);
1499 : :
1500 : 1308 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1501 : 1308 : R = mygprec(dblexp(-lrho), bit2);
1502 : 1308 : q = scalepol(p,R,bit2);
1503 : 1308 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1504 : : }
1505 : 13522 : bit += n;
1506 : 13522 : bit2 += n; R = invr(mygprec(R,bit2));
1507 : 13522 : *F = mygprec(scalepol(FF,R,bit2), bit);
1508 : 13522 : *G = mygprec(scalepol(GG,R,bit2), bit);
1509 : 13522 : }
1510 : :
1511 : : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1512 : : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1513 : : * where the maximum modulus of the roots of p is <=1.
1514 : : * Assume sum of roots is 0. */
1515 : : static void
1516 : 12214 : split_1(GEN p, long bit, GEN *F, GEN *G)
1517 : : {
1518 : 12214 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1519 : : GEN TWO, ctr, q, qq, FF, GG, v, gr, r, newq;
1520 : : double lrmin, lrmax, lthick;
1521 : 12214 : const double LOG3 = 1.098613;
1522 : :
1523 : 12214 : lrmax = logmax_modulus(p, 0.01);
1524 : 12214 : gr = mygprec(dblexp(-lrmax), bit2);
1525 : 12214 : q = scalepol(p,gr,bit2);
1526 : :
1527 : 12214 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1528 : 12214 : TWO = myreal_1(bit2); setexpo(TWO,1);
1529 : 12214 : v = cgetg(5,t_VEC);
1530 : 12214 : gel(v,1) = TWO;
1531 : 12214 : gel(v,2) = negr(TWO);
1532 : 12214 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1533 : 12214 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1534 : 12214 : q = mygprec(q,bit2); lthick = 0;
1535 : 12214 : newq = ctr = NULL; /* -Wall */
1536 [ + + ]: 12214 : imax = polreal? 3: 4;
1537 [ + + ]: 21931 : for (i=1; i<=imax; i++)
1538 : : {
1539 : 21651 : qq = RgX_translate(q, gel(v,i));
1540 : 21651 : lrmin = logmin_modulus(qq,0.05);
1541 [ + + ]: 21651 : if (LOG3 > lrmin + lthick)
1542 : : {
1543 : 21158 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1544 [ + + ]: 21158 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1545 : : }
1546 [ + + ]: 21651 : if (lthick > LOG2) break;
1547 [ + + ][ + + ]: 11421 : if (polreal && i==2 && lthick > LOG3 - LOG2) break;
[ + + ]
1548 : : }
1549 : 12214 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
1550 : 12214 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1551 : 12214 : r = gneg(mygprec(ctr,bit2));
1552 : 12214 : FF = RgX_translate(FF,r);
1553 : 12214 : GG = RgX_translate(GG,r);
1554 : :
1555 : 12214 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1556 : 12214 : *F = scalepol(FF,gr,bit2);
1557 : 12214 : *G = scalepol(GG,gr,bit2);
1558 : 12214 : }
1559 : :
1560 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1561 : : where the maximum modulus of the roots of p is < 0.5 */
1562 : : static int
1563 : 12214 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1564 : : {
1565 : : GEN q, b, FF, GG;
1566 : 12214 : long n = degpol(p), k, bit2, eq;
1567 : 12214 : double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
1568 : :
1569 : : /* beware double overflow */
1570 [ + + ][ + - ]: 12214 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
[ - + ]
1571 : :
1572 [ + + ]: 12214 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1573 : 12214 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1574 : :
1575 : 12214 : q = mygprec(p,bit2);
1576 : 12214 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1577 : 12214 : q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
1578 : 12214 : k = 0;
1579 [ + - ][ - + ]: 12214 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1580 [ - + ]: 12214 : || gequal0(gel(q,k+2)))) k++;
1581 [ - + ]: 12214 : if (k > 0)
1582 : : {
1583 [ # # ]: 0 : if (k > n/2) k = n/2;
1584 : 0 : bit2 += k<<1;
1585 : 0 : FF = monomial(myreal_1(bit2), k, 0);
1586 : 0 : GG = RgX_shift_shallow(q, -k);
1587 : : }
1588 : : else
1589 : : {
1590 : 12214 : split_1(q,bit2,&FF,&GG);
1591 : 12214 : bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
1592 : 12214 : FF = mygprec(FF,bit2);
1593 : : }
1594 : 12214 : GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
1595 : 12214 : *F = RgX_translate(FF, b);
1596 : 12214 : *G = RgX_translate(GG, b); return 1;
1597 : : }
1598 : :
1599 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1600 : : * Assume max_modulus(p) < 2 */
1601 : : static void
1602 : 12214 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1603 : : {
1604 : : GEN FF, GG;
1605 : : long n, bit2, normp;
1606 : :
1607 [ - + ]: 24428 : if (split_0_2(p,bit,F,G)) return;
1608 : :
1609 : 0 : normp = gexpo(p);
1610 : 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1611 : 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1612 : 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1613 : 0 : scalepol2n(FF,-2);
1614 : 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1615 : 0 : *F = mygprec(FF,bit2);
1616 : 12214 : *G = mygprec(GG,bit2);
1617 : : }
1618 : :
1619 : : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1620 : : static void
1621 : 13522 : split_0(GEN p, long bit, GEN *F, GEN *G)
1622 : : {
1623 : 13522 : const double LOG1_9 = 0.6418539;
1624 : 13522 : long n = degpol(p), k = 0;
1625 : : GEN q;
1626 : :
1627 [ - + ][ # # ]: 13522 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1628 [ - + ]: 13522 : if (k > 0)
1629 : : {
1630 [ # # ]: 0 : if (k > n/2) k = n/2;
1631 : 0 : *F = monomial(myreal_1(bit), k, 0);
1632 : 0 : *G = RgX_shift_shallow(p, -k);
1633 : : }
1634 : : else
1635 : : {
1636 : 13522 : double lr = logmax_modulus(p, 0.05);
1637 [ + + ]: 13522 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1638 : : else
1639 : : {
1640 : 10994 : q = RgX_recip_shallow(p);
1641 : 10994 : lr = logmax_modulus(q,0.05);
1642 [ + + ]: 10994 : if (lr < LOG1_9)
1643 : : {
1644 : 9686 : split_0_1(q, bit, F, G);
1645 : 9686 : *F = RgX_recip_shallow(*F);
1646 : 9686 : *G = RgX_recip_shallow(*G);
1647 : : }
1648 : : else
1649 : 1308 : split_2(p,bit,NULL, 1.2837,F,G);
1650 : : }
1651 : : }
1652 : 13522 : }
1653 : :
1654 : : /********************************************************************/
1655 : : /** **/
1656 : : /** ERROR ESTIMATE FOR THE ROOTS **/
1657 : : /** **/
1658 : : /********************************************************************/
1659 : :
1660 : : static GEN
1661 : 42096 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1662 : : {
1663 : 42096 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1664 : : long i, j;
1665 : :
1666 : 42096 : d = cgetg(n+1,t_VEC);
1667 [ + + ]: 661302 : for (i=1; i<=n; i++)
1668 : : {
1669 [ + + ]: 619206 : if (i!=k)
1670 : : {
1671 : 577110 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1672 : 577110 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1673 : : }
1674 : : }
1675 : 42096 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1676 [ + + ]: 42096 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1677 : 42096 : eps = mulrr(rho, shatzle);
1678 : 42096 : aux = shiftr(powru(rho,n), err);
1679 : :
1680 [ + + ][ + - ]: 134038 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
[ + + ]
1681 : : {
1682 : 91942 : GEN prod = NULL; /* 1. */
1683 : 91942 : long m = n;
1684 : 91942 : epsbis = mulrr(eps, dbltor(1.25));
1685 [ + + ]: 1672300 : for (i=1; i<=n; i++)
1686 : : {
1687 [ + + ][ + + ]: 1580358 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1688 : : {
1689 : 1473529 : GEN dif = subrr(gel(d,i),eps);
1690 [ + + ]: 1473529 : prod = prod? mulrr(prod, dif): dif;
1691 : 1473529 : m--;
1692 : : }
1693 : : }
1694 [ + + ]: 91942 : eps2 = prod? divrr(aux, prod): aux;
1695 [ + + ]: 91942 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1696 : 91942 : rap = divrr(eps,eps2); eps = eps2;
1697 : : }
1698 : 42096 : return eps;
1699 : : }
1700 : :
1701 : : /* round a complex or real number x to an absolute value of 2^(-bit) */
1702 : : static GEN
1703 : 91523 : mygprec_absolute(GEN x, long bit)
1704 : : {
1705 : : long e;
1706 : : GEN y;
1707 : :
1708 [ + + + ]: 91523 : switch(typ(x))
1709 : : {
1710 : : case t_REAL:
1711 : 64034 : e = expo(x) + bit;
1712 [ + + ][ - + ]: 64034 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1713 : : case t_COMPLEX:
1714 [ + + ]: 25057 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1715 : 24370 : y = cgetg(3,t_COMPLEX);
1716 : 24370 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1717 : 24370 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1718 : 24370 : return y;
1719 : 91523 : default: return x;
1720 : : }
1721 : : }
1722 : :
1723 : : static long
1724 : 6472 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1725 : : {
1726 : 6472 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1727 : : GEN sigma, shatzle;
1728 : :
1729 : 6472 : err += (long)log2((double)n) + 1;
1730 [ - + ]: 6472 : if (err > -2) return 0;
1731 : 6472 : sigma = real2n(-err, LOWDEFAULTPREC);
1732 : : /* 2 / ((s - 1)^(1/n) - 1) */
1733 : 6472 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1734 [ + + ]: 48568 : for (i=1; i<=n; i++)
1735 : : {
1736 : 42096 : pari_sp av = avma;
1737 : 42096 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1738 : 42096 : long e = gexpo(x);
1739 [ + + ]: 42096 : avma = av; if (e > e_max) e_max = e;
1740 : 42096 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1741 : : }
1742 : 6472 : return e_max;
1743 : : }
1744 : :
1745 : : /********************************************************************/
1746 : : /** **/
1747 : : /** MAIN **/
1748 : : /** **/
1749 : : /********************************************************************/
1750 : : static GEN
1751 : 31774 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1752 : :
1753 : : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1754 : : * returns prod (x-roots_pol[i]) */
1755 : : static GEN
1756 : 33516 : split_complete(GEN p, long bit, GEN roots_pol)
1757 : : {
1758 : 33516 : long n = degpol(p);
1759 : : pari_sp ltop;
1760 : : GEN p1, F, G, a, b, m1, m2;
1761 : :
1762 [ + + ]: 33516 : if (n == 1)
1763 : : {
1764 : 8214 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1765 : 8214 : (void)append_clone(roots_pol,a); return p;
1766 : : }
1767 : 25302 : ltop = avma;
1768 [ + + ]: 25302 : if (n == 2)
1769 : : {
1770 : 11780 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1771 : 11780 : F = gsqrt(F, nbits2prec(bit));
1772 : 11780 : p1 = ginv(gmul2n(gel(p,4),1));
1773 : 11780 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1774 : 11780 : b = gmul(gsub(F,gel(p,3)), p1);
1775 : 11780 : a = append_clone(roots_pol,a);
1776 : 11780 : b = append_clone(roots_pol,b); avma = ltop;
1777 : 11780 : a = mygprec(a, 3*bit);
1778 : 11780 : b = mygprec(b, 3*bit);
1779 : 11780 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1780 : : }
1781 : 13522 : split_0(p,bit,&F,&G);
1782 : 13522 : m1 = split_complete(F,bit,roots_pol);
1783 : 13522 : m2 = split_complete(G,bit,roots_pol);
1784 : 33516 : return gerepileupto(ltop, gmul(m1,m2));
1785 : : }
1786 : :
1787 : : static GEN
1788 : 168262 : quickabs(GEN x)
1789 : : {
1790 : 168262 : const long prec = DEFAULTPREC;
1791 : : GEN y;
1792 [ + + - + : 168262 : switch(typ(x))
- ]
1793 : : {
1794 : 167854 : case t_INT: y = itor(x, prec); setabssign(y); return y;
1795 : 396 : case t_REAL: y = rtor(x, prec); setabssign(y); return y;
1796 : 0 : case t_FRAC: y = fractor(x, prec); setabssign(y); return y;
1797 : : case t_COMPLEX: {
1798 : 12 : GEN a = gel(x,1), b = gel(x,2);
1799 : : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1800 [ - + ]: 12 : if (isintzero(a)) return cxcompotor(b, prec);
1801 [ - + ]: 12 : if (isintzero(b)) return cxcompotor(a, prec);
1802 : 12 : a = cxcompotor(a, prec);
1803 : 12 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1804 : : }
1805 : 0 : default: pari_err_TYPE("quickabs",x);
1806 : 168262 : return NULL;/*not reached*/
1807 : : }
1808 : :
1809 : : }
1810 : :
1811 : : /* bound ln |largest root of p| */
1812 : : double
1813 : 34522 : cauchy_bound(GEN p)
1814 : : {
1815 : 34522 : pari_sp av = avma;
1816 : 34522 : long i, n = degpol(p);
1817 : : GEN invlc;
1818 : 34522 : double Lmax = -pariINFINITY;
1819 : :
1820 [ - + ]: 34522 : if (n <= 0) pari_err_CONSTPOL("cauchy_bound");
1821 : 34522 : invlc = invr( quickabs(gel(p,n+2)) ); /* 1 / |lc(p)| */
1822 [ + + ]: 181991 : for (i = 0; i < n; i++)
1823 : : {
1824 : 147469 : GEN y = gel(p,i+2);
1825 : : double L;
1826 [ + + ]: 147469 : if (gequal0(y)) continue;
1827 : 133740 : L = mydbllogr(mulrr(quickabs(y), invlc)) / (n-i);
1828 [ + + ]: 133740 : if (L > Lmax) Lmax = L;
1829 : : }
1830 : 34522 : avma = av; return Lmax + LOG2;
1831 : : }
1832 : :
1833 : : static GEN
1834 : 48592 : mygprecrc_special(GEN x, long prec, long e)
1835 : : {
1836 : : GEN y;
1837 [ + + + ]: 48592 : switch(typ(x))
1838 : : {
1839 : : case t_REAL:
1840 [ - + ]: 420 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1841 [ + - ]: 420 : return (prec > realprec(x))? rtor(x, prec): x;
1842 : : case t_COMPLEX:
1843 : 12 : y = cgetg(3,t_COMPLEX);
1844 : 12 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1845 : 12 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1846 : 12 : return y;
1847 : 48592 : default: return x;
1848 : : }
1849 : : }
1850 : :
1851 : : /* like mygprec but keep at least the same precision as before */
1852 : : static GEN
1853 : 6472 : mygprec_special(GEN x, long bit)
1854 : : {
1855 : : long lx, i, e, prec;
1856 : : GEN y;
1857 : :
1858 [ - + ]: 6472 : if (bit < 0) bit = 0; /* should not happen */
1859 : 6472 : e = gexpo(x) - bit;
1860 : 6472 : prec = nbits2prec(bit);
1861 [ + - ]: 6472 : switch(typ(x))
1862 : : {
1863 : : case t_POL:
1864 : 6472 : y = cgetg_copy(x, &lx); y[1] = x[1];
1865 [ + + ]: 55040 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1866 : 6472 : break;
1867 : :
1868 : 0 : default: y = mygprecrc_special(x,prec,e);
1869 : : }
1870 : 6472 : return y;
1871 : : }
1872 : :
1873 : : static GEN
1874 : 3867 : fix_roots1(GEN r)
1875 : : {
1876 : 3867 : long i, l = lg(r);
1877 : 3867 : GEN allr = cgetg(l, t_VEC);
1878 [ + + ]: 28381 : for (i=1; i<l; i++)
1879 : : {
1880 : 24514 : GEN t = gel(r,i);
1881 : 24514 : gel(allr,i) = gcopy(t); gunclone(t);
1882 : : }
1883 : 3867 : return allr;
1884 : : }
1885 : : static GEN
1886 : 6472 : fix_roots(GEN r, GEN *m, long h, long bit)
1887 : : {
1888 : : long i, j, k, l, prec;
1889 : : GEN allr, ro1;
1890 [ + + ]: 6472 : if (h == 1) return fix_roots1(r);
1891 : 2605 : ro1 = initRUgen(h, bit);
1892 : 2605 : prec = nbits2prec(bit);
1893 : 2605 : l = lg(r)-1;
1894 : 2605 : allr = cgetg(h*l+1, t_VEC);
1895 [ + + ]: 9865 : for (k=1,i=1; i<=l; i++)
1896 : : {
1897 : 7260 : GEN p2, p1 = gel(r,i);
1898 [ + + ]: 7260 : p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1899 [ + + ]: 24842 : for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1900 : 7260 : gunclone(p1);
1901 : : }
1902 : 2605 : *m = roots_to_pol(allr, 0);
1903 : 6472 : return allr;
1904 : : }
1905 : :
1906 : : static GEN
1907 : 6378 : all_roots(GEN p, long bit)
1908 : : {
1909 : : GEN lc, pd, q, roots_pol, m;
1910 : 6378 : long bit0, bit2, n = degpol(p), i, e, h;
1911 : : pari_sp av;
1912 : :
1913 : 6378 : pd = RgX_deflate_max(p, &h); lc = leading_term(pd);
1914 [ + + ]: 6378 : e = (long)((2/LOG2) * cauchy_bound(pd)); if (e < 0) e = 0;
1915 : 6378 : bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
1916 : 6378 : bit2 = bit0; e = 0;
1917 : 6378 : for (av=avma,i=1;; i++,avma=av)
1918 : : {
1919 : 6472 : roots_pol = vectrunc_init(n+1);
1920 : 6472 : bit2 += e + (n << i);
1921 : 6472 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1922 : 6472 : q[1] = evalsigne(1)|evalvarn(0);
1923 : 6472 : m = split_complete(q,bit2,roots_pol);
1924 : 6472 : roots_pol = fix_roots(roots_pol, &m, h, bit2);
1925 : 6472 : q = mygprec_special(p,bit2); lc = leading_term(q);
1926 : 6472 : q[1] = evalsigne(1)|evalvarn(0);
1927 [ + + ]: 6472 : if (h > 1) m = gmul(m,lc);
1928 : :
1929 : 6472 : e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
1930 [ - + ]: 6472 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
1931 [ + - ]: 6472 : if (e < 0)
1932 : : {
1933 : 6472 : e = bit + a_posteriori_errors(p,roots_pol,e);
1934 [ + + ]: 6472 : if (e < 0) return roots_pol;
1935 : : }
1936 [ - + ]: 94 : if (DEBUGLEVEL > 7)
1937 : 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1938 : 94 : }
1939 : : }
1940 : :
1941 : : INLINE int
1942 : 1212 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1943 : :
1944 : : static int
1945 : 402 : isexactpol(GEN p)
1946 : : {
1947 : 402 : long i,n = degpol(p);
1948 [ + + ]: 1478 : for (i=0; i<=n; i++)
1949 [ + + ]: 1212 : if (!isexactscalar(gel(p,i+2))) return 0;
1950 : 402 : return 1;
1951 : : }
1952 : :
1953 : : static long
1954 : 1708 : isvalidcoeff(GEN x)
1955 : : {
1956 [ + + + ]: 1708 : switch (typ(x))
1957 : : {
1958 : 1692 : case t_INT: case t_REAL: case t_FRAC: return 1;
1959 [ + - ][ + - ]: 12 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
1960 : : }
1961 : 1708 : return 0;
1962 : : }
1963 : :
1964 : : static void
1965 : 406 : checkvalidpol(GEN p)
1966 : : {
1967 : 406 : long i,n = lg(p);
1968 [ + + ]: 2086 : for (i=2; i<n; i++)
1969 [ + + ]: 1684 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
1970 : 402 : }
1971 : :
1972 : : static GEN
1973 : 266 : solve_exact_pol(GEN p, long bit)
1974 : : {
1975 : 266 : long i, j, k, m, n = degpol(p), iroots = 0;
1976 : 266 : GEN ex, factors, v = zerovec(n);
1977 : :
1978 : 266 : factors = ZX_squff(Q_primpart(p), &ex);
1979 [ + + ]: 532 : for (i=1; i<lg(factors); i++)
1980 : : {
1981 : 266 : GEN roots_fact = all_roots(gel(factors,i), bit);
1982 : 266 : n = degpol(gel(factors,i));
1983 : 266 : m = ex[i];
1984 [ + + ]: 1064 : for (j=1; j<=n; j++)
1985 [ + + ]: 1596 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
1986 : : }
1987 : 266 : return v;
1988 : : }
1989 : :
1990 : : /* return the roots of p with absolute error bit */
1991 : : static GEN
1992 : 402 : roots_com(GEN q, long bit)
1993 : : {
1994 : : GEN L, p;
1995 : 402 : long v = RgX_valrem_inexact(q, &p);
1996 [ - + ]: 402 : if (lg(p) == 3) L = cgetg(1,t_VEC); /* constant polynomial */
1997 [ + + ]: 402 : else L = isexactpol(p)? solve_exact_pol(p,bit): all_roots(p,bit);
1998 [ + + ]: 402 : if (v)
1999 : : {
2000 : 56 : GEN M, z, t = gel(q,2);
2001 : : long i, x, y, l, n;
2002 : :
2003 [ + - ]: 56 : if (isrationalzero(t)) x = -bit;
2004 : : else
2005 : : {
2006 : 0 : n = gexpo(t);
2007 : 0 : x = n / v; l = degpol(q);
2008 [ # # ]: 0 : for (i = v; i <= l; i++)
2009 : : {
2010 : 0 : t = gel(q,i+2);
2011 [ # # ]: 0 : if (isrationalzero(t)) continue;
2012 : 0 : y = (n - gexpo(t)) / i;
2013 [ # # ]: 0 : if (y < x) x = y;
2014 : : }
2015 : : }
2016 : 56 : z = real_0_bit(x); l = v + lg(L);
2017 : 56 : M = cgetg(l, t_VEC); L -= v;
2018 [ + + ]: 112 : for (i = 1; i <= v; i++) gel(M,i) = z;
2019 [ + + ]: 168 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
2020 : 56 : L = M;
2021 : : }
2022 : 402 : return L;
2023 : : }
2024 : :
2025 : : static GEN
2026 : 24648 : tocomplex(GEN x, long l, long bit)
2027 : : {
2028 : : GEN y;
2029 [ + + ]: 24648 : if (typ(x) == t_COMPLEX)
2030 : : {
2031 [ + + ]: 24026 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
2032 : 2232 : x = gel(x,2);
2033 : 2232 : y = cgetg(3,t_COMPLEX);
2034 : 2232 : gel(y,1) = real_0_bit(-bit);
2035 : 2232 : gel(y,2) = mygprecrc(x, l, -bit);
2036 : : }
2037 : : else
2038 : : {
2039 : 622 : y = cgetg(3,t_COMPLEX);
2040 : 622 : gel(y,1) = mygprecrc(x, l, -bit);
2041 : 622 : gel(y,2) = real_0_bit(-bit);
2042 : : }
2043 : 24648 : return y;
2044 : : }
2045 : :
2046 : : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
2047 : : * up to 2^-e absolute error */
2048 : : static int
2049 : 84686 : cmp_complex_appr(void *E, GEN x, GEN y)
2050 : : {
2051 : 84686 : long e = (long)E;
2052 : : GEN z, xi, yi, xr, yr;
2053 : : long sxi, syi;
2054 [ + + ]: 84686 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
2055 : 38184 : else { xr = x; xi = NULL; sxi = 0; }
2056 [ + + ]: 84686 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
2057 : 35214 : else { yr = y; yi = NULL; syi = 0; }
2058 : : /* Compare absolute values of imaginary parts */
2059 [ + + ]: 84686 : if (!sxi)
2060 : : {
2061 [ + + ][ + - ]: 38874 : if (syi && expo(yi) >= e) return -1;
2062 : : /* |Im x| ~ |Im y| ~ 0 */
2063 : : }
2064 [ + + ]: 45812 : else if (!syi)
2065 : : {
2066 [ + - ][ + - ]: 1533 : if (sxi && expo(xi) >= e) return 1;
2067 : : /* |Im x| ~ |Im y| ~ 0 */
2068 : : }
2069 : : else
2070 : : {
2071 : : long sz;
2072 : 44279 : z = addrr_sign(xi, 1, yi, -1);
2073 : 44279 : sz = signe(z);
2074 [ + + ][ + - ]: 44279 : if (sz && expo(z) >= e) return (int)sz;
2075 : : }
2076 : : /* |Im x| ~ |Im y|, sort according to real parts */
2077 : 51688 : z = subrr(xr, yr);
2078 [ + + ]: 51688 : if (expo(z) >= e) return (int)signe(z);
2079 : : /* Re x ~ Re y. Place negative absolute value before positive */
2080 : 84686 : return (int) (sxi - syi);
2081 : : }
2082 : :
2083 : : static GEN
2084 : 6378 : clean_roots(GEN L, long l, long bit, long clean)
2085 : : {
2086 : 6378 : long i, n = lg(L), ex = 5 - bit;
2087 : 6378 : GEN res = cgetg(n,t_COL);
2088 [ + + ]: 47762 : for (i=1; i<n; i++)
2089 : : {
2090 : 41384 : GEN c = gel(L,i);
2091 [ + + ][ + + ]: 41384 : if (clean && isrealappr(c,ex))
2092 : : {
2093 [ - + ]: 16736 : if (typ(c) == t_COMPLEX) c = gel(c,1);
2094 : 16736 : c = mygprecrc(c, l, -bit);
2095 : : }
2096 : : else
2097 : 24648 : c = tocomplex(c, l, bit);
2098 : 41384 : gel(res,i) = c;
2099 : : }
2100 : 6378 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
2101 : 6378 : return res;
2102 : : }
2103 : :
2104 : : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
2105 : : static GEN
2106 : 406 : roots_aux(GEN p, long l, long clean)
2107 : : {
2108 : 406 : pari_sp av = avma;
2109 : : long bit;
2110 : : GEN L;
2111 : :
2112 [ - + ]: 406 : if (typ(p) != t_POL)
2113 : : {
2114 [ # # ]: 0 : if (gequal0(p)) pari_err_ROOTS0("roots");
2115 [ # # ]: 0 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
2116 : 0 : return cgetg(1,t_COL); /* constant polynomial */
2117 : : }
2118 [ - + ]: 406 : if (!signe(p)) pari_err_ROOTS0("roots");
2119 : 406 : checkvalidpol(p);
2120 [ - + ]: 402 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2121 [ - + ]: 402 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2122 : 402 : bit = prec2nbits(l);
2123 : 402 : L = roots_com(p, bit);
2124 : 402 : return gerepileupto(av, clean_roots(L, l, bit, clean));
2125 : : }
2126 : : GEN
2127 : 314 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
2128 : : /* clean up roots. If root is real replace it by its real part */
2129 : : GEN
2130 : 92 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
2131 : :
2132 : : /* private variant of conjvec. Allow non rational coefficients, shallow
2133 : : * function. */
2134 : : GEN
2135 : 44 : polmod_to_embed(GEN x, long prec)
2136 : : {
2137 : 44 : GEN v, T = gel(x,1), A = gel(x,2);
2138 : : long i, l;
2139 [ + - ][ - + ]: 44 : if (typ(A) != t_POL || varn(A) != varn(T))
2140 : : {
2141 : 0 : checkvalidpol(T);
2142 : 0 : return const_col(degpol(T), A);
2143 : : }
2144 : 44 : v = cleanroots(T,prec); l = lg(v);
2145 [ + + ]: 132 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
2146 : 44 : return v;
2147 : : }
2148 : :
2149 : : GEN
2150 : 5976 : QX_complex_roots(GEN p, long l)
2151 : : {
2152 : 5976 : pari_sp av = avma;
2153 : : long bit;
2154 : : GEN L;
2155 : :
2156 [ - + ]: 5976 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
2157 [ - + ]: 5976 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2158 [ - + ]: 5976 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2159 : 5976 : bit = prec2nbits(l);
2160 : 5976 : L = all_roots(Q_primpart(p), bit);
2161 : 5976 : return gerepileupto(av, clean_roots(L, l, bit, 1));
2162 : : }
|
