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 12762187 : addCC(GEN x, GEN y)
35 : {
36 : GEN z;
37 :
38 12762187 : if (typ(x) == t_INT)
39 : {
40 11383198 : if (typ(y) == t_INT) return addii(x,y);
41 : /* ty == t_COMPLEX */
42 1477 : z = cgetg(3,t_COMPLEX);
43 1477 : gel(z,1) = addii(x, gel(y,1));
44 1477 : gel(z,2) = icopy(gel(y,2)); return z;
45 : }
46 : /* tx == t_COMPLEX */
47 1378989 : z = cgetg(3,t_COMPLEX);
48 1378989 : if (typ(y) == t_INT)
49 : {
50 17568 : gel(z,1) = addii(gel(x,1),y);
51 17568 : gel(z,2) = icopy(gel(x,2)); return z;
52 : }
53 : /* ty == t_COMPLEX */
54 1361421 : gel(z,1) = addii(gel(x,1),gel(y,1));
55 1361421 : 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 23306487 : mulCC(GEN x, GEN y)
60 : {
61 : GEN z;
62 :
63 23306487 : if (typ(x) == t_INT)
64 : {
65 20189888 : if (typ(y) == t_INT) return mulii(x,y);
66 : /* ty == t_COMPLEX */
67 6928 : z = cgetg(3,t_COMPLEX);
68 6928 : gel(z,1) = mulii(x, gel(y,1));
69 6928 : gel(z,2) = mulii(x, gel(y,2)); return z;
70 : }
71 : /* tx == t_COMPLEX */
72 3116599 : z = cgetg(3,t_COMPLEX);
73 3116599 : if (typ(y) == t_INT)
74 : {
75 850165 : gel(z,1) = mulii(gel(x,1),y);
76 850165 : gel(z,2) = mulii(gel(x,2),y); return z;
77 : }
78 : /* ty == t_COMPLEX */
79 : {
80 2266434 : pari_sp av = avma, tetpil;
81 : GEN p1, p2;
82 :
83 2266434 : p1 = mulii(gel(x,1),gel(y,1));
84 2266434 : p2 = mulii(gel(x,2),gel(y,2));
85 4532868 : y = mulii(addii(gel(x,1),gel(x,2)),
86 4532868 : addii(gel(y,1),gel(y,2)));
87 2266434 : x = addii(p1,p2); tetpil = avma;
88 2266434 : gel(z,1) = subii(p1,p2);
89 2266434 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
90 2266434 : return z;
91 : }
92 : }
93 : /* fast squaring x: t_INT or t_COMPLEX(t_INT) */
94 : static GEN
95 17904149 : sqrCC(GEN x)
96 : {
97 : GEN z;
98 :
99 17904149 : if (typ(x) == t_INT) return sqri(x);
100 : /* tx == t_COMPLEX */
101 2845959 : z = cgetg(3,t_COMPLEX);
102 : {
103 2845959 : pari_sp av = avma, tetpil;
104 : GEN y, p1, p2;
105 :
106 2845959 : p1 = sqri(gel(x,1));
107 2845959 : p2 = sqri(gel(x,2));
108 2845959 : y = sqri(addii(gel(x,1),gel(x,2)));
109 2845959 : x = addii(p1,p2); tetpil = avma;
110 2845959 : gel(z,1) = subii(p1,p2);
111 2845959 : gel(z,2) = subii(y,x); gerepilecoeffssp(av,tetpil,z+1,2);
112 2845959 : return z;
113 : }
114 : }
115 :
116 : static void
117 2796891 : set_karasquare_limit(long bit)
118 : {
119 2796891 : if (bit<600) { KARASQUARE_LIMIT=8; COOKSQUARE_LIMIT=400; }
120 2471 : 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 2796891 : }
125 :
126 : /* assume lP > 0, lP = lgpol(P) */
127 : static GEN
128 6005627 : CX_square_spec(GEN P, long lP)
129 : {
130 : GEN s, t;
131 6005627 : long i, j, l, nn, n = lP - 1;
132 : pari_sp av;
133 :
134 6005627 : nn = n<<1; s = cgetg(nn+3,t_POL); s[1] = evalsigne(1)|evalvarn(0);
135 6005627 : gel(s,2) = sqrCC(gel(P,0)); /* i = 0 */
136 16549927 : for (i=1; i<=n; i++)
137 : {
138 10544300 : av = avma; l = (i+1)>>1;
139 10544300 : t = mulCC(gel(P,0), gel(P,i)); /* j = 0 */
140 10544300 : for (j=1; j<l; j++) t = addCC(t, mulCC(gel(P,j), gel(P,i-j)));
141 10544300 : t = gmul2n(t,1);
142 10544300 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
143 10544300 : gel(s,i+2) = gerepileupto(av, t);
144 : }
145 6005627 : gel(s,nn+2) = sqrCC(gel(P,n)); /* i = nn */
146 11898522 : for ( ; i<nn; i++)
147 : {
148 5892895 : av = avma; l = (i+1)>>1;
149 5892895 : t = mulCC(gel(P,i-n),gel(P,n)); /* j = i-n */
150 5892895 : for (j=i-n+1; j<l; j++) t = addCC(t, mulCC(gel(P,j),gel(P,i-j)));
151 5892895 : t = gmul2n(t,1);
152 5892895 : if ((i & 1) == 0) t = addCC(t, sqrCC(gel(P,i>>1)));
153 5892895 : gel(s,i+2) = gerepileupto(av, t);
154 : }
155 6005627 : return normalizepol_lg(s, nn+3);
156 : }
157 : /* not stack clean */
158 : static GEN
159 206584 : RgX_addspec(GEN x, long nx, GEN y, long ny)
160 : {
161 : GEN z, t;
162 : long i;
163 206584 : if (nx == ny) {
164 106817 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
165 106817 : for (i=0; i < nx; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
166 106817 : return normalizepol_lg(z, nx+2);
167 : }
168 99767 : if (ny < nx) {
169 99767 : z = cgetg(nx+2,t_POL); z[1] = evalsigne(1)|evalvarn(0); t = z+2;
170 99767 : for (i=0; i < ny; i++) gel(t,i) = gadd(gel(x,i),gel(y,i));
171 99767 : for ( ; i < nx; i++) gel(t,i) = gel(x,i);
172 99767 : 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 0 : 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 6213534 : karasquare(GEN P, long nP)
206 : {
207 : GEN Q, s0, s1, s2, a, t;
208 6213534 : long n0, n1, i, l, N, N0, N1, n = nP - 1; /* degree(P) */
209 : pari_sp av;
210 :
211 6213534 : if (n <= KARASQUARE_LIMIT) return nP? CX_square_spec(P, nP): pol_0(0);
212 206584 : av = avma;
213 206584 : n0 = (n>>1) + 1; n1 = nP - n0;
214 206584 : s0 = karasquare(P, n0); Q = P + n0;
215 206584 : s2 = karasquare(Q, n1);
216 206584 : s1 = RgX_addspec(P, n0, Q, n1);
217 206584 : s1 = RgX_sub(karasquare(s1+2, lgpol(s1)), RgX_add(s0,s2));
218 206584 : N = (n<<1) + 1;
219 206584 : a = cgetg(N + 2, t_POL); a[1] = evalsigne(1)|evalvarn(0);
220 206584 : t = a+2; l = lgpol(s0); s0 += 2; N0 = n0<<1;
221 206584 : for (i=0; i < l; i++) gel(t,i) = gel(s0,i);
222 206584 : for ( ; i < N0; i++) gel(t,i) = gen_0;
223 206584 : t = a+2 + N0; l = lgpol(s2); s2 += 2; N1 = N - N0;
224 206584 : for (i=0; i < l; i++) gel(t,i) = gel(s2,i);
225 206584 : for ( ; i < N1; i++) gel(t,i) = gen_0;
226 206584 : t = a+2 + n0; l = lgpol(s1); s1 += 2;
227 206584 : for (i=0; i < l; i++) gel(t,i) = gadd(gel(t,i), gel(s1,i));
228 206584 : return gerepilecopy(av, normalizepol_lg(a, N+2));
229 : }
230 : /* spec function. nP = lgpol(P) */
231 : static GEN
232 5593782 : cook_square(GEN P, long nP)
233 : {
234 : GEN Q, p0, p1, p2, p3, q, r, t, vp, vm;
235 5593782 : long n0, n3, i, j, n = nP - 1;
236 : pari_sp av;
237 :
238 5593782 : 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 0 : return gerepilecopy(av, normalizepol_lg(q, 2*n+3));
303 : }
304 :
305 : static GEN
306 2796891 : graeffe(GEN p)
307 : {
308 : GEN p0, p1, s0, s1;
309 2796891 : long n = degpol(p), n0, n1, i;
310 :
311 2796891 : if (!n) return gcopy(p);
312 2796891 : n0 = (n>>1)+1; n1 = n+1 - n0; /* n1 <= n0 <= n1+1 */
313 2796891 : p0 = new_chunk(n0);
314 2796891 : p1 = new_chunk(n1);
315 9673513 : for (i=0; i<n1; i++)
316 : {
317 6876622 : p0[i] = p[2+(i<<1)];
318 6876622 : p1[i] = p[3+(i<<1)];
319 : }
320 2796891 : if (n1 != n0)
321 1409583 : p0[i] = p[2+(i<<1)];
322 2796891 : s0 = cook_square(p0, n0);
323 2796891 : s1 = cook_square(p1, n1);
324 2796891 : return RgX_sub(s0, RgX_shift_shallow(s1,1));
325 : }
326 : GEN
327 4949 : ZX_graeffe(GEN p)
328 : {
329 4949 : pari_sp av = avma;
330 : GEN p0, p1, s0, s1;
331 4949 : long n = degpol(p);
332 :
333 4949 : if (!n) return ZX_copy(p);
334 4949 : RgX_even_odd(p, &p0, &p1);
335 : /* p = p0(x^2) + x p1(x^2) */
336 4949 : s0 = ZX_sqr(p0);
337 4949 : s1 = ZX_sqr(p1);
338 4949 : return gerepileupto(av, ZX_sub(s0, RgX_shift_shallow(s1,1)));
339 : }
340 : GEN
341 14 : polgraeffe(GEN p)
342 : {
343 14 : pari_sp av = avma;
344 : GEN p0, p1, s0, s1;
345 14 : long n = degpol(p);
346 :
347 14 : if (typ(p) != t_POL) pari_err_TYPE("polgraeffe",p);
348 14 : n = degpol(p);
349 14 : if (!n) return gcopy(p);
350 14 : RgX_even_odd(p, &p0, &p1);
351 : /* p = p0(x^2) + x p1(x^2) */
352 14 : s0 = RgX_sqr(p0);
353 14 : s1 = RgX_sqr(p1);
354 14 : return gerepileupto(av, RgX_sub(s0, RgX_shift_shallow(s1,1)));
355 : }
356 :
357 : /********************************************************************/
358 : /** **/
359 : /** MODULUS OF ROOTS **/
360 : /** **/
361 : /********************************************************************/
362 :
363 : /* Quick approximation to log2(|x|); first define y s.t. |y-x| < 2^-32 then
364 : * return y rounded to 2 ulp. In particular, if result < 2^21, absolute error
365 : * is bounded by 2^-31. If result > 2^21, it is correct to 2 ulp */
366 : static double
367 16134306 : mydbllog2i(GEN x)
368 : {
369 : #ifdef LONG_IS_64BIT
370 13862224 : const double W = 1/(4294967296. * 4294967296.); /* 2^-64 */
371 : #else
372 2272082 : const double W = 1/4294967296.; /*2^-32*/
373 : #endif
374 : GEN m;
375 16134306 : long lx = lgefint(x);
376 : double l;
377 16134306 : if (lx == 2) return -pariINFINITY;
378 16016496 : m = int_MSW(x);
379 16016496 : l = (double)(ulong)*m;
380 16016496 : if (lx == 3) return log2(l);
381 7457784 : 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 7457784 : return log2(l) + (double)(BITS_IN_LONG*(lx-3));
386 : }
387 :
388 : /* return log(|x|) or -pariINFINITY */
389 : static double
390 951796 : mydbllogr(GEN x) {
391 951796 : if (!signe(x)) return -pariINFINITY;
392 951796 : return LOG2*dbllog2r(x);
393 : }
394 :
395 : /* return log2(|x|) or -pariINFINITY */
396 : static double
397 9465841 : mydbllog2r(GEN x) {
398 9465841 : if (!signe(x)) return -pariINFINITY;
399 9428846 : return dbllog2r(x);
400 : }
401 : static double
402 4412434 : dbllog2mp(GEN x) { return typ(x) == t_INT? mydbllog2i(x): mydbllog2r(x); }
403 : double
404 22991518 : dbllog2(GEN z)
405 : {
406 : double x, y;
407 22991518 : switch(typ(z))
408 : {
409 12403166 : case t_INT: return mydbllog2i(z);
410 791 : case t_FRAC: return mydbllog2i(gel(z,1))-mydbllog2i(gel(z,2));
411 8381344 : case t_REAL: return mydbllog2r(z);
412 : default: /*t_COMPLEX*/
413 2206217 : x = dbllog2mp(gel(z,1));
414 2206217 : y = dbllog2mp(gel(z,2));
415 2206217 : if (fabs(x-y) > 10) return maxdd(x,y);
416 2089968 : return x + 0.5*log2(1 + exp2(2*(y-x)));
417 : }
418 : }
419 : static GEN /* beware overflow */
420 434788 : dblexp(double x) { return fabs(x) < 100.? dbltor(exp(x)): mpexp(dbltor(x)); }
421 :
422 : /* find s such that A_h <= 2^s <= 2 A_i for one h and all i < n = deg(p),
423 : * with A_i := (binom(n,i) lc(p) / p_i) ^ 1/(n-i), and p = sum p_i X^i */
424 : static long
425 2139318 : findpower(GEN p)
426 : {
427 2139318 : double x, L, mins = pariINFINITY;
428 2139318 : long n = degpol(p),i;
429 :
430 2139318 : L = dbllog2(gel(p,n+2)); /* log2(lc * binom(n,i)) */
431 11235679 : for (i=n-1; i>=0; i--)
432 : {
433 9096361 : L += log2((double)(i+1) / (double)(n-i));
434 9096361 : x = dbllog2(gel(p,i+2));
435 9096361 : if (x != -pariINFINITY)
436 : {
437 9020985 : double s = (L - x) / (double)(n-i);
438 9020985 : if (s < mins) mins = s;
439 : }
440 : }
441 2139318 : i = (long)ceil(mins);
442 2139318 : if (i - mins > 1 - 1e-12) i--;
443 2139318 : return i;
444 : }
445 :
446 : /* returns the exponent for logmodulus(), from the Newton diagram */
447 : static long
448 376817 : newton_polygon(GEN p, long k)
449 : {
450 376817 : pari_sp av = avma;
451 : double *logcoef, slope;
452 376817 : long n = degpol(p), i, j, h, l, *vertex;
453 :
454 376817 : init_dalloc();
455 376817 : logcoef = (double*)stack_malloc((n+1)*sizeof(double));
456 376817 : vertex = (long*)new_chunk(n+1);
457 :
458 : /* vertex[i] = 1 if i a vertex of convex hull, 0 otherwise */
459 376817 : for (i=0; i<=n; i++) { logcoef[i] = dbllog2(gel(p,2+i)); vertex[i] = 0; }
460 376817 : vertex[0] = 1; /* sentinel */
461 1598410 : for (i=0; i < n; i=h)
462 : {
463 1221593 : slope = logcoef[i+1]-logcoef[i];
464 5922110 : for (j = h = i+1; j<=n; j++)
465 : {
466 4700517 : double pij = (logcoef[j]-logcoef[i])/(double)(j-i);
467 4700517 : if (slope < pij) { slope = pij; h = j; }
468 : }
469 1221593 : vertex[h] = 1;
470 : }
471 376817 : h = k; while (!vertex[h]) h++;
472 376817 : l = k-1; while (!vertex[l]) l--;
473 376817 : avma = av;
474 376817 : return (long)floor((logcoef[h]-logcoef[l])/(double)(h-l) + 0.5);
475 : }
476 :
477 : /* change z into z*2^e, where z is real or complex of real */
478 : static void
479 2605595 : myshiftrc(GEN z, long e)
480 : {
481 2605595 : if (typ(z)==t_COMPLEX)
482 : {
483 516485 : if (signe(gel(z,1))) shiftr_inplace(gel(z,1), e);
484 516485 : if (signe(gel(z,2))) shiftr_inplace(gel(z,2), e);
485 : }
486 : else
487 2089110 : if (signe(z)) shiftr_inplace(z, e);
488 2605595 : }
489 :
490 : /* return z*2^e, where z is integer or complex of integer (destroy z) */
491 : static GEN
492 9704016 : myshiftic(GEN z, long e)
493 : {
494 9704016 : if (typ(z)==t_COMPLEX)
495 : {
496 1450142 : gel(z,1) = signe(gel(z,1))? mpshift(gel(z,1),e): gen_0;
497 1450142 : gel(z,2) = mpshift(gel(z,2),e);
498 1450142 : return z;
499 : }
500 8253874 : return signe(z)? mpshift(z,e): gen_0;
501 : }
502 :
503 : static GEN
504 476136 : RgX_gtofp_bit(GEN q, long bit)
505 : {
506 476136 : if (bit < 0) bit = 0;
507 476136 : return RgX_gtofp(q, nbits2prec(bit));
508 : }
509 :
510 : static GEN
511 17904534 : mygprecrc(GEN x, long prec, long e)
512 : {
513 : GEN y;
514 17904534 : switch(typ(x))
515 : {
516 13850864 : case t_REAL: return signe(x)? rtor(x, prec): real_0_bit(e);
517 : case t_COMPLEX:
518 3078201 : y = cgetg(3,t_COMPLEX);
519 3078201 : gel(y,1) = mygprecrc(gel(x,1),prec,e);
520 3078201 : gel(y,2) = mygprecrc(gel(x,2),prec,e);
521 3078201 : return y;
522 975469 : default: return gcopy(x);
523 : }
524 : }
525 :
526 : /* gprec behaves badly with the zero for polynomials.
527 : The second parameter in mygprec is the precision in base 2 */
528 : static GEN
529 4585439 : mygprec(GEN x, long bit)
530 : {
531 : long lx, i, e, prec;
532 : GEN y;
533 :
534 4585439 : if (bit < 0) bit = 0; /* should rarely happen */
535 4585439 : e = gexpo(x) - bit;
536 4585439 : prec = nbits2prec(bit);
537 4585439 : switch(typ(x))
538 : {
539 : case t_POL:
540 2720133 : y = cgetg_copy(x, &lx); y[1] = x[1];
541 2720133 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc(gel(x,i),prec,e);
542 2720133 : break;
543 :
544 1865306 : default: y = mygprecrc(x,prec,e);
545 : }
546 4585439 : return y;
547 : }
548 :
549 : /* normalize a polynomial p, that is change it with coefficients in Z[i],
550 : after making product by 2^shift */
551 : static GEN
552 1124968 : pol_to_gaussint(GEN p, long shift)
553 : {
554 1124968 : long i, l = lg(p);
555 1124968 : GEN q = cgetg(l, t_POL); q[1] = p[1];
556 1124968 : for (i=2; i<l; i++) gel(q,i) = gtrunc2n(gel(p,i), shift);
557 1124968 : return q;
558 : }
559 :
560 : /* returns a polynomial q in Z[i][x] keeping bit bits of p */
561 : static GEN
562 896492 : eval_rel_pol(GEN p, long bit)
563 : {
564 : long i;
565 6627809 : for (i = 2; i < lg(p); i++)
566 5731317 : if (gequal0(gel(p,i))) gel(p,i) = gen_0; /* bad behaviour of gexpo */
567 896492 : return pol_to_gaussint(p, bit-gexpo(p)+1);
568 : }
569 :
570 : /* returns p(R*x)/R^n (in R or R[i]), R = exp(lrho), bit bits of precision */
571 : static GEN
572 99383 : homothetie(GEN p, double lrho, long bit)
573 : {
574 : GEN q, r, t, iR;
575 99383 : long n = degpol(p), i;
576 :
577 99383 : iR = mygprec(dblexp(-lrho),bit);
578 99383 : q = mygprec(p, bit);
579 99383 : r = cgetg(n+3,t_POL); r[1] = p[1];
580 99383 : t = iR; r[n+2] = q[n+2];
581 665074 : for (i=n-1; i>0; i--)
582 : {
583 565691 : gel(r,i+2) = gmul(t, gel(q,i+2));
584 565691 : t = mulrr(t, iR);
585 : }
586 99383 : gel(r,2) = gmul(t, gel(q,2)); return r;
587 : }
588 :
589 : /* change q in 2^(n*e) p(x*2^(-e)), n=deg(q) [ ~as above with R = 2^-e ]*/
590 : static void
591 605293 : homothetie2n(GEN p, long e)
592 : {
593 605293 : if (e)
594 : {
595 455962 : long i,n = lg(p)-1;
596 455962 : for (i=2; i<=n; i++) myshiftrc(gel(p,i), (n-i)*e);
597 : }
598 605293 : }
599 :
600 : /* return 2^f * 2^(n*e) p(x*2^(-e)), n=deg(q) */
601 : static void
602 1910842 : homothetie_gauss(GEN p, long e, long f)
603 : {
604 1910842 : if (e || f)
605 : {
606 1748636 : long i, n = lg(p)-1;
607 1748636 : for (i=2; i<=n; i++) gel(p,i) = myshiftic(gel(p,i), f+(n-i)*e);
608 : }
609 1910842 : }
610 :
611 : /* Lower bound on the modulus of the largest root z_0
612 : * k is set to an upper bound for #{z roots, |z-z_0| < eps} */
613 : static double
614 2139318 : lower_bound(GEN p, long *k, double eps)
615 : {
616 2139318 : long n = degpol(p), i, j;
617 2139318 : pari_sp ltop = avma;
618 : GEN a, s, S, ilc;
619 : double r, R, rho;
620 :
621 2139318 : if (n < 4) { *k = n; return 0.; }
622 987181 : S = cgetg(5,t_VEC);
623 987181 : a = cgetg(5,t_VEC); ilc = gdiv(real_1(DEFAULTPREC), gel(p,n+2));
624 987181 : for (i=1; i<=4; i++) gel(a,i) = gmul(ilc,gel(p,n+2-i));
625 : /* i = 1 split out from next loop for efficiency and initialization */
626 987181 : s = gel(a,1);
627 987181 : gel(S,1) = gneg(s); /* Newton sum S_i */
628 987181 : rho = r = gtodouble(gabs(s,3));
629 987181 : R = r / n;
630 3948724 : for (i=2; i<=4; i++)
631 : {
632 2961543 : s = gmulsg(i,gel(a,i));
633 2961543 : for (j=1; j<i; j++) s = gadd(s, gmul(gel(S,j),gel(a,i-j)));
634 2961543 : gel(S,i) = gneg(s); /* Newton sum S_i */
635 2961543 : r = gtodouble(gabs(s,3));
636 2961543 : if (r > 0.)
637 : {
638 2951495 : r = exp(log(r/n) / (double)i);
639 2951495 : if (r > R) R = r;
640 : }
641 : }
642 987181 : if (R > 0. && eps < 1.2)
643 985331 : *k = (long)floor((rho/R + n) / (1 + exp(-eps)*cos(eps)));
644 : else
645 1850 : *k = n;
646 987181 : avma = ltop; return R;
647 : }
648 :
649 : /* log of modulus of the largest root of p with relative error tau. Assume
650 : * P(0) != 0 and P non constant */
651 : static double
652 228476 : logmax_modulus(GEN p, double tau)
653 : {
654 : GEN r, q, aux, gunr;
655 228476 : pari_sp av, ltop = avma;
656 228476 : long i,k,n=degpol(p),nn,bit,M,e;
657 228476 : double rho,eps, tau2 = (tau > 3.0)? 0.5: tau/6.;
658 :
659 228476 : r = cgeti(BIGDEFAULTPREC);
660 228476 : av = avma;
661 :
662 228476 : eps = - 1/log(1.5*tau2); /* > 0 */
663 228476 : bit = (long) ((double) n*log2(1./tau2)+3*log2((double) n))+1;
664 228476 : gunr = real_1_bit(bit+2*n);
665 228476 : aux = gdiv(gunr, gel(p,2+n));
666 228476 : q = RgX_Rg_mul(p, aux); gel(q,2+n) = gunr;
667 228476 : e = findpower(q);
668 228476 : homothetie2n(q,e);
669 228476 : affsi(e, r);
670 228476 : q = pol_to_gaussint(q, bit);
671 228476 : M = (long) (log2( log(4.*n) / (2*tau2) )) + 2;
672 228476 : nn = n;
673 228476 : for (i=0,e=0;;)
674 : { /* nn = deg(q) */
675 2139318 : rho = lower_bound(q, &k, eps);
676 2139318 : if (rho > exp2(-(double)e)) e = (long)-floor(log2(rho));
677 2139318 : affii(shifti(addis(r,e), 1), r);
678 2139318 : if (++i == M) break;
679 :
680 5732526 : bit = (long) ((double)k * log2(1./tau2) +
681 3821684 : (double)(nn-k)*log2(1./eps) + 3*log2((double)nn)) + 1;
682 1910842 : homothetie_gauss(q, e, bit-(long)floor(dbllog2(gel(q,2+nn))+0.5));
683 1910842 : nn -= RgX_valrem(q, &q);
684 1910842 : set_karasquare_limit(gexpo(q));
685 1910842 : q = gerepileupto(av, graeffe(q));
686 1910842 : tau2 *= 1.5; if (tau2 > 0.9) tau2 = 0.5;
687 1910842 : eps = -1/log(tau2); /* > 0 */
688 1910842 : e = findpower(q);
689 1910842 : }
690 228476 : if (!signe(r)) { avma = ltop; return 0.; }
691 211992 : r = itor(r, DEFAULTPREC); shiftr_inplace(r, -M);
692 211992 : avma = ltop; return -rtodbl(r) * LOG2; /* -log(2) sum e_i 2^-i */
693 : }
694 : /* assume P non constant */
695 : GEN
696 1842 : logmax_modulus_bound(GEN P)
697 : {
698 1842 : (void)RgX_valrem_inexact(P,&P);
699 1842 : return lg(P)==3? gen_0: dblexp(logmax_modulus(P, 0.01) + 0.01);
700 : }
701 :
702 : /* log of modulus of the smallest root of p, with relative error tau */
703 : static double
704 79955 : logmin_modulus(GEN p, double tau)
705 : {
706 79955 : pari_sp av = avma;
707 : double r;
708 :
709 79955 : if (gequal0(gel(p,2))) return -pariINFINITY;
710 79955 : r = - logmax_modulus(RgX_recip_shallow(p),tau);
711 79955 : avma = av; return r;
712 : }
713 :
714 : /* return the log of the k-th modulus (ascending order) of p, rel. error tau*/
715 : static double
716 41149 : logmodulus(GEN p, long k, double tau)
717 : {
718 : GEN q;
719 41149 : long i, kk = k, imax, n = degpol(p), nn, bit, e;
720 41149 : pari_sp av, ltop=avma;
721 41149 : double r, tau2 = tau/6;
722 :
723 41149 : bit = (long)(n * (2. + log2(3.*n/tau2)));
724 41149 : av = avma;
725 41149 : q = gprec_w(p, nbits2prec(bit));
726 41149 : q = RgX_gtofp_bit(q, bit);
727 41149 : e = newton_polygon(q,k);
728 41149 : r = (double)e;
729 41149 : homothetie2n(q,e);
730 41149 : imax = (long)(log2(3./tau) + log2(log(4.*n)))+1;
731 376817 : for (i=1; i<imax; i++)
732 : {
733 335668 : q = eval_rel_pol(q,bit);
734 335668 : kk -= RgX_valrem(q, &q);
735 335668 : nn = degpol(q);
736 :
737 335668 : set_karasquare_limit(bit);
738 335668 : q = gerepileupto(av, graeffe(q));
739 335668 : e = newton_polygon(q,kk);
740 335668 : r += e / exp2((double)i);
741 335668 : q = RgX_gtofp_bit(q, bit);
742 335668 : homothetie2n(q,e);
743 :
744 335668 : tau2 *= 1.5; if (tau2 > 1.) tau2 = 1.;
745 335668 : bit = 1 + (long)(nn*(2. + log2(3.*nn/tau2)));
746 : }
747 41149 : avma = ltop; return -r * LOG2;
748 : }
749 :
750 : /* return the log of the k-th modulus r_k of p, rel. error tau, knowing that
751 : * rmin < r_k < rmax. This information helps because we may reduce precision
752 : * quicker */
753 : static double
754 41149 : logpre_modulus(GEN p, long k, double tau, double lrmin, double lrmax)
755 : {
756 : GEN q;
757 41149 : long n = degpol(p), i, imax, imax2, bit;
758 41149 : pari_sp ltop = avma, av;
759 41149 : double lrho, aux, tau2 = tau/6.;
760 :
761 41149 : aux = (lrmax - lrmin) / 2. + 4*tau2;
762 41149 : imax = (long) log2(log((double)n)/ aux);
763 41149 : if (imax <= 0) return logmodulus(p,k,tau);
764 :
765 40321 : lrho = (lrmin + lrmax) / 2;
766 40321 : av = avma;
767 40321 : bit = (long)(n*(2. + aux / LOG2 - log2(tau2)));
768 40321 : q = homothetie(p, lrho, bit);
769 40321 : imax2 = (long)(log2(3./tau * log(4.*n))) + 1;
770 40321 : if (imax > imax2) imax = imax2;
771 :
772 126544 : for (i=0; i<imax; i++)
773 : {
774 86223 : q = eval_rel_pol(q,bit);
775 86223 : set_karasquare_limit(bit);
776 86223 : q = gerepileupto(av, graeffe(q));
777 86223 : aux = 2*aux + 2*tau2;
778 86223 : tau2 *= 1.5;
779 86223 : bit = (long)(n*(2. + aux / LOG2 - log2(1-exp(-tau2))));
780 86223 : q = RgX_gtofp_bit(q, bit);
781 : }
782 40321 : aux = exp2((double)imax);
783 40321 : aux = logmodulus(q,k, aux*tau/3.) / aux;
784 40321 : avma = ltop; return lrho + aux;
785 : }
786 :
787 : static double
788 48619 : ind_maxlog2(GEN q)
789 : {
790 48619 : long i, k = -1;
791 48619 : double L = - pariINFINITY;
792 134696 : for (i=0; i<=degpol(q); i++)
793 : {
794 86077 : double d = dbllog2(gel(q,2+i));
795 86077 : if (d > L) { L = d; k = i; }
796 : }
797 48619 : return k;
798 : }
799 :
800 : /* Returns k such that r_k e^(-tau) < R < r_{k+1} e^tau.
801 : * Assume that l <= k <= n-l */
802 : static long
803 59062 : dual_modulus(GEN p, double lrho, double tau, long l)
804 : {
805 59062 : long i, imax, delta_k = 0, n = degpol(p), nn, v2, v, bit, ll = l;
806 59062 : double tau2 = tau * 7./8.;
807 59062 : pari_sp av = avma;
808 : GEN q;
809 :
810 59062 : bit = 6*n - 5*l + (long)(n*(-log2(tau2) + tau2 * 8./7.));
811 59062 : q = homothetie(p, lrho, bit);
812 59062 : imax = (long)(log(log(2.*n)/tau2)/log(7./4.)+1);
813 :
814 523220 : for (i=0; i<imax; i++)
815 : {
816 474601 : q = eval_rel_pol(q,bit); v2 = n - degpol(q);
817 474601 : v = RgX_valrem(q, &q);
818 474601 : ll -= maxss(v, v2); if (ll < 0) ll = 0;
819 :
820 474601 : nn = degpol(q); delta_k += v;
821 474601 : if (!nn) return delta_k;
822 :
823 464158 : set_karasquare_limit(bit);
824 464158 : q = gerepileupto(av, graeffe(q));
825 464158 : tau2 *= 7./4.;
826 464158 : bit = 6*nn - 5*ll + (long)(nn*(-log2(tau2) + tau2 * 8./7.));
827 : }
828 48619 : avma = av; return delta_k + (long)ind_maxlog2(q);
829 : }
830 :
831 : /********************************************************************/
832 : /** **/
833 : /** FACTORS THROUGH CIRCLE INTEGRATION **/
834 : /** **/
835 : /********************************************************************/
836 : /* l power of 2 */
837 : static void
838 2118522 : fft(GEN Omega, GEN p, GEN f, long step, long l)
839 : {
840 : pari_sp ltop;
841 : long i, l1, l2, l3, rapi, step4;
842 : GEN f1, f2, f3, f02, f13, g02, g13, ff;
843 :
844 2118522 : if (l == 2)
845 : {
846 1237268 : gel(f,0) = gadd(gel(p,0),gel(p,step));
847 1237268 : gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
848 : }
849 881254 : if (l == 4)
850 : {
851 435946 : f1 = gadd(gel(p,0), gel(p,step<<1));
852 435946 : f2 = gsub(gel(p,0), gel(p,step<<1));
853 435946 : f3 = gadd(gel(p,step),gel(p,3*step));
854 435946 : f02= gsub(gel(p,step),gel(p,3*step));
855 435946 : f02 = mulcxI(f02);
856 435946 : gel(f,0) = gadd(f1, f3);
857 435946 : gel(f,1) = gadd(f2, f02);
858 435946 : gel(f,2) = gsub(f1, f3);
859 435946 : gel(f,3) = gsub(f2, f02); return;
860 : }
861 :
862 445308 : ltop = avma;
863 445308 : l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
864 445308 : fft(Omega,p, f, step4,l1);
865 445308 : fft(Omega,p+step, f+l1,step4,l1);
866 445308 : fft(Omega,p+(step<<1),f+l2,step4,l1);
867 445308 : fft(Omega,p+3*step, f+l3,step4,l1);
868 :
869 445308 : ff = cgetg(l+1,t_VEC);
870 1920238 : for (i=0; i<l1; i++)
871 : {
872 1474930 : rapi = step*i;
873 1474930 : f1 = gmul(gel(Omega,rapi), gel(f,i+l1));
874 1474930 : f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
875 1474930 : f3 = gmul(gel(Omega,3*rapi), gel(f,i+l3));
876 :
877 1474930 : f02 = gadd(gel(f,i),f2);
878 1474930 : g02 = gsub(gel(f,i),f2);
879 1474930 : f13 = gadd(f1,f3);
880 1474930 : g13 = mulcxI(gsub(f1,f3));
881 :
882 1474930 : gel(ff,i+1) = gadd(f02, f13);
883 1474930 : gel(ff,i+l1+1) = gadd(g02, g13);
884 1474930 : gel(ff,i+l2+1) = gsub(f02, f13);
885 1474930 : gel(ff,i+l3+1) = gsub(g02, g13);
886 : }
887 445308 : ff = gerepilecopy(ltop,ff);
888 445308 : for (i=0; i<l; i++) f[i] = ff[i+1];
889 : }
890 :
891 : GEN
892 0 : FFTinit(long k, long prec)
893 : {
894 0 : if (k <= 0) pari_err_DOMAIN("FFTinit", "k", "<=", gen_0, stoi(k));
895 0 : return grootsof1(1L << k, prec);
896 : }
897 :
898 : GEN
899 0 : FFT(GEN x, GEN Omega)
900 : {
901 0 : long i, l = lg(Omega), n = lg(x);
902 : GEN y, z;
903 0 : if (!is_vec_t(typ(x))) pari_err_TYPE("FFT",x);
904 0 : if (typ(Omega) != t_VEC) pari_err_TYPE("FFT",Omega);
905 0 : if (n > l) pari_err_DIM("FFT");
906 :
907 0 : if (n < l) {
908 0 : z = cgetg(l, t_VECSMALL); /* cf stackdummy */
909 0 : for (i = 1; i < n; i++) z[i] = x[i];
910 0 : for ( ; i < l; i++) gel(z,i) = gen_0;
911 : }
912 0 : else z = x;
913 0 : y = cgetg(l, t_VEC);
914 0 : fft(Omega+1, z+1, y+1, 1, l-1);
915 0 : return y;
916 : }
917 :
918 : /* returns 1 if p has only real coefficients, 0 else */
919 : static int
920 51401 : isreal(GEN p)
921 : {
922 : long i;
923 342362 : for (i = lg(p)-1; i > 1; i--)
924 303034 : if (typ(gel(p,i)) == t_COMPLEX) return 0;
925 39328 : return 1;
926 : }
927 :
928 : /* x non complex */
929 : static GEN
930 38642 : abs_update_r(GEN x, double *mu) {
931 38642 : GEN y = gtofp(x, DEFAULTPREC);
932 38642 : double ly = mydbllogr(y); if (ly < *mu) *mu = ly;
933 38642 : setabssign(y); return y;
934 : }
935 : /* return |x|, low accuracy. Set *mu = min(log(y), *mu) */
936 : static GEN
937 858656 : abs_update(GEN x, double *mu) {
938 : GEN y, xr, yr;
939 : double ly;
940 858656 : if (typ(x) != t_COMPLEX) return abs_update_r(x, mu);
941 821022 : xr = gel(x,1);
942 821022 : yr = gel(x,2);
943 821022 : if (gequal0(xr)) return abs_update_r(yr,mu);
944 820896 : if (gequal0(yr)) return abs_update_r(xr,mu);
945 : /* have to treat 0 specially: 0E-10 + 1e-20 = 0E-10 */
946 820014 : xr = gtofp(xr, DEFAULTPREC);
947 820014 : yr = gtofp(yr, DEFAULTPREC);
948 820014 : y = sqrtr(addrr(sqrr(xr), sqrr(yr)));
949 820014 : ly = mydbllogr(y); if (ly < *mu) *mu = ly;
950 820014 : return y;
951 : }
952 :
953 : static void
954 55696 : initdft(GEN *Omega, GEN *prim, long N, long Lmax, long bit)
955 : {
956 55696 : long prec = nbits2prec(bit);
957 55696 : *Omega = grootsof1(Lmax, prec) + 1;
958 55696 : *prim = rootsof1u_cx(N, prec);
959 55696 : }
960 :
961 : static void
962 26737 : parameters(GEN p, long *LMAX, double *mu, double *gamma,
963 : int polreal, double param, double param2)
964 : {
965 : GEN q, pc, Omega, A, RU, prim, g, TWO;
966 26737 : long n = degpol(p), bit, NN, K, i, j, Lmax;
967 26737 : pari_sp av2, av = avma;
968 :
969 26737 : bit = gexpo(p) + (long)param2+8;
970 26737 : Lmax = 4; while (Lmax <= n) Lmax <<= 1;
971 26737 : NN = (long)(param*3.14)+1; if (NN < Lmax) NN = Lmax;
972 26737 : K = NN/Lmax; if (K & 1) K++;
973 26737 : NN = Lmax*K;
974 26737 : if (polreal) K = K/2+1;
975 :
976 26737 : initdft(&Omega, &prim, NN, Lmax, bit);
977 26737 : q = mygprec(p,bit) + 2;
978 26737 : A = cgetg(Lmax+1,t_VEC); A++;
979 26737 : pc= cgetg(Lmax+1,t_VEC); pc++;
980 26737 : for (i=0; i <= n; i++) gel(pc,i)= gel(q,i);
981 26737 : for ( ; i<Lmax; i++) gel(pc,i) = gen_0;
982 :
983 26737 : *mu = pariINFINITY;
984 26737 : g = real_0_bit(-bit);
985 26737 : TWO = real2n(1, DEFAULTPREC);
986 26737 : av2 = avma;
987 26737 : RU = gen_1;
988 104683 : for (i=0; i<K; i++)
989 : {
990 77946 : if (i) {
991 51209 : GEN z = RU;
992 395567 : for (j=1; j<n; j++)
993 : {
994 344358 : gel(pc,j) = gmul(gel(q,j),z);
995 344358 : z = gmul(z,RU); /* RU = prim^i, z=prim^(ij) */
996 : }
997 51209 : gel(pc,n) = gmul(gel(q,n),z);
998 : }
999 :
1000 77946 : fft(Omega,pc,A,1,Lmax);
1001 91108 : if (polreal && i>0 && i<K-1)
1002 13162 : for (j=0; j<Lmax; j++) g = addrr(g, divrr(TWO, abs_update(gel(A,j),mu)));
1003 : else
1004 64784 : for (j=0; j<Lmax; j++) g = addrr(g, invr(abs_update(gel(A,j),mu)));
1005 77946 : RU = gmul(RU, prim);
1006 77946 : if (gc_needed(av,1))
1007 : {
1008 0 : if(DEBUGMEM>1) pari_warn(warnmem,"parameters");
1009 0 : gerepileall(av2,2, &g,&RU);
1010 : }
1011 : }
1012 26737 : *gamma = mydbllog2r(divru(g,NN));
1013 26737 : *LMAX = Lmax; avma = av;
1014 26737 : }
1015 :
1016 : /* NN is a multiple of Lmax */
1017 : static void
1018 28959 : dft(GEN p, long k, long NN, long Lmax, long bit, GEN F, GEN H, long polreal)
1019 : {
1020 : GEN Omega, q, qd, pc, pd, A, B, C, RU, aux, U, W, prim, prim2;
1021 28959 : long n = degpol(p), i, j, K;
1022 : pari_sp ltop;
1023 :
1024 28959 : initdft(&Omega, &prim, NN, Lmax, bit);
1025 28959 : RU = cgetg(n+2,t_VEC) + 1;
1026 :
1027 28959 : K = NN/Lmax; if (polreal) K = K/2+1;
1028 28959 : q = mygprec(p,bit);
1029 28959 : qd = RgX_deriv(q);
1030 :
1031 28959 : A = cgetg(Lmax+1,t_VEC); A++;
1032 28959 : B = cgetg(Lmax+1,t_VEC); B++;
1033 28959 : C = cgetg(Lmax+1,t_VEC); C++;
1034 28959 : pc = cgetg(Lmax+1,t_VEC); pc++;
1035 28959 : pd = cgetg(Lmax+1,t_VEC); pd++;
1036 28959 : pc[0] = q[2]; for (i=n+1; i<Lmax; i++) gel(pc,i) = gen_0;
1037 28959 : pd[0] = qd[2]; for (i=n; i<Lmax; i++) gel(pd,i) = gen_0;
1038 :
1039 28959 : ltop = avma;
1040 28959 : W = cgetg(k+1,t_VEC);
1041 28959 : U = cgetg(k+1,t_VEC);
1042 28959 : for (i=1; i<=k; i++) gel(W,i) = gel(U,i) = gen_0;
1043 :
1044 28959 : gel(RU,0) = gen_1;
1045 28959 : prim2 = gen_1;
1046 93795 : for (i=0; i<K; i++)
1047 : {
1048 64836 : gel(RU,1) = prim2;
1049 64836 : for (j=1; j<n; j++) gel(RU,j+1) = gmul(gel(RU,j),prim2);
1050 : /* RU[j] = prim^(ij)= prim2^j */
1051 :
1052 64836 : for (j=1; j<n; j++) gel(pd,j) = gmul(gel(qd,j+2),gel(RU,j));
1053 64836 : fft(Omega,pd,A,1,Lmax);
1054 64836 : for (j=1; j<=n; j++) gel(pc,j) = gmul(gel(q,j+2),gel(RU,j));
1055 64836 : fft(Omega,pc,B,1,Lmax);
1056 64836 : for (j=0; j<Lmax; j++) gel(C,j) = ginv(gel(B,j));
1057 64836 : for (j=0; j<Lmax; j++) gel(B,j) = gmul(gel(A,j),gel(C,j));
1058 64836 : fft(Omega,B,A,1,Lmax);
1059 64836 : fft(Omega,C,B,1,Lmax);
1060 :
1061 64836 : if (polreal) /* p has real coefficients */
1062 : {
1063 54008 : if (i>0 && i<K-1)
1064 : {
1065 41278 : for (j=1; j<=k; j++)
1066 : {
1067 35215 : gel(W,j) = gadd(gel(W,j), gshift(mulreal(gel(A,j+1),gel(RU,j+1)),1));
1068 35215 : gel(U,j) = gadd(gel(U,j), gshift(mulreal(gel(B,j),gel(RU,j)),1));
1069 : }
1070 : }
1071 : else
1072 : {
1073 136410 : for (j=1; j<=k; j++)
1074 : {
1075 94528 : gel(W,j) = gadd(gel(W,j), mulreal(gel(A,j+1),gel(RU,j+1)));
1076 94528 : gel(U,j) = gadd(gel(U,j), mulreal(gel(B,j),gel(RU,j)));
1077 : }
1078 : }
1079 : }
1080 : else
1081 : {
1082 50504 : for (j=1; j<=k; j++)
1083 : {
1084 33613 : gel(W,j) = gadd(gel(W,j), gmul(gel(A,j+1),gel(RU,j+1)));
1085 33613 : gel(U,j) = gadd(gel(U,j), gmul(gel(B,j),gel(RU,j)));
1086 : }
1087 : }
1088 64836 : prim2 = gmul(prim2,prim);
1089 64836 : gerepileall(ltop,3, &W,&U,&prim2);
1090 : }
1091 :
1092 90480 : for (i=1; i<=k; i++)
1093 : {
1094 61521 : aux=gel(W,i);
1095 61521 : for (j=1; j<i; j++) aux = gadd(aux, gmul(gel(W,i-j),gel(F,k+2-j)));
1096 61521 : gel(F,k+2-i) = gdivgs(aux,-i*NN);
1097 : }
1098 90480 : for (i=0; i<k; i++)
1099 : {
1100 61521 : aux=gel(U,k-i);
1101 61521 : for (j=1+i; j<k; j++) aux = gadd(aux,gmul(gel(F,2+j),gel(U,j-i)));
1102 61521 : gel(H,i+2) = gdivgs(aux,NN);
1103 : }
1104 28959 : }
1105 :
1106 : #define NEWTON_MAX 10
1107 : static GEN
1108 146124 : refine_H(GEN F, GEN G, GEN HH, long bit, long Sbit)
1109 : {
1110 146124 : GEN H = HH, D, aux;
1111 146124 : pari_sp ltop = avma;
1112 : long error, i, bit1, bit2;
1113 :
1114 146124 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F)); error = gexpo(D);
1115 146124 : bit2 = bit + Sbit;
1116 273430 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
1117 : {
1118 127306 : if (gc_needed(ltop,1))
1119 : {
1120 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_H");
1121 0 : gerepileall(ltop,2, &D,&H);
1122 : }
1123 127306 : bit1 = -error + Sbit;
1124 127306 : aux = RgX_mul(mygprec(H,bit1), mygprec(D,bit1));
1125 127306 : aux = RgX_rem(mygprec(aux,bit1), mygprec(F,bit1));
1126 :
1127 127306 : bit1 = -error*2 + Sbit; if (bit1 > bit2) bit1 = bit2;
1128 127306 : H = RgX_add(mygprec(H,bit1), aux);
1129 127306 : D = Rg_RgX_sub(gen_1, RgX_rem(RgX_mul(H,G),F));
1130 127306 : error = gexpo(D); if (error < -bit1) error = -bit1;
1131 : }
1132 146124 : if (error > -bit/2) return NULL; /* FAIL */
1133 146067 : return gerepilecopy(ltop,H);
1134 : }
1135 :
1136 : /* return 0 if fails, 1 else */
1137 : static long
1138 28959 : refine_F(GEN p, GEN *F, GEN *G, GEN H, long bit, double gamma)
1139 : {
1140 28959 : GEN f0, FF, GG, r, HH = H;
1141 28959 : long error, i, bit1 = 0, bit2, Sbit, Sbit2, enh, normF, normG, n = degpol(p);
1142 28959 : pari_sp av = avma;
1143 :
1144 28959 : FF = *F; GG = RgX_divrem(p, FF, &r);
1145 28959 : error = gexpo(r); if (error <= -bit) error = 1-bit;
1146 28959 : normF = gexpo(FF);
1147 28959 : normG = gexpo(GG);
1148 28959 : enh = gexpo(H); if (enh < 0) enh = 0;
1149 28959 : Sbit = normF + 2*normG + enh + (long)(4.*log2((double)n)+gamma) + 1;
1150 28959 : Sbit2 = enh + 2*(normF+normG) + (long)(2.*gamma+5.*log2((double)n)) + 1;
1151 28959 : bit2 = bit + Sbit;
1152 175026 : for (i=0; error>-bit && i<NEWTON_MAX && error<=0; i++)
1153 : {
1154 146124 : if (bit1 == bit2 && i >= 2) { Sbit += n; Sbit2 += n; bit2 += n; }
1155 146124 : if (gc_needed(av,1))
1156 : {
1157 0 : if(DEBUGMEM>1) pari_warn(warnmem,"refine_F");
1158 0 : gerepileall(av,4, &FF,&GG,&r,&HH);
1159 : }
1160 :
1161 146124 : bit1 = -error + Sbit2;
1162 146124 : HH = refine_H(mygprec(FF,bit1), mygprec(GG,bit1), mygprec(HH,bit1),
1163 : 1-error, Sbit2);
1164 146124 : if (!HH) return 0; /* FAIL */
1165 :
1166 146067 : bit1 = -error + Sbit;
1167 146067 : r = RgX_mul(mygprec(HH,bit1), mygprec(r,bit1));
1168 146067 : f0 = RgX_rem(mygprec(r,bit1), mygprec(FF,bit1));
1169 :
1170 146067 : bit1 = -2*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1171 146067 : FF = gadd(mygprec(FF,bit1),f0);
1172 :
1173 146067 : bit1 = -3*error + Sbit; if (bit1 > bit2) bit1 = bit2;
1174 146067 : GG = RgX_divrem(mygprec(p,bit1), mygprec(FF,bit1), &r);
1175 146067 : error = gexpo(r); if (error < -bit1) error = -bit1;
1176 : }
1177 28902 : if (error>-bit) return 0; /* FAIL */
1178 26737 : *F = FF; *G = GG; return 1;
1179 : }
1180 :
1181 : /* returns F and G from the unit circle U such that |p-FG|<2^(-bit) |cd|,
1182 : where cd is the leading coefficient of p */
1183 : static void
1184 26737 : split_fromU(GEN p, long k, double delta, long bit,
1185 : GEN *F, GEN *G, double param, double param2)
1186 : {
1187 : GEN pp, FF, GG, H;
1188 26737 : long n = degpol(p), NN, bit2, Lmax;
1189 26737 : int polreal = isreal(p);
1190 : pari_sp ltop;
1191 : double mu, gamma;
1192 :
1193 26737 : pp = gdiv(p, gel(p,2+n));
1194 26737 : parameters(pp, &Lmax,&mu,&gamma, polreal,param,param2);
1195 :
1196 26737 : H = cgetg(k+2,t_POL); H[1] = p[1];
1197 26737 : FF = cgetg(k+3,t_POL); FF[1]= p[1];
1198 26737 : gel(FF,k+2) = gen_1;
1199 :
1200 26737 : NN = (long)(0.5/delta); NN |= 1; if (NN < 2) NN = 2;
1201 26737 : NN *= Lmax; ltop = avma;
1202 : for(;;)
1203 : {
1204 28959 : bit2 = (long)(((double)NN*delta-mu)/LOG2) + gexpo(pp) + 8;
1205 28959 : dft(pp, k, NN, Lmax, bit2, FF, H, polreal);
1206 28959 : if (refine_F(pp,&FF,&GG,H,bit,gamma)) break;
1207 2222 : NN <<= 1; avma = ltop;
1208 2222 : }
1209 26737 : *G = gmul(GG,gel(p,2+n)); *F = FF;
1210 26737 : }
1211 :
1212 : static void
1213 26737 : optimize_split(GEN p, long k, double delta, long bit,
1214 : GEN *F, GEN *G, double param, double param2)
1215 : {
1216 26737 : long n = degpol(p);
1217 : GEN FF, GG;
1218 :
1219 26737 : if (k <= n/2)
1220 20543 : split_fromU(p,k,delta,bit,F,G,param,param2);
1221 : else
1222 : {
1223 6194 : split_fromU(RgX_recip_shallow(p),n-k,delta,bit,&FF,&GG,param,param2);
1224 6194 : *F = RgX_recip_shallow(GG);
1225 6194 : *G = RgX_recip_shallow(FF);
1226 : }
1227 26737 : }
1228 :
1229 : /********************************************************************/
1230 : /** **/
1231 : /** SEARCH FOR SEPARATING CIRCLE **/
1232 : /** **/
1233 : /********************************************************************/
1234 :
1235 : /* return p(2^e*x) *2^(-n*e) */
1236 : static void
1237 0 : scalepol2n(GEN p, long e)
1238 : {
1239 0 : long i,n=lg(p)-1;
1240 0 : for (i=2; i<=n; i++) gel(p,i) = gmul2n(gel(p,i),(i-n)*e);
1241 0 : }
1242 :
1243 : /* returns p(x/R)*R^n */
1244 : static GEN
1245 228195 : scalepol(GEN p, GEN R, long bit)
1246 : {
1247 : GEN q,aux,gR;
1248 : long i;
1249 :
1250 228195 : aux = gR = mygprec(R,bit); q = mygprec(p,bit);
1251 1127335 : for (i=lg(p)-2; i>=2; i--)
1252 : {
1253 899140 : gel(q,i) = gmul(aux,gel(q,i));
1254 899140 : aux = gmul(aux,gR);
1255 : }
1256 228195 : return q;
1257 : }
1258 :
1259 : /* return (conj(a)X-1)^n * p[ (X-a) / (conj(a)X-1) ] */
1260 : static GEN
1261 73992 : conformal_pol(GEN p, GEN a)
1262 : {
1263 73992 : GEN z, r, ma = gneg(a), ca = gconj(a);
1264 73992 : long n = degpol(p), i;
1265 73992 : pari_sp av = avma;
1266 :
1267 73992 : z = mkpoln(2, ca, gen_m1);
1268 73992 : r = scalarpol(gel(p,2+n), 0);
1269 285198 : for (i=n-1; ; i--)
1270 : {
1271 285198 : r = addmulXn(r, gmul(ma,r), 1); /* r *= (X - a) */
1272 285198 : r = gadd(r, gmul(z, gel(p,2+i)));
1273 359190 : if (i == 0) return gerepileupto(av, r);
1274 211206 : z = addmulXn(gmul(z,ca), gneg(z), 1); /* z *= conj(a)X - 1 */
1275 211206 : if (gc_needed(av,2))
1276 : {
1277 0 : if(DEBUGMEM>1) pari_warn(warnmem,"conformal_pol");
1278 0 : gerepileall(av,2, &r,&z);
1279 : }
1280 211206 : }
1281 : }
1282 :
1283 : static const double UNDEF = -100000.;
1284 :
1285 : static double
1286 26737 : logradius(double *radii, GEN p, long k, double aux, double *delta)
1287 : {
1288 26737 : long i, n = degpol(p);
1289 : double lrho, lrmin, lrmax;
1290 26737 : if (k > 1)
1291 : {
1292 18321 : i = k-1; while (i>0 && radii[i] == UNDEF) i--;
1293 18321 : lrmin = logpre_modulus(p,k,aux, radii[i], radii[k]);
1294 : }
1295 : else /* k=1 */
1296 8416 : lrmin = logmin_modulus(p,aux);
1297 26737 : radii[k] = lrmin;
1298 :
1299 26737 : if (k+1<n)
1300 : {
1301 22828 : i = k+2; while (i<=n && radii[i] == UNDEF) i++;
1302 22828 : lrmax = logpre_modulus(p,k+1,aux, radii[k+1], radii[i]);
1303 : }
1304 : else /* k+1=n */
1305 3909 : lrmax = logmax_modulus(p,aux);
1306 26737 : radii[k+1] = lrmax;
1307 :
1308 26737 : lrho = radii[k];
1309 65275 : for (i=k-1; i>=1; i--)
1310 : {
1311 38538 : if (radii[i] == UNDEF || radii[i] > lrho)
1312 26880 : radii[i] = lrho;
1313 : else
1314 11658 : lrho = radii[i];
1315 : }
1316 26737 : lrho = radii[k+1];
1317 125834 : for (i=k+1; i<=n; i++)
1318 : {
1319 99097 : if (radii[i] == UNDEF || radii[i] < lrho)
1320 62159 : radii[i] = lrho;
1321 : else
1322 36938 : lrho = radii[i];
1323 : }
1324 26737 : *delta = (lrmax - lrmin) / 2;
1325 26737 : if (*delta > 1.) *delta = 1.;
1326 26737 : return (lrmin + lrmax) / 2;
1327 : }
1328 :
1329 : static void
1330 26737 : update_radius(long n, double *radii, double lrho, double *par, double *par2)
1331 : {
1332 26737 : double t, param = 0., param2 = 0.;
1333 : long i;
1334 191109 : for (i=1; i<=n; i++)
1335 : {
1336 164372 : radii[i] -= lrho;
1337 164372 : t = fabs(rtodbl( invr(subsr(1, dblexp(radii[i]))) ));
1338 164372 : param += t; if (t > 1.) param2 += log2(t);
1339 : }
1340 26737 : *par = param; *par2 = param2;
1341 26737 : }
1342 :
1343 : /* apply the conformal mapping then split from U */
1344 : static void
1345 24664 : conformal_mapping(double *radii, GEN ctr, GEN p, long k, long bit,
1346 : double aux, GEN *F,GEN *G)
1347 : {
1348 24664 : long bit2, n = degpol(p), i;
1349 24664 : pari_sp ltop = avma, av;
1350 : GEN q, FF, GG, a, R;
1351 : double lrho, delta, param, param2;
1352 : /* n * (2.*log2(2.732)+log2(1.5)) + 1 */
1353 24664 : bit2 = bit + (long)(n*3.4848775) + 1;
1354 24664 : a = sqrtr_abs( stor(3, 2*MEDDEFAULTPREC - 2) );
1355 24664 : a = divrs(a, -6);
1356 24664 : a = gmul(mygprec(a,bit2), mygprec(ctr,bit2)); /* a = -ctr/2sqrt(3) */
1357 :
1358 24664 : av = avma;
1359 24664 : q = conformal_pol(mygprec(p,bit2), a);
1360 167263 : for (i=1; i<=n; i++)
1361 142599 : if (radii[i] != UNDEF) /* update array radii */
1362 : {
1363 93140 : pari_sp av2 = avma;
1364 93140 : GEN t, r = dblexp(radii[i]), r2 = sqrr(r);
1365 : /* 2(r^2 - 1) / (r^2 - 3(r-1)) */
1366 93140 : t = divrr(shiftr((subrs(r2,1)),1), subrr(r2, mulur(3,subrs(r,1))));
1367 93140 : radii[i] = mydbllogr(addsr(1,t)) / 2;
1368 93140 : avma = av2;
1369 : }
1370 24664 : lrho = logradius(radii, q,k,aux/10., &delta);
1371 24664 : update_radius(n, radii, lrho, ¶m, ¶m2);
1372 :
1373 24664 : bit2 += (long)(n * fabs(lrho)/LOG2 + 1.);
1374 24664 : R = mygprec(dblexp(-lrho), bit2);
1375 24664 : q = scalepol(q,R,bit2);
1376 24664 : gerepileall(av,2, &q,&R);
1377 :
1378 24664 : optimize_split(q,k,delta,bit2,&FF,&GG,param,param2);
1379 24664 : bit2 += n; R = invr(R);
1380 24664 : FF = scalepol(FF,R,bit2);
1381 24664 : GG = scalepol(GG,R,bit2);
1382 :
1383 24664 : a = mygprec(a,bit2);
1384 24664 : FF = conformal_pol(FF,a);
1385 24664 : GG = conformal_pol(GG,a);
1386 :
1387 24664 : a = invr(subsr(1, gnorm(a)));
1388 24664 : FF = RgX_Rg_mul(FF, powru(a,k));
1389 24664 : GG = RgX_Rg_mul(GG, powru(a,n-k));
1390 :
1391 24664 : *F = mygprec(FF,bit+n);
1392 24664 : *G = mygprec(GG,bit+n); gerepileall(ltop,2, F,G);
1393 24664 : }
1394 :
1395 : /* split p, this time without scaling. returns in F and G two polynomials
1396 : * such that |p-FG|< 2^(-bit)|p| */
1397 : static void
1398 26737 : split_2(GEN p, long bit, GEN ctr, double thickness, GEN *F, GEN *G)
1399 : {
1400 : GEN q, FF, GG, R;
1401 : double aux, delta, param, param2;
1402 26737 : long n = degpol(p), i, j, k, bit2;
1403 : double lrmin, lrmax, lrho, *radii;
1404 :
1405 26737 : init_dalloc();
1406 26737 : radii = (double*) stack_malloc((n+1) * sizeof(double));
1407 :
1408 26737 : for (i=2; i<n; i++) radii[i] = UNDEF;
1409 26737 : aux = thickness/(double)(4 * n);
1410 26737 : lrmin = logmin_modulus(p, aux);
1411 26737 : lrmax = logmax_modulus(p, aux);
1412 26737 : radii[1] = lrmin;
1413 26737 : radii[n] = lrmax;
1414 26737 : i = 1; j = n;
1415 26737 : lrho = (lrmin + lrmax) / 2;
1416 26737 : k = dual_modulus(p, lrho, aux, 1);
1417 26737 : if (5*k < n || (n < 2*k && 5*k < 4*n))
1418 5217 : { lrmax = lrho; j=k+1; radii[j] = lrho; }
1419 : else
1420 21520 : { lrmin = lrho; i=k; radii[i] = lrho; }
1421 85799 : while (j > i+1)
1422 : {
1423 32325 : if (i+j == n+1)
1424 11536 : lrho = (lrmin + lrmax) / 2;
1425 : else
1426 : {
1427 20789 : double kappa = 2. - log(1. + minss(i,n-j)) / log(1. + minss(j,n-i));
1428 20789 : if (i+j < n+1) lrho = lrmax * kappa + lrmin;
1429 14878 : else lrho = lrmin * kappa + lrmax;
1430 20789 : lrho /= 1+kappa;
1431 : }
1432 32325 : aux = (lrmax - lrmin) / (4*(j-i));
1433 32325 : k = dual_modulus(p, lrho, aux, minss(i,n+1-j));
1434 32325 : if (k-i < j-k-1 || (k-i == j-k-1 && 2*k > n))
1435 22263 : { lrmax = lrho; j=k+1; radii[j] = lrho - aux; }
1436 : else
1437 10062 : { lrmin = lrho; i=k; radii[i] = lrho + aux; }
1438 : }
1439 26737 : aux = lrmax - lrmin;
1440 :
1441 26737 : if (ctr)
1442 : {
1443 24664 : lrho = (lrmax + lrmin) / 2;
1444 167263 : for (i=1; i<=n; i++)
1445 142599 : if (radii[i] != UNDEF) radii[i] -= lrho;
1446 :
1447 24664 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1448 24664 : R = mygprec(dblexp(-lrho), bit2);
1449 24664 : q = scalepol(p,R,bit2);
1450 24664 : conformal_mapping(radii, ctr, q, k, bit2, aux, &FF, &GG);
1451 : }
1452 : else
1453 : {
1454 2073 : lrho = logradius(radii, p, k, aux/10., &delta);
1455 2073 : update_radius(n, radii, lrho, ¶m, ¶m2);
1456 :
1457 2073 : bit2 = bit + (long)(n * fabs(lrho)/LOG2 + 1.);
1458 2073 : R = mygprec(dblexp(-lrho), bit2);
1459 2073 : q = scalepol(p,R,bit2);
1460 2073 : optimize_split(q, k, delta, bit2, &FF, &GG, param, param2);
1461 : }
1462 26737 : bit += n;
1463 26737 : bit2 += n; R = invr(mygprec(R,bit2));
1464 26737 : *F = mygprec(scalepol(FF,R,bit2), bit);
1465 26737 : *G = mygprec(scalepol(GG,R,bit2), bit);
1466 26737 : }
1467 :
1468 : /* procedure corresponding to steps 5,6,.. page 44 in RR n. 1852 */
1469 : /* put in F and G two polynomial such that |p-FG|<2^(-bit)|p|
1470 : * where the maximum modulus of the roots of p is <=1.
1471 : * Assume sum of roots is 0. */
1472 : static void
1473 24664 : split_1(GEN p, long bit, GEN *F, GEN *G)
1474 : {
1475 24664 : long i, imax, n = degpol(p), polreal = isreal(p), ep = gexpo(p), bit2 = bit+n;
1476 : GEN ctr, q, qq, FF, GG, v, gr, r, newq;
1477 : double lrmin, lrmax, lthick;
1478 24664 : const double LOG3 = 1.098613;
1479 :
1480 24664 : lrmax = logmax_modulus(p, 0.01);
1481 24664 : gr = mygprec(dblexp(-lrmax), bit2);
1482 24664 : q = scalepol(p,gr,bit2);
1483 :
1484 24664 : bit2 = bit + gexpo(q) - ep + (long)((double)n*2.*log2(3.)+1);
1485 24664 : v = cgetg(5,t_VEC);
1486 24664 : gel(v,1) = gen_2;
1487 24664 : gel(v,2) = gen_m2;
1488 24664 : gel(v,3) = mkcomplex(gen_0, gel(v,1));
1489 24664 : gel(v,4) = mkcomplex(gen_0, gel(v,2));
1490 24664 : q = mygprec(q,bit2); lthick = 0;
1491 24664 : newq = ctr = NULL; /* -Wall */
1492 24664 : imax = polreal? 3: 4;
1493 45236 : for (i=1; i<=imax; i++)
1494 : {
1495 44802 : qq = RgX_translate(q, gel(v,i));
1496 44802 : lrmin = logmin_modulus(qq,0.05);
1497 44802 : if (LOG3 > lrmin + lthick)
1498 : {
1499 43767 : double lquo = logmax_modulus(qq,0.05) - lrmin;
1500 43767 : if (lquo > lthick) { lthick = lquo; newq = qq; ctr = gel(v,i); }
1501 : }
1502 44802 : if (lthick > LOG2) break;
1503 24392 : if (polreal && i==2 && lthick > LOG3 - LOG2) break;
1504 : }
1505 24664 : bit2 = bit + gexpo(newq) - ep + (long)(n*LOG3/LOG2 + 1);
1506 24664 : split_2(newq, bit2, ctr, lthick, &FF, &GG);
1507 24664 : r = gneg(mygprec(ctr,bit2));
1508 24664 : FF = RgX_translate(FF,r);
1509 24664 : GG = RgX_translate(GG,r);
1510 :
1511 24664 : gr = invr(gr); bit2 = bit - ep + gexpo(FF)+gexpo(GG);
1512 24664 : *F = scalepol(FF,gr,bit2);
1513 24664 : *G = scalepol(GG,gr,bit2);
1514 24664 : }
1515 :
1516 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|,
1517 : where the maximum modulus of the roots of p is < 0.5 */
1518 : static int
1519 24799 : split_0_2(GEN p, long bit, GEN *F, GEN *G)
1520 : {
1521 : GEN q, b, FF, GG;
1522 24799 : long n = degpol(p), k, bit2, eq;
1523 24799 : double aux = dbllog2(gel(p,n+1)) - dbllog2(gel(p,n+2));
1524 :
1525 : /* beware double overflow */
1526 24799 : if (aux >= 0 && (aux > 1e4 || exp2(aux) > 2.5*n)) return 0;
1527 :
1528 24799 : aux = (aux < -300)? 0.: n*log2(1 + exp2(aux)/(double)n);
1529 24799 : bit2 = bit+1 + (long)(log2((double)n) + aux);
1530 :
1531 24799 : q = mygprec(p,bit2);
1532 24799 : b = gdivgs(gdiv(gel(q,n+1),gel(q,n+2)),-n);
1533 24799 : q = RgX_translate(q,b); gel(q,n+1) = gen_0; eq = gexpo(q);
1534 24799 : k = 0;
1535 49963 : while (k <= n/2 && (- gexpo(gel(q,k+2)) > bit2 + 2*(n-k) + eq
1536 25042 : || gequal0(gel(q,k+2)))) k++;
1537 24799 : if (k > 0)
1538 : {
1539 135 : if (k > n/2) k = n/2;
1540 135 : bit2 += k<<1;
1541 135 : FF = pol_xn(k, 0);
1542 135 : GG = RgX_shift_shallow(q, -k);
1543 : }
1544 : else
1545 : {
1546 24664 : split_1(q,bit2,&FF,&GG);
1547 24664 : bit2 = bit + gexpo(FF) + gexpo(GG) - gexpo(p) + (long)aux+1;
1548 24664 : FF = mygprec(FF,bit2);
1549 : }
1550 24799 : GG = mygprec(GG,bit2); b = mygprec(gneg(b),bit2);
1551 24799 : *F = RgX_translate(FF, b);
1552 24799 : *G = RgX_translate(GG, b); return 1;
1553 : }
1554 :
1555 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P|.
1556 : * Assume max_modulus(p) < 2 */
1557 : static void
1558 24799 : split_0_1(GEN p, long bit, GEN *F, GEN *G)
1559 : {
1560 : GEN FF, GG;
1561 : long n, bit2, normp;
1562 :
1563 49598 : if (split_0_2(p,bit,F,G)) return;
1564 :
1565 0 : normp = gexpo(p);
1566 0 : scalepol2n(p,2); /* p := 4^(-n) p(4*x) */
1567 0 : n = degpol(p); bit2 = bit + 2*n + gexpo(p) - normp;
1568 0 : split_1(mygprec(p,bit2), bit2,&FF,&GG);
1569 0 : scalepol2n(FF,-2);
1570 0 : scalepol2n(GG,-2); bit2 = bit + gexpo(FF) + gexpo(GG) - normp;
1571 0 : *F = mygprec(FF,bit2);
1572 0 : *G = mygprec(GG,bit2);
1573 : }
1574 :
1575 : /* put in F and G two polynomials such that |P-FG|<2^(-bit)|P| */
1576 : static void
1577 26872 : split_0(GEN p, long bit, GEN *F, GEN *G)
1578 : {
1579 26872 : const double LOG1_9 = 0.6418539;
1580 26872 : long n = degpol(p), k = 0;
1581 : GEN q;
1582 :
1583 26872 : while (gexpo(gel(p,k+2)) < -bit && k <= n/2) k++;
1584 26872 : if (k > 0)
1585 : {
1586 0 : if (k > n/2) k = n/2;
1587 0 : *F = pol_xn(k, 0);
1588 0 : *G = RgX_shift_shallow(p, -k);
1589 : }
1590 : else
1591 : {
1592 26872 : double lr = logmax_modulus(p, 0.05);
1593 26872 : if (lr < LOG1_9) split_0_1(p, bit, F, G);
1594 : else
1595 : {
1596 20744 : q = RgX_recip_shallow(p);
1597 20744 : lr = logmax_modulus(q,0.05);
1598 20744 : if (lr < LOG1_9)
1599 : {
1600 18671 : split_0_1(q, bit, F, G);
1601 18671 : *F = RgX_recip_shallow(*F);
1602 18671 : *G = RgX_recip_shallow(*G);
1603 : }
1604 : else
1605 2073 : split_2(p,bit,NULL, 1.2837,F,G);
1606 : }
1607 : }
1608 26872 : }
1609 :
1610 : /********************************************************************/
1611 : /** **/
1612 : /** ERROR ESTIMATE FOR THE ROOTS **/
1613 : /** **/
1614 : /********************************************************************/
1615 :
1616 : static GEN
1617 86117 : root_error(long n, long k, GEN roots_pol, long err, GEN shatzle)
1618 : {
1619 86117 : GEN rho, d, eps, epsbis, eps2, aux, rap = NULL;
1620 : long i, j;
1621 :
1622 86117 : d = cgetg(n+1,t_VEC);
1623 1348320 : for (i=1; i<=n; i++)
1624 : {
1625 1262203 : if (i!=k)
1626 : {
1627 1176086 : aux = gsub(gel(roots_pol,i), gel(roots_pol,k));
1628 1176086 : gel(d,i) = gabs(mygprec(aux,31), DEFAULTPREC);
1629 : }
1630 : }
1631 86117 : rho = gabs(mygprec(gel(roots_pol,k),31), DEFAULTPREC);
1632 86117 : if (expo(rho) < 0) rho = real_1(DEFAULTPREC);
1633 86117 : eps = mulrr(rho, shatzle);
1634 86117 : aux = shiftr(powru(rho,n), err);
1635 :
1636 271891 : for (j=1; j<=2 || (j<=5 && cmprr(rap, dbltor(1.2)) > 0); j++)
1637 : {
1638 185774 : GEN prod = NULL; /* 1. */
1639 185774 : long m = n;
1640 185774 : epsbis = mulrr(eps, dbltor(1.25));
1641 3320879 : for (i=1; i<=n; i++)
1642 : {
1643 3135105 : if (i != k && cmprr(gel(d,i),epsbis) > 0)
1644 : {
1645 2918895 : GEN dif = subrr(gel(d,i),eps);
1646 2918895 : prod = prod? mulrr(prod, dif): dif;
1647 2918895 : m--;
1648 : }
1649 : }
1650 185774 : eps2 = prod? divrr(aux, prod): aux;
1651 185774 : if (m > 1) eps2 = sqrtnr(shiftr(eps2, 2*m-2), m);
1652 185774 : rap = divrr(eps,eps2); eps = eps2;
1653 : }
1654 86117 : return eps;
1655 : }
1656 :
1657 : /* round a complex or real number x to an absolute value of 2^(-bit) */
1658 : static GEN
1659 204021 : mygprec_absolute(GEN x, long bit)
1660 : {
1661 : long e;
1662 : GEN y;
1663 :
1664 204021 : switch(typ(x))
1665 : {
1666 : case t_REAL:
1667 137943 : e = expo(x) + bit;
1668 137943 : return (e <= 0 || !signe(x))? real_0_bit(-bit): rtor(x, nbits2prec(e));
1669 : case t_COMPLEX:
1670 59742 : if (gexpo(gel(x,2)) < -bit) return mygprec_absolute(gel(x,1),bit);
1671 58162 : y = cgetg(3,t_COMPLEX);
1672 58162 : gel(y,1) = mygprec_absolute(gel(x,1),bit);
1673 58162 : gel(y,2) = mygprec_absolute(gel(x,2),bit);
1674 58162 : return y;
1675 6336 : default: return x;
1676 : }
1677 : }
1678 :
1679 : static long
1680 13096 : a_posteriori_errors(GEN p, GEN roots_pol, long err)
1681 : {
1682 13096 : long i, n = degpol(p), e_max = -(long)EXPOBITS;
1683 : GEN sigma, shatzle;
1684 :
1685 13096 : err += (long)log2((double)n) + 1;
1686 13096 : if (err > -2) return 0;
1687 13096 : sigma = real2n(-err, LOWDEFAULTPREC);
1688 : /* 2 / ((s - 1)^(1/n) - 1) */
1689 13096 : shatzle = divur(2, subrs(sqrtnr(subrs(sigma,1),n), 1));
1690 99213 : for (i=1; i<=n; i++)
1691 : {
1692 86117 : pari_sp av = avma;
1693 86117 : GEN x = root_error(n,i,roots_pol,err,shatzle);
1694 86117 : long e = gexpo(x);
1695 86117 : avma = av; if (e > e_max) e_max = e;
1696 86117 : gel(roots_pol,i) = mygprec_absolute(gel(roots_pol,i), -e);
1697 : }
1698 13096 : return e_max;
1699 : }
1700 :
1701 : /********************************************************************/
1702 : /** **/
1703 : /** MAIN **/
1704 : /** **/
1705 : /********************************************************************/
1706 : static GEN
1707 64602 : append_clone(GEN r, GEN a) { a = gclone(a); vectrunc_append(r, a); return a; }
1708 :
1709 : /* put roots in placeholder roots_pol so that |P - L_1...L_n| < 2^(-bit)|P|
1710 : * returns prod (x-roots_pol[i]) */
1711 : static GEN
1712 66840 : split_complete(GEN p, long bit, GEN roots_pol)
1713 : {
1714 66840 : long n = degpol(p);
1715 : pari_sp ltop;
1716 : GEN p1, F, G, a, b, m1, m2;
1717 :
1718 66840 : if (n == 1)
1719 : {
1720 15334 : a = gneg_i(gdiv(gel(p,2), gel(p,3)));
1721 15334 : (void)append_clone(roots_pol,a); return p;
1722 : }
1723 51506 : ltop = avma;
1724 51506 : if (n == 2)
1725 : {
1726 24634 : F = gsub(gsqr(gel(p,3)), gmul2n(gmul(gel(p,2),gel(p,4)), 2));
1727 24634 : F = gsqrt(F, nbits2prec(bit));
1728 24634 : p1 = ginv(gmul2n(gel(p,4),1));
1729 24634 : a = gneg_i(gmul(gadd(F,gel(p,3)), p1));
1730 24634 : b = gmul(gsub(F,gel(p,3)), p1);
1731 24634 : a = append_clone(roots_pol,a);
1732 24634 : b = append_clone(roots_pol,b); avma = ltop;
1733 24634 : a = mygprec(a, 3*bit);
1734 24634 : b = mygprec(b, 3*bit);
1735 24634 : return gmul(gel(p,4), mkpoln(3, gen_1, gneg(gadd(a,b)), gmul(a,b)));
1736 : }
1737 26872 : split_0(p,bit,&F,&G);
1738 26872 : m1 = split_complete(F,bit,roots_pol);
1739 26872 : m2 = split_complete(G,bit,roots_pol);
1740 26872 : return gerepileupto(ltop, gmul(m1,m2));
1741 : }
1742 :
1743 : static GEN
1744 374884 : quicktofp(GEN x)
1745 : {
1746 374884 : const long prec = DEFAULTPREC;
1747 374884 : switch(typ(x))
1748 : {
1749 371531 : case t_INT: return itor(x, prec);
1750 3332 : case t_REAL: return rtor(x, prec);
1751 0 : case t_FRAC: return fractor(x, prec);
1752 : case t_COMPLEX: {
1753 21 : GEN a = gel(x,1), b = gel(x,2);
1754 : /* avoid problem with 0, e.g. x = 0 + I*1e-100. We don't want |x| = 0. */
1755 21 : if (isintzero(a)) return cxcompotor(b, prec);
1756 21 : if (isintzero(b)) return cxcompotor(a, prec);
1757 21 : a = cxcompotor(a, prec);
1758 21 : b = cxcompotor(b, prec); return sqrtr(addrr(sqrr(a), sqrr(b)));
1759 : }
1760 0 : default: pari_err_TYPE("quicktofp",x);
1761 0 : return NULL;/*not reached*/
1762 : }
1763 :
1764 : }
1765 :
1766 : /* bound log_2 |largest root of p| (Fujiwara's bound) */
1767 : double
1768 77515 : fujiwara_bound(GEN p)
1769 : {
1770 77515 : pari_sp av = avma;
1771 77515 : long i, n = degpol(p);
1772 : GEN cc;
1773 : double loglc, Lmax;
1774 :
1775 77515 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1776 77515 : loglc = mydbllog2r( quicktofp(gel(p,n+2)) ); /* log_2 |lc(p)| */
1777 77515 : cc = gel(p, 2);
1778 77515 : if (gequal0(cc))
1779 7264 : Lmax = -pariINFINITY-1;
1780 : else
1781 70251 : Lmax = (mydbllog2r(quicktofp(cc)) - loglc - 1) / n;
1782 360916 : for (i = 1; i < n; i++)
1783 : {
1784 283401 : GEN y = gel(p,i+2);
1785 : double L;
1786 283401 : if (gequal0(y)) continue;
1787 227118 : L = (mydbllog2r(quicktofp(y)) - loglc) / (n-i);
1788 227118 : if (L > Lmax) Lmax = L;
1789 : }
1790 77515 : avma = av; return Lmax + 1;
1791 : }
1792 :
1793 : /* Fujiwara's bound, real roots. Based on the following remark: if
1794 : * p = x^n + sum a_i x^i and q = x^n + sum min(a_i,0)x^i
1795 : * then for all x >= 0, p(x) >= q(x). Thus any bound for the (positive) roots
1796 : * of q is a bound for the positive roots of p. */
1797 : double
1798 9821 : fujiwara_bound_real(GEN p, long sign)
1799 : {
1800 9821 : pari_sp av = avma;
1801 : GEN x;
1802 9821 : long n = degpol(p), i, signodd, signeven;
1803 : double fb;
1804 9821 : if (n <= 0) pari_err_CONSTPOL("fujiwara_bound");
1805 9821 : x = shallowcopy(p);
1806 9821 : if (gsigne(gel(x, n+2)) > 0)
1807 : {
1808 9821 : signeven = 1;
1809 9821 : signodd = sign;
1810 : }
1811 : else
1812 : {
1813 0 : signeven = -1;
1814 0 : signodd = -sign;
1815 : }
1816 57718 : for (i = 0; i < n; i++)
1817 : {
1818 47897 : if ((n - i) % 2)
1819 : {
1820 25938 : if (gsigne(gel(x, i+2)) == signodd ) gel(x, i+2) = gen_0;
1821 : }
1822 : else
1823 : {
1824 21959 : if (gsigne(gel(x, i+2)) == signeven) gel(x, i+2) = gen_0;
1825 : }
1826 : }
1827 9821 : fb = fujiwara_bound(x);
1828 9821 : avma = av; return fb;
1829 : }
1830 :
1831 : static GEN
1832 99255 : mygprecrc_special(GEN x, long prec, long e)
1833 : {
1834 : GEN y;
1835 99255 : switch(typ(x))
1836 : {
1837 : case t_REAL:
1838 3905 : if (!signe(x)) return real_0_bit(minss(e, expo(x)));
1839 3905 : return (prec > realprec(x))? rtor(x, prec): x;
1840 : case t_COMPLEX:
1841 21 : y = cgetg(3,t_COMPLEX);
1842 21 : gel(y,1) = mygprecrc_special(gel(x,1),prec,e);
1843 21 : gel(y,2) = mygprecrc_special(gel(x,2),prec,e);
1844 21 : return y;
1845 95329 : default: return x;
1846 : }
1847 : }
1848 :
1849 : /* like mygprec but keep at least the same precision as before */
1850 : static GEN
1851 13096 : mygprec_special(GEN x, long bit)
1852 : {
1853 : long lx, i, e, prec;
1854 : GEN y;
1855 :
1856 13096 : if (bit < 0) bit = 0; /* should not happen */
1857 13096 : e = gexpo(x) - bit;
1858 13096 : prec = nbits2prec(bit);
1859 13096 : switch(typ(x))
1860 : {
1861 : case t_POL:
1862 13096 : y = cgetg_copy(x, &lx); y[1] = x[1];
1863 13096 : for (i=2; i<lx; i++) gel(y,i) = mygprecrc_special(gel(x,i),prec,e);
1864 13096 : break;
1865 :
1866 0 : default: y = mygprecrc_special(x,prec,e);
1867 : }
1868 13096 : return y;
1869 : }
1870 :
1871 : static GEN
1872 7568 : fix_roots1(GEN r)
1873 : {
1874 7568 : long i, l = lg(r);
1875 7568 : GEN allr = cgetg(l, t_VEC);
1876 57657 : for (i=1; i<l; i++)
1877 : {
1878 50089 : GEN t = gel(r,i);
1879 50089 : gel(allr,i) = gcopy(t); gunclone(t);
1880 : }
1881 7568 : return allr;
1882 : }
1883 : static GEN
1884 13096 : fix_roots(GEN r, GEN *m, long h, long bit)
1885 : {
1886 : long i, j, k, l, prec;
1887 : GEN allr, ro1;
1888 13096 : if (h == 1) return fix_roots1(r);
1889 5528 : prec = nbits2prec(bit);
1890 5528 : ro1 = grootsof1(h, prec) + 1;
1891 5528 : l = lg(r)-1;
1892 5528 : allr = cgetg(h*l+1, t_VEC);
1893 20041 : for (k=1,i=1; i<=l; i++)
1894 : {
1895 14513 : GEN p2, p1 = gel(r,i);
1896 14513 : p2 = (h == 2)? gsqrt(p1, prec): gsqrtn(p1, utoipos(h), NULL, prec);
1897 14513 : for (j=0; j<h; j++) gel(allr,k++) = gmul(p2, gel(ro1,j));
1898 14513 : gunclone(p1);
1899 : }
1900 5528 : *m = roots_to_pol(allr, 0);
1901 5528 : return allr;
1902 : }
1903 :
1904 : static GEN
1905 364 : RgX_normalize1(GEN x)
1906 : {
1907 364 : long i, n = lg(x)-1;
1908 : GEN y;
1909 371 : for (i = n; i > 1; i--)
1910 371 : if (!gequal0( gel(x,i) )) break;
1911 364 : if (i == n) return x;
1912 7 : pari_warn(warner,"normalizing a polynomial with 0 leading term");
1913 7 : if (i == 1) pari_err_ROOTS0("roots");
1914 7 : y = cgetg(i+1, t_POL); y[1] = x[1];
1915 7 : for (; i > 1; i--) gel(y,i) = gel(x,i);
1916 7 : return y;
1917 : }
1918 :
1919 : static GEN
1920 12903 : all_roots(GEN p, long bit)
1921 : {
1922 : GEN lc, pd, q, roots_pol, m;
1923 12903 : long bit0, bit2, i, e, h, n = degpol(p);
1924 : pari_sp av;
1925 :
1926 12903 : pd = RgX_deflate_max(p, &h); lc = leading_coeff(pd);
1927 12903 : e = (long)(2 * fujiwara_bound(pd)); if (e < 0) e = 0;
1928 12903 : bit0 = bit + gexpo(pd) - gexpo(lc) + (long)log2(n/h)+1+e;
1929 12903 : bit2 = bit0; e = 0;
1930 13096 : for (av=avma,i=1;; i++,avma=av)
1931 : {
1932 13096 : roots_pol = vectrunc_init(n+1);
1933 13096 : bit2 += e + (n << i);
1934 13096 : q = RgX_gtofp_bit(mygprec(pd,bit2), bit2);
1935 13096 : q[1] = evalsigne(1)|evalvarn(0);
1936 13096 : m = split_complete(q,bit2,roots_pol);
1937 13096 : roots_pol = fix_roots(roots_pol, &m, h, bit2);
1938 13096 : q = mygprec_special(p,bit2); lc = leading_coeff(q);
1939 13096 : q[1] = evalsigne(1)|evalvarn(0);
1940 13096 : if (h > 1) m = gmul(m,lc);
1941 :
1942 13096 : e = gexpo(gsub(q, m)) - gexpo(lc) + (long)log2((double)n) + 1;
1943 13096 : if (e < -2*bit2) e = -2*bit2; /* avoid e = -pariINFINITY */
1944 13096 : if (e < 0)
1945 : {
1946 13096 : e = bit + a_posteriori_errors(p,roots_pol,e);
1947 25999 : if (e < 0) return roots_pol;
1948 : }
1949 193 : if (DEBUGLEVEL > 7)
1950 0 : err_printf("all_roots: restarting, i = %ld, e = %ld\n", i,e);
1951 193 : }
1952 : }
1953 :
1954 : INLINE int
1955 3644 : isexactscalar(GEN x) { long tx = typ(x); return is_rational_t(tx); }
1956 :
1957 : static int
1958 1177 : isexactpol(GEN p)
1959 : {
1960 1177 : long i,n = degpol(p);
1961 4457 : for (i=0; i<=n; i++)
1962 3644 : if (!isexactscalar(gel(p,i+2))) return 0;
1963 813 : return 1;
1964 : }
1965 :
1966 : static long
1967 7102 : isvalidcoeff(GEN x)
1968 : {
1969 7102 : switch (typ(x))
1970 : {
1971 7046 : case t_INT: case t_REAL: case t_FRAC: return 1;
1972 42 : case t_COMPLEX: return isvalidcoeff(gel(x,1)) && isvalidcoeff(gel(x,2));
1973 : }
1974 14 : return 0;
1975 : }
1976 :
1977 : static void
1978 1184 : checkvalidpol(GEN p)
1979 : {
1980 1184 : long i,n = lg(p);
1981 8181 : for (i=2; i<n; i++)
1982 7004 : if (!isvalidcoeff(gel(p,i))) pari_err_TYPE("roots", gel(p,i));
1983 1177 : }
1984 :
1985 : static GEN
1986 813 : solve_exact_pol(GEN p, long bit)
1987 : {
1988 813 : long i, j, k, m, n = degpol(p), iroots = 0;
1989 813 : GEN ex, factors, v = zerovec(n);
1990 :
1991 813 : factors = ZX_squff(Q_primpart(p), &ex);
1992 1626 : for (i=1; i<lg(factors); i++)
1993 : {
1994 813 : GEN roots_fact = all_roots(gel(factors,i), bit);
1995 813 : n = degpol(gel(factors,i));
1996 813 : m = ex[i];
1997 3238 : for (j=1; j<=n; j++)
1998 2425 : for (k=1; k<=m; k++) v[++iroots] = roots_fact[j];
1999 : }
2000 813 : return v;
2001 : }
2002 :
2003 : /* return the roots of p with absolute error bit */
2004 : static GEN
2005 1177 : roots_com(GEN q, long bit)
2006 : {
2007 : GEN L, p;
2008 1177 : long v = RgX_valrem_inexact(q, &p);
2009 1177 : int ex = isexactpol(p);
2010 1177 : if (!ex) p = RgX_normalize1(p);
2011 1177 : if (lg(p) == 3)
2012 7 : L = cgetg(1,t_VEC); /* constant polynomial */
2013 : else
2014 1170 : L = ex? solve_exact_pol(p,bit): all_roots(p,bit);
2015 1177 : if (v)
2016 : {
2017 105 : GEN M, z, t = gel(q,2);
2018 : long i, x, y, l, n;
2019 :
2020 105 : if (isrationalzero(t)) x = -bit;
2021 : else
2022 : {
2023 14 : n = gexpo(t);
2024 14 : x = n / v; l = degpol(q);
2025 56 : for (i = v; i <= l; i++)
2026 : {
2027 42 : t = gel(q,i+2);
2028 42 : if (isrationalzero(t)) continue;
2029 42 : y = (n - gexpo(t)) / i;
2030 42 : if (y < x) x = y;
2031 : }
2032 : }
2033 105 : z = real_0_bit(x); l = v + lg(L);
2034 105 : M = cgetg(l, t_VEC); L -= v;
2035 105 : for (i = 1; i <= v; i++) gel(M,i) = z;
2036 105 : for ( ; i < l; i++) gel(M,i) = gel(L,i);
2037 105 : L = M;
2038 : }
2039 1177 : return L;
2040 : }
2041 :
2042 : static GEN
2043 61272 : tocomplex(GEN x, long l, long bit)
2044 : {
2045 : GEN y;
2046 61272 : if (typ(x) == t_COMPLEX)
2047 : {
2048 57272 : if (signe(gel(x,1))) return mygprecrc(x, l, -bit);
2049 6046 : x = gel(x,2);
2050 6046 : y = cgetg(3,t_COMPLEX);
2051 6046 : gel(y,1) = real_0_bit(-bit);
2052 6046 : gel(y,2) = mygprecrc(x, l, -bit);
2053 : }
2054 : else
2055 : {
2056 4000 : y = cgetg(3,t_COMPLEX);
2057 4000 : gel(y,1) = mygprecrc(x, l, -bit);
2058 4000 : gel(y,2) = real_0_bit(-bit);
2059 : }
2060 10046 : return y;
2061 : }
2062 :
2063 : /* x,y are t_COMPLEX of t_REALs or t_REAL, compare lexicographically,
2064 : * up to 2^-e absolute error */
2065 : static int
2066 170283 : cmp_complex_appr(void *E, GEN x, GEN y)
2067 : {
2068 170283 : long e = (long)E;
2069 : GEN z, xi, yi, xr, yr;
2070 : long sxi, syi;
2071 170283 : if (typ(x) == t_COMPLEX) { xr = gel(x,1); xi = gel(x,2); sxi = signe(xi); }
2072 59241 : else { xr = x; xi = NULL; sxi = 0; }
2073 170283 : if (typ(y) == t_COMPLEX) { yr = gel(y,1); yi = gel(y,2); syi = signe(yi); }
2074 52861 : else { yr = y; yi = NULL; syi = 0; }
2075 : /* Compare absolute values of imaginary parts */
2076 170283 : if (!sxi)
2077 : {
2078 66664 : if (syi && expo(yi) >= e) return -1;
2079 : /* |Im x| ~ |Im y| ~ 0 */
2080 : }
2081 103619 : else if (!syi)
2082 : {
2083 2904 : if (sxi && expo(xi) >= e) return 1;
2084 : /* |Im x| ~ |Im y| ~ 0 */
2085 : }
2086 : else
2087 : {
2088 : long sz;
2089 100715 : z = addrr_sign(xi, 1, yi, -1);
2090 100715 : sz = signe(z);
2091 100715 : if (sz && expo(z) >= e) return (int)sz;
2092 : }
2093 : /* |Im x| ~ |Im y|, sort according to real parts */
2094 99181 : z = subrr(xr, yr);
2095 99181 : if (expo(z) >= e) return (int)signe(z);
2096 : /* Re x ~ Re y. Place negative absolute value before positive */
2097 25274 : return (int) (sxi - syi);
2098 : }
2099 :
2100 : static GEN
2101 12910 : clean_roots(GEN L, long l, long bit, long clean)
2102 : {
2103 12910 : long i, n = lg(L), ex = 5 - bit;
2104 12910 : GEN res = cgetg(n,t_COL);
2105 97245 : for (i=1; i<n; i++)
2106 : {
2107 84335 : GEN c = gel(L,i);
2108 84335 : if (clean && isrealappr(c,ex))
2109 : {
2110 23063 : if (typ(c) == t_COMPLEX) c = gel(c,1);
2111 23063 : c = mygprecrc(c, l, -bit);
2112 : }
2113 : else
2114 61272 : c = tocomplex(c, l, bit);
2115 84335 : gel(res,i) = c;
2116 : }
2117 12910 : gen_sort_inplace(res, (void*)ex, &cmp_complex_appr, NULL);
2118 12910 : return res;
2119 : }
2120 :
2121 : /* the vector of roots of p, with absolute error 2^(- prec2nbits(l)) */
2122 : static GEN
2123 1205 : roots_aux(GEN p, long l, long clean)
2124 : {
2125 1205 : pari_sp av = avma;
2126 : long bit;
2127 : GEN L;
2128 :
2129 1205 : if (typ(p) != t_POL)
2130 : {
2131 21 : if (gequal0(p)) pari_err_ROOTS0("roots");
2132 14 : if (!isvalidcoeff(p)) pari_err_TYPE("roots",p);
2133 7 : return cgetg(1,t_COL); /* constant polynomial */
2134 : }
2135 1184 : if (!signe(p)) pari_err_ROOTS0("roots");
2136 1184 : checkvalidpol(p);
2137 1177 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2138 1177 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2139 1177 : bit = prec2nbits(l);
2140 1177 : L = roots_com(p, bit);
2141 1177 : return gerepileupto(av, clean_roots(L, l, bit, clean));
2142 : }
2143 : GEN
2144 1058 : roots(GEN p, long l) { return roots_aux(p,l, 0); }
2145 : /* clean up roots. If root is real replace it by its real part */
2146 : GEN
2147 147 : cleanroots(GEN p, long l) { return roots_aux(p,l, 1); }
2148 :
2149 : /* private variant of conjvec. Allow non rational coefficients, shallow
2150 : * function. */
2151 : GEN
2152 77 : polmod_to_embed(GEN x, long prec)
2153 : {
2154 77 : GEN v, T = gel(x,1), A = gel(x,2);
2155 : long i, l;
2156 77 : if (typ(A) != t_POL || varn(A) != varn(T))
2157 : {
2158 0 : checkvalidpol(T);
2159 0 : return const_col(degpol(T), A);
2160 : }
2161 77 : v = cleanroots(T,prec); l = lg(v);
2162 77 : for (i=1; i<l; i++) gel(v,i) = poleval(A,gel(v,i));
2163 77 : return v;
2164 : }
2165 :
2166 : GEN
2167 11733 : QX_complex_roots(GEN p, long l)
2168 : {
2169 11733 : pari_sp av = avma;
2170 : long bit;
2171 : GEN L;
2172 :
2173 11733 : if (!signe(p)) pari_err_ROOTS0("QX_complex_roots");
2174 11733 : if (lg(p) == 3) return cgetg(1,t_COL); /* constant polynomial */
2175 11733 : if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
2176 11733 : bit = prec2nbits(l);
2177 11733 : L = all_roots(Q_primpart(p), bit);
2178 11733 : return gerepileupto(av, clean_roots(L, l, bit, 1));
2179 : }
2180 :
2181 : /********************************************************************/
2182 : /** **/
2183 : /** REAL ROOTS OF INTEGER POLYNOMIAL **/
2184 : /** **/
2185 : /********************************************************************/
2186 :
2187 : /* Count sign changes in the coefficients of (x+1)^deg(P)*P(1/(x+1))
2188 : * The inversion is implicit (we take coefficients backwards). Roots of P
2189 : * at 0 and 1 (mapped to oo and 0) are ignored here and must be dealt with
2190 : * by the caller */
2191 : static long
2192 62864 : X2XP1(GEN P, long deg, int *root1, GEN *Premapped)
2193 : {
2194 62864 : const pari_sp av = avma;
2195 62864 : GEN v = shallowcopy(P);
2196 : long i, j, vlim, nb, s;
2197 :
2198 62864 : for (i = 0, vlim = deg+2;;)
2199 : {
2200 62934 : for (j = 2; j < vlim; j++) gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2201 62934 : s = -signe(gel(v, vlim));
2202 62934 : vlim--; i++; if (s) break;
2203 70 : }
2204 62864 : if (vlim == deg+1) *root1 = 0;
2205 : else
2206 : {
2207 70 : *root1 = 1;
2208 70 : if (Premapped) setlg(v, vlim + 2);
2209 : }
2210 :
2211 62864 : nb = 0;
2212 297419 : for (; i < deg; i++)
2213 : {
2214 281001 : long s2 = -signe(gel(v, 2));
2215 281001 : int flag = (s2 == s);
2216 3662344 : for (j = 2; j < vlim; j++)
2217 : {
2218 3381343 : gel(v, j+1) = addii(gel(v, j), gel(v, j+1));
2219 3381343 : if (flag) flag = (s2 != signe(gel(v, j+1)));
2220 : }
2221 281001 : if (s == signe(gel(v, vlim)))
2222 : {
2223 60280 : if (++nb >= 2) { avma = av; return 2; }
2224 39290 : s = -s;
2225 : }
2226 : /* if flag is set there will be no further sign changes */
2227 260011 : if (flag && (!Premapped || !nb)) goto END;
2228 234555 : vlim--;
2229 234555 : if (gc_needed(av, 3))
2230 : {
2231 0 : if (!Premapped) setlg(v, vlim + 2);
2232 0 : v = gerepileupto(av, v);
2233 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "X2XP1");
2234 : }
2235 : }
2236 16418 : if (s == signe(gel(v, vlim))) nb++;
2237 : END:
2238 41874 : if (Premapped && nb == 1) *Premapped = v; else avma = av;
2239 41874 : return nb;
2240 : }
2241 :
2242 : static long
2243 0 : _intervalcmp(GEN x, GEN y)
2244 : {
2245 0 : if (typ(x) == t_VEC) x = gel(x, 1);
2246 0 : if (typ(y) == t_VEC) y = gel(y, 1);
2247 0 : return gcmp(x, y);
2248 : }
2249 :
2250 : static GEN
2251 632087 : _gen_nored(void *E, GEN x) { (void)E; return x; }
2252 : static GEN
2253 3092137 : _mp_add(void *E, GEN x, GEN y) { (void)E; return mpadd(x, y); }
2254 : static GEN
2255 0 : _mp_sub(void *E, GEN x, GEN y) { (void)E; return mpsub(x, y); }
2256 : static GEN
2257 653339 : _mp_mul(void *E, GEN x, GEN y) { (void)E; return mpmul(x, y); }
2258 : static GEN
2259 646892 : _mp_sqr(void *E, GEN x) { (void)E; return mpsqr(x); }
2260 : static GEN
2261 736660 : _gen_one(void *E) { (void)E; return gen_1; }
2262 : static GEN
2263 9965 : _gen_zero(void *E) { (void)E; return gen_0; }
2264 :
2265 : static struct bb_algebra mp_algebra = { _gen_nored, _mp_add, _mp_sub,
2266 : _mp_mul, _mp_sqr, _gen_one, _gen_zero };
2267 :
2268 : static GEN
2269 3575282 : _mp_cmul(void *E, GEN P, long a, GEN x) {(void)E; return mpmul(gel(P,a+2), x);}
2270 :
2271 : /* Split the polynom P in two parts, whose coeffs have constant sign:
2272 : * P(X) = X^D*Pp + Pm. Also compute the two parts of the derivative of P,
2273 : * Pprimem = Pm', Pprimep = X*Pp'+ D*Pp => P' = X^(D-1)*Pprimep + Pprimem;
2274 : * Pprimep[i] = (i+D) Pp[i]. Return D */
2275 : static long
2276 7411 : split_pols(GEN P, GEN *pPp, GEN *pPm, GEN *pPprimep, GEN *pPprimem)
2277 : {
2278 7411 : long i, D, dP = degpol(P), s0 = signe(gel(P,2));
2279 : GEN Pp, Pm, Pprimep, Pprimem;
2280 49645 : for(i=1; i <= dP; i++)
2281 49645 : if (signe(gel(P, i+2)) == -s0) break;
2282 7411 : D = i;
2283 7411 : Pm = cgetg(D + 2, t_POL);
2284 7411 : Pprimem = cgetg(D + 1, t_POL);
2285 7411 : Pp = cgetg(dP-D + 3, t_POL);
2286 7411 : Pprimep = cgetg(dP-D + 3, t_POL);
2287 7411 : Pm[1] = Pp[1] = Pprimem[1] = Pprimep[1] = P[1];
2288 57056 : for(i=0; i < D; i++)
2289 : {
2290 49645 : GEN c = gel(P, i+2);
2291 49645 : gel(Pm, i+2) = c;
2292 49645 : if (i) gel(Pprimem, i+1) = mului(i, c);
2293 : }
2294 55866 : for(; i <= dP; i++)
2295 : {
2296 48455 : GEN c = gel(P, i+2);
2297 48455 : gel(Pp, i+2-D) = c;
2298 48455 : gel(Pprimep, i+2-D) = mului(i, c);
2299 : }
2300 7411 : *pPm = normalizepol_lg(Pm, D+2);
2301 7411 : *pPprimem = normalizepol_lg(Pprimem, D+1);
2302 7411 : *pPp = normalizepol_lg(Pp, dP-D+3);
2303 7411 : *pPprimep = normalizepol_lg(Pprimep, dP-D+3);
2304 7411 : return dP - degpol(*pPp);
2305 : }
2306 :
2307 : static GEN
2308 246555 : bkeval_single_power(long d, GEN V)
2309 : {
2310 246555 : long mp = lg(V) - 2;
2311 246555 : if (d > mp) return gmul(gpowgs(gel(V, mp+1), d/mp), gel(V, (d%mp)+1));
2312 246555 : return gel(V, d+1);
2313 : }
2314 :
2315 : static GEN
2316 246555 : splitpoleval(GEN Pp, GEN Pm, GEN pows, long D, long bitprec)
2317 : {
2318 246555 : GEN vp = gen_bkeval_powers(Pp, degpol(Pp), pows, NULL, &mp_algebra, _mp_cmul);
2319 246555 : GEN vm = gen_bkeval_powers(Pm, degpol(Pm), pows, NULL, &mp_algebra, _mp_cmul);
2320 246555 : GEN xa = bkeval_single_power(D, pows);
2321 : GEN r;
2322 246555 : if (!signe(vp)) return vm;
2323 246555 : vp = gmul(vp, xa);
2324 246555 : r = gadd(vp, vm);
2325 246555 : if (gexpo(vp) - (signe(r)? gexpo(r): 0) > prec2nbits(realprec(vp)) - bitprec)
2326 17500 : return NULL;
2327 229055 : return r;
2328 : }
2329 :
2330 : /* optimized Cauchy bound for P = X^D*Pp + Pm, D > deg(Pm) */
2331 : static GEN
2332 7411 : splitcauchy(GEN Pp, GEN Pm, long prec)
2333 : {
2334 7411 : GEN S = gel(Pp,2), A = gel(Pm,2);
2335 7411 : long i, lPm = lg(Pm), lPp = lg(Pp);
2336 7411 : for (i=3; i < lPm; i++) { GEN c = gel(Pm,i); if (abscmpii(A, c) < 0) A = c; }
2337 7411 : for (i=3; i < lPp; i++) S = addii(S, gel(Pp, i));
2338 7411 : return subsr(1, rdivii(A, S, prec)); /* 1 + |Pm|_oo / |Pp|_1 */
2339 : }
2340 :
2341 : /* Newton for polynom P, P(0)!=0, with unique sign change => one root in ]0,oo[
2342 : * P' has also at most one zero there */
2343 : static GEN
2344 7411 : polsolve(GEN P, long bitprec)
2345 : {
2346 7411 : pari_sp av = avma;
2347 : GEN Pp, Pm, Pprimep, Pprimem, pows;
2348 : GEN Pprime, Pprime2;
2349 7411 : GEN ra, rb, rc, rcold = NULL, Pc, Ppc;
2350 7411 : long deg = degpol(P), D, rt;
2351 7411 : long s0, cprec = DEFAULTPREC, prec = nbits2prec(bitprec);
2352 : long bitaddprec, addprec;
2353 7411 : long expoold = LONG_MAX, iter;
2354 7411 : if (deg == 1)
2355 0 : return gerepileuptoleaf(av, rdivii(negi(gel(P,2)), gel(P,3), prec));
2356 7411 : Pprime = ZX_deriv(P);
2357 7411 : Pprime2 = ZX_deriv(Pprime);
2358 7411 : bitaddprec = 1 + 2*expu(deg); addprec = nbits2prec(bitaddprec);
2359 7411 : D = split_pols(P, &Pp, &Pm, &Pprimep, &Pprimem); /* P = X^D*Pp + Pm */
2360 7411 : s0 = signe(gel(P, 2));
2361 7411 : rt = maxss(D, brent_kung_optpow(maxss(degpol(Pp), degpol(Pm)), 2, 1));
2362 7411 : rb = splitcauchy(Pp, Pm, DEFAULTPREC);
2363 : for(;;)
2364 : {
2365 7411 : pows = gen_powers(rb, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2366 7411 : Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec);
2367 7411 : if (!Pc) { cprec++; rb = rtor(rb, cprec); continue; }
2368 7411 : if (signe(Pc) != s0) break;
2369 0 : shiftr_inplace(rb,1);
2370 0 : }
2371 7411 : ra = NULL;
2372 7411 : iter = 0;
2373 : for(;;)
2374 : {
2375 : GEN wdth;
2376 69536 : iter++;
2377 69536 : if (ra)
2378 46653 : rc = shiftr(addrr(ra, rb), -1);
2379 : else
2380 22883 : rc = shiftr(rb, -1);
2381 : do
2382 : {
2383 69536 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2384 69536 : Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+2);
2385 69536 : if (Pc) break;
2386 0 : cprec++; rc = rtor(rc, cprec);
2387 0 : } while (1);
2388 69536 : if (signe(Pc) == s0)
2389 30147 : ra = rc;
2390 : else
2391 39389 : rb = rc;
2392 69536 : if (!ra) continue;
2393 54064 : wdth = subrr(rb, ra);
2394 54064 : if (!(iter % 8))
2395 : {
2396 7992 : GEN m1 = poleval(Pprime, ra), M2;
2397 7992 : if (signe(m1) == s0) continue;
2398 7775 : M2 = poleval(Pprime2, rb);
2399 7775 : if (abscmprr(gmul(M2, wdth), shiftr(m1, 1)) > 0) continue;
2400 7411 : break;
2401 : }
2402 46072 : else if (gexpo(wdth) <= -bitprec)
2403 0 : break;
2404 62125 : }
2405 7411 : rc = rb;
2406 7411 : iter = 0;
2407 : for(;;)
2408 : {
2409 : long exponew;
2410 : GEN dist;
2411 84804 : iter++;
2412 84804 : rcold = rc;
2413 84804 : pows = gen_powers(rc, rt, 1, NULL, _mp_sqr, _mp_mul, _gen_one);
2414 84804 : Ppc = splitpoleval(Pprimep, Pprimem, pows, D-1, bitaddprec+4);
2415 84804 : if (Ppc)
2416 84804 : Pc = splitpoleval(Pp, Pm, pows, D, bitaddprec+4);
2417 84804 : if (!Ppc || !Pc)
2418 : {
2419 17500 : if (cprec >= prec+addprec)
2420 1878 : cprec++;
2421 : else
2422 15622 : cprec = minss(2*cprec, prec+addprec);
2423 17500 : rc = rtor(rc, cprec); /* Backtrack one step */
2424 17500 : continue;
2425 : }
2426 67304 : dist = typ(Ppc) == t_REAL? divrr(Pc, Ppc): divri(Pc, Ppc);
2427 67304 : rc = subrr(rc, dist);
2428 67304 : if (cmprr(ra, rc) > 0 || cmprr(rb, rc) < 0)
2429 : {
2430 0 : if (cprec >= prec+addprec) break;
2431 0 : cprec = minss(2*cprec, prec+addprec);
2432 0 : rc = rtor(rcold, cprec); /* Backtrack one step */
2433 0 : continue;
2434 : }
2435 67304 : if (expoold == LONG_MAX) { expoold = expo(dist); continue; }
2436 57711 : exponew = expo(dist);
2437 57711 : if (exponew < -bitprec - 1)
2438 : {
2439 8828 : if (cprec >= prec+addprec) break;
2440 1417 : cprec = minss(2*cprec, prec+addprec);
2441 1417 : rc = rtor(rc, cprec); continue;
2442 : }
2443 48883 : if (exponew > expoold - 2)
2444 : {
2445 2182 : if (cprec >= prec+addprec) break;
2446 2182 : cprec = minss(2*cprec, prec+addprec);
2447 2182 : rc = rtor(rc, cprec);
2448 2182 : expoold = LONG_MAX; continue;
2449 : }
2450 46701 : expoold = exponew;
2451 77393 : }
2452 7411 : return gerepileuptoleaf(av, rtor(rc, prec));
2453 : }
2454 :
2455 : static GEN
2456 8316 : usp(GEN Q0, long deg, long flag, long bitprec)
2457 : {
2458 8316 : const pari_sp av = avma;
2459 : GEN Q, sol, c, Lc, Lk;
2460 8316 : long listsize = 64, nbr = 0, nb_todo, ind, deg0, indf, i, k, nb;
2461 :
2462 :
2463 8316 : sol = const_col(deg, gen_0);
2464 8316 : deg0 = deg;
2465 8316 : Lc = const_vec(listsize, gen_0);
2466 8316 : Lk = cgetg(listsize+1, t_VECSMALL);
2467 8316 : c = gen_0;
2468 8316 : k = Lk[1] = 0;
2469 8316 : ind = 1; indf = 2;
2470 8316 : Q = leafcopy(Q0);
2471 :
2472 8316 : nb_todo = 1;
2473 79496 : while (nb_todo)
2474 : {
2475 62864 : GEN nc = gel(Lc, ind), Qremapped;
2476 : pari_sp av2;
2477 : int root1;
2478 62864 : if (Lk[ind] == k + 1)
2479 : {
2480 24313 : deg0 = deg;
2481 24313 : setlg(Q0, deg + 3);
2482 24313 : Q0 = ZX_rescale2n(Q0, 1);
2483 24313 : Q = Q_primpart(Q0);
2484 24313 : c = gen_0;
2485 : }
2486 :
2487 62864 : if (!equalii(nc, c)) Q = ZX_translate(Q, subii(nc, c));
2488 :
2489 62864 : k = Lk[ind];
2490 62864 : c = nc;
2491 62864 : ind++;
2492 62864 : nb_todo--;
2493 :
2494 62864 : if (equalii(gel(Q, 2), gen_0))
2495 : { /* Q(0) = 0 */
2496 154 : GEN s = gmul2n(c, -k);
2497 : long j;
2498 196 : for (j = 1; j <= nbr; j++)
2499 126 : if (gequal(gel(sol, j), s)) break;
2500 154 : if (j > nbr) gel(sol, ++nbr) = s;
2501 :
2502 154 : deg0--;
2503 154 : for (j = 2; j <= deg0 + 2; j++) gel(Q, j) = gel(Q, j+1);
2504 154 : setlg(Q, j);
2505 : }
2506 :
2507 62864 : av2 = avma;
2508 62864 : nb = X2XP1(Q, deg0, &root1, flag == 1 ? &Qremapped : NULL);
2509 :
2510 62864 : if (nb == 0) /* no root in this open interval */;
2511 41405 : else if (nb == 1) /* exactly one root */
2512 : {
2513 14131 : GEN s = gen_0;
2514 14131 : if (flag == 0)
2515 0 : s = mkvec2(gmul2n(c,-k), gmul2n(addiu(c,1),-k));
2516 14131 : else if (flag == 1) /* Caveat: Qremapped is the reciprocal polynomial */
2517 : {
2518 7411 : s = polsolve(Qremapped, bitprec+1);
2519 7411 : s = divrr(s, addsr(1, s));
2520 7411 : s = gmul2n(addir(c, s), -k);
2521 7411 : s = rtor(s, nbits2prec(bitprec));
2522 : }
2523 14131 : gel(sol, ++nbr) = gerepileupto(av2, s);
2524 : }
2525 : else
2526 : { /* unknown, add two nodes to refine */
2527 27274 : if (indf + 2 > listsize)
2528 : {
2529 154 : if (ind>1)
2530 : {
2531 805 : for (i = ind; i < indf; i++)
2532 : {
2533 651 : gel(Lc, i-ind+1) = gel(Lc, i);
2534 651 : Lk[i-ind+1] = Lk[i];
2535 : }
2536 154 : indf -= ind-1;
2537 154 : ind = 1;
2538 : }
2539 154 : if (indf + 2 > listsize)
2540 : {
2541 0 : listsize *= 2;
2542 0 : Lc = vec_lengthen(Lc, listsize);
2543 0 : Lk = vecsmall_lengthen(Lk, listsize);
2544 : }
2545 154 : for (i = indf; i <= listsize; i++) gel(Lc, i) = gen_0;
2546 : }
2547 27274 : nc = shifti(c, 1);
2548 27274 : gel(Lc, indf) = nc;
2549 27274 : gel(Lc, indf + 1) = addiu(nc, 1);
2550 27274 : Lk[indf] = Lk[indf + 1] = k + 1;
2551 27274 : indf += 2;
2552 27274 : nb_todo += 2;
2553 : }
2554 62864 : if (root1)
2555 : { /* Q(1) = 0 */
2556 70 : GEN s = gmul2n(addiu(c,1), -k);
2557 : long j;
2558 112 : for (j = 1; j <= nbr; j++)
2559 49 : if (gequal(gel(sol, j), s)) break;
2560 70 : if (j > nbr) gel(sol, ++nbr) = s;
2561 : }
2562 :
2563 62864 : if (gc_needed(av, 2))
2564 : {
2565 0 : gerepileall(av, 6, &Q0, &Q, &c, &Lc, &Lk, &sol);
2566 0 : if (DEBUGMEM > 1) pari_warn(warnmem, "ZX_Uspensky", avma);
2567 : }
2568 : }
2569 :
2570 8316 : setlg(sol, nbr+1);
2571 8316 : return gerepilecopy(av, sol);
2572 : }
2573 :
2574 : static GEN
2575 2226 : ZX_Uspensky_cst_pol(long nbz, long flag, long bitprec)
2576 : {
2577 2226 : switch(flag)
2578 : {
2579 0 : case 0: return zerocol(nbz);
2580 0 : case 1: retconst_col(nbz, real_0_bit(bitprec));
2581 2226 : default: return utoi(nbz);
2582 : }
2583 : }
2584 :
2585 : /* ZX_Uspensky(P, [a,a], flag) */
2586 : static GEN
2587 28 : ZX_Uspensky_equal(GEN P, GEN a, long flag)
2588 : {
2589 28 : if (typ(a) != t_INFINITY && gequal0(poleval(P, a)))
2590 21 : return flag <= 1 ? mkcol(a): gen_1;
2591 : else
2592 7 : return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2593 : }
2594 :
2595 : GEN
2596 9900 : ZX_Uspensky(GEN P, GEN ab, long flag, long bitprec)
2597 : {
2598 9900 : pari_sp av = avma;
2599 9900 : GEN a, b, sol = NULL, Pcur;
2600 : double fb;
2601 : long nbz, deg;
2602 :
2603 9900 : deg = degpol(P);
2604 9900 : if (deg == 0) return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2605 9900 : if (ab)
2606 : {
2607 3887 : if (typ(ab) == t_VEC)
2608 : {
2609 175 : if (lg(ab) != 3) pari_err_DIM("ZX_Uspensky");
2610 175 : a = gel(ab, 1);
2611 175 : b = gel(ab, 2);
2612 : }
2613 : else
2614 : {
2615 3712 : a = ab;
2616 3712 : b = mkoo();
2617 : }
2618 : }
2619 : else
2620 : {
2621 6013 : a = mkmoo();
2622 6013 : b = mkoo();
2623 : }
2624 9900 : switch (gcmp(a, b))
2625 : {
2626 7 : case 1: avma = av; return flag <= 1 ? cgetg(1, t_COL) : gen_0;
2627 21 : case 0: return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag));
2628 : }
2629 9872 : nbz = ZX_valrem(P, &Pcur);
2630 9872 : deg -= nbz;
2631 9872 : if (!nbz) Pcur = P;
2632 9872 : if (nbz && (gsigne(a) > 0 || gsigne(b) < 0)) nbz = 0;
2633 9872 : if (deg == 0) { avma = av; return ZX_Uspensky_cst_pol(nbz, flag, bitprec); }
2634 8759 : if (deg == 1)
2635 : {
2636 3222 : sol = gdiv(gneg(gel(Pcur, 2)), pollead(Pcur, -1));
2637 3222 : if (gcmp(a, sol) > 0 || gcmp(sol, b) > 0)
2638 : {
2639 1113 : avma = av;
2640 1113 : return ZX_Uspensky_cst_pol(nbz, flag, bitprec);
2641 : }
2642 2109 : if (flag >= 2) { avma = av; return utoi(nbz+1); }
2643 954 : sol = gconcat(zerocol(nbz), mkcol(sol));
2644 954 : if (flag == 1) sol = RgC_gtofp(sol, nbits2prec(bitprec));
2645 954 : return gerepilecopy(av, sol);
2646 : }
2647 5537 : switch(flag)
2648 : {
2649 : case 0:
2650 0 : sol = zerocol(nbz);
2651 0 : break;
2652 : case 1:
2653 1575 : sol = const_col(nbz, real_0_bit(bitprec));
2654 1575 : break;
2655 : /* case 2: nothing */
2656 : }
2657 :
2658 5537 : if (typ(a) == t_INFINITY && typ(b) != t_INFINITY && gsigne(b))
2659 : {
2660 35 : fb = fujiwara_bound_real(Pcur, -1);
2661 35 : if (fb > -pariINFINITY) a = negi(int2n((long)ceil(fb))); else a = gen_0;
2662 : }
2663 5537 : if (typ(b) == t_INFINITY && typ(a) != t_INFINITY && gsigne(a))
2664 : {
2665 21 : fb = fujiwara_bound_real(Pcur, 1);
2666 21 : if (fb > -pariINFINITY) b = int2n((long)ceil(fb)); else b = gen_0;
2667 : }
2668 :
2669 5537 : if (typ(a) != t_INFINITY && typ(b) != t_INFINITY)
2670 : {
2671 : GEN den, diff, unscaledres, co, Pdiv, ascaled;
2672 : pari_sp av1;
2673 : long i;
2674 140 : if (gequal(a,b)) /* can occur if one of a,b was initially a t_INFINITY */
2675 7 : return gerepilecopy(av, ZX_Uspensky_equal(P, a, flag));
2676 133 : den = lcmii(Q_denom(a), Q_denom(b));
2677 133 : if (!is_pm1(den))
2678 : {
2679 28 : Pcur = ZX_rescale(Pcur, den);
2680 28 : ascaled = gmul(a, den);
2681 : }
2682 : else
2683 : {
2684 105 : den = NULL;
2685 105 : ascaled = a;
2686 : }
2687 133 : diff = subii(den ? gmul(b,den) : b, ascaled);
2688 133 : Pcur = ZX_unscale(ZX_translate(Pcur, ascaled), diff);
2689 133 : av1 = avma;
2690 133 : Pdiv = cgetg(deg+2, t_POL);
2691 133 : Pdiv[1] = Pcur[1];
2692 133 : co = gel(Pcur, deg+2);
2693 665 : for (i = deg; --i >= 0; )
2694 : {
2695 399 : gel(Pdiv, i+2) = co;
2696 399 : co = addii(co, gel(Pcur, i+2));
2697 : }
2698 133 : if (!signe(co))
2699 : {
2700 56 : Pcur = Pdiv;
2701 56 : deg--;
2702 56 : if (flag <= 1)
2703 14 : sol = gconcat(sol, b);
2704 : else
2705 42 : nbz++;
2706 : }
2707 : else
2708 77 : avma = av1;
2709 133 : unscaledres = usp(Pcur, deg, flag, bitprec);
2710 133 : if (flag <= 1)
2711 : {
2712 63 : for (i = 1; i < lg(unscaledres); i++)
2713 : {
2714 35 : GEN z = gmul(diff, gel(unscaledres, i));
2715 35 : if (typ(z) == t_VEC)
2716 : {
2717 0 : gel(z, 1) = gadd(ascaled, gel(z, 1));
2718 0 : gel(z, 2) = gadd(ascaled, gel(z, 2));
2719 : }
2720 : else
2721 35 : z = gadd(ascaled, z);
2722 35 : if (den) z = gdiv(z, den);
2723 35 : gel(unscaledres, i) = z;
2724 : }
2725 28 : sol = gconcat(sol, unscaledres);
2726 : }
2727 : else
2728 105 : nbz += lg(unscaledres) - 1;
2729 : }
2730 5530 : if (typ(b) == t_INFINITY && (fb=fujiwara_bound_real(Pcur, 1)) > -pariINFINITY)
2731 : {
2732 : GEN Pcurp, unscaledres;
2733 4452 : long bp = (long)ceil(fb);
2734 4452 : if (bp < 0) bp = 0;
2735 4452 : Pcurp = ZX_unscale2n(Pcur, bp);
2736 4452 : unscaledres = usp(Pcurp, deg, flag, bitprec);
2737 4452 : if (flag <= 1)
2738 1540 : sol = shallowconcat(sol, gmul2n(unscaledres, bp));
2739 : else
2740 2912 : nbz += lg(unscaledres)-1;
2741 : }
2742 5530 : if (typ(a) == t_INFINITY && (fb=fujiwara_bound_real(Pcur,-1)) > -pariINFINITY)
2743 : {
2744 : GEN Pcurm, unscaledres;
2745 3731 : long i, bm = (long)ceil(fb);
2746 3731 : if (bm < 0) bm = 0;
2747 3731 : Pcurm = ZX_unscale2n(ZX_unscale(Pcur, gen_m1), bm);
2748 3731 : unscaledres = usp(Pcurm, deg, flag, bitprec);
2749 3731 : if (flag <= 1)
2750 : {
2751 3688 : for (i = 1; i < lg(unscaledres); i++)
2752 : {
2753 2470 : GEN z = gneg(gmul2n(gel(unscaledres, i), bm));
2754 2470 : if (typ(z) == t_VEC) swap(gel(z, 1), gel(z, 2));
2755 2470 : gel(unscaledres, i) = z;
2756 : }
2757 1218 : sol = shallowconcat(unscaledres, sol);
2758 : }
2759 : else
2760 2513 : nbz += lg(unscaledres)-1;
2761 : }
2762 5530 : if (flag >= 2) return utoi(nbz);
2763 1575 : if (flag)
2764 1575 : sol = sort(sol);
2765 : else
2766 0 : sol = gen_sort(sol, (void *)_intervalcmp, cmp_nodata);
2767 1575 : return gerepileupto(av, sol);
2768 : }
2769 :
2770 : /* x a scalar */
2771 : static GEN
2772 35 : rootsdeg0(GEN x)
2773 : {
2774 35 : if (!is_rational_t(typ(x))) pari_err_TYPE("realroots",x);
2775 28 : if (gequal0(x)) pari_err_ROOTS0("realroots");
2776 14 : return cgetg(1,t_COL); /* constant polynomial */
2777 : }
2778 : static void
2779 1708 : checkbound(GEN a)
2780 : {
2781 1708 : switch(typ(a))
2782 : {
2783 1708 : case t_INT: case t_FRAC: case t_INFINITY: break;
2784 0 : default: pari_err_TYPE("polrealroots", a);
2785 : }
2786 1708 : }
2787 : static GEN
2788 4328 : check_ab(GEN ab)
2789 : {
2790 : GEN a, b;
2791 4328 : if (!ab) return NULL;
2792 854 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("polrootsreal",ab);
2793 854 : a = gel(ab,1); checkbound(a);
2794 854 : b = gel(ab,2); checkbound(b);
2795 1547 : if (typ(a) == t_INFINITY && inf_get_sign(a) < 0 &&
2796 1351 : typ(b) == t_INFINITY && inf_get_sign(b) > 0) ab = NULL;
2797 854 : return ab;
2798 : }
2799 : GEN
2800 3530 : realroots(GEN P, GEN ab, long prec)
2801 : {
2802 3530 : pari_sp av = avma;
2803 3530 : long nrr = 0;
2804 3530 : GEN sol = NULL, fa, ex;
2805 : long i, j, v;
2806 :
2807 3530 : ab = check_ab(ab);
2808 3530 : if (typ(P) != t_POL) return rootsdeg0(P);
2809 3509 : switch(degpol(P))
2810 : {
2811 7 : case -1: return rootsdeg0(gen_0);
2812 7 : case 0: return rootsdeg0(gel(P,2));
2813 : }
2814 3495 : P = Q_primpart(P);
2815 3495 : RgX_check_ZX(P,"realroots");
2816 3495 : v = ZX_valrem(P,&P);
2817 3495 : if (v && (!ab || (gsigne(gel(ab,1)) <= 0 && gsigne(gel(ab,2)) >= 0)))
2818 1008 : sol = const_col(v, real_0(prec));
2819 3495 : fa = ZX_squff(P, &ex);
2820 6003 : for (i = 1; i < lg(fa); i++)
2821 : {
2822 2508 : GEN Pi = gel(fa, i), soli, soli2 = NULL;
2823 2508 : long n, nrri = 0, h;
2824 2508 : if (ab)
2825 42 : h = 1;
2826 : else
2827 2466 : Pi = ZX_deflate_max(Pi, &h);
2828 2508 : soli = ZX_Uspensky(Pi, h%2 ? ab: gen_0, 1, prec2nbits(prec));
2829 2508 : n = lg(soli);
2830 2508 : if (!(h % 2)) soli2 = cgetg(n, t_COL);
2831 10005 : for (j = 1; j < n; j++)
2832 : {
2833 7497 : GEN elt = gel(soli, j); /* != 0 */
2834 7497 : if (typ(elt) != t_REAL)
2835 : {
2836 70 : nrri++; if (h > 1 && !(h % 2)) nrri++;
2837 70 : elt = gtofp(elt, prec);
2838 70 : gel(soli, j) = elt;
2839 : }
2840 7497 : if (h > 1)
2841 : {
2842 : GEN r;
2843 1472 : if (h == 2)
2844 1458 : r = sqrtr(elt);
2845 : else
2846 : {
2847 14 : if (signe(elt) < 0)
2848 7 : r = negr(sqrtnr(negr(elt), h));
2849 : else
2850 7 : r = sqrtnr(elt, h);
2851 : }
2852 1472 : gel(soli, j) = r;
2853 1472 : if (!(h % 2)) gel(soli2, j) = negr(r);
2854 : }
2855 : }
2856 2508 : if (!(h % 2)) soli = shallowconcat(soli, soli2);
2857 2508 : if (ex[i] > 1) soli = shallowconcat1( const_vec(ex[i], soli) );
2858 2508 : sol = sol? shallowconcat(sol, soli): soli;
2859 2508 : nrr += ex[i]*nrri;
2860 : }
2861 3495 : if (!sol) { avma = av; return cgetg(1,t_COL); }
2862 :
2863 3481 : if (DEBUGLEVEL > 4)
2864 : {
2865 0 : err_printf("Number of real roots: %d\n", lg(sol)-1);
2866 0 : err_printf(" -- of which 2-integral: %ld\n", nrr);
2867 : }
2868 3481 : return gerepileupto(av, sort(sol));
2869 : }
2870 :
2871 : /* P non-constant, squarefree ZX */
2872 : long
2873 7224 : ZX_sturm(GEN P)
2874 : {
2875 7224 : pari_sp av = avma;
2876 : long h, r;
2877 7224 : P = ZX_deflate_max(P, &h);
2878 7224 : if (odd(h))
2879 4487 : r = itos(ZX_Uspensky(P, NULL, 2, 0));
2880 : else
2881 2737 : r = 2*itos(ZX_Uspensky(P, gen_0, 2, 0));
2882 7224 : avma = av; return r;
2883 : }
2884 : /* P non-constant, squarefree ZX */
2885 : long
2886 798 : ZX_sturmpart(GEN P, GEN ab)
2887 : {
2888 798 : pari_sp av = avma;
2889 : long r;
2890 798 : if (!check_ab(ab)) return ZX_sturm(P);
2891 133 : r = itos(ZX_Uspensky(P, ab, 2, 0));
2892 133 : avma = av; return r;
2893 : }
|
