Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - arith2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 496 558 88.9 %
Date: 2022-07-03 07:33:15 Functions: 68 76 89.5 %
Legend: Lines: hit not hit

          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          49 : boundfact(GEN n, ulong lim)
      26             : {
      27          49 :   switch(typ(n))
      28             :   {
      29          35 :     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, merge_factor(a,b,(void*)&cmpii,cmp_nodata));
      36             :     }
      37             :   }
      38           0 :   pari_err_TYPE("boundfact",n);
      39             :   return NULL; /* LCOV_EXCL_LINE */
      40             : }
      41             : 
      42             : /* NOT memory clean */
      43             : GEN
      44       37104 : Z_lsmoothen(GEN N, GEN L, GEN *pP, GEN *pE)
      45             : {
      46       37104 :   long i, j, l = lg(L);
      47       37104 :   GEN E = new_chunk(l), P = new_chunk(l);
      48      391918 :   for (i = j = 1; i < l; i++)
      49             :   {
      50      373589 :     ulong p = uel(L,i);
      51      373589 :     long v = Z_lvalrem(N, p, &N);
      52      373589 :     if (v) { P[j] = p; E[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
      53             :   }
      54       37104 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); if (pP) *pP = P;
      55       37104 :   E[0] = evaltyp(t_VECSMALL) | evallg(j); if (pE) *pE = E; return N;
      56             : }
      57             : GEN
      58       27105 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pE)
      59             : {
      60       27105 :   long i, j, l = lg(L);
      61       27105 :   GEN E = new_chunk(l), P = new_chunk(l);
      62       45064 :   for (i = j = 1; i < l; i++)
      63             :   {
      64       44974 :     GEN p = gel(L,i);
      65       44974 :     long v = Z_pvalrem(N, p, &N);
      66       44977 :     if (v)
      67             :     {
      68       28096 :       gel(P,j) = p; gel(E,j) = utoipos(v); j++;
      69       28100 :      if (is_pm1(N)) { N = NULL; break; }
      70             :     }
      71             :   }
      72       27160 :   P[0] = evaltyp(t_COL) | evallg(j); if (pP) *pP = P;
      73       27133 :   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      122100 : factoru_pow(ulong n)
      82             : {
      83      122100 :   GEN f = cgetg(4,t_VEC);
      84      122100 :   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      122100 :   (void)new_chunk((15 + 1)*3);
      89      122101 :   F = factoru(n);
      90      122101 :   P = gel(F,1);
      91      122101 :   E = gel(F,2); l = lg(P);
      92      122101 :   set_avma(av);
      93      122101 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
      94      122103 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
      95      122103 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
      96      280556 :   for(i = 1; i < l; i++)
      97             :   {
      98      158454 :     p[i] = P[i];
      99      158454 :     e[i] = E[i];
     100      158454 :     c[i] = upowuu(p[i], e[i]);
     101             :   }
     102      122102 :   return f;
     103             : }
     104             : 
     105             : static GEN
     106       16254 : factorlim(GEN n, ulong lim)
     107       16254 : { 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        7385 : factor_pn_1_limit(GEN p, long n, ulong lim)
     112             : {
     113        7385 :   pari_sp av = avma;
     114        7385 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
     115        7385 :   long i, pp = itos_or_0(p);
     116       15911 :   for(i=2; i<lg(d); i++)
     117             :   {
     118             :     GEN B;
     119        8526 :     if (pp && d[i]%pp==0 && (
     120        3444 :        ((pp&3)==1 && (d[i]&1)) ||
     121        3423 :        ((pp&3)==3 && (d[i]&3)==2) ||
     122        3304 :        (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 = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     127         343 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     128             :     }
     129             :     else
     130        8183 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     131        8526 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     132             :   }
     133        7385 :   return gerepilecopy(av, A);
     134             : }
     135             : GEN
     136        7385 : factor_pn_1(GEN p, ulong n)
     137        7385 : { 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? merge_factor(B, fa, (void*)&cmpii, cmp_nodata): 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    12962313 : RgV_is_ZVpos(GEN v)
     201             : {
     202    12962313 :   long i, l = lg(v);
     203    39362512 :   for (i = 1; i < l; i++)
     204             :   {
     205    26406058 :     GEN c = gel(v,i);
     206    26406058 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     207             :   }
     208    12956454 :   return 1;
     209             : }
     210             : /* check whether v is a ZV with nonzero entries */
     211             : int
     212       42789 : RgV_is_ZVnon0(GEN v)
     213             : {
     214       42789 :   long i, l = lg(v);
     215       99474 :   for (i = 1; i < l; i++)
     216             :   {
     217       56741 :     GEN c = gel(v,i);
     218       56741 :     if (typ(c) != t_INT || !signe(c)) return 0;
     219             :   }
     220       42733 :   return 1;
     221             : }
     222             : /* check whether v is a ZV with nonzero entries OR exactly [0] */
     223             : static int
     224      663834 : RgV_is_ZV0(GEN v)
     225             : {
     226      663834 :   long i, l = lg(v);
     227     2326252 :   for (i = 1; i < l; i++)
     228             :   {
     229     1662530 :     GEN c = gel(v,i);
     230             :     long s;
     231     1662530 :     if (typ(c) != t_INT) return 0;
     232     1662530 :     s = signe(c);
     233     1662530 :     if (!s) return (l == 2);
     234             :   }
     235      663722 :   return 1;
     236             : }
     237             : 
     238             : static int
     239      214236 : RgV_is_prV(GEN v)
     240             : {
     241      214236 :   long l = lg(v), i;
     242      251917 :   for (i = 1; i < l; i++)
     243      212121 :     if (!checkprid_i(gel(v,i))) return 0;
     244       39796 :   return 1;
     245             : }
     246             : int
     247      257937 : is_nf_factor(GEN F)
     248             : {
     249      239296 :   return typ(F) == t_MAT && lg(F) == 3
     250      497233 :     && 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     6811365 : is_Z_factor_i(GEN f)
     261     6811365 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
     262             : int
     263     6104750 : is_Z_factorpos(GEN f)
     264     6104750 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     265             : int
     266      663826 : is_Z_factor(GEN f)
     267      663826 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     268             : /* as is_Z_factorpos, also allow factor(0) */
     269             : int
     270       42791 : is_Z_factornon0(GEN f)
     271       42791 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     272             : GEN
     273       52381 : clean_Z_factor(GEN f)
     274             : {
     275       52381 :   GEN P = gel(f,1);
     276       52381 :   long n = lg(P)-1;
     277       52381 :   if (n && equalim1(gel(P,1)))
     278        3066 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     279       49315 :   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         602 : check_arith_pos(GEN n, const char *f) {
     302         602 :   switch(typ(n))
     303             :   {
     304         553 :     case t_INT:
     305         553 :       if (signe(n) <= 0) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
     306         553 :       return NULL;
     307          20 :     case t_VEC:
     308          20 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT || signe(gel(n,1)) <= 0) break;
     309          20 :       n = gel(n,2); /* fall through */
     310          42 :     case t_MAT:
     311          42 :       if (!is_Z_factorpos(n)) break;
     312          42 :       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     2844811 : check_arith_non0(GEN n, const char *f) {
     320     2844811 :   switch(typ(n))
     321             :   {
     322     2802411 :     case t_INT:
     323     2802411 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     324     2802369 :       return NULL;
     325       27685 :     case t_VEC:
     326       27685 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT || !signe(gel(n,1))) break;
     327       27685 :       n = gel(n,2); /* fall through */
     328       42553 :     case t_MAT:
     329       42553 :       if (!is_Z_factornon0(n)) break;
     330       42504 :       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      980581 : check_arith_all(GEN n, const char *f) {
     338      980581 :   switch(typ(n))
     339             :   {
     340      316743 :     case t_INT:
     341      316743 :       return NULL;
     342      657986 :     case t_VEC:
     343      657986 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     344      657985 :       n = gel(n,2); /* fall through */
     345      663830 :     case t_MAT:
     346      663830 :       if (!is_Z_factor(n)) break;
     347      663843 :       return n;
     348             :   }
     349           8 :   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      362299 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     360             : {
     361             :   GEN E, P;
     362      362299 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     363      362299 :   P = gel(F,1);
     364      362299 :   E = gel(F,2);
     365      362299 :   RgV_check_ZV(E, "divisors");
     366      362299 :   *isint = RgV_is_ZV(P);
     367      362299 :   if (*isint)
     368             :   {
     369      362285 :     long i, l = lg(P);
     370             :     /* skip -1 */
     371      362285 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     372             :     /* test for 0 */
     373     1174383 :     for (i = 1; i < l; i++)
     374      812105 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     375           7 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     376             :   }
     377      362292 :   *pP = P;
     378      362292 :   *pE = E;
     379      362292 : }
     380             : static void
     381       16219 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     382             : 
     383             : int
     384      378525 : divisors_init(GEN n, GEN *pP, GEN *pE)
     385             : {
     386             :   long i,l;
     387             :   GEN E, P, e;
     388             :   int isint;
     389             : 
     390      378525 :   switch(typ(n))
     391             :   {
     392       16198 :     case t_INT:
     393       16198 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     394       16198 :       set_fact(absZ_factor(n), &P,&E);
     395       16198 :       isint = 1; break;
     396        1330 :     case t_VEC:
     397        1330 :       if (lg(n) != 3 || typ(gel(n,2)) !=t_MAT) pari_err_TYPE("divisors",n);
     398        1323 :       set_fact_check(gel(n,2), &P,&E, &isint);
     399        1323 :       break;
     400      360976 :     case t_MAT:
     401      360976 :       set_fact_check(n, &P,&E, &isint);
     402      360969 :       break;
     403          21 :     default:
     404          21 :       set_fact(factor(n), &P,&E);
     405          21 :       isint = 0; break;
     406             :   }
     407      378511 :   l = lg(P);
     408      378511 :   e = cgetg(l, t_VECSMALL);
     409     1231461 :   for (i=1; i<l; i++)
     410             :   {
     411      852957 :     e[i] = itos(gel(E,i));
     412      852957 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     413             :   }
     414      378504 :   *pP = P; *pE = e; return isint;
     415             : }
     416             : 
     417             : static long
     418      378455 : ndiv(GEN E)
     419             : {
     420             :   long i, l;
     421      378455 :   GEN e = cgetg_copy(E, &l); /* left on stack */
     422             :   ulong n;
     423     1231230 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     424      378455 :   n = itou_or_0( zv_prod_Z(e) );
     425      378455 :   if (!n || n & ~LGBITS) pari_err_OVERFLOW("divisors");
     426      378455 :   return n;
     427             : }
     428             : static int
     429       15904 : 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       23219 : fa_clean(GEN P, GEN E)
     434             : {
     435       23219 :   long i, j, l = lg(E);
     436       23219 :   GEN Q = cgetg(l, t_COL);
     437       66220 :   for (i = j = 1; i < l; i++)
     438       43001 :     if (E[i]) { gel(Q,j) = gel(P,i); E[j] = E[i]; j++; }
     439       23219 :   setlg(Q,j);
     440       23219 :   setlg(E,j); return mkmat2(Q,Flc_to_ZC(E));
     441             : }
     442             : GEN
     443       12110 : divisors_factored(GEN N)
     444             : {
     445       12110 :   pari_sp av = avma;
     446             :   GEN *d, *t1, *t2, *t3, D, P, E;
     447       12110 :   int isint = divisors_init(N, &P, &E);
     448       12110 :   GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
     449       12110 :   long i, j, l, n = ndiv(E);
     450             : 
     451       12110 :   D = cgetg(n+1,t_VEC); d = (GEN*)D;
     452       12110 :   l = lg(E);
     453       12110 :   *++d = mkvec2(gen_1, const_vecsmall(l-1,0));
     454       33208 :   for (i=1; i<l; i++)
     455       29778 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     456       19789 :       for (t2=d, t3=t1; t3<t2; )
     457             :       {
     458             :         GEN a, b;
     459       11109 :         a = mul(gel(*++t3,1), gel(P,i));
     460       11109 :         b = leafcopy(gel(*t3,2)); b[i]++;
     461       11109 :         *++d = mkvec2(a,b);
     462             :       }
     463       12110 :   if (isint) gen_sort_inplace(D,NULL,&cmpi1,NULL);
     464       35329 :   for (i = 1; i <= n; i++) gmael(D,i,2) = fa_clean(P, gmael(D,i,2));
     465       12110 :   return gerepilecopy(av, D);
     466             : }
     467             : static int
     468           0 : cmpu1(void *E, GEN va, GEN vb)
     469           0 : { long a = va[1], b = vb[1]; (void)E; return a>b? 1: (a<b? -1: 0); }
     470             : static GEN
     471           0 : fa_clean_u(GEN P, GEN E)
     472             : {
     473           0 :   long i, j, l = lg(E);
     474           0 :   GEN Q = cgetg(l, t_VECSMALL);
     475           0 :   for (i = j = 1; i < l; i++)
     476           0 :     if (E[i]) { Q[j] = P[i]; E[j] = E[i]; j++; }
     477           0 :   setlg(Q,j);
     478           0 :   setlg(E,j); return mkmat2(Q,E);
     479             : }
     480             : GEN
     481           0 : divisorsu_fact_factored(GEN fa)
     482             : {
     483           0 :   pari_sp av = avma;
     484           0 :   GEN *d, *t1, *t2, *t3, vD, D, P = gel(fa,1), E = gel(fa,2);
     485           0 :   long i, j, l, n = ndiv(E);
     486             : 
     487           0 :   D = cgetg(n+1,t_VEC); d = (GEN*)D;
     488           0 :   l = lg(E);
     489           0 :   *++d = mkvec2((GEN)1, const_vecsmall(l-1,0));
     490           0 :   for (i=1; i<l; i++)
     491           0 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     492           0 :       for (t2=d, t3=t1; t3<t2; )
     493             :       {
     494             :         ulong a;
     495             :         GEN b;
     496           0 :         a = (*++t3)[1] * P[i];
     497           0 :         b = leafcopy(gel(*t3,2)); b[i]++;
     498           0 :         *++d = mkvec2((GEN)a,b);
     499             :       }
     500           0 :   gen_sort_inplace(D,NULL,&cmpu1,NULL);
     501           0 :   vD = cgetg(n+1, t_VECSMALL);
     502           0 :   for (i = 1; i <= n; i++)
     503             :   {
     504           0 :     vD[i] = umael(D,i,1);
     505           0 :     gel(D,i) = fa_clean_u(P, gmael(D,i,2));
     506             :   }
     507           0 :   return gerepilecopy(av, mkvec2(vD,D));
     508             : }
     509             : GEN
     510      366366 : divisors(GEN N)
     511             : {
     512             :   long i, j, l;
     513             :   GEN *d, *t1, *t2, *t3, D, P, E;
     514      366366 :   int isint = divisors_init(N, &P, &E);
     515      366345 :   GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
     516             : 
     517      366345 :   D = cgetg(ndiv(E)+1,t_VEC); d = (GEN*)D;
     518      366345 :   l = lg(E);
     519      366345 :   *++d = gen_1;
     520     1198022 :   for (i=1; i<l; i++)
     521     2231803 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     522     5402194 :       for (t2=d, t3=t1; t3<t2; ) *++d = mul(*++t3, gel(P,i));
     523      366345 :   if (isint) ZV_sort_inplace(D);
     524      366345 :   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       12670 : divisorsu_moebius(GEN P)
     535             : {
     536             :   GEN d, t, t2, t3;
     537       12670 :   long i, l = lg(P);
     538       12670 :   d = t = cgetg((1 << (l-1)) + 1, t_VECSMALL);
     539       12670 :   *++d = 1;
     540       29316 :   for (i=1; i<l; i++)
     541       39662 :     for (t2=d, t3=t; t3<t2; ) *(++d) = *(++t3) * -P[i];
     542       12670 :   return t;
     543             : }
     544             : GEN
     545    25949670 : divisorsu_fact(GEN fa)
     546             : {
     547    25949670 :   GEN d, t, t1, t2, t3, P = gel(fa,1), E = gel(fa,2);
     548    25949670 :   long i, j, l = lg(P);
     549    25949670 :   d = t = cgetg(numdivu_fact(fa) + 1,t_VECSMALL);
     550    28934065 :   *++d = 1;
     551    95746150 :   for (i=1; i<l; i++)
     552   145709592 :     for (t1=t,j=E[i]; j; j--,t1=t2)
     553   283947570 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     554    28934065 :   vecsmall_sort(t); return t;
     555             : }
     556             : GEN
     557      144374 : divisorsu(ulong N)
     558             : {
     559      144374 :   pari_sp av = avma;
     560      144374 :   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       56077 : arith_n(GEN A)
     602             : {
     603       56077 :   switch(typ(A))
     604             :   {
     605       52710 :     case t_INT: return A;
     606        1190 :     case t_VEC: return gel(A,1);
     607        2177 :     default: return factorback(A);
     608             :   }
     609             : }
     610             : static GEN
     611       54089 : core2_i(GEN n)
     612             : {
     613       54089 :   GEN f = core(n);
     614       54089 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     615       54054 :   return mkvec2(f, sqrtint(diviiexact(arith_n(n), f)));
     616             : }
     617             : GEN
     618       54047 : 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        5796 : _mod4(GEN c) {
     625        5796 :   long r, s = signe(c);
     626        5796 :   if (!s) return 0;
     627        5796 :   r = mod4(c); if (s < 0) r = 4-r;
     628        5796 :   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        5754 : coredisc(GEN n)
     642             : {
     643        5754 :   pari_sp av = avma;
     644        5754 :   GEN c = core(n);
     645        5754 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     646        2065 :   return gerepileuptoint(av, shifti(c,2));
     647             : }
     648             : 
     649             : GEN
     650          42 : coredisc2(GEN n)
     651             : {
     652          42 :   pari_sp av = avma;
     653          42 :   GEN y = core2_i(n);
     654          42 :   GEN c = gel(y,1), f = gel(y,2);
     655          42 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
     656           7 :   y = cgetg(3,t_VEC);
     657           7 :   gel(y,1) = shifti(c,2);
     658           7 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
     659             : }
     660             : 
     661             : GEN
     662          14 : 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       60776 : coredisc2_fact(GEN fa, long s, GEN *pP, GEN *pE)
     668             : {
     669       60776 :   GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2), D = s > 0? gen_1: gen_m1;
     670       60776 :   long l = lg(P0), i, j;
     671             : 
     672       60776 :   E = cgetg(l, t_VECSMALL);
     673       60776 :   P = cgetg(l, t_VEC);
     674      252287 :   for (i = j = 1; i < l; i++)
     675             :   {
     676      191516 :     long e = itos(gel(E0,i));
     677      191531 :     GEN p = gel(P0,i);
     678      191531 :     if (odd(e)) D = mulii(D, p);
     679      191493 :     e >>= 1; if (e) { gel(P, j) = p; E[j] = e; j++; }
     680             :   }
     681       60771 :   if (Mod4(D) != 1)
     682             :   {
     683       17003 :     D = shifti(D, 2);
     684       17003 :     if (!--E[1])
     685             :     {
     686       16996 :       P[1] = P[0]; P++;
     687       16996 :       E[1] = E[0]; E++; j--;
     688             :     }
     689             :   }
     690       60775 :   setlg(P,j); *pP = P;
     691       60775 :   setlg(E,j); *pE = E; return D;
     692             : }
     693             : ulong
     694        4358 : coredisc2u_fact(GEN fa, long s, GEN *pP, GEN *pE)
     695             : {
     696        4358 :   GEN P, E, P0 = gel(fa,1), E0 = gel(fa,2);
     697        4358 :   ulong D = 1;
     698        4358 :   long i, j, l = lg(P0);
     699             : 
     700        4358 :   E = cgetg(l, t_VECSMALL);
     701        4358 :   P = cgetg(l, t_VECSMALL);
     702       11702 :   for (i = j = 1; i < l; i++)
     703             :   {
     704        7344 :     long e = E0[i], p = P0[i];
     705        7344 :     if (odd(e)) D *= p;
     706        7344 :     e >>= 1; if (e) { P[j] = p; E[j] = e; j++; }
     707             :   }
     708        4358 :   if ((D & 3) != (s > 0? 1: 3))
     709             :   {
     710         541 :     D *= 4;
     711         541 :     if (!--E[1])
     712             :     {
     713         506 :       P[1] = P[0]; P++;
     714         506 :       E[1] = E[0]; E++; j--;
     715             :     }
     716             :   }
     717        4358 :   setlg(P,j); *pP = P;
     718        4358 :   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         828 : bigomegau(ulong n)
     747             : {
     748             :   pari_sp av;
     749         828 :   if (n == 1) return 0;
     750         814 :   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      272541 : eulerphiu_fact(GEN f)
     775             : {
     776      272541 :   GEN P = gel(f,1), E = gel(f,2);
     777      272541 :   long i, m = 1, l = lg(P);
     778      894180 :   for (i = 1; i < l; i++)
     779             :   {
     780      621639 :     ulong p = P[i], e = E[i];
     781      621639 :     if (!e) continue;
     782      621625 :     if (p == 2)
     783      207413 :     { if (e > 1) m <<= e-1; }
     784             :     else
     785             :     {
     786      414212 :       m *= (p-1);
     787      414212 :       if (e > 1) m *= upowuu(p, e-1);
     788             :     }
     789             :   }
     790      272541 :   return m;
     791             : }
     792             : ulong
     793      174030 : eulerphiu(ulong n)
     794             : {
     795             :   pari_sp av;
     796      174030 :   if (!n) return 2;
     797      174030 :   av = avma; return gc_long(av, eulerphiu_fact(factoru(n)));
     798             : }
     799             : GEN
     800      146622 : eulerphi(GEN n)
     801             : {
     802      146622 :   pari_sp av = avma;
     803             :   GEN Q, F, P, E;
     804             :   long i, l;
     805             : 
     806      146622 :   if ((F = check_arith_all(n,"eulerphi")))
     807             :   {
     808        1295 :     F = clean_Z_factor(F);
     809        1295 :     n = arith_n(n);
     810        1295 :     if (lgefint(n) == 3)
     811             :     {
     812             :       ulong e;
     813        1281 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
     814        1281 :       e = eulerphiu_fact(F);
     815        1281 :       set_avma(av); return utoipos(e);
     816             :     }
     817             :   }
     818      145327 :   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    25994363 : numdivu_fact(GEN fa)
     838             : {
     839    25994363 :   GEN E = gel(fa,2);
     840    25994363 :   long n = 1, i, l = lg(E);
     841    88569187 :   for (i = 1; i < l; i++) n *= E[i]+1;
     842    25994363 :   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       54954 : u_euler_sumdiv(ulong p, long v)
     880             : {
     881       54954 :   GEN u = utoipos(1 + p); /* can't overflow */
     882       75184 :   for (; v > 1; v--) u = addui(1, mului(p, u));
     883       54954 :   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       33848 : usumdiv_fact(GEN f)
     969             : {
     970       33848 :   GEN P = gel(f,1), E = gel(f,2);
     971       33848 :   long i, l = lg(P);
     972       33848 :   GEN v = cgetg(l, t_VEC);
     973       88802 :   for (i=1; i<l; i++) gel(v,i) = u_euler_sumdiv(P[i],E[i]);
     974       33848 :   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        2842 : uissquarefree_fact(GEN f)
     988             : {
     989        2842 :   GEN E = gel(f,2);
     990        2842 :   long i, l = lg(E);
     991        2842 :   if (l == 2) return umael(f,1,1)? E[1] == 1: 0; /* handle factor(0) */
     992        5607 :   for (i = 1; i < l; i++)
     993        3913 :     if (E[i] > 1) return 0;
     994        1694 :   return 1;
     995             : }
     996             : long
     997     2451954 : uissquarefree(ulong n)
     998             : {
     999     2451954 :   if (!n) return 0;
    1000     2451954 :   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        4235 : issquarefree(GEN x)
    1025             : {
    1026             :   pari_sp av;
    1027             :   GEN d;
    1028        4235 :   switch(typ(x))
    1029             :   {
    1030        1421 :     case t_INT: return Z_issquarefree(x);
    1031        1407 :     case t_POL:
    1032        1407 :       if (!signe(x)) return 0;
    1033        1407 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
    1034        1407 :       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             : }

Generated by: LCOV version 1.13