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