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