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 : /*********************************************************************/
16 : /** ARITHMETIC FUNCTIONS **/
17 : /** (second part) **/
18 : /*********************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : #define DEBUGLEVEL DEBUGLEVEL_arith
23 :
24 : GEN
25 84 : boundfact(GEN n, ulong lim)
26 : {
27 84 : switch(typ(n))
28 : {
29 70 : case t_INT: return Z_factor_limit(n,lim);
30 14 : case t_FRAC: {
31 14 : pari_sp av = avma;
32 14 : GEN a = Z_factor_limit(gel(n,1),lim);
33 14 : GEN b = Z_factor_limit(gel(n,2),lim);
34 14 : gel(b,2) = ZC_neg(gel(b,2));
35 14 : return gerepilecopy(av, ZM_merge_factor(a,b));
36 : }
37 : }
38 0 : pari_err_TYPE("boundfact",n);
39 : return NULL; /* LCOV_EXCL_LINE */
40 : }
41 :
42 : /* NOT memory clean */
43 : GEN
44 19302 : Z_lsmoothen(GEN N, GEN L, GEN *pP, GEN *pE)
45 : {
46 19302 : long i, j, l = lg(L);
47 19302 : GEN E = new_chunk(l), P = new_chunk(l);
48 181040 : for (i = j = 1; i < l; i++)
49 : {
50 172632 : ulong p = uel(L,i);
51 172632 : long v = Z_lvalrem(N, p, &N);
52 172632 : if (v) { P[j] = p; E[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
53 : }
54 19302 : P[0] = evaltyp(t_VECSMALL) | _evallg(j); if (pP) *pP = P;
55 19302 : E[0] = evaltyp(t_VECSMALL) | _evallg(j); if (pE) *pE = E; return N;
56 : }
57 : GEN
58 82157 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pE)
59 : {
60 82157 : long i, j, l = lg(L);
61 82157 : GEN E = new_chunk(l), P = new_chunk(l);
62 289492 : for (i = j = 1; i < l; i++)
63 : {
64 264340 : GEN p = gel(L,i);
65 264340 : long v = Z_pvalrem(N, p, &N);
66 264398 : if (v)
67 : {
68 195657 : gel(P,j) = p; gel(E,j) = utoipos(v); j++;
69 195658 : if (is_pm1(N)) { N = NULL; break; }
70 : }
71 : }
72 82235 : P[0] = evaltyp(t_COL) | _evallg(j); if (pP) *pP = P;
73 82235 : E[0] = evaltyp(t_COL) | _evallg(j); if (pE) *pE = E; return N;
74 : }
75 : /***********************************************************************/
76 : /** SIMPLE FACTORISATIONS **/
77 : /***********************************************************************/
78 : /* Factor n and output [p,e,c] where
79 : * p, e and c are vecsmall with n = prod{p[i]^e[i]} and c[i] = p[i]^e[i] */
80 : GEN
81 141222 : factoru_pow(ulong n)
82 : {
83 141222 : GEN f = cgetg(4,t_VEC);
84 141222 : pari_sp av = avma;
85 : GEN F, P, E, p, e, c;
86 : long i, l;
87 : /* enough room to store <= 15 * [p,e,p^e] (OK if n < 2^64) */
88 141222 : (void)new_chunk((15 + 1)*3);
89 141222 : F = factoru(n);
90 141221 : P = gel(F,1);
91 141221 : E = gel(F,2); l = lg(P);
92 141221 : set_avma(av);
93 141221 : gel(f,1) = p = cgetg(l,t_VECSMALL);
94 141221 : gel(f,2) = e = cgetg(l,t_VECSMALL);
95 141221 : gel(f,3) = c = cgetg(l,t_VECSMALL);
96 315446 : for(i = 1; i < l; i++)
97 : {
98 174225 : p[i] = P[i];
99 174225 : e[i] = E[i];
100 174225 : c[i] = upowuu(p[i], e[i]);
101 : }
102 141221 : return f;
103 : }
104 :
105 : static GEN
106 16443 : factorlim(GEN n, ulong lim)
107 16443 : { return lim? Z_factor_limit(n, lim): Z_factor(n); }
108 : /* factor p^n - 1, assuming p prime. If lim != 0, limit factorization to
109 : * primes <= lim */
110 : GEN
111 7476 : factor_pn_1_limit(GEN p, long n, ulong lim)
112 : {
113 7476 : pari_sp av = avma;
114 7476 : GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
115 7476 : long i, pp = itos_or_0(p);
116 16100 : for(i=2; i<lg(d); i++)
117 : {
118 : GEN B;
119 8624 : if (pp && d[i]%pp==0 && (
120 3479 : ((pp&3)==1 && (d[i]&1)) ||
121 3458 : ((pp&3)==3 && (d[i]&3)==2) ||
122 3339 : (pp==2 && (d[i]&7)==4)))
123 343 : {
124 343 : GEN f=factor_Aurifeuille_prime(p,d[i]);
125 343 : B = factorlim(f, lim);
126 343 : A = ZM_merge_factor(A, B);
127 343 : B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
128 : }
129 : else
130 8281 : B = factorlim(polcyclo_eval(d[i],p), lim);
131 8624 : A = ZM_merge_factor(A, B);
132 : }
133 7476 : return gerepilecopy(av, A);
134 : }
135 : GEN
136 7476 : factor_pn_1(GEN p, ulong n)
137 7476 : { return factor_pn_1_limit(p, n, 0); }
138 :
139 : #if 0
140 : static GEN
141 : to_mat(GEN p, long e) {
142 : GEN B = cgetg(3, t_MAT);
143 : gel(B,1) = mkcol(p);
144 : gel(B,2) = mkcol(utoipos(e)); return B;
145 : }
146 : /* factor phi(n) */
147 : GEN
148 : factor_eulerphi(GEN n)
149 : {
150 : pari_sp av = avma;
151 : GEN B = NULL, A, P, E, AP, AE;
152 : long i, l, v = vali(n);
153 :
154 : l = lg(n);
155 : /* result requires less than this: at most expi(n) primes */
156 : (void)new_chunk(bit_accuracy(l) * (l /*p*/ + 3 /*e*/ + 2 /*vectors*/) + 3+2);
157 : if (v) { n = shifti(n, -v); v--; }
158 : A = Z_factor(n); P = gel(A,1); E = gel(A,2); l = lg(P);
159 : for(i = 1; i < l; i++)
160 : {
161 : GEN p = gel(P,i), q = subiu(p,1), fa;
162 : long e = itos(gel(E,i)), w;
163 :
164 : w = vali(q); v += w; q = shifti(q,-w);
165 : if (! is_pm1(q))
166 : {
167 : fa = Z_factor(q);
168 : B = B? ZM_merge_factor(B, fa): fa;
169 : }
170 : if (e > 1) {
171 : if (B) {
172 : gel(B,1) = vec_append(gel(B,1), p);
173 : gel(B,2) = vec_append(gel(B,2), utoipos(e-1));
174 : } else
175 : B = to_mat(p, e-1);
176 : }
177 : }
178 : set_avma(av);
179 : if (!B) return v? to_mat(gen_2, v): trivial_fact();
180 : A = cgetg(3, t_MAT);
181 : P = gel(B,1); E = gel(B,2); l = lg(P);
182 : AP = cgetg(l+1, t_COL); gel(A,1) = AP; AP++;
183 : AE = cgetg(l+1, t_COL); gel(A,2) = AE; AE++;
184 : /* prepend "2^v" */
185 : gel(AP,0) = gen_2;
186 : gel(AE,0) = utoipos(v);
187 : for (i = 1; i < l; i++)
188 : {
189 : gel(AP,i) = icopy(gel(P,i));
190 : gel(AE,i) = icopy(gel(E,i));
191 : }
192 : return A;
193 : }
194 : #endif
195 :
196 : /***********************************************************************/
197 : /** CHECK FACTORIZATION FOR ARITHMETIC FUNCTIONS **/
198 : /***********************************************************************/
199 : int
200 15636400 : RgV_is_ZVpos(GEN v)
201 : {
202 15636400 : long i, l = lg(v);
203 45496445 : for (i = 1; i < l; i++)
204 : {
205 29870188 : GEN c = gel(v,i);
206 29870188 : if (typ(c) != t_INT || signe(c) <= 0) return 0;
207 : }
208 15626257 : return 1;
209 : }
210 : /* check whether v is a ZV with nonzero entries */
211 : int
212 22868 : RgV_is_ZVnon0(GEN v)
213 : {
214 22868 : long i, l = lg(v);
215 164394 : for (i = 1; i < l; i++)
216 : {
217 141582 : GEN c = gel(v,i);
218 141582 : if (typ(c) != t_INT || !signe(c)) return 0;
219 : }
220 22812 : return 1;
221 : }
222 : /* check whether v is a ZV with nonzero entries OR exactly [0] */
223 : static int
224 707684 : RgV_is_ZV0(GEN v)
225 : {
226 707684 : long i, l = lg(v);
227 2429932 : for (i = 1; i < l; i++)
228 : {
229 1722388 : GEN c = gel(v,i);
230 : long s;
231 1722388 : if (typ(c) != t_INT) return 0;
232 1722388 : s = signe(c);
233 1722388 : if (!s) return (l == 2);
234 : }
235 707544 : return 1;
236 : }
237 :
238 : int
239 302920 : RgV_is_prV(GEN v)
240 : {
241 302920 : long l = lg(v), i;
242 658622 : for (i = 1; i < l; i++)
243 550840 : if (!checkprid_i(gel(v,i))) return 0;
244 107782 : return 1;
245 : }
246 : int
247 260939 : is_nf_factor(GEN F)
248 : {
249 242109 : return typ(F) == t_MAT && lg(F) == 3
250 503048 : && RgV_is_prV(gel(F,1)) && RgV_is_ZVpos(gel(F,2));
251 : }
252 : int
253 14 : is_nf_extfactor(GEN F)
254 : {
255 14 : return typ(F) == t_MAT && lg(F) == 3
256 28 : && RgV_is_prV(gel(F,1)) && RgV_is_ZV(gel(F,2));
257 : }
258 :
259 : static int
260 8152721 : is_Z_factor_i(GEN f)
261 8152721 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
262 : int
263 7430578 : is_Z_factorpos(GEN f)
264 7430578 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
265 : int
266 707686 : is_Z_factor(GEN f)
267 707686 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
268 : /* as is_Z_factorpos, also allow factor(0) */
269 : int
270 14462 : is_Z_factornon0(GEN f)
271 14462 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
272 : GEN
273 26376 : clean_Z_factor(GEN f)
274 : {
275 26376 : GEN P = gel(f,1);
276 26376 : long n = lg(P)-1;
277 26376 : if (n && equalim1(gel(P,1)))
278 3066 : return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
279 23310 : return f;
280 : }
281 : GEN
282 0 : fuse_Z_factor(GEN f, GEN B)
283 : {
284 0 : GEN P = gel(f,1), E = gel(f,2), P2,E2;
285 0 : long i, l = lg(P);
286 0 : if (l == 1) return f;
287 0 : for (i = 1; i < l; i++)
288 0 : if (abscmpii(gel(P,i), B) > 0) break;
289 0 : if (i == l) return f;
290 : /* tail / initial segment */
291 0 : P2 = vecslice(P, i, l-1); P = vecslice(P, 1, i-1);
292 0 : E2 = vecslice(E, i, l-1); E = vecslice(E, 1, i-1);
293 0 : P = vec_append(P, factorback2(P2,E2));
294 0 : E = vec_append(E, gen_1);
295 0 : return mkmat2(P, E);
296 : }
297 :
298 : /* n attached to a factorization of a positive integer: either N (t_INT)
299 : * a factorization matrix faN, or a t_VEC: [N, faN] */
300 : GEN
301 630 : check_arith_pos(GEN n, const char *f) {
302 630 : switch(typ(n))
303 : {
304 602 : case t_INT:
305 602 : if (signe(n) <= 0) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
306 602 : return NULL;
307 7 : case t_VEC:
308 7 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT || signe(gel(n,1)) <= 0) break;
309 7 : n = gel(n,2); /* fall through */
310 21 : case t_MAT:
311 21 : if (!is_Z_factorpos(n)) break;
312 21 : return n;
313 : }
314 7 : pari_err_TYPE(f,n);
315 : return NULL;/*LCOV_EXCL_LINE*/
316 : }
317 : /* n attached to a factorization of a nonzero integer */
318 : GEN
319 2621052 : check_arith_non0(GEN n, const char *f) {
320 2621052 : switch(typ(n))
321 : {
322 2606835 : case t_INT:
323 2606835 : if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
324 2606785 : return NULL;
325 1232 : case t_VEC:
326 1232 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT || !signe(gel(n,1))) break;
327 1232 : n = gel(n,2); /* fall through */
328 14231 : case t_MAT:
329 14231 : if (!is_Z_factornon0(n)) break;
330 14175 : return n;
331 : }
332 48 : pari_err_TYPE(f,n);
333 : return NULL;/*LCOV_EXCL_LINE*/
334 : }
335 : /* n attached to a factorization of an integer */
336 : GEN
337 2898622 : check_arith_all(GEN n, const char *f) {
338 2898622 : switch(typ(n))
339 : {
340 2190927 : case t_INT:
341 2190927 : return NULL;
342 687395 : case t_VEC:
343 687395 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
344 687393 : n = gel(n,2); /* fall through */
345 707686 : case t_MAT:
346 707686 : if (!is_Z_factor(n)) break;
347 707684 : return n;
348 : }
349 9 : pari_err_TYPE(f,n);
350 : return NULL;/*LCOV_EXCL_LINE*/
351 : }
352 :
353 : /***********************************************************************/
354 : /** MISCELLANEOUS ARITHMETIC FUNCTIONS **/
355 : /** (ultimately depend on Z_factor()) **/
356 : /***********************************************************************/
357 : /* set P,E from F. Check whether F is an integer and kill "factor" -1 */
358 : static void
359 362649 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
360 : {
361 : GEN E, P;
362 362649 : if (lg(F) != 3) pari_err_TYPE("divisors",F);
363 362649 : P = gel(F,1);
364 362649 : E = gel(F,2);
365 362649 : RgV_check_ZV(E, "divisors");
366 362649 : *isint = RgV_is_ZV(P);
367 362649 : if (*isint)
368 : {
369 362635 : long i, l = lg(P);
370 : /* skip -1 */
371 362635 : if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
372 : /* test for 0 */
373 1174964 : for (i = 1; i < l; i++)
374 812336 : if (!signe(gel(P,i)) && signe(gel(E,i)))
375 7 : pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
376 : }
377 362642 : *pP = P;
378 362642 : *pE = E;
379 362642 : }
380 : static void
381 16247 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
382 :
383 : int
384 378903 : divisors_init(GEN n, GEN *pP, GEN *pE)
385 : {
386 : long i,l;
387 : GEN E, P, e;
388 : int isint;
389 :
390 378903 : switch(typ(n))
391 : {
392 16226 : case t_INT:
393 16226 : if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
394 16226 : set_fact(absZ_factor(n), &P,&E);
395 16226 : isint = 1; break;
396 1666 : case t_VEC:
397 1666 : if (lg(n) != 3 || typ(gel(n,2)) !=t_MAT) pari_err_TYPE("divisors",n);
398 1659 : set_fact_check(gel(n,2), &P,&E, &isint);
399 1659 : break;
400 360990 : case t_MAT:
401 360990 : set_fact_check(n, &P,&E, &isint);
402 360983 : break;
403 21 : default:
404 21 : set_fact(factor(n), &P,&E);
405 21 : isint = 0; break;
406 : }
407 378889 : l = lg(P);
408 378889 : e = cgetg(l, t_VECSMALL);
409 1232119 : for (i=1; i<l; i++)
410 : {
411 853237 : e[i] = itos(gel(E,i));
412 853237 : if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
413 : }
414 378882 : *pP = P; *pE = e; return isint;
415 : }
416 :
417 : static long
418 379183 : ndiv(GEN E)
419 : {
420 : long i, l;
421 379183 : GEN e = cgetg_copy(E, &l); /* left on stack */
422 : ulong n;
423 1232812 : for (i=1; i<l; i++) e[i] = E[i]+1;
424 379183 : n = itou_or_0( zv_prod_Z(e) );
425 379183 : if (!n || n & ~LGBITS) pari_err_OVERFLOW("divisors");
426 379183 : return n;
427 : }
428 : static int
429 16646 : cmpi1(void *E, GEN a, GEN b) { (void)E; return cmpii(gel(a,1), gel(b,1)); }
430 : /* P a t_COL of objects, E a t_VECSMALL of exponents, return cleaned-up
431 : * factorization (removing 0 exponents) as a t_MAT with 2 cols. */
432 : static GEN
433 24171 : fa_clean(GEN P, GEN E)
434 : {
435 24171 : long i, j, l = lg(E);
436 24171 : GEN Q = cgetg(l, t_COL);
437 68054 : for (i = j = 1; i < l; i++)
438 43883 : if (E[i]) { gel(Q,j) = gel(P,i); E[j] = E[i]; j++; }
439 24171 : setlg(Q,j);
440 24171 : setlg(E,j); return mkmat2(Q,Flc_to_ZC(E));
441 : }
442 : GEN
443 12453 : divisors_factored(GEN N)
444 : {
445 12453 : pari_sp av = avma;
446 : GEN *d, *t1, *t2, *t3, D, P, E;
447 12453 : int isint = divisors_init(N, &P, &E);
448 12453 : GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
449 12453 : long i, j, l, n = ndiv(E);
450 :
451 12453 : D = cgetg(n+1,t_VEC); d = (GEN*)D;
452 12453 : l = lg(E);
453 12453 : *++d = mkvec2(gen_1, const_vecsmall(l-1,0));
454 33789 : for (i=1; i<l; i++)
455 30618 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
456 21000 : for (t2=d, t3=t1; t3<t2; )
457 : {
458 : GEN a, b;
459 11718 : a = mul(gel(*++t3,1), gel(P,i));
460 11718 : b = leafcopy(gel(*t3,2)); b[i]++;
461 11718 : *++d = mkvec2(a,b);
462 : }
463 12453 : if (isint) gen_sort_inplace(D,NULL,&cmpi1,NULL);
464 36624 : for (i = 1; i <= n; i++) gmael(D,i,2) = fa_clean(P, gmael(D,i,2));
465 12453 : return gerepilecopy(av, D);
466 : }
467 : static int
468 1701 : cmpu1(void *E, GEN va, GEN vb)
469 1701 : { long a = va[1], b = vb[1]; (void)E; return a>b? 1: (a<b? -1: 0); }
470 : static GEN
471 1498 : fa_clean_u(GEN P, GEN E)
472 : {
473 1498 : long i, j, l = lg(E);
474 1498 : GEN Q = cgetg(l, t_VECSMALL);
475 4417 : for (i = j = 1; i < l; i++)
476 2919 : if (E[i]) { Q[j] = P[i]; E[j] = E[i]; j++; }
477 1498 : setlg(Q,j);
478 1498 : setlg(E,j); return mkmat2(Q,E);
479 : }
480 : GEN
481 350 : divisorsu_fact_factored(GEN fa)
482 : {
483 350 : pari_sp av = avma;
484 350 : GEN *d, *t1, *t2, *t3, vD, D, P = gel(fa,1), E = gel(fa,2);
485 350 : long i, j, l, n = ndiv(E);
486 :
487 350 : D = cgetg(n+1,t_VEC); d = (GEN*)D;
488 350 : l = lg(E);
489 350 : *++d = mkvec2((GEN)1, const_vecsmall(l-1,0));
490 924 : for (i=1; i<l; i++)
491 1330 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
492 1904 : for (t2=d, t3=t1; t3<t2; )
493 : {
494 : ulong a;
495 : GEN b;
496 1148 : a = (*++t3)[1] * P[i];
497 1148 : b = leafcopy(gel(*t3,2)); b[i]++;
498 1148 : *++d = mkvec2((GEN)a,b);
499 : }
500 350 : gen_sort_inplace(D,NULL,&cmpu1,NULL);
501 350 : vD = cgetg(n+1, t_VECSMALL);
502 1848 : for (i = 1; i <= n; i++)
503 : {
504 1498 : vD[i] = umael(D,i,1);
505 1498 : gel(D,i) = fa_clean_u(P, gmael(D,i,2));
506 : }
507 350 : return gerepilecopy(av, mkvec2(vD,D));
508 : }
509 : GEN
510 366401 : divisors(GEN N)
511 : {
512 : long i, j, l;
513 : GEN *d, *t1, *t2, *t3, D, P, E;
514 366401 : int isint = divisors_init(N, &P, &E);
515 366380 : GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
516 :
517 366380 : D = cgetg(ndiv(E)+1,t_VEC); d = (GEN*)D;
518 366380 : l = lg(E);
519 366380 : *++d = gen_1;
520 1198099 : for (i=1; i<l; i++)
521 2232013 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
522 5402873 : for (t2=d, t3=t1; t3<t2; ) *++d = mul(*++t3, gel(P,i));
523 366380 : if (isint) ZV_sort_inplace(D);
524 366380 : return D;
525 : }
526 : GEN
527 784 : divisors0(GEN N, long flag)
528 : {
529 784 : if (flag && flag != 1) pari_err_FLAG("divisors");
530 784 : return flag? divisors_factored(N): divisors(N);
531 : }
532 :
533 : GEN
534 13034 : divisorsu_moebius(GEN P)
535 : {
536 : GEN d, t, t2, t3;
537 13034 : long i, l = lg(P);
538 13034 : d = t = cgetg((1 << (l-1)) + 1, t_VECSMALL);
539 13034 : *++d = 1;
540 30177 : for (i=1; i<l; i++)
541 40845 : for (t2=d, t3=t; t3<t2; ) *(++d) = *(++t3) * -P[i];
542 13034 : return t;
543 : }
544 : GEN
545 18340432 : divisorsu_fact(GEN fa)
546 : {
547 18340432 : GEN d, t, t1, t2, t3, P = gel(fa,1), E = gel(fa,2);
548 18340432 : long i, j, l = lg(P);
549 18340432 : d = t = cgetg(numdivu_fact(fa) + 1,t_VECSMALL);
550 19760392 : *++d = 1;
551 64934439 : for (i=1; i<l; i++)
552 99349019 : for (t1=t,j=E[i]; j; j--,t1=t2)
553 196250375 : for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
554 19760392 : vecsmall_sort(t); return t;
555 : }
556 : GEN
557 161780 : divisorsu(ulong N)
558 : {
559 161780 : pari_sp av = avma;
560 161780 : return gerepileupto(av, divisorsu_fact(factoru(N)));
561 : }
562 :
563 : static GEN
564 0 : corefa(GEN fa)
565 : {
566 0 : GEN P = gel(fa,1), E = gel(fa,2), c = gen_1;
567 : long i;
568 0 : for (i=1; i<lg(P); i++)
569 0 : if (mod2(gel(E,i))) c = mulii(c,gel(P,i));
570 0 : return c;
571 : }
572 : static GEN
573 0 : core2fa(GEN fa)
574 : {
575 0 : GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
576 : long i;
577 0 : for (i=1; i<lg(P); i++)
578 : {
579 0 : long e = itos(gel(E,i));
580 0 : if (e & 1) c = mulii(c, gel(P,i));
581 0 : if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
582 : }
583 0 : return mkvec2(c,f);
584 : }
585 : GEN
586 0 : corepartial(GEN n, long all)
587 : {
588 0 : pari_sp av = avma;
589 0 : if (typ(n) != t_INT) pari_err_TYPE("corepartial",n);
590 0 : return gerepileuptoint(av, corefa(Z_factor_limit(n,all)));
591 : }
592 : GEN
593 0 : core2partial(GEN n, long all)
594 : {
595 0 : pari_sp av = avma;
596 0 : if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
597 0 : return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
598 : }
599 : /* given an arithmetic function argument, return the underlying integer */
600 : static GEN
601 58050 : arith_n(GEN A)
602 : {
603 58050 : switch(typ(A))
604 : {
605 53010 : case t_INT: return A;
606 2863 : case t_VEC: return gel(A,1);
607 2177 : default: return factorback(A);
608 : }
609 : }
610 : static GEN
611 53766 : core2_i(GEN n)
612 : {
613 53766 : GEN f = core(n);
614 53766 : if (!signe(f)) return mkvec2(gen_0, gen_1);
615 53731 : return mkvec2(f, sqrtint(diviiexact(arith_n(n), f)));
616 : }
617 : GEN
618 53682 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
619 :
620 : GEN
621 3101 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
622 :
623 : static long
624 8596 : _mod4(GEN c) {
625 8596 : long r, s = signe(c);
626 8596 : if (!s) return 0;
627 8596 : r = mod4(c); if (s < 0) r = 4-r;
628 8596 : return r;
629 : }
630 :
631 : long
632 56 : corediscs(long D, ulong *f)
633 : { /* D = f^2 d */
634 56 : long d = D >= 0? (long)coreu(D) : -(long)coreu(-(ulong)D);
635 56 : if ((((ulong)d)&3UL) != 1) d *= 4;
636 56 : if (f) *f = usqrt((ulong)(D/d));
637 56 : return d;
638 : }
639 :
640 : GEN
641 8512 : coredisc(GEN n)
642 : {
643 8512 : pari_sp av = avma;
644 8512 : GEN c = core(n);
645 8512 : if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
646 3045 : return gerepileuptoint(av, shifti(c,2));
647 : }
648 :
649 : GEN
650 84 : coredisc2(GEN n)
651 : {
652 84 : pari_sp av = avma;
653 84 : GEN y = core2_i(n);
654 84 : GEN c = gel(y,1), f = gel(y,2);
655 84 : if (_mod4(c)<=1) return gerepilecopy(av, y);
656 14 : y = cgetg(3,t_VEC);
657 14 : gel(y,1) = shifti(c,2);
658 14 : gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
659 : }
660 :
661 : GEN
662 56 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
663 :
664 : /* Write x = Df^2, where D = fundamental discriminant,
665 : * P^E = factorisation of conductor f */
666 : GEN
667 129817 : coredisc2_fact(GEN fa, long s, GEN *pP, GEN *pE)
668 : {
669 129817 : GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2), D = s > 0? gen_1: gen_m1;
670 129817 : long l = lg(P0), i, j;
671 :
672 129817 : E = cgetg(l, t_VECSMALL);
673 129817 : P = cgetg(l, t_VEC);
674 492140 : for (i = j = 1; i < l; i++)
675 : {
676 362330 : long e = itos(gel(E0,i));
677 362330 : GEN p = gel(P0,i);
678 362330 : if (odd(e)) D = mulii(D, p);
679 362314 : e >>= 1; if (e) { gel(P, j) = p; E[j] = e; j++; }
680 : }
681 129810 : if (Mod4(D) != 1)
682 : {
683 40327 : D = shifti(D, 2);
684 40327 : if (!--E[1])
685 : {
686 34496 : P[1] = P[0]; P++;
687 34496 : E[1] = E[0]; E++; j--;
688 : }
689 : }
690 129810 : setlg(P,j); *pP = P;
691 129811 : setlg(E,j); *pE = E; return D;
692 : }
693 : ulong
694 7137 : coredisc2u_fact(GEN fa, long s, GEN *pP, GEN *pE)
695 : {
696 7137 : GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2);
697 7137 : ulong D = 1;
698 7137 : long i, j, l = lg(P0);
699 :
700 7137 : E = cgetg(l, t_VECSMALL);
701 7137 : P = cgetg(l, t_VECSMALL);
702 21303 : for (i = j = 1; i < l; i++)
703 : {
704 14166 : long e = E0[i], p = P0[i];
705 14166 : if (odd(e)) D *= p;
706 14166 : e >>= 1; if (e) { P[j] = p; E[j] = e; j++; }
707 : }
708 7137 : if ((D & 3) != (s > 0? 1: 3))
709 : {
710 1612 : D *= 4;
711 1612 : if (!--E[1])
712 : {
713 1199 : P[1] = P[0]; P++;
714 1199 : E[1] = E[0]; E++; j--;
715 : }
716 : }
717 7137 : setlg(P,j); *pP = P;
718 7137 : setlg(E,j); *pE = E; return D;
719 : }
720 :
721 : long
722 815 : omegau(ulong n)
723 : {
724 : pari_sp av;
725 815 : if (n == 1UL) return 0;
726 801 : av = avma; return gc_long(av, nbrows(factoru(n)));
727 : }
728 : long
729 1589 : omega(GEN n)
730 : {
731 : pari_sp av;
732 : GEN F, P;
733 1589 : if ((F = check_arith_non0(n,"omega"))) {
734 : long n;
735 721 : P = gel(F,1); n = lg(P)-1;
736 721 : if (n && equalim1(gel(P,1))) n--;
737 721 : return n;
738 : }
739 854 : if (lgefint(n) == 3) return omegau(n[2]);
740 39 : av = avma;
741 39 : F = absZ_factor(n);
742 39 : return gc_long(av, nbrows(F));
743 : }
744 :
745 : long
746 835 : bigomegau(ulong n)
747 : {
748 : pari_sp av;
749 835 : if (n == 1) return 0;
750 821 : av = avma; return gc_long(av, zv_sum(gel(factoru(n),2)));
751 : }
752 : long
753 1589 : bigomega(GEN n)
754 : {
755 1589 : pari_sp av = avma;
756 : GEN F, E;
757 1589 : if ((F = check_arith_non0(n,"bigomega")))
758 : {
759 721 : GEN P = gel(F,1);
760 721 : long n = lg(P)-1;
761 721 : E = gel(F,2);
762 721 : if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
763 : }
764 854 : else if (lgefint(n) == 3)
765 821 : return bigomegau(n[2]);
766 : else
767 33 : E = gel(absZ_factor(n), 2);
768 754 : E = ZV_to_zv(E);
769 754 : return gc_long(av, zv_sum(E));
770 : }
771 :
772 : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
773 : ulong
774 634791 : eulerphiu_fact(GEN f)
775 : {
776 634791 : GEN P = gel(f,1), E = gel(f,2);
777 634791 : long i, m = 1, l = lg(P);
778 1830941 : for (i = 1; i < l; i++)
779 : {
780 1196150 : ulong p = P[i], e = E[i];
781 1196150 : if (!e) continue;
782 1196136 : if (p == 2)
783 341568 : { if (e > 1) m <<= e-1; }
784 : else
785 : {
786 854568 : m *= (p-1);
787 854568 : if (e > 1) m *= upowuu(p, e-1);
788 : }
789 : }
790 634791 : return m;
791 : }
792 : ulong
793 223703 : eulerphiu(ulong n)
794 : {
795 : pari_sp av;
796 223703 : if (!n) return 2;
797 223703 : av = avma; return gc_long(av, eulerphiu_fact(factoru(n)));
798 : }
799 : GEN
800 150906 : eulerphi(GEN n)
801 : {
802 150906 : pari_sp av = avma;
803 : GEN Q, F, P, E;
804 : long i, l;
805 :
806 150906 : if ((F = check_arith_all(n,"eulerphi")))
807 : {
808 3591 : F = clean_Z_factor(F);
809 3591 : n = arith_n(n);
810 3591 : if (lgefint(n) == 3)
811 : {
812 : ulong e;
813 3577 : F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
814 3577 : e = eulerphiu_fact(F);
815 3577 : return gc_utoipos(av, e);
816 : }
817 : }
818 147315 : else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
819 : else
820 53 : F = absZ_factor(n);
821 67 : if (!signe(n)) return gen_2;
822 39 : P = gel(F,1);
823 39 : E = gel(F,2); l = lg(P);
824 39 : Q = cgetg(l, t_VEC);
825 252 : for (i = 1; i < l; i++)
826 : {
827 213 : GEN p = gel(P,i), q;
828 213 : ulong v = itou(gel(E,i));
829 213 : q = subiu(p,1);
830 213 : if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
831 213 : gel(Q,i) = q;
832 : }
833 39 : return gerepileuptoint(av, ZV_prod(Q));
834 : }
835 :
836 : long
837 18470888 : numdivu_fact(GEN fa)
838 : {
839 18470888 : GEN E = gel(fa,2);
840 18470888 : long n = 1, i, l = lg(E);
841 63002412 : for (i = 1; i < l; i++) n *= E[i]+1;
842 18470888 : return n;
843 : }
844 : long
845 815 : numdivu(long N)
846 : {
847 : pari_sp av;
848 815 : if (N == 1) return 1;
849 801 : av = avma; return gc_long(av, numdivu_fact(factoru(N)));
850 : }
851 : static GEN
852 760 : numdiv_aux(GEN F)
853 : {
854 760 : GEN x, E = gel(F,2);
855 760 : long i, l = lg(E);
856 760 : x = cgetg(l, t_VECSMALL);
857 2079 : for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
858 760 : return x;
859 : }
860 : GEN
861 1589 : numdiv(GEN n)
862 : {
863 1589 : pari_sp av = avma;
864 : GEN F, E;
865 1589 : if ((F = check_arith_non0(n,"numdiv")))
866 : {
867 721 : F = clean_Z_factor(F);
868 721 : E = numdiv_aux(F);
869 : }
870 854 : else if (lgefint(n) == 3)
871 815 : return utoipos(numdivu(n[2]));
872 : else
873 39 : E = numdiv_aux(absZ_factor(n));
874 760 : return gerepileuptoint(av, zv_prod_Z(E));
875 : }
876 :
877 : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
878 : static GEN
879 55598 : u_euler_sumdiv(ulong p, long v)
880 : {
881 55598 : GEN u = utoipos(1 + p); /* can't overflow */
882 76080 : for (; v > 1; v--) u = addui(1, mului(p, u));
883 55598 : return u;
884 : }
885 : /* 1 + q + ... + q^v */
886 : static GEN
887 18168 : euler_sumdiv(GEN q, long v)
888 : {
889 18168 : GEN u = addui(1, q);
890 24734 : for (; v > 1; v--) u = addui(1, mulii(q, u));
891 18168 : return u;
892 : }
893 :
894 : static GEN
895 1474 : sumdiv_aux(GEN F)
896 : {
897 1474 : GEN x, P = gel(F,1), E = gel(F,2);
898 1474 : long i, l = lg(P);
899 1474 : x = cgetg(l, t_VEC);
900 3892 : for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
901 1474 : return ZV_prod(x);
902 : }
903 : GEN
904 3017 : sumdiv(GEN n)
905 : {
906 3017 : pari_sp av = avma;
907 : GEN F, v;
908 :
909 3017 : if ((F = check_arith_non0(n,"sumdiv")))
910 : {
911 1442 : F = clean_Z_factor(F);
912 1442 : v = sumdiv_aux(F);
913 : }
914 1561 : else if (lgefint(n) == 3)
915 : {
916 1529 : if (n[2] == 1) return gen_1;
917 1501 : F = factoru(n[2]);
918 1501 : v = usumdiv_fact(F);
919 : }
920 : else
921 32 : v = sumdiv_aux(absZ_factor(n));
922 2975 : return gerepileuptoint(av, v);
923 : }
924 :
925 : static GEN
926 1506 : sumdivk_aux(GEN F, long k)
927 : {
928 1506 : GEN x, P = gel(F,1), E = gel(F,2);
929 1506 : long i, l = lg(P);
930 1506 : x = cgetg(l, t_VEC);
931 4116 : for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(powiu(gel(P,i),k), gel(E,i)[2]);
932 1506 : return ZV_prod(x);
933 : }
934 : GEN
935 7987 : sumdivk(GEN n, long k)
936 : {
937 7987 : pari_sp av = avma;
938 : GEN F, v;
939 : long k1;
940 :
941 7987 : if (!k) return numdiv(n);
942 7987 : if (k == 1) return sumdiv(n);
943 6398 : if ((F = check_arith_non0(n,"sumdivk"))) F = clean_Z_factor(F);
944 6356 : k1 = k; if (k < 0) k = -k;
945 6356 : if (k == 1)
946 1428 : v = sumdiv(F? F: n);
947 : else
948 : {
949 4928 : if (F)
950 1442 : v = sumdivk_aux(F,k);
951 3486 : else if (lgefint(n) == 3)
952 : {
953 3422 : if (n[2] == 1) return gen_1;
954 3324 : F = factoru(n[2]);
955 3324 : v = usumdivk_fact(F,k);
956 : }
957 : else
958 64 : v = sumdivk_aux(absZ_factor(n), k);
959 4830 : if (k1 > 0) return gerepileuptoint(av, v);
960 : }
961 :
962 1435 : if (F) n = arith_n(n);
963 1435 : if (k != 1) n = powiu(n,k);
964 1435 : return gerepileupto(av, gdiv(v, n));
965 : }
966 :
967 : GEN
968 34352 : usumdiv_fact(GEN f)
969 : {
970 34352 : GEN P = gel(f,1), E = gel(f,2);
971 34352 : long i, l = lg(P);
972 34352 : GEN v = cgetg(l, t_VEC);
973 89950 : for (i=1; i<l; i++) gel(v,i) = u_euler_sumdiv(P[i],E[i]);
974 34352 : return ZV_prod(v);
975 : }
976 : GEN
977 9372 : usumdivk_fact(GEN f, ulong k)
978 : {
979 9372 : GEN P = gel(f,1), E = gel(f,2);
980 9372 : long i, l = lg(P);
981 9372 : GEN v = cgetg(l, t_VEC);
982 22512 : for (i=1; i<l; i++) gel(v,i) = euler_sumdiv(powuu(P[i],k),E[i]);
983 9372 : return ZV_prod(v);
984 : }
985 :
986 : long
987 2919 : uissquarefree_fact(GEN f)
988 : {
989 2919 : GEN E = gel(f,2);
990 2919 : long i, l = lg(E);
991 2919 : if (l == 2) return umael(f,1,1)? E[1] == 1: 0; /* handle factor(0) */
992 5824 : for (i = 1; i < l; i++)
993 4060 : if (E[i] > 1) return 0;
994 1764 : return 1;
995 : }
996 : long
997 2451618 : uissquarefree(ulong n)
998 : {
999 2451618 : if (!n) return 0;
1000 2451618 : return moebiusu(n)? 1: 0;
1001 : }
1002 : long
1003 2214 : Z_issquarefree(GEN n)
1004 : {
1005 2214 : switch(lgefint(n))
1006 : {
1007 14 : case 2: return 0;
1008 1424 : case 3: return uissquarefree(n[2]);
1009 : }
1010 776 : return moebius(n)? 1: 0;
1011 : }
1012 :
1013 : long
1014 2814 : Z_issquarefree_fact(GEN F)
1015 : {
1016 2814 : GEN E = gel(F,2);
1017 2814 : long i, l = lg(E);
1018 2814 : if (l == 2) return signe(gcoeff(F,1,1))? equali1(gel(E,1)): 0;
1019 6300 : for(i = 1; i < l; i++)
1020 4956 : if (!equali1(gel(E,i))) return 0;
1021 1344 : return 1;
1022 : }
1023 : long
1024 6454 : issquarefree(GEN x)
1025 : {
1026 : pari_sp av;
1027 : GEN d;
1028 6454 : switch(typ(x))
1029 : {
1030 1421 : case t_INT: return Z_issquarefree(x);
1031 3626 : case t_POL:
1032 3626 : if (!signe(x)) return 0;
1033 3626 : av = avma; d = RgX_gcd(x, RgX_deriv(x));
1034 3626 : return gc_long(av, lg(d)==3);
1035 1407 : case t_VEC:
1036 1407 : case t_MAT: return Z_issquarefree_fact(check_arith_all(x,"issquarefree"));
1037 0 : default: pari_err_TYPE("issquarefree",x);
1038 : return 0; /* LCOV_EXCL_LINE */
1039 : }
1040 : }
|
