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; either version 2 of the License, or (at your option) any later
8 version. It is distributed in the hope that it will be useful, but WITHOUT
9 ANY WARRANTY WHATSOEVER.
10
11 Check the License for details. You should have received a copy of it, along
12 with the package; see the file 'COPYING'. If not, write to the Free Software
13 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14
15 #include "pari.h"
16 #include "paripriv.h"
17
18 /*******************************************************************/
19 /** **/
20 /** SPECIAL POLYNOMIALS **/
21 /** **/
22 /*******************************************************************/
23 /* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
24 * T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
25 * where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
26 * and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
27 GEN
polchebyshev1(long n,long v)28 polchebyshev1(long n, long v) /* Assume 4*n < LONG_MAX */
29 {
30 long k, l;
31 pari_sp av;
32 GEN q,a,r;
33
34 if (v<0) v = 0;
35 /* polchebyshev(-n,1) = polchebyshev(n,1) */
36 if (n < 0) n = -n;
37 if (n==0) return pol_1(v);
38 if (n==1) return pol_x(v);
39
40 q = cgetg(n+3, t_POL); r = q + n+2;
41 a = int2n(n-1);
42 gel(r--,0) = a;
43 gel(r--,0) = gen_0;
44 for (k=1,l=n; l>1; k++,l-=2)
45 {
46 av = avma;
47 a = diviuuexact(muluui(l, l-1, a), 4*k, n-k);
48 togglesign(a); a = gerepileuptoint(av, a);
49 gel(r--,0) = a;
50 gel(r--,0) = gen_0;
51 }
52 q[1] = evalsigne(1) | evalvarn(v);
53 return q;
54 }
55 static void
polchebyshev1_eval_aux(long n,GEN x,GEN * pt1,GEN * pt2)56 polchebyshev1_eval_aux(long n, GEN x, GEN *pt1, GEN *pt2)
57 {
58 GEN t1, t2, b;
59 if (n == 1) { *pt1 = gen_1; *pt2 = x; return; }
60 if (n == 0) { *pt1 = x; *pt2 = gen_1; return; }
61 polchebyshev1_eval_aux((n+1) >> 1, x, &t1, &t2);
62 b = gsub(gmul(gmul2n(t1,1), t2), x);
63 if (odd(n)) { *pt1 = gadd(gmul2n(gsqr(t1), 1), gen_m1); *pt2 = b; }
64 else { *pt1 = b; *pt2 = gadd(gmul2n(gsqr(t2), 1), gen_m1); }
65 }
66 static GEN
polchebyshev1_eval(long n,GEN x)67 polchebyshev1_eval(long n, GEN x)
68 {
69 GEN t1, t2;
70 long i, v;
71 pari_sp av;
72
73 if (n < 0) n = -n;
74 if (n==0) return gen_1;
75 if (n==1) return gcopy(x);
76 av = avma;
77 v = u_lvalrem(n, 2, (ulong*)&n);
78 polchebyshev1_eval_aux((n+1)>>1, x, &t1, &t2);
79 if (n != 1) t2 = gsub(gmul(gmul2n(t1,1), t2), x);
80 for (i = 1; i <= v; i++) t2 = gadd(gmul2n(gsqr(t2), 1), gen_m1);
81 return gerepileupto(av, t2);
82 }
83
84 /* Chebychev polynomial of the second kind U(n,x): the coefficient in front of
85 * x^(n-2*m) is (-1)^m * 2^(n-2m)*(n-m)!/m!/(n-2m)! for m=0,1,...,n/2 */
86 GEN
polchebyshev2(long n,long v)87 polchebyshev2(long n, long v)
88 {
89 pari_sp av;
90 GEN q, a, r;
91 long m;
92 int neg = 0;
93
94 if (v<0) v = 0;
95 /* polchebyshev(-n,2) = -polchebyshev(n-2,2) */
96 if (n < 0) {
97 if (n == -1) return zeropol(v);
98 neg = 1; n = -n-2;
99 }
100 if (n==0) return neg ? scalar_ZX_shallow(gen_m1, v): pol_1(v);
101
102 q = cgetg(n+3, t_POL); r = q + n+2;
103 a = int2n(n);
104 if (neg) togglesign(a);
105 gel(r--,0) = a;
106 gel(r--,0) = gen_0;
107 for (m=1; 2*m<= n; m++)
108 {
109 av = avma;
110 a = diviuuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m, n-m+1);
111 togglesign(a); a = gerepileuptoint(av, a);
112 gel(r--,0) = a;
113 gel(r--,0) = gen_0;
114 }
115 q[1] = evalsigne(1) | evalvarn(v);
116 return q;
117 }
118 static void
polchebyshev2_eval_aux(long n,GEN x,GEN * pu1,GEN * pu2)119 polchebyshev2_eval_aux(long n, GEN x, GEN *pu1, GEN *pu2)
120 {
121 GEN u1, u2, u, mu1;
122 if (n == 1) { *pu1 = gen_1; *pu2 = gmul2n(x,1); return; }
123 if (n == 0) { *pu1 = gen_0; *pu2 = gen_1; return; }
124 polchebyshev2_eval_aux(n >> 1, x, &u1, &u2);
125 mu1 = gneg(u1);
126 u = gmul(gadd(u2,u1), gadd(u2,mu1));
127 if (odd(n)) { *pu1 = u; *pu2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1)); }
128 else { *pu2 = u; *pu1 = gmul(gmul2n(u1,1), gadd(u2, gmul(x,mu1))); }
129 }
130 static GEN
polchebyshev2_eval(long n,GEN x)131 polchebyshev2_eval(long n, GEN x)
132 {
133 GEN u1, u2, mu1;
134 long neg = 0;
135 pari_sp av;
136
137 if (n < 0) {
138 if (n == -1) return gen_0;
139 neg = 1; n = -n-2;
140 }
141 if (n==0) return neg ? gen_m1: gen_1;
142 av = avma;
143 polchebyshev2_eval_aux(n>>1, x, &u1, &u2);
144 mu1 = gneg(u1);
145 if (odd(n)) u2 = gmul(gmul2n(u2,1), gadd(gmul(x,u2), mu1));
146 else u2 = gmul(gadd(u2,u1), gadd(u2,mu1));
147 if (neg) u2 = gneg(u2);
148 return gerepileupto(av, u2);
149 }
150
151 GEN
polchebyshev(long n,long kind,long v)152 polchebyshev(long n, long kind, long v)
153 {
154 switch (kind)
155 {
156 case 1: return polchebyshev1(n, v);
157 case 2: return polchebyshev2(n, v);
158 default: pari_err_FLAG("polchebyshev");
159 }
160 return NULL; /* LCOV_EXCL_LINE */
161 }
162 GEN
polchebyshev_eval(long n,long kind,GEN x)163 polchebyshev_eval(long n, long kind, GEN x)
164 {
165 if (!x) return polchebyshev(n, kind, 0);
166 if (gequalX(x)) return polchebyshev(n, kind, varn(x));
167 switch (kind)
168 {
169 case 1: return polchebyshev1_eval(n, x);
170 case 2: return polchebyshev2_eval(n, x);
171 default: pari_err_FLAG("polchebyshev");
172 }
173 return NULL; /* LCOV_EXCL_LINE */
174 }
175
176 /* Hermite polynomial H(n,x): H(n+1) = 2x H(n) - 2n H(n-1)
177 * The coefficient in front of x^(n-2*m) is
178 * (-1)^m * n! * 2^(n-2m)/m!/(n-2m)! for m=0,1,...,n/2.. */
179 GEN
polhermite(long n,long v)180 polhermite(long n, long v)
181 {
182 long m;
183 pari_sp av;
184 GEN q,a,r;
185
186 if (v<0) v = 0;
187 if (n==0) return pol_1(v);
188
189 q = cgetg(n+3, t_POL); r = q + n+2;
190 a = int2n(n);
191 gel(r--,0) = a;
192 gel(r--,0) = gen_0;
193 for (m=1; 2*m<= n; m++)
194 {
195 av = avma;
196 a = diviuexact(muluui(n-2*m+2, n-2*m+1, a), 4*m);
197 togglesign(a);
198 gel(r--,0) = a = gerepileuptoint(av, a);
199 gel(r--,0) = gen_0;
200 }
201 q[1] = evalsigne(1) | evalvarn(v);
202 return q;
203 }
204 static void
err_hermite(long n)205 err_hermite(long n)
206 { pari_err_DOMAIN("polhermite", "degree", "<", gen_0, stoi(n)); }
207 GEN
polhermite_eval0(long n,GEN x,long flag)208 polhermite_eval0(long n, GEN x, long flag)
209 {
210 long i;
211 pari_sp av, av2;
212 GEN x2, u, v;
213
214 if (n < 0) err_hermite(n);
215 if (!x || gequalX(x))
216 {
217 long v = x? varn(x): 0;
218 if (flag)
219 {
220 if (!n) err_hermite(-1);
221 retmkvec2(polhermite(n-1,v),polhermite(n,v));
222 }
223 return polhermite(n, v);
224 }
225 if (n==0)
226 {
227 if (flag) err_hermite(-1);
228 return gen_1;
229 }
230 if (n==1)
231 {
232 if (flag) retmkvec2(gen_1, gmul2n(x,1));
233 return gmul2n(x,1);
234 }
235 av = avma; x2 = gmul2n(x,1); v = gen_1; u = x2;
236 av2= avma;
237 for (i=1; i<n; i++)
238 { /* u = H_i(x), v = H_{i-1}(x), compute t = H_{i+1}(x) */
239 GEN t;
240 if ((i & 0xff) == 0) gerepileall(av2,2,&u, &v);
241 t = gsub(gmul(x2, u), gmulsg(2*i,v));
242 v = u; u = t;
243 }
244 if (flag) return gerepilecopy(av, mkvec2(v, u));
245 return gerepileupto(av, u);
246 }
247 GEN
polhermite_eval(long n,GEN x)248 polhermite_eval(long n, GEN x) { return polhermite_eval0(n, x, 0); }
249
250 /* Legendre polynomial
251 * L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1)
252 * L(n) = 2^-n sum_{k=0}^{n/2} a_k x^(n-2k)
253 * where a_k = (-1)^k (2n-2k)! / k! (n-k)! (n-2k)! is an integer
254 * and a_0 = binom(2n,n), a_k / a_{k-1} = - (n-2k+1)(n-2k+2) / 2k (2n-2k+1) */
255 GEN
pollegendre(long n,long v)256 pollegendre(long n, long v)
257 {
258 long k, l;
259 pari_sp av;
260 GEN a, r, q;
261
262 if (v<0) v = 0;
263 /* pollegendre(-n) = pollegendre(n-1) */
264 if (n < 0) n = -n-1;
265 if (n==0) return pol_1(v);
266 if (n==1) return pol_x(v);
267
268 av = avma;
269 q = cgetg(n+3, t_POL); r = q + n+2;
270 gel(r--,0) = a = binomialuu(n<<1,n);
271 gel(r--,0) = gen_0;
272 for (k=1,l=n; l>1; k++,l-=2)
273 { /* l = n-2*k+2 */
274 av = avma;
275 a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
276 togglesign(a); a = gerepileuptoint(av, a);
277 gel(r--,0) = a;
278 gel(r--,0) = gen_0;
279 }
280 q[1] = evalsigne(1) | evalvarn(v);
281 return gerepileupto(av, gmul2n(q,-n));
282 }
283 /* q such that Ln * 2^n = q(x^2) [n even] or x q(x^2) [n odd] */
284 GEN
pollegendre_reduced(long n,long v)285 pollegendre_reduced(long n, long v)
286 {
287 long k, l, N;
288 pari_sp av;
289 GEN a, r, q;
290
291 if (v<0) v = 0;
292 /* pollegendre(-n) = pollegendre(n-1) */
293 if (n < 0) n = -n-1;
294 if (n<=1) return n? scalarpol_shallow(gen_2,v): pol_1(v);
295
296 N = n >> 1;
297 q = cgetg(N+3, t_POL); r = q + N+2;
298 gel(r--,0) = a = binomialuu(n<<1,n);
299 for (k=1,l=n; l>1; k++,l-=2)
300 { /* l = n-2*k+2 */
301 av = avma;
302 a = diviuuexact(muluui(l, l-1, a), 2*k, n+l-1);
303 togglesign(a);
304 gel(r--,0) = a = gerepileuptoint(av, a);
305 }
306 q[1] = evalsigne(1) | evalvarn(v);
307 return q;
308 }
309
310 GEN
pollegendre_eval0(long n,GEN x,long flag)311 pollegendre_eval0(long n, GEN x, long flag)
312 {
313 pari_sp av;
314 GEN u, v;
315 long i;
316
317 if (n < 0) n = -n-1; /* L(-n) = L(n-1) */
318 /* n >= 0 */
319 if (flag && flag != 1) pari_err_FLAG("pollegendre");
320 if (!x || gequalX(x))
321 {
322 long v = x? varn(x): 0;
323 if (flag) retmkvec2(pollegendre(n-1,v), pollegendre(n,v));
324 return pollegendre(n, v);
325 }
326 if (n==0)
327 {
328 if (flag) retmkvec2(gen_1, gcopy(x));
329 return gen_1;
330 }
331 if (n==1)
332 {
333 if (flag) retmkvec2(gcopy(x), gen_1);
334 return gcopy(x);
335 }
336 av = avma; v = gen_1; u = x;
337 for (i=1; i<n; i++)
338 { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
339 GEN t;
340 if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
341 t = gdivgs(gsub(gmul(gmulsg(2*i+1,x), u), gmulsg(i,v)), i+1);
342 v = u; u = t;
343 }
344 if (flag) return gerepilecopy(av, mkvec2(v, u));
345 return gerepileupto(av, u);
346 }
347 GEN
pollegendre_eval(long n,GEN x)348 pollegendre_eval(long n, GEN x) { return pollegendre_eval0(n, x, 0); }
349
350 /* Laguerre polynomial
351 * L0^a = 1; L1^a = -X+a+1;
352 * (n+1)*L^a(n+1) = (-X+(2*n+a+1))*L^a(n) - (n+a)*L^a(n-1)
353 * L^a(n) = sum_{k=0}^n (-1)^k * binom(n+a,n-k) * x^k/k! */
354 GEN
pollaguerre(long n,GEN a,long v)355 pollaguerre(long n, GEN a, long v)
356 {
357 pari_sp av = avma;
358 GEN L = cgetg(n+3, t_POL), c1 = gen_1, c2 = mpfact(n);
359 long i;
360
361 L[1] = evalsigne(1) | evalvarn(v);
362 if (odd(n)) togglesign_safe(&c2);
363 for (i = n; i >= 0; i--)
364 {
365 gel(L, i+2) = gdiv(c1, c2);
366 if (i)
367 {
368 c2 = divis(c2,-i);
369 c1 = gdivgs(gmul(c1, gaddsg(i,a)), n+1-i);
370 }
371 }
372 return gerepilecopy(av, L);
373 }
374 static void
err_lag(long n)375 err_lag(long n)
376 { pari_err_DOMAIN("pollaguerre", "degree", "<", gen_0, stoi(n)); }
377 GEN
pollaguerre_eval0(long n,GEN a,GEN x,long flag)378 pollaguerre_eval0(long n, GEN a, GEN x, long flag)
379 {
380 pari_sp av = avma;
381 long i;
382 GEN v, u;
383
384 if (n < 0) err_lag(n);
385 if (flag && flag != 1) pari_err_FLAG("pollaguerre");
386 if (!a) a = gen_0;
387 if (!x || gequalX(x))
388 {
389 long v = x? varn(x): 0;
390 if (flag)
391 {
392 if (!n) err_lag(-1);
393 retmkvec2(pollaguerre(n-1,a,v), pollaguerre(n,a,v));
394 }
395 return pollaguerre(n,a,v);
396 }
397 if (n==0)
398 {
399 if (flag) err_lag(-1);
400 return gen_1;
401 }
402 if (n==1)
403 {
404 if (flag) retmkvec2(gsub(gaddgs(a,1),x), gen_1);
405 return gsub(gaddgs(a,1),x);
406 }
407 av = avma; v = gen_1; u = gsub(gaddgs(a,1),x);
408 for (i=1; i<n; i++)
409 { /* u = P_i(x), v = P_{i-1}(x), compute t = P_{i+1}(x) */
410 GEN t;
411 if ((i & 0xff) == 0) gerepileall(av,2,&u, &v);
412 t = gdivgs(gsub(gmul(gsub(gaddsg(2*i+1,a),x), u), gmul(gaddsg(i,a),v)), i+1);
413 v = u; u = t;
414 }
415 if (flag) return gerepilecopy(av, mkvec2(v, u));
416 return gerepileupto(av, u);
417 }
418 GEN
pollaguerre_eval(long n,GEN x,GEN a)419 pollaguerre_eval(long n, GEN x, GEN a) { return pollaguerre_eval0(n, x, a, 0); }
420
421 /* polcyclo(p) = X^(p-1) + ... + 1 */
422 static GEN
polcyclo_prime(long p,long v)423 polcyclo_prime(long p, long v)
424 {
425 GEN T = cgetg(p+2, t_POL);
426 long i;
427 T[1] = evalsigne(1) | evalvarn(v);
428 for (i = 2; i < p+2; i++) gel(T,i) = gen_1;
429 return T;
430 }
431
432 /* cyclotomic polynomial */
433 GEN
polcyclo(long n,long v)434 polcyclo(long n, long v)
435 {
436 long s, q, i, l;
437 pari_sp av=avma;
438 GEN T, P;
439
440 if (v<0) v = 0;
441 if (n < 3)
442 switch(n)
443 {
444 case 1: return deg1pol_shallow(gen_1, gen_m1, v);
445 case 2: return deg1pol_shallow(gen_1, gen_1, v);
446 default: pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
447 }
448 P = gel(factoru(n), 1); l = lg(P);
449 s = P[1]; T = polcyclo_prime(s, v);
450 for (i = 2; i < l; i++)
451 { /* Phi_{np}(X) = Phi_n(X^p) / Phi_n(X) */
452 s *= P[i];
453 T = RgX_div(RgX_inflate(T, P[i]), T);
454 }
455 /* s = squarefree part of n */
456 q = n / s;
457 if (q == 1) return gerepileupto(av, T);
458 return gerepilecopy(av, RgX_inflate(T,q));
459 }
460
461 /* cyclotomic polynomial */
462 GEN
polcyclo_eval(long n,GEN x)463 polcyclo_eval(long n, GEN x)
464 {
465 pari_sp av= avma;
466 GEN P, md, xd, yneg, ypos;
467 long l, s, i, j, q, tx;
468 long root_of_1 = 0;
469
470 if (!x) return polcyclo(n, 0);
471 tx = typ(x);
472 if (gequalX(x)) return polcyclo(n, varn(x));
473 if (n <= 0) pari_err_DOMAIN("polcyclo", "index", "<=", gen_0, stoi(n));
474 if (n == 1) return gsubgs(x, 1);
475 if (tx == t_INT && !signe(x)) return gen_1;
476 while ((n & 3) == 0) { n >>= 1; x = gsqr(x); } /* Phi_4n(x) = Phi_2n(x^2) */
477 /* n not divisible by 4 */
478 if (n == 2) return gerepileupto(av, gaddgs(x,1));
479 if (!odd(n)) { n >>= 1; x = gneg(x); } /* Phi_2n(x) = Phi_n(-x) for n>1 odd */
480 /* n odd > 2. s largest squarefree divisor of n */
481 P = gel(factoru(n), 1); s = zv_prod(P);
482 /* replace n by largest squarefree divisor */
483 q = n/s; if (q != 1) { x = gpowgs(x, q); n = s; }
484 l = lg(P)-1;
485 /* n squarefree odd > 2, l distinct prime divisors. Now handle x = 1 or -1 */
486 if (tx == t_INT) { /* shortcut */
487 if (is_pm1(x))
488 {
489 set_avma(av);
490 if (signe(x) > 0 && l == 1) return utoipos(P[1]);
491 return gen_1;
492 }
493 } else {
494 if (gequal1(x))
495 { /* n is prime, return n; multiply by x to keep the type */
496 if (l == 1) return gerepileupto(av, gmulgs(x,n));
497 return gerepilecopy(av, x); /* else 1 */
498 }
499 if (gequalm1(x)) return gerepileupto(av, gneg(x)); /* -1 */
500 }
501 /* Heuristic: evaluation will probably not improve things */
502 if (tx == t_POL || tx == t_MAT || lg(x) > n)
503 return gerepileupto(av, poleval(polcyclo(n,0), x));
504
505 xd = cgetg((1L<<l) + 1, t_VEC); /* the x^d, where d | n */
506 md = cgetg((1L<<l) + 1, t_VECSMALL); /* the mu(d), where d | n */
507 gel(xd, 1) = x;
508 md[1] = 1;
509 /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d).
510 * If x has exact order D, n = Dq, then the result is 0 if q = 1. Otherwise
511 * the factors with x^d-1, D|d are omitted and we multiply at the end by
512 * prod_{d | q} d^mu(q/d) = q if prime, 1 otherwise */
513 /* We store the factors with mu(d)= 1 (resp.-1) in ypos (resp yneg).
514 * At the end we return ypos/yneg if mu(n)=1 and yneg/ypos if mu(n)=-1 */
515 ypos = gsubgs(x,1);
516 yneg = gen_1;
517 for (i = 1; i <= l; i++)
518 {
519 long ti = 1L<<(i-1), p = P[i];
520 for (j = 1; j <= ti; j++) {
521 GEN X = gpowgs(gel(xd,j), p), t = gsubgs(X,1);
522 gel(xd,ti+j) = X;
523 md[ti+j] = -md[j];
524 if (gequal0(t))
525 { /* x^d = 1; root_of_1 := the smallest index ti+j such that X == 1
526 * (whose bits code d: bit i-1 is set iff P[i] | d). If no such index
527 * exists, then root_of_1 remains 0. Do not multiply with X-1 if X = 1,
528 * we handle these factors at the end */
529 if (!root_of_1) root_of_1 = ti+j;
530 }
531 else
532 {
533 if (md[ti+j] == 1) ypos = gmul(ypos, t);
534 else yneg = gmul(yneg, t);
535 }
536 }
537 }
538 ypos = odd(l)? gdiv(yneg,ypos): gdiv(ypos,yneg);
539 if (root_of_1)
540 {
541 GEN X = gel(xd,(1<<l)); /* = x^n = 1 */
542 long bitmask_q = (1<<l) - root_of_1;
543 /* bitmask_q encodes q = n/d: bit (i-1) is 1 iff P[i] | q */
544
545 /* x is a root of unity. If bitmask_q = 0, then x was a primitive n-th
546 * root of 1 and the result is zero. Return X - 1 to preserve type. */
547 if (!bitmask_q) return gerepileupto(av, gsubgs(X, 1));
548 /* x is a primitive d-th root of unity, where d|n and d<n: we
549 * must multiply ypos by if(isprime(n/d), n/d, 1) */
550 ypos = gmul(ypos, X); /* multiply by X = 1 to preserve type */
551 /* If bitmask_q = 1<<(i-1) for some i <= l, then q == P[i] and we multiply
552 * by P[i]; otherwise q is composite and nothing more needs to be done */
553 if (!(bitmask_q & (bitmask_q-1))) /* detects power of 2, since bitmask!=0 */
554 {
555 i = vals(bitmask_q)+1; /* q = P[i] */
556 ypos = gmulgs(ypos, P[i]);
557 }
558 }
559 return gerepileupto(av, ypos);
560 }
561 /********************************************************************/
562 /** **/
563 /** HILBERT & PASCAL MATRICES **/
564 /** **/
565 /********************************************************************/
566 GEN
mathilbert(long n)567 mathilbert(long n) /* Hilbert matrix of order n */
568 {
569 long i,j;
570 GEN p;
571
572 if (n < 0) pari_err_DOMAIN("mathilbert", "dimension", "<", gen_0, stoi(n));
573 p = cgetg(n+1,t_MAT);
574 for (j=1; j<=n; j++)
575 {
576 gel(p,j) = cgetg(n+1,t_COL);
577 for (i=1+(j==1); i<=n; i++)
578 gcoeff(p,i,j) = mkfrac(gen_1, utoipos(i+j-1));
579 }
580 if (n) gcoeff(p,1,1) = gen_1;
581 return p;
582 }
583
584 /* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
585 GEN
matqpascal(long n,GEN q)586 matqpascal(long n, GEN q)
587 {
588 long i, j, I;
589 pari_sp av = avma;
590 GEN m, qpow = NULL; /* gcc -Wall */
591
592 if (n < -1) pari_err_DOMAIN("matpascal", "n", "<", gen_m1, stoi(n));
593 n++; m = cgetg(n+1,t_MAT);
594 for (j=1; j<=n; j++) gel(m,j) = cgetg(n+1,t_COL);
595 if (q)
596 {
597 I = (n+1)/2;
598 if (I > 1) { qpow = new_chunk(I+1); gel(qpow,2)=q; }
599 for (j=3; j<=I; j++) gel(qpow,j) = gmul(q, gel(qpow,j-1));
600 }
601 for (i=1; i<=n; i++)
602 {
603 I = (i+1)/2; gcoeff(m,i,1)= gen_1;
604 if (q)
605 {
606 for (j=2; j<=I; j++)
607 gcoeff(m,i,j) = gadd(gmul(gel(qpow,j),gcoeff(m,i-1,j)),
608 gcoeff(m,i-1,j-1));
609 }
610 else
611 {
612 for (j=2; j<=I; j++)
613 gcoeff(m,i,j) = addii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
614 }
615 for ( ; j<=i; j++) gcoeff(m,i,j) = gcoeff(m,i,i+1-j);
616 for ( ; j<=n; j++) gcoeff(m,i,j) = gen_0;
617 }
618 return gerepilecopy(av, m);
619 }
620
621 GEN
eulerianpol(long N,long v)622 eulerianpol(long N, long v)
623 {
624 pari_sp av = avma;
625 long n, n2, k = 0;
626 GEN A;
627 if (v < 0) v = 0;
628 if (N <= 0) pari_err_DOMAIN("eulerianpol", "index", "<=", gen_0, stoi(N));
629 if (N == 1) return pol_1(v);
630 if (N == 2) return deg1pol_shallow(gen_1, gen_1, v);
631 A = cgetg(N+1, t_VEC);
632 gel(A,1) = gen_1; gel(A,2) = gen_1; /* A_2 = x+1 */
633 for (n = 3; n <= N; n++)
634 { /* A(n,k) = (n-k)A(n-1,k-1) + (k+1)A(n-1,k) */
635 n2 = n >> 1;
636 if (odd(n)) gel(A,n2+1) = mului(n+1, gel(A,n2));
637 for (k = n2-1; k; k--)
638 gel(A,k+1) = addii(mului(n-k, gel(A,k)), mului(k+1, gel(A,k+1)));
639 if (gc_needed(av,1))
640 {
641 if (DEBUGMEM>1) pari_warn(warnmem,"eulerianpol, %ld/%ld",n,N);
642 for (k = odd(n)? n2+1: n2; k < N; k++) gel(A,k+1) = gen_0;
643 A = gerepilecopy(av, A);
644 }
645 }
646 k = N >> 1; if (odd(N)) k++;
647 for (; k < N; k++) gel(A,k+1) = gel(A, N-k);
648 return gerepilecopy(av, RgV_to_RgX(A, v));
649 }
650
651 /******************************************************************/
652 /** **/
653 /** PRECISION CHANGES **/
654 /** **/
655 /******************************************************************/
656
657 GEN
gprec(GEN x,long d)658 gprec(GEN x, long d)
659 {
660 pari_sp av = avma;
661 if (d <= 0) pari_err_DOMAIN("gprec", "precision", "<=", gen_0, stoi(d));
662 return gerepilecopy(av, gprec_w(x, ndec2prec(d)));
663 }
664
665 /* not GC-safe; precision given in word length (including codewords) */
666 GEN
gprec_w(GEN x,long pr)667 gprec_w(GEN x, long pr)
668 {
669 long lx, i;
670 GEN y;
671
672 switch(typ(x))
673 {
674 case t_REAL:
675 if (signe(x)) return realprec(x) != pr? rtor(x,pr): x;
676 i = -prec2nbits(pr);
677 if (i < expo(x)) return real_0_bit(i);
678 y = cgetr(2); y[1] = x[1]; return y;
679 case t_COMPLEX:
680 y = cgetg(3, t_COMPLEX);
681 gel(y,1) = gprec_w(gel(x,1),pr);
682 gel(y,2) = gprec_w(gel(x,2),pr);
683 break;
684 case t_POL: case t_SER:
685 y = cgetg_copy(x, &lx); y[1] = x[1];
686 for (i=2; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
687 break;
688 case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
689 y = cgetg_copy(x, &lx);
690 for (i=1; i<lx; i++) gel(y,i) = gprec_w(gel(x,i),pr);
691 break;
692 default: return x;
693 }
694 return y;
695 }
696 /* not GC-safe; precision given in word length (including codewords) */
697 GEN
gprec_wensure(GEN x,long pr)698 gprec_wensure(GEN x, long pr)
699 {
700 long lx, i;
701 GEN y;
702
703 switch(typ(x))
704 {
705 case t_REAL:
706 if (signe(x)) return realprec(x) < pr? rtor(x,pr): x;
707 i = -prec2nbits(pr);
708 if (i < expo(x)) return real_0_bit(i);
709 y = cgetr(2); y[1] = x[1]; return y;
710 case t_COMPLEX:
711 y = cgetg(3, t_COMPLEX);
712 gel(y,1) = gprec_wensure(gel(x,1),pr);
713 gel(y,2) = gprec_wensure(gel(x,2),pr);
714 break;
715 case t_POL: case t_SER:
716 y = cgetg_copy(x, &lx); y[1] = x[1];
717 for (i=2; i<lx; i++) gel(y,i) = gprec_wensure(gel(x,i),pr);
718 break;
719 case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
720 y = cgetg_copy(x, &lx);
721 for (i=1; i<lx; i++) gel(y,i) = gprec_wensure(gel(x,i),pr);
722 break;
723 default: return x;
724 }
725 return y;
726 }
727
728 /* not GC-safe; precision given in word length (including codewords),
729 * truncate mantissa to precision 'pr' but never increase it */
730 GEN
gprec_wtrunc(GEN x,long pr)731 gprec_wtrunc(GEN x, long pr)
732 {
733 long lx, i;
734 GEN y;
735
736 switch(typ(x))
737 {
738 case t_REAL:
739 return (signe(x) && realprec(x) > pr)? rtor(x,pr): x;
740 case t_COMPLEX:
741 y = cgetg(3, t_COMPLEX);
742 gel(y,1) = gprec_wtrunc(gel(x,1),pr);
743 gel(y,2) = gprec_wtrunc(gel(x,2),pr);
744 break;
745 case t_POL:
746 case t_SER:
747 y = cgetg_copy(x, &lx); y[1] = x[1];
748 for (i=2; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
749 break;
750 case t_POLMOD: case t_RFRAC: case t_VEC: case t_COL: case t_MAT:
751 y = cgetg_copy(x, &lx);
752 for (i=1; i<lx; i++) gel(y,i) = gprec_wtrunc(gel(x,i),pr);
753 break;
754 default: return x;
755 }
756 return y;
757 }
758
759 /********************************************************************/
760 /** **/
761 /** SERIES TRANSFORMS **/
762 /** **/
763 /********************************************************************/
764 /** LAPLACE TRANSFORM (OF A SERIES) **/
765 /********************************************************************/
766 static GEN
serlaplace(GEN x)767 serlaplace(GEN x)
768 {
769 long i, l = lg(x), e = valp(x);
770 GEN t, y = cgetg(l,t_SER);
771 if (e < 0) pari_err_DOMAIN("laplace","valuation","<",gen_0,stoi(e));
772 t = mpfact(e); y[1] = x[1];
773 for (i=2; i<l; i++)
774 {
775 gel(y,i) = gmul(t, gel(x,i));
776 e++; t = mului(e,t);
777 }
778 return y;
779 }
780 static GEN
pollaplace(GEN x)781 pollaplace(GEN x)
782 {
783 long i, e = 0, l = lg(x);
784 GEN t = gen_1, y = cgetg(l,t_POL);
785 y[1] = x[1];
786 for (i=2; i<l; i++)
787 {
788 gel(y,i) = gmul(t, gel(x,i));
789 e++; t = mului(e,t);
790 }
791 return y;
792 }
793 GEN
laplace(GEN x)794 laplace(GEN x)
795 {
796 pari_sp av = avma;
797 switch(typ(x))
798 {
799 case t_POL: x = pollaplace(x); break;
800 case t_SER: x = serlaplace(x); break;
801 default: if (is_scalar_t(typ(x))) return gcopy(x);
802 pari_err_TYPE("laplace",x);
803 }
804 return gerepilecopy(av, x);
805 }
806
807 /********************************************************************/
808 /** CONVOLUTION PRODUCT (OF TWO SERIES) **/
809 /********************************************************************/
810 GEN
convol(GEN x,GEN y)811 convol(GEN x, GEN y)
812 {
813 long j, lx, ly, ex, ey, vx = varn(x);
814 GEN z;
815
816 if (typ(x) != t_SER) pari_err_TYPE("convol",x);
817 if (typ(y) != t_SER) pari_err_TYPE("convol",y);
818 if (varn(y) != vx) pari_err_VAR("convol", x,y);
819 ex = valp(x);
820 ey = valp(y);
821 if (ser_isexactzero(x))
822 return scalarser(gadd(Rg_get_0(x), Rg_get_0(y)), varn(x), maxss(ex,ey));
823 lx = lg(x) + ex; x -= ex;
824 ly = lg(y) + ey; y -= ey;
825 /* inputs shifted: x[i] and y[i] now correspond to monomials of same degree */
826 if (ly < lx) lx = ly; /* min length */
827 if (ex < ey) ex = ey; /* max valuation */
828 if (lx - ex < 3) return zeroser(vx, lx-2);
829
830 z = cgetg(lx - ex, t_SER);
831 z[1] = evalvalp(ex) | evalvarn(vx);
832 for (j = ex+2; j<lx; j++) gel(z,j-ex) = gmul(gel(x,j),gel(y,j));
833 return normalize(z);
834 }
835
836 /***********************************************************************/
837 /* OPERATIONS ON DIRICHLET SERIES: *, / */
838 /* (+, -, scalar multiplication are done on the corresponding vectors) */
839 /***********************************************************************/
840 static long
dirval(GEN x)841 dirval(GEN x)
842 {
843 long i = 1, lx = lg(x);
844 while (i < lx && gequal0(gel(x,i))) i++;
845 return i;
846 }
847
848 GEN
dirmul(GEN x,GEN y)849 dirmul(GEN x, GEN y)
850 {
851 pari_sp av = avma, av2;
852 long nx, ny, nz, dx, dy, i, j, k;
853 GEN z;
854
855 if (typ(x)!=t_VEC) pari_err_TYPE("dirmul",x);
856 if (typ(y)!=t_VEC) pari_err_TYPE("dirmul",y);
857 dx = dirval(x); nx = lg(x)-1;
858 dy = dirval(y); ny = lg(y)-1;
859 if (ny-dy < nx-dx) { swap(x,y); lswap(nx,ny); lswap(dx,dy); }
860 nz = minss(nx*dy,ny*dx);
861 y = RgV_kill0(y);
862 av2 = avma;
863 z = zerovec(nz);
864 for (j=dx; j<=nx; j++)
865 {
866 GEN c = gel(x,j);
867 if (gequal0(c)) continue;
868 if (gequal1(c))
869 {
870 for (k=dy,i=j*dy; i<=nz; i+=j,k++)
871 if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gel(y,k));
872 }
873 else if (gequalm1(c))
874 {
875 for (k=dy,i=j*dy; i<=nz; i+=j,k++)
876 if (gel(y,k)) gel(z,i) = gsub(gel(z,i),gel(y,k));
877 }
878 else
879 {
880 for (k=dy,i=j*dy; i<=nz; i+=j,k++)
881 if (gel(y,k)) gel(z,i) = gadd(gel(z,i),gmul(c,gel(y,k)));
882 }
883 if (gc_needed(av2,3))
884 {
885 if (DEBUGMEM>1) pari_warn(warnmem,"dirmul, %ld/%ld",j,nx);
886 z = gerepilecopy(av2,z);
887 }
888 }
889 return gerepilecopy(av,z);
890 }
891
892 GEN
dirdiv(GEN x,GEN y)893 dirdiv(GEN x, GEN y)
894 {
895 pari_sp av = avma, av2;
896 long nx,ny,nz, dx,dy, i,j,k;
897 GEN p1;
898
899 if (typ(x)!=t_VEC) pari_err_TYPE("dirdiv",x);
900 if (typ(y)!=t_VEC) pari_err_TYPE("dirdiv",y);
901 dx = dirval(x); nx = lg(x)-1;
902 dy = dirval(y); ny = lg(y)-1;
903 if (dy != 1 || !ny) pari_err_INV("dirdiv",y);
904 nz = minss(nx,ny*dx);
905 p1 = gel(y,1);
906 if (gequal1(p1)) p1 = NULL; else y = gdiv(y,p1);
907 y = RgV_kill0(y);
908 av2 = avma;
909 x = p1 ? gdiv(x,p1): leafcopy(x);
910 for (j=1; j<dx; j++) gel(x,j) = gen_0;
911 setlg(x,nz+1);
912 for (j=dx; j<=nz; j++)
913 {
914 GEN c = gel(x,j);
915 if (gequal0(c)) continue;
916 if (gequal1(c))
917 {
918 for (i=j+j,k=2; i<=nz; i+=j,k++)
919 if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gel(y,k));
920 }
921 else if (gequalm1(c))
922 {
923 for (i=j+j,k=2; i<=nz; i+=j,k++)
924 if (gel(y,k)) gel(x,i) = gadd(gel(x,i),gel(y,k));
925 }
926 else
927 {
928 for (i=j+j,k=2; i<=nz; i+=j,k++)
929 if (gel(y,k)) gel(x,i) = gsub(gel(x,i),gmul(c,gel(y,k)));
930 }
931 if (gc_needed(av2,3))
932 {
933 if (DEBUGMEM>1) pari_warn(warnmem,"dirdiv, %ld/%ld",j,nz);
934 x = gerepilecopy(av2,x);
935 }
936 }
937 return gerepilecopy(av,x);
938 }
939
940 /*******************************************************************/
941 /** **/
942 /** COMBINATORICS **/
943 /** **/
944 /*******************************************************************/
945 /** BINOMIAL COEFFICIENTS **/
946 /*******************************************************************/
947 GEN
binomialuu(ulong n,ulong k)948 binomialuu(ulong n, ulong k)
949 {
950 pari_sp ltop = avma;
951 GEN z;
952 if (k > n) return gen_0;
953 k = minuu(k,n-k);
954 if (!k) return gen_1;
955 if (k == 1) return utoipos(n);
956 z = diviiexact(mulu_interval(n-k+1, n), mulu_interval(2UL, k));
957 return gerepileuptoint(ltop,z);
958 }
959
960 GEN
binomial(GEN n,long k)961 binomial(GEN n, long k)
962 {
963 long i, prec;
964 pari_sp av;
965 GEN y;
966
967 if (k <= 1)
968 {
969 if (is_noncalc_t(typ(n))) pari_err_TYPE("binomial",n);
970 if (k < 0) return gen_0;
971 if (k == 0) return gen_1;
972 return gcopy(n);
973 }
974 av = avma;
975 if (typ(n) == t_INT)
976 {
977 if (signe(n) > 0)
978 {
979 GEN z = subiu(n,k);
980 if (cmpis(z,k) < 0)
981 {
982 k = itos(z); set_avma(av);
983 if (k <= 1)
984 {
985 if (k < 0) return gen_0;
986 if (k == 0) return gen_1;
987 return icopy(n);
988 }
989 }
990 }
991 /* k > 1 */
992 if (lgefint(n) == 3 && signe(n) > 0)
993 {
994 y = binomialuu(itou(n),(ulong)k);
995 return gerepileupto(av, y);
996 }
997 else
998 {
999 y = cgetg(k+1,t_VEC);
1000 for (i=1; i<=k; i++) gel(y,i) = subiu(n,i-1);
1001 y = ZV_prod(y);
1002 }
1003 y = diviiexact(y, mpfact(k));
1004 return gerepileuptoint(av, y);
1005 }
1006
1007 prec = precision(n);
1008 if (prec && k > 200 + 0.8*prec2nbits(prec)) {
1009 GEN A = mpfactr(k, prec), B = ggamma(gsubgs(n,k-1), prec);
1010 return gerepileupto(av, gdiv(ggamma(gaddgs(n,1), prec), gmul(A,B)));
1011 }
1012
1013 y = cgetg(k+1,t_VEC);
1014 for (i=1; i<=k; i++) gel(y,i) = gsubgs(n,i-1);
1015 return gerepileupto(av, gdiv(RgV_prod(y), mpfact(k)));
1016 }
1017
1018 GEN
binomial0(GEN x,GEN k)1019 binomial0(GEN x, GEN k)
1020 {
1021 if (!k)
1022 {
1023 if (typ(x) != t_INT || signe(x) < 0) pari_err_TYPE("binomial", x);
1024 return vecbinomial(itos(x));
1025 }
1026 if (typ(k) != t_INT) pari_err_TYPE("binomial", k);
1027 return binomial(x, itos(k));
1028 }
1029
1030 /* Assume n >= 0, return bin, bin[k+1] = binomial(n, k) */
1031 GEN
vecbinomial(long n)1032 vecbinomial(long n)
1033 {
1034 long d, k;
1035 GEN C;
1036 if (!n) return mkvec(gen_1);
1037 C = cgetg(n+2, t_VEC) + 1; /* C[k] = binomial(n, k) */
1038 gel(C,0) = gen_1;
1039 gel(C,1) = utoipos(n); d = (n + 1) >> 1;
1040 for (k=2; k <= d; k++)
1041 {
1042 pari_sp av = avma;
1043 gel(C,k) = gerepileuptoint(av, diviuexact(mului(n-k+1, gel(C,k-1)), k));
1044 }
1045 for ( ; k <= n; k++) gel(C,k) = gel(C,n-k);
1046 return C - 1;
1047 }
1048
1049 /********************************************************************/
1050 /** STIRLING NUMBERS **/
1051 /********************************************************************/
1052 /* Stirling number of the 2nd kind. The number of ways of partitioning
1053 a set of n elements into m nonempty subsets. */
1054 GEN
stirling2(ulong n,ulong m)1055 stirling2(ulong n, ulong m)
1056 {
1057 pari_sp av = avma;
1058 GEN s, bmk;
1059 ulong k;
1060 if (n==0) return (m == 0)? gen_1: gen_0;
1061 if (m > n || m == 0) return gen_0;
1062 if (m==n) return gen_1;
1063 /* k = 0 */
1064 bmk = gen_1; s = powuu(m, n);
1065 for (k = 1; k <= ((m-1)>>1); ++k)
1066 { /* bmk = binomial(m, k) */
1067 GEN c, kn, mkn;
1068 bmk = diviuexact(mului(m-k+1, bmk), k);
1069 kn = powuu(k, n); mkn = powuu(m-k, n);
1070 c = odd(m)? subii(mkn,kn): addii(mkn,kn);
1071 c = mulii(bmk, c);
1072 s = odd(k)? subii(s, c): addii(s, c);
1073 if (gc_needed(av,2))
1074 {
1075 if(DEBUGMEM>1) pari_warn(warnmem,"stirling2");
1076 gerepileall(av, 2, &s, &bmk);
1077 }
1078 }
1079 /* k = m/2 */
1080 if (!odd(m))
1081 {
1082 GEN c;
1083 bmk = diviuexact(mului(k+1, bmk), k);
1084 c = mulii(bmk, powuu(k,n));
1085 s = odd(k)? subii(s, c): addii(s, c);
1086 }
1087 return gerepileuptoint(av, diviiexact(s, mpfact(m)));
1088 }
1089
1090 /* Stirling number of the first kind. Up to the sign, the number of
1091 permutations of n symbols which have exactly m cycles. */
1092 GEN
stirling1(ulong n,ulong m)1093 stirling1(ulong n, ulong m)
1094 {
1095 pari_sp ltop=avma;
1096 ulong k;
1097 GEN s, t;
1098 if (n < m) return gen_0;
1099 else if (n==m) return gen_1;
1100 /* t = binomial(n-1+k, m-1) * binomial(2n-m, n-m-k) */
1101 /* k = n-m > 0 */
1102 t = binomialuu(2*n-m-1, m-1);
1103 s = mulii(t, stirling2(2*(n-m), n-m));
1104 if (odd(n-m)) togglesign(s);
1105 for (k = n-m-1; k > 0; --k)
1106 {
1107 GEN c;
1108 t = diviuuexact(muluui(n-m+k+1, n+k+1, t), n+k, n-m-k);
1109 c = mulii(t, stirling2(n-m+k, k));
1110 s = odd(k)? subii(s, c): addii(s, c);
1111 if ((k & 0x1f) == 0) {
1112 t = gerepileuptoint(ltop, t);
1113 s = gerepileuptoint(avma, s);
1114 }
1115 }
1116 return gerepileuptoint(ltop, s);
1117 }
1118
1119 GEN
stirling(long n,long m,long flag)1120 stirling(long n, long m, long flag)
1121 {
1122 if (n < 0) pari_err_DOMAIN("stirling", "n", "<", gen_0, stoi(n));
1123 if (m < 0) pari_err_DOMAIN("stirling", "m", "<", gen_0, stoi(m));
1124 switch (flag)
1125 {
1126 case 1: return stirling1((ulong)n,(ulong)m);
1127 case 2: return stirling2((ulong)n,(ulong)m);
1128 default: pari_err_FLAG("stirling");
1129 }
1130 return NULL; /*LCOV_EXCL_LINE*/
1131 }
1132
1133 /*******************************************************************/
1134 /** **/
1135 /** RECIPROCAL POLYNOMIAL **/
1136 /** **/
1137 /*******************************************************************/
1138 /* return coefficients s.t x = x_0 X^n + ... + x_n */
1139 GEN
polrecip(GEN x)1140 polrecip(GEN x)
1141 {
1142 long tx = typ(x);
1143 if (is_scalar_t(tx)) return gcopy(x);
1144 if (tx != t_POL) pari_err_TYPE("polrecip",x);
1145 return RgX_recip(x);
1146 }
1147
1148 /********************************************************************/
1149 /** **/
1150 /** POLYNOMIAL INTERPOLATION **/
1151 /** **/
1152 /********************************************************************/
1153 static GEN
RgV_polint_fast(GEN X,GEN Y,long v)1154 RgV_polint_fast(GEN X, GEN Y, long v)
1155 {
1156 GEN p, pol;
1157 long t, pa;
1158 if (X) t = RgV_type2(X,Y, &p, &pol, &pa);
1159 else t = Rg_type(Y, &p, &pol, &pa);
1160 if (t != t_INTMOD) return NULL;
1161 Y = RgC_to_FpC(Y, p);
1162 X = X? RgC_to_FpC(X, p): identity_ZV(lg(Y)-1);
1163 return FpX_to_mod(FpV_polint(X, Y, p, v), p);
1164 }
1165 /* allow X = NULL for [1,...,n] */
1166 GEN
RgV_polint(GEN X,GEN Y,long v)1167 RgV_polint(GEN X, GEN Y, long v)
1168 {
1169 pari_sp av0 = avma, av;
1170 GEN Q, P = NULL;
1171 long i, l = lg(Y);
1172 if ((Q = RgV_polint_fast(X,Y,v))) return Q;
1173 if (!X) X = identity_ZV(l-1);
1174 Q = roots_to_pol(X, v); av = avma;
1175 for (i=1; i<l; i++)
1176 {
1177 GEN inv, T, dP;
1178 if (gequal0(gel(Y,i))) continue;
1179 T = RgX_div_by_X_x(Q, gel(X,i), NULL);
1180 inv = ginv(poleval(T,gel(X,i)));
1181 dP = RgX_Rg_mul(T, gmul(gel(Y,i),inv));
1182 P = P? RgX_add(P, dP): dP;
1183 if (gc_needed(av,2))
1184 {
1185 if (DEBUGMEM>1) pari_warn(warnmem,"RgV_polint i = %ld/%ld", i, l-1);
1186 P = gerepileupto(av, P);
1187 }
1188 }
1189 if (!P) { set_avma(av); return zeropol(v); }
1190 return gerepileupto(av0, P);
1191 }
1192 static int
inC(GEN x)1193 inC(GEN x)
1194 {
1195 switch(typ(x)) {
1196 case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD: return 1;
1197 default: return 0;
1198 }
1199 }
1200 static long
check_dy(GEN X,GEN x,long n)1201 check_dy(GEN X, GEN x, long n)
1202 {
1203 GEN D = NULL;
1204 long i, ns = 0;
1205 if (!inC(x)) return -1;
1206 for (i = 0; i < n; i++)
1207 {
1208 GEN t = gsub(x, gel(X,i));
1209 if (!inC(t)) return -1;
1210 t = gabs(t, DEFAULTPREC);
1211 if (!D || gcmp(t,D) < 0) { ns = i; D = t; }
1212 }
1213 /* X[ns] is closest to x */
1214 return ns;
1215 }
1216 /* X,Y are "spec" GEN vectors with n > 0 components ( at X[0], ... X[n-1] ) */
1217 GEN
polintspec(GEN X,GEN Y,GEN x,long n,long * pe)1218 polintspec(GEN X, GEN Y, GEN x, long n, long *pe)
1219 {
1220 long i, m, ns;
1221 pari_sp av = avma, av2;
1222 GEN y, c, d, dy = NULL; /* gcc -Wall */
1223
1224 if (pe) *pe = -HIGHEXPOBIT;
1225 if (n == 1) return gmul(gel(Y,0), Rg_get_1(x));
1226 if (!X) X = identity_ZV(n) + 1;
1227 av2 = avma;
1228 ns = check_dy(X, x, n); if (ns < 0) { pe = NULL; ns = 0; }
1229 c = cgetg(n+1, t_VEC);
1230 d = cgetg(n+1, t_VEC); for (i=0; i<n; i++) gel(c,i+1) = gel(d,i+1) = gel(Y,i);
1231 y = gel(d,ns+1);
1232 /* divided differences */
1233 for (m = 1; m < n; m++)
1234 {
1235 for (i = 0; i < n-m; i++)
1236 {
1237 GEN ho = gsub(gel(X,i),x), hp = gsub(gel(X,i+m),x), den = gsub(ho,hp);
1238 if (gequal0(den))
1239 {
1240 char *x1 = stack_sprintf("X[%ld]", i+1);
1241 char *x2 = stack_sprintf("X[%ld]", i+m+1);
1242 pari_err_DOMAIN("polinterpolate",x1,"=",strtoGENstr(x2), X);
1243 }
1244 den = gdiv(gsub(gel(c,i+2),gel(d,i+1)), den);
1245 gel(c,i+1) = gmul(ho,den);
1246 gel(d,i+1) = gmul(hp,den);
1247 }
1248 dy = (2*ns < n-m)? gel(c,ns+1): gel(d,ns--);
1249 y = gadd(y,dy);
1250 if (gc_needed(av2,2))
1251 {
1252 if (DEBUGMEM>1) pari_warn(warnmem,"polint, %ld/%ld",m,n-1);
1253 gerepileall(av2, 4, &y, &c, &d, &dy);
1254 }
1255 }
1256 if (pe && inC(dy)) *pe = gexpo(dy);
1257 return gerepileupto(av, y);
1258 }
1259
1260 GEN
polint_i(GEN X,GEN Y,GEN t,long * pe)1261 polint_i(GEN X, GEN Y, GEN t, long *pe)
1262 {
1263 long lx = lg(X), vt;
1264
1265 if (! is_vec_t(typ(X))) pari_err_TYPE("polinterpolate",X);
1266 if (Y)
1267 {
1268 if (! is_vec_t(typ(Y))) pari_err_TYPE("polinterpolate",Y);
1269 if (lx != lg(Y)) pari_err_DIM("polinterpolate");
1270 }
1271 else
1272 {
1273 Y = X;
1274 X = NULL;
1275 }
1276 if (pe) *pe = -HIGHEXPOBIT;
1277 vt = t? gvar(t): 0;
1278 if (vt != NO_VARIABLE)
1279 { /* formal interpolation */
1280 pari_sp av;
1281 long v0, vY = gvar(Y);
1282 GEN P;
1283 if (X) vY = varnmax(vY, gvar(X));
1284 /* shortcut */
1285 if (varncmp(vY, vt) > 0 && (!t || gequalX(t))) return RgV_polint(X, Y, vt);
1286 av = avma;
1287 /* first interpolate in high priority variable, then substitute t */
1288 v0 = fetch_var_higher();
1289 P = RgV_polint(X, Y, v0);
1290 P = gsubst(P, v0, t? t: pol_x(0));
1291 (void)delete_var();
1292 return gerepileupto(av, P);
1293 }
1294 /* numerical interpolation */
1295 if (lx == 1) return Rg_get_0(t);
1296 return polintspec(X? X+1: NULL,Y+1,t,lx-1, pe);
1297 }
1298 GEN
polint(GEN X,GEN Y,GEN t,GEN * pe)1299 polint(GEN X, GEN Y, GEN t, GEN *pe)
1300 {
1301 long e;
1302 GEN p = polint_i(X, Y, t, &e);
1303 if (pe) *pe = stoi(e);
1304 return p;
1305 }
1306
1307 /********************************************************************/
1308 /** **/
1309 /** MODREVERSE **/
1310 /** **/
1311 /********************************************************************/
1312 static void
err_reverse(GEN x,GEN T)1313 err_reverse(GEN x, GEN T)
1314 {
1315 pari_err_DOMAIN("modreverse","deg(minpoly(z))", "<", stoi(degpol(T)),
1316 mkpolmod(x,T));
1317 }
1318
1319 /* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
1320 GEN
RgXQ_reverse(GEN a,GEN T)1321 RgXQ_reverse(GEN a, GEN T)
1322 {
1323 pari_sp av = avma;
1324 long n = degpol(T);
1325 GEN y;
1326
1327 if (n <= 1) {
1328 if (n <= 0) return gcopy(a);
1329 return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
1330 }
1331 if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
1332 y = RgXV_to_RgM(RgXQ_powers(a,n-1,T), n);
1333 y = RgM_solve(y, col_ei(n, 2));
1334 if (!y) err_reverse(a,T);
1335 return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
1336 }
1337 GEN
QXQ_reverse(GEN a,GEN T)1338 QXQ_reverse(GEN a, GEN T)
1339 {
1340 pari_sp av = avma;
1341 long n = degpol(T);
1342 GEN y;
1343
1344 if (n <= 1) {
1345 if (n <= 0) return gcopy(a);
1346 return gerepileupto(av, gneg(gdiv(gel(T,2), gel(T,3))));
1347 }
1348 if (typ(a) != t_POL || !signe(a)) err_reverse(a,T);
1349 if (gequalX(a)) return gcopy(a);
1350 y = RgXV_to_RgM(QXQ_powers(a,n-1,T), n);
1351 y = QM_gauss(y, col_ei(n, 2));
1352 if (!y) err_reverse(a,T);
1353 return gerepilecopy(av, RgV_to_RgX(y, varn(T)));
1354 }
1355
1356 GEN
modreverse(GEN x)1357 modreverse(GEN x)
1358 {
1359 long v, n;
1360 GEN T, a;
1361
1362 if (typ(x)!=t_POLMOD) pari_err_TYPE("modreverse",x);
1363 T = gel(x,1); n = degpol(T); if (n <= 0) return gcopy(x);
1364 a = gel(x,2);
1365 v = varn(T);
1366 retmkpolmod(RgXQ_reverse(a, T),
1367 (n==1)? gsub(pol_x(v), a): RgXQ_charpoly(a, T, v));
1368 }
1369
1370 /********************************************************************/
1371 /** **/
1372 /** MERGESORT **/
1373 /** **/
1374 /********************************************************************/
1375 static int
cmp_small(GEN x,GEN y)1376 cmp_small(GEN x, GEN y) {
1377 long a = (long)x, b = (long)y;
1378 return a>b? 1: (a<b? -1: 0);
1379 }
1380
1381 static int
veccmp(void * data,GEN x,GEN y)1382 veccmp(void *data, GEN x, GEN y)
1383 {
1384 GEN k = (GEN)data;
1385 long i, s, lk = lg(k), lx = minss(lg(x), lg(y));
1386
1387 if (!is_vec_t(typ(x))) pari_err_TYPE("lexicographic vecsort",x);
1388 if (!is_vec_t(typ(y))) pari_err_TYPE("lexicographic vecsort",y);
1389 for (i=1; i<lk; i++)
1390 {
1391 long c = k[i];
1392 if (c >= lx)
1393 pari_err_TYPE("lexicographic vecsort, index too large", stoi(c));
1394 s = lexcmp(gel(x,c), gel(y,c));
1395 if (s) return s;
1396 }
1397 return 0;
1398 }
1399
1400 /* return permutation sorting v[1..n], removing duplicates. Assume n > 0 */
1401 static GEN
gen_sortspec_uniq(GEN v,long n,void * E,int (* cmp)(void *,GEN,GEN))1402 gen_sortspec_uniq(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
1403 {
1404 pari_sp av;
1405 long NX, nx, ny, m, ix, iy, i;
1406 GEN x, y, w, W;
1407 int s;
1408 switch(n)
1409 {
1410 case 1: return mkvecsmall(1);
1411 case 2:
1412 s = cmp(E,gel(v,1),gel(v,2));
1413 if (s < 0) return mkvecsmall2(1,2);
1414 else if (s > 0) return mkvecsmall2(2,1);
1415 return mkvecsmall(1);
1416 case 3:
1417 s = cmp(E,gel(v,1),gel(v,2));
1418 if (s < 0) {
1419 s = cmp(E,gel(v,2),gel(v,3));
1420 if (s < 0) return mkvecsmall3(1,2,3);
1421 else if (s == 0) return mkvecsmall2(1,2);
1422 s = cmp(E,gel(v,1),gel(v,3));
1423 if (s < 0) return mkvecsmall3(1,3,2);
1424 else if (s > 0) return mkvecsmall3(3,1,2);
1425 return mkvecsmall2(1,2);
1426 } else if (s > 0) {
1427 s = cmp(E,gel(v,1),gel(v,3));
1428 if (s < 0) return mkvecsmall3(2,1,3);
1429 else if (s == 0) return mkvecsmall2(2,1);
1430 s = cmp(E,gel(v,2),gel(v,3));
1431 if (s < 0) return mkvecsmall3(2,3,1);
1432 else if (s > 0) return mkvecsmall3(3,2,1);
1433 return mkvecsmall2(2,1);
1434 } else {
1435 s = cmp(E,gel(v,1),gel(v,3));
1436 if (s < 0) return mkvecsmall2(1,3);
1437 else if (s == 0) return mkvecsmall(1);
1438 return mkvecsmall2(3,1);
1439 }
1440 }
1441 NX = nx = n>>1; ny = n-nx;
1442 av = avma;
1443 x = gen_sortspec_uniq(v, nx,E,cmp); nx = lg(x)-1;
1444 y = gen_sortspec_uniq(v+NX,ny,E,cmp); ny = lg(y)-1;
1445 w = cgetg(n+1, t_VECSMALL);
1446 m = ix = iy = 1;
1447 while (ix<=nx && iy<=ny)
1448 {
1449 s = cmp(E, gel(v,x[ix]), gel(v,y[iy]+NX));
1450 if (s < 0)
1451 w[m++] = x[ix++];
1452 else if (s > 0)
1453 w[m++] = y[iy++]+NX;
1454 else {
1455 w[m++] = x[ix++];
1456 iy++;
1457 }
1458 }
1459 while (ix<=nx) w[m++] = x[ix++];
1460 while (iy<=ny) w[m++] = y[iy++]+NX;
1461 set_avma(av);
1462 W = cgetg(m, t_VECSMALL);
1463 for (i = 1; i < m; i++) W[i] = w[i];
1464 return W;
1465 }
1466
1467 /* return permutation sorting v[1..n]. Assume n > 0 */
1468 static GEN
gen_sortspec(GEN v,long n,void * E,int (* cmp)(void *,GEN,GEN))1469 gen_sortspec(GEN v, long n, void *E, int (*cmp)(void*,GEN,GEN))
1470 {
1471 long nx, ny, m, ix, iy;
1472 GEN x, y, w;
1473 switch(n)
1474 {
1475 case 1:
1476 (void)cmp(E,gel(v,1),gel(v,1)); /* check for type error */
1477 return mkvecsmall(1);
1478 case 2:
1479 return cmp(E,gel(v,1),gel(v,2)) <= 0? mkvecsmall2(1,2)
1480 : mkvecsmall2(2,1);
1481 case 3:
1482 if (cmp(E,gel(v,1),gel(v,2)) <= 0) {
1483 if (cmp(E,gel(v,2),gel(v,3)) <= 0) return mkvecsmall3(1,2,3);
1484 return (cmp(E,gel(v,1),gel(v,3)) <= 0)? mkvecsmall3(1,3,2)
1485 : mkvecsmall3(3,1,2);
1486 } else {
1487 if (cmp(E,gel(v,1),gel(v,3)) <= 0) return mkvecsmall3(2,1,3);
1488 return (cmp(E,gel(v,2),gel(v,3)) <= 0)? mkvecsmall3(2,3,1)
1489 : mkvecsmall3(3,2,1);
1490 }
1491 }
1492 nx = n>>1; ny = n-nx;
1493 w = cgetg(n+1,t_VECSMALL);
1494 x = gen_sortspec(v, nx,E,cmp);
1495 y = gen_sortspec(v+nx,ny,E,cmp);
1496 m = ix = iy = 1;
1497 while (ix<=nx && iy<=ny)
1498 if (cmp(E, gel(v,x[ix]), gel(v,y[iy]+nx))<=0)
1499 w[m++] = x[ix++];
1500 else
1501 w[m++] = y[iy++]+nx;
1502 while (ix<=nx) w[m++] = x[ix++];
1503 while (iy<=ny) w[m++] = y[iy++]+nx;
1504 set_avma((pari_sp)w); return w;
1505 }
1506
1507 static void
init_sort(GEN * x,long * tx,long * lx)1508 init_sort(GEN *x, long *tx, long *lx)
1509 {
1510 *tx = typ(*x);
1511 if (*tx == t_LIST)
1512 {
1513 if (list_typ(*x)!=t_LIST_RAW) pari_err_TYPE("sort",*x);
1514 *x = list_data(*x);
1515 *lx = *x? lg(*x): 1;
1516 } else {
1517 if (!is_matvec_t(*tx) && *tx != t_VECSMALL) pari_err_TYPE("gen_sort",*x);
1518 *lx = lg(*x);
1519 }
1520 }
1521
1522 /* (x o y)[1..lx-1], destroy y */
1523 INLINE GEN
sort_extract(GEN x,GEN y,long tx,long lx)1524 sort_extract(GEN x, GEN y, long tx, long lx)
1525 {
1526 long i;
1527 switch(tx)
1528 {
1529 case t_VECSMALL:
1530 for (i=1; i<lx; i++) y[i] = x[y[i]];
1531 break;
1532 case t_LIST:
1533 settyp(y,t_VEC);
1534 for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
1535 return gtolist(y);
1536 default:
1537 settyp(y,tx);
1538 for (i=1; i<lx; i++) gel(y,i) = gcopy(gel(x,y[i]));
1539 }
1540 return y;
1541 }
1542
1543 static GEN
triv_sort(long tx)1544 triv_sort(long tx) { return tx == t_LIST? mklist(): cgetg(1, tx); }
1545 /* Sort the vector x, using cmp to compare entries. */
1546 GEN
gen_sort_uniq(GEN x,void * E,int (* cmp)(void *,GEN,GEN))1547 gen_sort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
1548 {
1549 long tx, lx;
1550 GEN y;
1551
1552 init_sort(&x, &tx, &lx);
1553 if (lx==1) return triv_sort(tx);
1554 y = gen_sortspec_uniq(x,lx-1,E,cmp);
1555 return sort_extract(x, y, tx, lg(y)); /* lg(y) <= lx */
1556 }
1557 /* Sort the vector x, using cmp to compare entries. */
1558 GEN
gen_sort(GEN x,void * E,int (* cmp)(void *,GEN,GEN))1559 gen_sort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
1560 {
1561 long tx, lx;
1562 GEN y;
1563
1564 init_sort(&x, &tx, &lx);
1565 if (lx==1) return triv_sort(tx);
1566 y = gen_sortspec(x,lx-1,E,cmp);
1567 return sort_extract(x, y, tx, lx);
1568 }
1569 /* indirect sort: return the permutation that would sort x */
1570 GEN
gen_indexsort_uniq(GEN x,void * E,int (* cmp)(void *,GEN,GEN))1571 gen_indexsort_uniq(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
1572 {
1573 long tx, lx;
1574 init_sort(&x, &tx, &lx);
1575 if (lx==1) return cgetg(1, t_VECSMALL);
1576 return gen_sortspec_uniq(x,lx-1,E,cmp);
1577 }
1578 /* indirect sort: return the permutation that would sort x */
1579 GEN
gen_indexsort(GEN x,void * E,int (* cmp)(void *,GEN,GEN))1580 gen_indexsort(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
1581 {
1582 long tx, lx;
1583 init_sort(&x, &tx, &lx);
1584 if (lx==1) return cgetg(1, t_VECSMALL);
1585 return gen_sortspec(x,lx-1,E,cmp);
1586 }
1587
1588 /* Sort the vector x in place, using cmp to compare entries */
1589 void
gen_sort_inplace(GEN x,void * E,int (* cmp)(void *,GEN,GEN),GEN * perm)1590 gen_sort_inplace(GEN x, void *E, int (*cmp)(void*,GEN,GEN), GEN *perm)
1591 {
1592 long tx, lx, i;
1593 pari_sp av = avma;
1594 GEN y;
1595
1596 init_sort(&x, &tx, &lx);
1597 if (lx<=2)
1598 {
1599 if (perm) *perm = lx == 1? cgetg(1, t_VECSMALL): mkvecsmall(1);
1600 return;
1601 }
1602 y = gen_sortspec(x,lx-1, E, cmp);
1603 if (perm)
1604 {
1605 GEN z = new_chunk(lx);
1606 for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
1607 for (i=1; i<lx; i++) gel(x,i) = gel(z,i);
1608 *perm = y;
1609 set_avma((pari_sp)y);
1610 } else {
1611 for (i=1; i<lx; i++) gel(y,i) = gel(x,y[i]);
1612 for (i=1; i<lx; i++) gel(x,i) = gel(y,i);
1613 set_avma(av);
1614 }
1615 }
1616 GEN
gen_sort_shallow(GEN x,void * E,int (* cmp)(void *,GEN,GEN))1617 gen_sort_shallow(GEN x, void *E, int (*cmp)(void*,GEN,GEN))
1618 {
1619 long tx, lx, i;
1620 pari_sp av;
1621 GEN y, z;
1622
1623 init_sort(&x, &tx, &lx);
1624 if (lx<=2) return x;
1625 z = cgetg(lx, tx); av = avma;
1626 y = gen_sortspec(x,lx-1, E, cmp);
1627 for (i=1; i<lx; i++) gel(z,i) = gel(x,y[i]);
1628 return gc_const(av, z);
1629 }
1630
1631 static int
closurecmp(void * data,GEN x,GEN y)1632 closurecmp(void *data, GEN x, GEN y)
1633 {
1634 pari_sp av = avma;
1635 long s = gsigne(closure_callgen2((GEN)data, x,y));
1636 set_avma(av); return s;
1637 }
1638 static void
check_positive_entries(GEN k)1639 check_positive_entries(GEN k)
1640 {
1641 long i, l = lg(k);
1642 for (i=1; i<l; i++)
1643 if (k[i] <= 0) pari_err_DOMAIN("sort_function", "index", "<", gen_0, stoi(k[i]));
1644 }
1645
1646 typedef int (*CMP_FUN)(void*,GEN,GEN);
1647 /* return NULL if t_CLOSURE k is a "key" (arity 1) and not a sorting func */
1648 static CMP_FUN
sort_function(void ** E,GEN x,GEN k)1649 sort_function(void **E, GEN x, GEN k)
1650 {
1651 int (*cmp)(GEN,GEN) = &lexcmp;
1652 long tx = typ(x);
1653 if (!k)
1654 {
1655 *E = (void*)((typ(x) == t_VECSMALL)? cmp_small: cmp);
1656 return &cmp_nodata;
1657 }
1658 if (tx == t_VECSMALL) pari_err_TYPE("sort_function", x);
1659 switch(typ(k))
1660 {
1661 case t_INT: k = mkvecsmall(itos(k)); break;
1662 case t_VEC: case t_COL: k = ZV_to_zv(k); break;
1663 case t_VECSMALL: break;
1664 case t_CLOSURE:
1665 if (closure_is_variadic(k))
1666 pari_err_TYPE("sort_function, variadic cmpf",k);
1667 *E = (void*)k;
1668 switch(closure_arity(k))
1669 {
1670 case 1: return NULL; /* wrt key */
1671 case 2: return &closurecmp;
1672 default: pari_err_TYPE("sort_function, cmpf arity != 1, 2",k);
1673 }
1674 default: pari_err_TYPE("sort_function",k);
1675 }
1676 check_positive_entries(k);
1677 *E = (void*)k; return &veccmp;
1678 }
1679
1680 #define cmp_IND 1
1681 #define cmp_LEX 2 /* FIXME: backward compatibility, ignored */
1682 #define cmp_REV 4
1683 #define cmp_UNIQ 8
1684 GEN
vecsort0(GEN x,GEN k,long flag)1685 vecsort0(GEN x, GEN k, long flag)
1686 {
1687 void *E;
1688 int (*CMP)(void*,GEN,GEN) = sort_function(&E, x, k);
1689
1690 if (flag < 0 || flag > (cmp_REV|cmp_LEX|cmp_IND|cmp_UNIQ))
1691 pari_err_FLAG("vecsort");
1692 if (!CMP)
1693 { /* wrt key: precompute all values, O(n) calls instead of O(n log n) */
1694 pari_sp av = avma;
1695 GEN v, y;
1696 long i, tx, lx;
1697 init_sort(&x, &tx, &lx);
1698 if (lx == 1) return flag&cmp_IND? cgetg(1,t_VECSMALL): triv_sort(tx);
1699 v = cgetg(lx, t_VEC);
1700 for (i = 1; i < lx; i++) gel(v,i) = closure_callgen1(k, gel(x,i));
1701 y = vecsort0(v, NULL, flag | cmp_IND);
1702 y = flag&cmp_IND? y: sort_extract(x, y, tx, lg(y));
1703 return gerepileupto(av, y);
1704 }
1705 if (flag&cmp_UNIQ)
1706 x = flag&cmp_IND? gen_indexsort_uniq(x, E, CMP): gen_sort_uniq(x, E, CMP);
1707 else
1708 x = flag&cmp_IND? gen_indexsort(x, E, CMP): gen_sort(x, E, CMP);
1709 if (flag & cmp_REV)
1710 { /* reverse order */
1711 GEN y = x;
1712 if (typ(x)==t_LIST) { y = list_data(x); if (!y) return x; }
1713 vecreverse_inplace(y);
1714 }
1715 return x;
1716 }
1717
1718 GEN
indexsort(GEN x)1719 indexsort(GEN x) { return gen_indexsort(x, (void*)&gcmp, cmp_nodata); }
1720 GEN
indexlexsort(GEN x)1721 indexlexsort(GEN x) { return gen_indexsort(x, (void*)&lexcmp, cmp_nodata); }
1722 GEN
indexvecsort(GEN x,GEN k)1723 indexvecsort(GEN x, GEN k)
1724 {
1725 if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
1726 return gen_indexsort(x, (void*)k, &veccmp);
1727 }
1728
1729 GEN
sort(GEN x)1730 sort(GEN x) { return gen_sort(x, (void*)gcmp, cmp_nodata); }
1731 GEN
lexsort(GEN x)1732 lexsort(GEN x) { return gen_sort(x, (void*)lexcmp, cmp_nodata); }
1733 GEN
vecsort(GEN x,GEN k)1734 vecsort(GEN x, GEN k)
1735 {
1736 if (typ(k) != t_VECSMALL) pari_err_TYPE("vecsort",k);
1737 return gen_sort(x, (void*)k, &veccmp);
1738 }
1739 /* adapted from gen_search; don't export: keys of T[i] should be precomputed */
1740 static long
key_search(GEN T,GEN x,GEN code)1741 key_search(GEN T, GEN x, GEN code)
1742 {
1743 long u = lg(T)-1, i, l, s;
1744
1745 if (!u) return 0;
1746 l = 1; x = closure_callgen1(code, x);
1747 do
1748 {
1749 i = (l+u)>>1; s = lexcmp(x, closure_callgen1(code, gel(T,i)));
1750 if (!s) return i;
1751 if (s<0) u=i-1; else l=i+1;
1752 } while (u>=l);
1753 return 0;
1754 }
1755 long
vecsearch(GEN v,GEN x,GEN k)1756 vecsearch(GEN v, GEN x, GEN k)
1757 {
1758 pari_sp av = avma;
1759 void *E;
1760 int (*CMP)(void*,GEN,GEN) = sort_function(&E, v, k);
1761 long r, tv = typ(v);
1762 if (tv == t_VECSMALL)
1763 x = (GEN)itos(x);
1764 else if (!is_matvec_t(tv)) pari_err_TYPE("vecsearch", v);
1765 r = CMP? gen_search(v, x, 0, E, CMP): key_search(v, x, k);
1766 return gc_long(av, r);
1767 }
1768
1769 GEN
ZV_indexsort(GEN L)1770 ZV_indexsort(GEN L) { return gen_indexsort(L, (void*)&cmpii, &cmp_nodata); }
1771 GEN
ZV_sort(GEN L)1772 ZV_sort(GEN L) { return gen_sort(L, (void*)&cmpii, &cmp_nodata); }
1773 GEN
ZV_sort_uniq(GEN L)1774 ZV_sort_uniq(GEN L) { return gen_sort_uniq(L, (void*)&cmpii, &cmp_nodata); }
1775 void
ZV_sort_inplace(GEN L)1776 ZV_sort_inplace(GEN L) { gen_sort_inplace(L, (void*)&cmpii, &cmp_nodata,NULL); }
1777
1778 GEN
vec_equiv(GEN F)1779 vec_equiv(GEN F)
1780 {
1781 pari_sp av = avma;
1782 long j, k, L = lg(F);
1783 GEN w = cgetg(L, t_VEC);
1784 GEN perm = gen_indexsort(F, (void*)&cmp_universal, cmp_nodata);
1785 for (j = k = 1; j < L;)
1786 {
1787 GEN v = cgetg(L, t_VECSMALL);
1788 long l = 1, o = perm[j];
1789 v[l++] = o;
1790 for (j++; j < L; v[l++] = perm[j++])
1791 if (!gequal(gel(F,o), gel(F, perm[j]))) break;
1792 setlg(v, l); gel(w, k++) = v;
1793 }
1794 setlg(w, k); return gerepilecopy(av,w);
1795 }
1796
1797 GEN
vec_reduce(GEN v,GEN * pE)1798 vec_reduce(GEN v, GEN *pE)
1799 {
1800 GEN E, F, P = gen_indexsort(v, (void*)cmp_universal, cmp_nodata);
1801 long i, m, l;
1802 F = cgetg_copy(v, &l);
1803 *pE = E = cgetg(l, t_VECSMALL);
1804 for (i = m = 1; i < l;)
1805 {
1806 GEN u = gel(v, P[i]);
1807 long k;
1808 for(k = i + 1; k < l; k++)
1809 if (cmp_universal(gel(v, P[k]), u)) break;
1810 E[m] = k - i; gel(F, m) = u; i = k; m++;
1811 }
1812 setlg(F, m);
1813 setlg(E, m); return F;
1814 }
1815
1816 /********************************************************************/
1817 /** SEARCH IN SORTED VECTOR **/
1818 /********************************************************************/
1819 /* index of x in table T, 0 otherwise */
1820 long
tablesearch(GEN T,GEN x,int (* cmp)(GEN,GEN))1821 tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
1822 {
1823 long l = 1, u = lg(T)-1, i, s;
1824
1825 while (u>=l)
1826 {
1827 i = (l+u)>>1; s = cmp(x, gel(T,i));
1828 if (!s) return i;
1829 if (s<0) u=i-1; else l=i+1;
1830 }
1831 return 0;
1832 }
1833
1834 /* looks if x belongs to the set T and returns the index if yes, 0 if no */
1835 long
gen_search(GEN T,GEN x,long flag,void * data,int (* cmp)(void *,GEN,GEN))1836 gen_search(GEN T, GEN x, long flag, void *data, int (*cmp)(void*,GEN,GEN))
1837 {
1838 long u = lg(T)-1, i, l, s;
1839
1840 if (!u) return flag? 1: 0;
1841 l = 1;
1842 do
1843 {
1844 i = (l+u)>>1; s = cmp(data, x, gel(T,i));
1845 if (!s) return flag? 0: i;
1846 if (s<0) u=i-1; else l=i+1;
1847 } while (u>=l);
1848 if (!flag) return 0;
1849 return (s<0)? i: i+1;
1850 }
1851
1852 long
ZV_search(GEN x,GEN y)1853 ZV_search(GEN x, GEN y) { return tablesearch(x, y, cmpii); }
1854
1855 long
zv_search(GEN T,long x)1856 zv_search(GEN T, long x)
1857 {
1858 long l = 1, u = lg(T)-1;
1859 while (u>=l)
1860 {
1861 long i = (l+u)>>1;
1862 if (x < T[i]) u = i-1;
1863 else if (x > T[i]) l = i+1;
1864 else return i;
1865 }
1866 return 0;
1867 }
1868
1869 /********************************************************************/
1870 /** COMPARISON FUNCTIONS **/
1871 /********************************************************************/
1872 int
cmp_nodata(void * data,GEN x,GEN y)1873 cmp_nodata(void *data, GEN x, GEN y)
1874 {
1875 int (*cmp)(GEN,GEN)=(int (*)(GEN,GEN)) data;
1876 return cmp(x,y);
1877 }
1878
1879 /* assume x and y come from the same idealprimedec call (uniformizer unique) */
1880 int
cmp_prime_over_p(GEN x,GEN y)1881 cmp_prime_over_p(GEN x, GEN y)
1882 {
1883 long k = pr_get_f(x) - pr_get_f(y); /* diff. between residue degree */
1884 return k? ((k > 0)? 1: -1)
1885 : ZV_cmp(pr_get_gen(x), pr_get_gen(y));
1886 }
1887
1888 int
cmp_prime_ideal(GEN x,GEN y)1889 cmp_prime_ideal(GEN x, GEN y)
1890 {
1891 int k = cmpii(pr_get_p(x), pr_get_p(y));
1892 return k? k: cmp_prime_over_p(x,y);
1893 }
1894
1895 /* assume x and y are t_POL in the same variable whose coeffs can be
1896 * compared (used to sort polynomial factorizations) */
1897 int
gen_cmp_RgX(void * data,GEN x,GEN y)1898 gen_cmp_RgX(void *data, GEN x, GEN y)
1899 {
1900 int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
1901 long i, lx = lg(x), ly = lg(y);
1902 int fl;
1903 if (lx > ly) return 1;
1904 if (lx < ly) return -1;
1905 for (i=lx-1; i>1; i--)
1906 if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
1907 return 0;
1908 }
1909
1910 static int
cmp_RgX_Rg(GEN x,GEN y)1911 cmp_RgX_Rg(GEN x, GEN y)
1912 {
1913 long lx = lgpol(x), ly;
1914 if (lx > 1) return 1;
1915 ly = gequal0(y) ? 0:1;
1916 if (lx > ly) return 1;
1917 if (lx < ly) return -1;
1918 if (lx==0) return 0;
1919 return gcmp(gel(x,2), y);
1920 }
1921 int
cmp_RgX(GEN x,GEN y)1922 cmp_RgX(GEN x, GEN y)
1923 {
1924 if (typ(x) == t_POLMOD) x = gel(x,2);
1925 if (typ(y) == t_POLMOD) y = gel(y,2);
1926 if (typ(x) == t_POL) {
1927 if (typ(y) != t_POL) return cmp_RgX_Rg(x, y);
1928 } else {
1929 if (typ(y) != t_POL) return gcmp(x,y);
1930 return - cmp_RgX_Rg(y,x);
1931 }
1932 return gen_cmp_RgX((void*)&gcmp,x,y);
1933 }
1934
1935 int
cmp_Flx(GEN x,GEN y)1936 cmp_Flx(GEN x, GEN y)
1937 {
1938 long i, lx = lg(x), ly = lg(y);
1939 if (lx > ly) return 1;
1940 if (lx < ly) return -1;
1941 for (i=lx-1; i>1; i--)
1942 if (uel(x,i) != uel(y,i)) return uel(x,i)<uel(y,i)? -1: 1;
1943 return 0;
1944 }
1945 /********************************************************************/
1946 /** MERGE & SORT FACTORIZATIONS **/
1947 /********************************************************************/
1948 /* merge fx, fy two factorizations, whose 1st column is sorted in strictly
1949 * increasing order wrt cmp */
1950 GEN
merge_factor(GEN fx,GEN fy,void * data,int (* cmp)(void *,GEN,GEN))1951 merge_factor(GEN fx, GEN fy, void *data, int (*cmp)(void *,GEN,GEN))
1952 {
1953 GEN x = gel(fx,1), e = gel(fx,2), M, E;
1954 GEN y = gel(fy,1), f = gel(fy,2);
1955 long ix, iy, m, lx = lg(x), ly = lg(y), l = lx+ly-1;
1956
1957 M = cgetg(l, t_COL);
1958 E = cgetg(l, t_COL);
1959
1960 m = ix = iy = 1;
1961 while (ix<lx && iy<ly)
1962 {
1963 int s = cmp(data, gel(x,ix), gel(y,iy));
1964 if (s < 0)
1965 { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; }
1966 else if (s == 0)
1967 {
1968 GEN z = gel(x,ix), g = addii(gel(e,ix), gel(f,iy));
1969 iy++; ix++; if (!signe(g)) continue;
1970 gel(M,m) = z; gel(E,m) = g;
1971 }
1972 else
1973 { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; }
1974 m++;
1975 }
1976 while (ix<lx) { gel(M,m) = gel(x,ix); gel(E,m) = gel(e,ix); ix++; m++; }
1977 while (iy<ly) { gel(M,m) = gel(y,iy); gel(E,m) = gel(f,iy); iy++; m++; }
1978 setlg(M, m);
1979 setlg(E, m); return mkmat2(M, E);
1980 }
1981 /* merge two sorted vectors, removing duplicates. Shallow */
1982 GEN
merge_sort_uniq(GEN x,GEN y,void * data,int (* cmp)(void *,GEN,GEN))1983 merge_sort_uniq(GEN x, GEN y, void *data, int (*cmp)(void *,GEN,GEN))
1984 {
1985 long i, j, k, lx = lg(x), ly = lg(y);
1986 GEN z = cgetg(lx + ly - 1, typ(x));
1987 i = j = k = 1;
1988 while (i<lx && j<ly)
1989 {
1990 int s = cmp(data, gel(x,i), gel(y,j));
1991 if (s < 0)
1992 gel(z,k++) = gel(x,i++);
1993 else if (s > 0)
1994 gel(z,k++) = gel(y,j++);
1995 else
1996 { gel(z,k++) = gel(x,i++); j++; }
1997 }
1998 while (i<lx) gel(z,k++) = gel(x,i++);
1999 while (j<ly) gel(z,k++) = gel(y,j++);
2000 setlg(z, k); return z;
2001 }
2002 /* in case of equal keys in x,y, take the key from x */
2003 static GEN
ZV_union_shallow_t(GEN x,GEN y,long t)2004 ZV_union_shallow_t(GEN x, GEN y, long t)
2005 {
2006 long i, j, k, lx = lg(x), ly = lg(y);
2007 GEN z = cgetg(lx + ly - 1, t);
2008 i = j = k = 1;
2009 while (i<lx && j<ly)
2010 {
2011 int s = cmpii(gel(x,i), gel(y,j));
2012 if (s < 0)
2013 gel(z,k++) = gel(x,i++);
2014 else if (s > 0)
2015 gel(z,k++) = gel(y,j++);
2016 else
2017 { gel(z,k++) = gel(x,i++); j++; }
2018 }
2019 while (i < lx) gel(z,k++) = gel(x,i++);
2020 while (j < ly) gel(z,k++) = gel(y,j++);
2021 setlg(z, k); return z;
2022 }
2023 GEN
ZV_union_shallow(GEN x,GEN y)2024 ZV_union_shallow(GEN x, GEN y)
2025 { return ZV_union_shallow_t(x, y, t_VEC); }
2026 GEN
ZC_union_shallow(GEN x,GEN y)2027 ZC_union_shallow(GEN x, GEN y)
2028 { return ZV_union_shallow_t(x, y, t_COL); }
2029
2030 /* sort generic factorization, in place */
2031 GEN
sort_factor(GEN y,void * data,int (* cmp)(void *,GEN,GEN))2032 sort_factor(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
2033 {
2034 GEN a, b, A, B, w;
2035 pari_sp av;
2036 long n, i;
2037
2038 a = gel(y,1); n = lg(a); if (n == 1) return y;
2039 b = gel(y,2); av = avma;
2040 A = new_chunk(n);
2041 B = new_chunk(n);
2042 w = gen_sortspec(a, n-1, data, cmp);
2043 for (i=1; i<n; i++) { long k=w[i]; gel(A,i) = gel(a,k); gel(B,i) = gel(b,k); }
2044 for (i=1; i<n; i++) { gel(a,i) = gel(A,i); gel(b,i) = gel(B,i); }
2045 set_avma(av); return y;
2046 }
2047 /* sort polynomial factorization, in place */
2048 GEN
sort_factor_pol(GEN y,int (* cmp)(GEN,GEN))2049 sort_factor_pol(GEN y,int (*cmp)(GEN,GEN))
2050 {
2051 (void)sort_factor(y,(void*)cmp, &gen_cmp_RgX);
2052 return y;
2053 }
2054
2055 /***********************************************************************/
2056 /* */
2057 /* SET OPERATIONS */
2058 /* */
2059 /***********************************************************************/
2060 GEN
gtoset(GEN x)2061 gtoset(GEN x)
2062 {
2063 long lx;
2064 if (!x) return cgetg(1, t_VEC);
2065 switch(typ(x))
2066 {
2067 case t_VEC:
2068 case t_COL: lx = lg(x); break;
2069 case t_LIST:
2070 if (list_typ(x)==t_LIST_MAP) return mapdomain(x);
2071 x = list_data(x); lx = x? lg(x): 1; break;
2072 case t_VECSMALL: lx = lg(x); x = zv_to_ZV(x); break;
2073 default: return mkveccopy(x);
2074 }
2075 if (lx==1) return cgetg(1,t_VEC);
2076 x = gen_sort_uniq(x, (void*)&cmp_universal, cmp_nodata);
2077 settyp(x, t_VEC); /* it may be t_COL */
2078 return x;
2079 }
2080
2081 long
setisset(GEN x)2082 setisset(GEN x)
2083 {
2084 long i, lx = lg(x);
2085
2086 if (typ(x) != t_VEC) return 0;
2087 if (lx == 1) return 1;
2088 for (i=1; i<lx-1; i++)
2089 if (cmp_universal(gel(x,i+1), gel(x,i)) <= 0) return 0;
2090 return 1;
2091 }
2092
2093 long
setsearch(GEN T,GEN y,long flag)2094 setsearch(GEN T, GEN y, long flag)
2095 {
2096 long lx;
2097 switch(typ(T))
2098 {
2099 case t_VEC: lx = lg(T); break;
2100 case t_LIST:
2101 if (list_typ(T) != t_LIST_RAW) pari_err_TYPE("setsearch",T);
2102 T = list_data(T); lx = T? lg(T): 1; break;
2103 default: pari_err_TYPE("setsearch",T);
2104 return 0; /*LCOV_EXCL_LINE*/
2105 }
2106 if (lx==1) return flag? 1: 0;
2107 return gen_search(T,y,flag,(void*)cmp_universal,cmp_nodata);
2108 }
2109
2110 GEN
setunion_i(GEN x,GEN y)2111 setunion_i(GEN x, GEN y)
2112 { return merge_sort_uniq(x,y, (void*)cmp_universal, cmp_nodata); }
2113
2114 GEN
setunion(GEN x,GEN y)2115 setunion(GEN x, GEN y)
2116 {
2117 pari_sp av = avma;
2118 if (typ(x) != t_VEC) pari_err_TYPE("setunion",x);
2119 if (typ(y) != t_VEC) pari_err_TYPE("setunion",y);
2120 return gerepilecopy(av, setunion_i(x, y));
2121 }
2122
2123 GEN
setintersect(GEN x,GEN y)2124 setintersect(GEN x, GEN y)
2125 {
2126 long ix = 1, iy = 1, iz = 1, lx = lg(x), ly = lg(y);
2127 pari_sp av = avma;
2128 GEN z = cgetg(lx,t_VEC);
2129 if (typ(x) != t_VEC) pari_err_TYPE("setintersect",x);
2130 if (typ(y) != t_VEC) pari_err_TYPE("setintersect",y);
2131 while (ix < lx && iy < ly)
2132 {
2133 int c = cmp_universal(gel(x,ix), gel(y,iy));
2134 if (c < 0) ix++;
2135 else if (c > 0) iy++;
2136 else { gel(z, iz++) = gel(x,ix); ix++; iy++; }
2137 }
2138 setlg(z,iz); return gerepilecopy(av,z);
2139 }
2140
2141 GEN
gen_setminus(GEN A,GEN B,int (* cmp)(GEN,GEN))2142 gen_setminus(GEN A, GEN B, int (*cmp)(GEN,GEN))
2143 {
2144 pari_sp ltop = avma;
2145 long i = 1, j = 1, k = 1, lx = lg(A), ly = lg(B);
2146 GEN diff = cgetg(lx,t_VEC);
2147 while (i < lx && j < ly)
2148 switch ( cmp(gel(A,i),gel(B,j)) )
2149 {
2150 case -1: gel(diff,k++) = gel(A,i++); break;
2151 case 1: j++; break;
2152 case 0: i++; break;
2153 }
2154 while (i < lx) gel(diff,k++) = gel(A,i++);
2155 setlg(diff,k);
2156 return gerepilecopy(ltop,diff);
2157 }
2158
2159 GEN
setminus(GEN x,GEN y)2160 setminus(GEN x, GEN y)
2161 {
2162 if (typ(x) != t_VEC) pari_err_TYPE("setminus",x);
2163 if (typ(y) != t_VEC) pari_err_TYPE("setminus",y);
2164 return gen_setminus(x,y,cmp_universal);
2165 }
2166
2167 GEN
setbinop(GEN f,GEN x,GEN y)2168 setbinop(GEN f, GEN x, GEN y)
2169 {
2170 pari_sp av = avma;
2171 long i, j, lx, ly, k = 1;
2172 GEN z;
2173 if (typ(f) != t_CLOSURE || closure_arity(f) != 2 || closure_is_variadic(f))
2174 pari_err_TYPE("setbinop [function needs exactly 2 arguments]",f);
2175 lx = lg(x);
2176 if (typ(x) != t_VEC) pari_err_TYPE("setbinop", x);
2177 if (y == NULL) { /* assume x = y and f symmetric */
2178 z = cgetg((((lx-1)*lx) >> 1) + 1, t_VEC);
2179 for (i = 1; i < lx; i++)
2180 for (j = i; j < lx; j++)
2181 gel(z, k++) = closure_callgen2(f, gel(x,i),gel(x,j));
2182 } else {
2183 ly = lg(y);
2184 if (typ(y) != t_VEC) pari_err_TYPE("setbinop", y);
2185 z = cgetg((lx-1)*(ly-1) + 1, t_VEC);
2186 for (i = 1; i < lx; i++)
2187 for (j = 1; j < ly; j++)
2188 gel(z, k++) = closure_callgen2(f, gel(x,i),gel(y,j));
2189 }
2190 return gerepileupto(av, gtoset(z));
2191 }
2192