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 70 : boundfact(GEN n, ulong lim)
26 : {
27 70 : switch(typ(n))
28 : {
29 58 : case t_INT: return Z_factor_limit(n,lim);
30 12 : case t_FRAC: {
31 12 : pari_sp av = avma;
32 12 : GEN a = Z_factor_limit(gel(n,1),lim);
33 12 : GEN b = Z_factor_limit(gel(n,2),lim);
34 12 : gel(b,2) = ZC_neg(gel(b,2));
35 12 : return gc_GEN(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 16301 : Z_lsmoothen(GEN N, GEN L, GEN *pP, GEN *pE)
45 : {
46 16301 : long i, j, l = lg(L);
47 16301 : GEN E = new_chunk(l), P = new_chunk(l);
48 151905 : for (i = j = 1; i < l; i++)
49 : {
50 144874 : ulong p = uel(L,i);
51 144874 : long v = Z_lvalrem(N, p, &N);
52 144874 : if (v) { P[j] = p; E[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
53 : }
54 16301 : P[0] = evaltyp(t_VECSMALL) | _evallg(j); if (pP) *pP = P;
55 16301 : E[0] = evaltyp(t_VECSMALL) | _evallg(j); if (pE) *pE = E; return N;
56 : }
57 : GEN
58 66725 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pE)
59 : {
60 66725 : long i, j, l = lg(L);
61 66725 : GEN E = new_chunk(l), P = new_chunk(l);
62 241464 : for (i = j = 1; i < l; i++)
63 : {
64 222565 : GEN p = gel(L,i);
65 222565 : long v = Z_pvalrem(N, p, &N);
66 222565 : if (v)
67 : {
68 165625 : gel(P,j) = p; gel(E,j) = utoipos(v); j++;
69 165625 : if (is_pm1(N)) { N = NULL; break; }
70 : }
71 : }
72 66725 : P[0] = evaltyp(t_COL) | _evallg(j); if (pP) *pP = P;
73 66725 : 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 118977 : factoru_pow(ulong n)
82 : {
83 118977 : GEN f = cgetg(4,t_VEC);
84 118977 : 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 118977 : (void)new_chunk((15 + 1)*3);
89 118977 : F = factoru(n);
90 118977 : P = gel(F,1);
91 118977 : E = gel(F,2); l = lg(P);
92 118977 : set_avma(av);
93 118977 : gel(f,1) = p = cgetg(l,t_VECSMALL);
94 118977 : gel(f,2) = e = cgetg(l,t_VECSMALL);
95 118977 : gel(f,3) = c = cgetg(l,t_VECSMALL);
96 263287 : for(i = 1; i < l; i++)
97 : {
98 144310 : p[i] = P[i];
99 144310 : e[i] = E[i];
100 144310 : c[i] = upowuu(p[i], e[i]);
101 : }
102 118977 : return f;
103 : }
104 :
105 : static GEN
106 13931 : factorlim(GEN n, ulong lim)
107 13931 : { 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 6330 : factor_pn_1_limit(GEN p, long n, ulong lim)
112 : {
113 6330 : pari_sp av = avma;
114 6330 : GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
115 6330 : long i, pp = itos_or_0(p);
116 13639 : for(i=2; i<lg(d); i++)
117 : {
118 : GEN B;
119 7309 : if (pp && d[i]%pp==0 && (
120 2945 : ((pp&3)==1 && (d[i]&1)) ||
121 2927 : ((pp&3)==3 && (d[i]&3)==2) ||
122 2825 : (pp==2 && (d[i]&7)==4)))
123 292 : {
124 292 : GEN f=factor_Aurifeuille_prime(p,d[i]);
125 292 : B = factorlim(f, lim);
126 292 : A = ZM_merge_factor(A, B);
127 292 : B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
128 : }
129 : else
130 7017 : B = factorlim(polcyclo_eval(d[i],p), lim);
131 7309 : A = ZM_merge_factor(A, B);
132 : }
133 6330 : return gc_GEN(av, A);
134 : }
135 : GEN
136 6330 : factor_pn_1(GEN p, ulong n)
137 6330 : { 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 13678183 : RgV_is_ZVpos(GEN v)
201 : {
202 13678183 : long i, l = lg(v);
203 39505101 : for (i = 1; i < l; i++)
204 : {
205 25835613 : GEN c = gel(v,i);
206 25835613 : if (typ(c) != t_INT || signe(c) <= 0) return 0;
207 : }
208 13669488 : return 1;
209 : }
210 : /* check whether v is a ZV with nonzero entries */
211 : int
212 17717 : RgV_is_ZVnon0(GEN v)
213 : {
214 17717 : long i, l = lg(v);
215 135352 : for (i = 1; i < l; i++)
216 : {
217 117675 : GEN c = gel(v,i);
218 117675 : if (typ(c) != t_INT || !signe(c)) return 0;
219 : }
220 17677 : return 1;
221 : }
222 : /* check whether v is a ZV with nonzero entries OR exactly [0] */
223 : static int
224 603862 : RgV_is_ZV0(GEN v)
225 : {
226 603862 : long i, l = lg(v);
227 2074137 : for (i = 1; i < l; i++)
228 : {
229 1470377 : GEN c = gel(v,i);
230 : long s;
231 1470377 : if (typ(c) != t_INT) return 0;
232 1470377 : s = signe(c);
233 1470377 : if (!s) return (l == 2);
234 : }
235 603760 : return 1;
236 : }
237 :
238 : int
239 277298 : RgV_is_prV(GEN v)
240 : {
241 277298 : long l = lg(v), i;
242 624168 : for (i = 1; i < l; i++)
243 513895 : if (!checkprid_i(gel(v,i))) return 0;
244 110273 : return 1;
245 : }
246 : int
247 234079 : is_nf_factor(GEN F)
248 : {
249 218026 : return typ(F) == t_MAT && lg(F) == 3
250 452105 : && RgV_is_prV(gel(F,1)) && RgV_is_ZVpos(gel(F,2));
251 : }
252 : int
253 12 : is_nf_extfactor(GEN F)
254 : {
255 12 : return typ(F) == t_MAT && lg(F) == 3
256 24 : && RgV_is_prV(gel(F,1)) && RgV_is_ZV(gel(F,2));
257 : }
258 :
259 : static int
260 7118077 : is_Z_factor_i(GEN f)
261 7118077 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
262 : int
263 6503703 : is_Z_factorpos(GEN f)
264 6503703 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
265 : int
266 603862 : is_Z_factor(GEN f)
267 603862 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
268 : /* as is_Z_factorpos, also allow factor(0) */
269 : int
270 10512 : is_Z_factornon0(GEN f)
271 10512 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
272 : GEN
273 20849 : clean_Z_factor(GEN f)
274 : {
275 20849 : GEN P = gel(f,1);
276 20849 : long n = lg(P)-1;
277 20849 : if (n && equalim1(gel(P,1)))
278 2190 : return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
279 18659 : 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 476 : check_arith_pos(GEN n, const char *f) {
302 476 : switch(typ(n))
303 : {
304 455 : case t_INT:
305 455 : if (signe(n) <= 0) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
306 455 : return NULL;
307 5 : case t_VEC:
308 5 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT || signe(gel(n,1)) <= 0) break;
309 5 : n = gel(n,2); /* fall through */
310 16 : case t_MAT:
311 16 : if (!is_Z_factorpos(n)) break;
312 16 : return n;
313 : }
314 5 : 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 1841895 : check_arith_non0(GEN n, const char *f) {
320 1841895 : switch(typ(n))
321 : {
322 1831577 : case t_INT:
323 1831577 : if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
324 1831531 : return NULL;
325 1018 : case t_VEC:
326 1018 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT || !signe(gel(n,1))) break;
327 1018 : n = gel(n,2); /* fall through */
328 10318 : case t_MAT:
329 10318 : if (!is_Z_factornon0(n)) break;
330 10278 : return n;
331 : }
332 40 : pari_err_TYPE(f,n);
333 : return NULL;/*LCOV_EXCL_LINE*/
334 : }
335 : /* n attached to a factorization of an integer */
336 : GEN
337 2200196 : check_arith_all(GEN n, const char *f) {
338 2200196 : switch(typ(n))
339 : {
340 1596329 : case t_INT:
341 1596329 : return NULL;
342 588942 : case t_VEC:
343 588942 : if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
344 588942 : n = gel(n,2); /* fall through */
345 603862 : case t_MAT:
346 603862 : if (!is_Z_factor(n)) break;
347 603862 : return n;
348 : }
349 5 : 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 259250 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
360 : {
361 : GEN E, P;
362 259250 : if (lg(F) != 3) pari_err_TYPE("divisors",F);
363 259250 : P = gel(F,1);
364 259250 : E = gel(F,2);
365 259250 : RgV_check_ZV(E, "divisors");
366 259250 : *isint = RgV_is_ZV(P);
367 259250 : if (*isint)
368 : {
369 259240 : long i, l = lg(P);
370 : /* skip -1 */
371 259240 : if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
372 : /* test for 0 */
373 839615 : for (i = 1; i < l; i++)
374 580380 : if (!signe(gel(P,i)) && signe(gel(E,i)))
375 5 : pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
376 : }
377 259245 : *pP = P;
378 259245 : *pE = E;
379 259245 : }
380 : static void
381 3092 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
382 :
383 : int
384 262347 : divisors_init(GEN n, GEN *pP, GEN *pE)
385 : {
386 : long i,l;
387 : GEN E, P, e;
388 : int isint;
389 :
390 262347 : switch(typ(n))
391 : {
392 3076 : case t_INT:
393 3076 : if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
394 3076 : set_fact(absZ_factor(n), &P,&E);
395 3076 : isint = 1; break;
396 1384 : case t_VEC:
397 1384 : if (lg(n) != 3 || typ(gel(n,2)) !=t_MAT) pari_err_TYPE("divisors",n);
398 1379 : set_fact_check(gel(n,2), &P,&E, &isint);
399 1379 : break;
400 257871 : case t_MAT:
401 257871 : set_fact_check(n, &P,&E, &isint);
402 257866 : break;
403 16 : default:
404 16 : set_fact(factor(n), &P,&E);
405 16 : isint = 0; break;
406 : }
407 262337 : l = lg(P);
408 262337 : e = cgetg(l, t_VECSMALL);
409 847374 : for (i=1; i<l; i++)
410 : {
411 585043 : e[i] = itos(gel(E,i));
412 585043 : if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
413 : }
414 262331 : *pP = P; *pE = e; return isint;
415 : }
416 :
417 : static long
418 262560 : ndiv(GEN E)
419 : {
420 : long i, l;
421 262560 : GEN e = cgetg_copy(E, &l); /* left on stack */
422 : ulong n;
423 847907 : for (i=1; i<l; i++) e[i] = E[i]+1;
424 262560 : n = itou_or_0( zv_prod_Z(e) );
425 262560 : if (!n || n & ~LGBITS) pari_err_OVERFLOW("divisors");
426 262560 : return n;
427 : }
428 : static int
429 13149 : 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 18272 : fa_clean(GEN P, GEN E)
434 : {
435 18272 : long i, j, l = lg(E);
436 18272 : GEN Q = cgetg(l, t_COL);
437 51235 : for (i = j = 1; i < l; i++)
438 32963 : if (E[i]) { gel(Q,j) = gel(P,i); E[j] = E[i]; j++; }
439 18272 : setlg(Q,j);
440 18272 : setlg(E,j); return mkmat2(Q,Flc_to_ZC(E));
441 : }
442 : GEN
443 9170 : divisors_factored(GEN N)
444 : {
445 9170 : pari_sp av = avma;
446 : GEN *d, *t1, *t2, *t3, D, P, E;
447 9170 : int isint = divisors_init(N, &P, &E);
448 9170 : GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
449 9170 : long i, j, l, n = ndiv(E);
450 :
451 9170 : D = cgetg(n+1,t_VEC); d = (GEN*)D;
452 9170 : l = lg(E);
453 9170 : *++d = mkvec2(gen_1, const_vecsmall(l-1,0));
454 24707 : for (i=1; i<l; i++)
455 22675 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
456 16240 : for (t2=d, t3=t1; t3<t2; )
457 : {
458 : GEN a, b;
459 9102 : a = mul(gel(*++t3,1), gel(P,i));
460 9102 : b = leafcopy(gel(*t3,2)); b[i]++;
461 9102 : *++d = mkvec2(a,b);
462 : }
463 9170 : if (isint) gen_sort_inplace(D,NULL,&cmpi1,NULL);
464 27442 : for (i = 1; i <= n; i++) gmael(D,i,2) = fa_clean(P, gmael(D,i,2));
465 9170 : return gc_GEN(av, D);
466 : }
467 : static int
468 1290 : cmpu1(void *E, GEN va, GEN vb)
469 1290 : { long a = va[1], b = vb[1]; (void)E; return a>b? 1: (a<b? -1: 0); }
470 : static GEN
471 1137 : fa_clean_u(GEN P, GEN E)
472 : {
473 1137 : long i, j, l = lg(E);
474 1137 : GEN Q = cgetg(l, t_VECSMALL);
475 3355 : for (i = j = 1; i < l; i++)
476 2218 : if (E[i]) { Q[j] = P[i]; E[j] = E[i]; j++; }
477 1137 : setlg(Q,j);
478 1137 : setlg(E,j); return mkmat2(Q,E);
479 : }
480 : GEN
481 264 : divisorsu_fact_factored(GEN fa)
482 : {
483 264 : pari_sp av = avma;
484 264 : GEN *d, *t1, *t2, *t3, vD, D, P = gel(fa,1), E = gel(fa,2);
485 264 : long i, j, l, n = ndiv(E);
486 :
487 264 : D = cgetg(n+1,t_VEC); d = (GEN*)D;
488 264 : l = lg(E);
489 264 : *++d = mkvec2((GEN)1, const_vecsmall(l-1,0));
490 699 : for (i=1; i<l; i++)
491 1008 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
492 1446 : for (t2=d, t3=t1; t3<t2; )
493 : {
494 : ulong a;
495 : GEN b;
496 873 : a = (*++t3)[1] * P[i];
497 873 : b = leafcopy(gel(*t3,2)); b[i]++;
498 873 : *++d = mkvec2((GEN)a,b);
499 : }
500 264 : gen_sort_inplace(D,NULL,&cmpu1,NULL);
501 264 : vD = cgetg(n+1, t_VECSMALL);
502 1401 : for (i = 1; i <= n; i++)
503 : {
504 1137 : vD[i] = umael(D,i,1);
505 1137 : gel(D,i) = fa_clean_u(P, gmael(D,i,2));
506 : }
507 264 : return gc_GEN(av, mkvec2(vD,D));
508 : }
509 : GEN
510 253142 : divisors(GEN N)
511 : {
512 : long i, j, l;
513 : GEN *d, *t1, *t2, *t3, D, P, E;
514 253142 : int isint = divisors_init(N, &P, &E);
515 253126 : GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
516 :
517 253126 : D = cgetg(ndiv(E)+1,t_VEC); d = (GEN*)D;
518 253126 : l = lg(E);
519 253126 : *++d = gen_1;
520 822501 : for (i=1; i<l; i++)
521 1539300 : for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
522 3720240 : for (t2=d, t3=t1; t3<t2; ) *++d = mul(*++t3, gel(P,i));
523 253126 : if (isint) ZV_sort_inplace(D);
524 253126 : return D;
525 : }
526 : GEN
527 666 : divisors0(GEN N, long flag)
528 : {
529 666 : if (flag && flag != 1) pari_err_FLAG("divisors");
530 666 : return flag? divisors_factored(N): divisors(N);
531 : }
532 :
533 : GEN
534 19715 : divisorsu_moebius(GEN P)
535 : {
536 : GEN d, t, t2, t3;
537 19715 : long i, l = lg(P);
538 19715 : d = t = cgetg((1 << (l-1)) + 1, t_VECSMALL);
539 19715 : *++d = 1;
540 51350 : for (i=1; i<l; i++)
541 81333 : for (t2=d, t3=t; t3<t2; ) *(++d) = *(++t3) * -P[i];
542 19715 : return t;
543 : }
544 : GEN
545 8245160 : divisorsu_fact(GEN fa)
546 : {
547 8245160 : GEN d, t, t1, t2, t3, P = gel(fa,1), E = gel(fa,2);
548 8245160 : long i, j, l = lg(P);
549 8245160 : d = t = cgetg(numdivu_fact(fa) + 1,t_VECSMALL);
550 8245160 : *++d = 1;
551 27860770 : for (i=1; i<l; i++)
552 43726213 : for (t1=t,j=E[i]; j; j--,t1=t2)
553 88512698 : for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
554 8245160 : vecsmall_sort(t); return t;
555 : }
556 : GEN
557 135275 : divisorsu(ulong N)
558 : {
559 135275 : pari_sp av = avma;
560 135275 : return gc_upto(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 gc_INT(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 gc_GEN(av, core2fa(Z_factor_limit(n,all)));
598 : }
599 : /* given an arithmetic function argument, return the underlying integer */
600 : static GEN
601 50106 : arith_n(GEN A)
602 : {
603 50106 : switch(typ(A))
604 : {
605 46101 : case t_INT: return A;
606 2450 : case t_VEC: return gel(A,1);
607 1555 : default: return factorback(A);
608 : }
609 : }
610 : static GEN
611 46641 : core2_i(GEN n)
612 : {
613 46641 : GEN f = core(n);
614 46641 : if (!signe(f)) return mkvec2(gen_0, gen_1);
615 46616 : return mkvec2(f, sqrtint(diviiexact(arith_n(n), f)));
616 : }
617 : GEN
618 46570 : core2(GEN n) { pari_sp av = avma; return gc_GEN(av, core2_i(n)); }
619 :
620 : GEN
621 2236 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
622 :
623 : static long
624 7366 : _mod4(GEN c) {
625 7366 : long r, s = signe(c);
626 7366 : if (!s) return 0;
627 7366 : r = mod4(c); if (s < 0) r = 4-r;
628 7366 : return r;
629 : }
630 :
631 : long
632 48 : corediscs(long D, ulong *f)
633 : { /* D = f^2 d */
634 48 : long d = D >= 0? (long)coreu(D) : -(long)coreu(-(ulong)D);
635 48 : if ((((ulong)d)&3UL) != 1) d *= 4;
636 48 : if (f) *f = usqrt((ulong)(D/d));
637 48 : return d;
638 : }
639 :
640 : GEN
641 7295 : coredisc(GEN n)
642 : {
643 7295 : pari_sp av = avma;
644 7295 : GEN c = core(n);
645 7295 : if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
646 2617 : return gc_INT(av, shifti(c,2));
647 : }
648 :
649 : GEN
650 71 : coredisc2(GEN n)
651 : {
652 71 : pari_sp av = avma;
653 71 : GEN y = core2_i(n);
654 71 : GEN c = gel(y,1), f = gel(y,2);
655 71 : if (_mod4(c)<=1) return gc_GEN(av, y);
656 11 : y = cgetg(3,t_VEC);
657 11 : gel(y,1) = shifti(c,2);
658 11 : gel(y,2) = gmul2n(f,-1); return gc_upto(av, y);
659 : }
660 :
661 : GEN
662 46 : 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 92702 : coredisc2_fact(GEN fa, long s, GEN *pP, GEN *pE)
668 : {
669 92702 : GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2), D = s > 0? gen_1: gen_m1;
670 92702 : long l = lg(P0), i, j;
671 :
672 92702 : E = cgetg(l, t_VECSMALL);
673 92702 : P = cgetg(l, t_VEC);
674 351428 : for (i = j = 1; i < l; i++)
675 : {
676 258726 : long e = itos(gel(E0,i));
677 258726 : GEN p = gel(P0,i);
678 258726 : if (odd(e)) D = mulii(D, p);
679 258726 : e >>= 1; if (e) { gel(P, j) = p; E[j] = e; j++; }
680 : }
681 92702 : if (Mod4(D) != 1)
682 : {
683 28795 : D = shifti(D, 2);
684 28795 : if (!--E[1])
685 : {
686 24630 : P[1] = P[0]; P++;
687 24630 : E[1] = E[0]; E++; j--;
688 : }
689 : }
690 92702 : setlg(P,j); *pP = P;
691 92702 : setlg(E,j); *pE = E; return D;
692 : }
693 : ulong
694 5608 : coredisc2u_fact(GEN fa, long s, GEN *pP, GEN *pE)
695 : {
696 5608 : GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2);
697 5608 : ulong D = 1;
698 5608 : long i, j, l = lg(P0);
699 :
700 5608 : E = cgetg(l, t_VECSMALL);
701 5608 : P = cgetg(l, t_VECSMALL);
702 16860 : for (i = j = 1; i < l; i++)
703 : {
704 11252 : long e = E0[i], p = P0[i];
705 11252 : if (odd(e)) D *= p;
706 11252 : e >>= 1; if (e) { P[j] = p; E[j] = e; j++; }
707 : }
708 5608 : if ((D & 3) != (s > 0? 1: 3))
709 : {
710 1223 : D *= 4;
711 1223 : if (!--E[1])
712 : {
713 887 : P[1] = P[0]; P++;
714 887 : E[1] = E[0]; E++; j--;
715 : }
716 : }
717 5608 : setlg(P,j); *pP = P;
718 5608 : setlg(E,j); *pE = E; return D;
719 : }
720 :
721 : long
722 597 : omegau(ulong n)
723 : {
724 : pari_sp av;
725 597 : if (n == 1UL) return 0;
726 587 : av = avma; return gc_long(av, nbrows(factoru(n)));
727 : }
728 : long
729 1155 : omega(GEN n)
730 : {
731 : pari_sp av;
732 : GEN F, P;
733 1155 : if ((F = check_arith_non0(n,"omega"))) {
734 : long n;
735 515 : P = gel(F,1); n = lg(P)-1;
736 515 : if (n && equalim1(gel(P,1))) n--;
737 515 : return n;
738 : }
739 630 : if (lgefint(n) == 3) return omegau(n[2]);
740 33 : av = avma;
741 33 : F = absZ_factor(n);
742 33 : return gc_long(av, nbrows(F));
743 : }
744 :
745 : long
746 613 : bigomegau(ulong n)
747 : {
748 : pari_sp av;
749 613 : if (n == 1) return 0;
750 603 : av = avma; return gc_long(av, zv_sum(gel(factoru(n),2)));
751 : }
752 : long
753 1155 : bigomega(GEN n)
754 : {
755 1155 : pari_sp av = avma;
756 : GEN F, E;
757 1155 : if ((F = check_arith_non0(n,"bigomega")))
758 : {
759 515 : GEN P = gel(F,1);
760 515 : long n = lg(P)-1;
761 515 : E = gel(F,2);
762 515 : if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
763 : }
764 630 : else if (lgefint(n) == 3)
765 601 : return bigomegau(n[2]);
766 : else
767 29 : E = gel(absZ_factor(n), 2);
768 544 : E = ZV_to_zv(E);
769 544 : return gc_long(av, zv_sum(E));
770 : }
771 :
772 : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
773 : ulong
774 520350 : eulerphiu_fact(GEN f)
775 : {
776 520350 : GEN P = gel(f,1), E = gel(f,2);
777 520350 : long i, m = 1, l = lg(P);
778 1474846 : for (i = 1; i < l; i++)
779 : {
780 954496 : ulong p = P[i], e = E[i];
781 954496 : if (!e) continue;
782 954485 : if (p == 2)
783 275342 : { if (e > 1) m <<= e-1; }
784 : else
785 : {
786 679143 : m *= (p-1);
787 679143 : if (e > 1) m *= upowuu(p, e-1);
788 : }
789 : }
790 520350 : return m;
791 : }
792 : ulong
793 167857 : eulerphiu(ulong n)
794 : {
795 : pari_sp av;
796 167857 : if (!n) return 2;
797 167857 : av = avma; return gc_long(av, eulerphiu_fact(factoru(n)));
798 : }
799 : GEN
800 108566 : eulerphi(GEN n)
801 : {
802 108566 : pari_sp av = avma;
803 : GEN Q, F, P, E;
804 : long i, l;
805 :
806 108566 : if ((F = check_arith_all(n,"eulerphi")))
807 : {
808 2970 : F = clean_Z_factor(F);
809 2970 : n = arith_n(n);
810 2970 : if (lgefint(n) == 3)
811 : {
812 : ulong e;
813 2960 : F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
814 2960 : e = eulerphiu_fact(F);
815 2960 : return gc_utoipos(av, e);
816 : }
817 : }
818 105596 : else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
819 : else
820 43 : F = absZ_factor(n);
821 53 : if (!signe(n)) return gen_2;
822 33 : P = gel(F,1);
823 33 : E = gel(F,2); l = lg(P);
824 33 : Q = cgetg(l, t_VEC);
825 216 : for (i = 1; i < l; i++)
826 : {
827 183 : GEN p = gel(P,i), q;
828 183 : ulong v = itou(gel(E,i));
829 183 : q = subiu(p,1);
830 183 : if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
831 183 : gel(Q,i) = q;
832 : }
833 33 : return gc_INT(av, ZV_prod(Q));
834 : }
835 :
836 : long
837 8248231 : numdivu_fact(GEN fa)
838 : {
839 8248231 : GEN E = gel(fa,2);
840 8248231 : long n = 1, i, l = lg(E);
841 27869206 : for (i = 1; i < l; i++) n *= E[i]+1;
842 8248231 : return n;
843 : }
844 : long
845 597 : numdivu(long N)
846 : {
847 : pari_sp av;
848 597 : if (N == 1) return 1;
849 587 : av = avma; return gc_long(av, numdivu_fact(factoru(N)));
850 : }
851 : static GEN
852 548 : numdiv_aux(GEN F)
853 : {
854 548 : GEN x, E = gel(F,2);
855 548 : long i, l = lg(E);
856 548 : x = cgetg(l, t_VECSMALL);
857 1521 : for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
858 548 : return x;
859 : }
860 : GEN
861 1155 : numdiv(GEN n)
862 : {
863 1155 : pari_sp av = avma;
864 : GEN F, E;
865 1155 : if ((F = check_arith_non0(n,"numdiv")))
866 : {
867 515 : F = clean_Z_factor(F);
868 515 : E = numdiv_aux(F);
869 : }
870 630 : else if (lgefint(n) == 3)
871 597 : return utoipos(numdivu(n[2]));
872 : else
873 33 : E = numdiv_aux(absZ_factor(n));
874 548 : return gc_INT(av, zv_prod_Z(E));
875 : }
876 :
877 : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
878 : static GEN
879 61080 : u_euler_sumdiv(ulong p, long v)
880 : {
881 61080 : GEN u = utoipos(1 + p); /* can't overflow */
882 83760 : for (; v > 1; v--) u = addui(1, mului(p, u));
883 61080 : return u;
884 : }
885 : /* 1 + q + ... + q^v */
886 : static GEN
887 15055 : euler_sumdiv(GEN q, long v)
888 : {
889 15055 : GEN u = addui(1, q);
890 20416 : for (; v > 1; v--) u = addui(1, mulii(q, u));
891 15055 : return u;
892 : }
893 :
894 : static GEN
895 1058 : sumdiv_aux(GEN F)
896 : {
897 1058 : GEN x, P = gel(F,1), E = gel(F,2);
898 1058 : long i, l = lg(P);
899 1058 : x = cgetg(l, t_VEC);
900 2816 : for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
901 1058 : return ZV_prod(x);
902 : }
903 : GEN
904 2175 : sumdiv(GEN n)
905 : {
906 2175 : pari_sp av = avma;
907 : GEN F, v;
908 :
909 2175 : if ((F = check_arith_non0(n,"sumdiv")))
910 : {
911 1030 : F = clean_Z_factor(F);
912 1030 : v = sumdiv_aux(F);
913 : }
914 1135 : else if (lgefint(n) == 3)
915 : {
916 1107 : if (n[2] == 1) return gen_1;
917 1087 : F = factoru(n[2]);
918 1087 : v = usumdiv_fact(F);
919 : }
920 : else
921 28 : v = sumdiv_aux(absZ_factor(n));
922 2145 : return gc_INT(av, v);
923 : }
924 :
925 : static GEN
926 1086 : sumdivk_aux(GEN F, long k)
927 : {
928 1086 : GEN x, P = gel(F,1), E = gel(F,2);
929 1086 : long i, l = lg(P);
930 1086 : x = cgetg(l, t_VEC);
931 3012 : for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(powiu(gel(P,i),k), gel(E,i)[2]);
932 1086 : return ZV_prod(x);
933 : }
934 : GEN
935 6019 : sumdivk(GEN n, long k)
936 : {
937 6019 : pari_sp av = avma;
938 : GEN F, v;
939 : long k1;
940 :
941 6019 : if (!k) return numdiv(n);
942 6019 : if (k == 1) return sumdiv(n);
943 4864 : if ((F = check_arith_non0(n,"sumdivk"))) F = clean_Z_factor(F);
944 4834 : k1 = k; if (k < 0) k = -k;
945 4834 : if (k == 1)
946 1020 : v = sumdiv(F? F: n);
947 : else
948 : {
949 3814 : if (F)
950 1030 : v = sumdivk_aux(F,k);
951 2784 : else if (lgefint(n) == 3)
952 : {
953 2728 : if (n[2] == 1) return gen_1;
954 2648 : F = factoru(n[2]);
955 2648 : v = usumdivk_fact(F,k);
956 : }
957 : else
958 56 : v = sumdivk_aux(absZ_factor(n), k);
959 3734 : if (k1 > 0) return gc_INT(av, v);
960 : }
961 :
962 1025 : if (F) n = arith_n(n);
963 1025 : if (k != 1) n = powiu(n,k);
964 1025 : return gc_upto(av, gdiv(v, n));
965 : }
966 :
967 : GEN
968 37207 : usumdiv_fact(GEN f)
969 : {
970 37207 : GEN P = gel(f,1), E = gel(f,2);
971 37207 : long i, l = lg(P);
972 37207 : GEN v = cgetg(l, t_VEC);
973 98287 : for (i=1; i<l; i++) gel(v,i) = u_euler_sumdiv(P[i],E[i]);
974 37207 : return ZV_prod(v);
975 : }
976 : GEN
977 8182 : usumdivk_fact(GEN f, ulong k)
978 : {
979 8182 : GEN P = gel(f,1), E = gel(f,2);
980 8182 : long i, l = lg(P);
981 8182 : GEN v = cgetg(l, t_VEC);
982 19553 : for (i=1; i<l; i++) gel(v,i) = euler_sumdiv(powuu(P[i],k),E[i]);
983 8182 : return ZV_prod(v);
984 : }
985 :
986 : long
987 2493 : uissquarefree_fact(GEN f)
988 : {
989 2493 : GEN E = gel(f,2);
990 2493 : long i, l = lg(E);
991 2493 : if (l == 2) return umael(f,1,1)? E[1] == 1: 0; /* handle factor(0) */
992 4977 : for (i = 1; i < l; i++)
993 3469 : if (E[i] > 1) return 0;
994 1508 : return 1;
995 : }
996 : long
997 2162572 : uissquarefree(ulong n)
998 : {
999 2162572 : if (!n) return 0;
1000 2162572 : return moebiusu(n)? 1: 0;
1001 : }
1002 : long
1003 2004 : Z_issquarefree(GEN n)
1004 : {
1005 2004 : switch(lgefint(n))
1006 : {
1007 12 : case 2: return 0;
1008 1219 : case 3: return uissquarefree(n[2]);
1009 : }
1010 773 : return moebius(n)? 1: 0;
1011 : }
1012 :
1013 : long
1014 2211 : Z_issquarefree_fact(GEN F)
1015 : {
1016 2211 : GEN E = gel(F,2);
1017 2211 : long i, l = lg(E);
1018 2211 : if (l == 2) return signe(gcoeff(F,1,1))? equali1(gel(E,1)): 0;
1019 4950 : for(i = 1; i < l; i++)
1020 3894 : if (!equali1(gel(E,i))) return 0;
1021 1056 : return 1;
1022 : }
1023 : long
1024 5468 : issquarefree(GEN x)
1025 : {
1026 : pari_sp av;
1027 : GEN d;
1028 5468 : switch(typ(x))
1029 : {
1030 1217 : case t_INT: return Z_issquarefree(x);
1031 3045 : case t_POL:
1032 3045 : if (!signe(x)) return 0;
1033 3045 : av = avma; d = RgX_gcd(x, RgX_deriv(x));
1034 3045 : return gc_long(av, lg(d)==3);
1035 1206 : case t_VEC:
1036 1206 : 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 : }
|
