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 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.8.0 lcov report (development 19214-1621e44) Lines: 543 589 92.2 %
Date: 2016-07-26 07:10:39 Functions: 76 79 96.2 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*********************************************************************/
      15             : /**                                                                 **/
      16             : /**                     ARITHMETIC FUNCTIONS                        **/
      17             : /**                        (second part)                            **/
      18             : /**                                                                 **/
      19             : /*********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : GEN
      24          35 : boundfact(GEN n, ulong lim)
      25             : {
      26          35 :   switch(typ(n))
      27             :   {
      28          21 :     case t_INT: return Z_factor_limit(n,lim);
      29             :     case t_FRAC: {
      30          14 :       pari_sp av = avma;
      31          14 :       GEN a = Z_factor_limit(gel(n,1),lim);
      32          14 :       GEN b = Z_factor_limit(gel(n,2),lim);
      33          14 :       gel(b,2) = ZC_neg(gel(b,2));
      34          14 :       return gerepilecopy(av, merge_factor_i(a,b));
      35             :     }
      36             :   }
      37           0 :   pari_err_TYPE("boundfact",n);
      38           0 :   return NULL; /* not reached */
      39             : }
      40             : 
      41             : /* NOT memory clean */
      42             : GEN
      43       12460 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
      44             : {
      45       12460 :   long i, j, l = lg(L);
      46       12460 :   GEN e = new_chunk(l), P = new_chunk(l);
      47       87500 :   for (i = j = 1; i < l; i++)
      48             :   {
      49       83825 :     ulong p = uel(L,i);
      50       83825 :     long v = Z_lvalrem(N, p, &N);
      51       83825 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
      52             :   }
      53       12460 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
      54       12460 :   e[0] = evaltyp(t_VECSMALL) | evallg(j); *pe = e; return N;
      55             : }
      56             : /***********************************************************************/
      57             : /**                                                                   **/
      58             : /**                    SIMPLE FACTORISATIONS                          **/
      59             : /**                                                                   **/
      60             : /***********************************************************************/
      61             : /* Factor n and output [p,e,c] where
      62             :  * p, e and c are vecsmall with n = prod{p[i]^e[i]} and c[i] = p[i]^e[i] */
      63             : GEN
      64       28290 : factoru_pow(ulong n)
      65             : {
      66       28290 :   GEN f = cgetg(4,t_VEC);
      67       28289 :   pari_sp av = avma;
      68             :   GEN F, P, E, p, e, c;
      69             :   long i, l;
      70             :   /* enough room to store <= 15 * [p,e,p^e] (OK if n < 2^64) */
      71       28289 :   (void)new_chunk((15 + 1)*3);
      72       28289 :   F = factoru(n);
      73       28288 :   P = gel(F,1);
      74       28288 :   E = gel(F,2); l = lg(P);
      75       28288 :   avma = av;
      76       28288 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
      77       28289 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
      78       28289 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
      79       89605 :   for(i = 1; i < l; i++)
      80             :   {
      81       61315 :     p[i] = P[i];
      82       61315 :     e[i] = E[i];
      83       61315 :     c[i] = upowuu(p[i], e[i]);
      84             :   }
      85       28290 :   return f;
      86             : }
      87             : 
      88             : static GEN
      89       80517 : factorlim(GEN n, ulong lim)
      90       80517 : { return lim? Z_factor_limit(n, lim): Z_factor(n); }
      91             : /* factor p^n - 1, assuming p prime. If lim != 0, limit factorization to
      92             :  * primes <= lim */
      93             : GEN
      94       62835 : factor_pn_1_limit(GEN p, long n, ulong lim)
      95             : {
      96       62835 :   pari_sp av = avma;
      97       62835 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
      98       62835 :   long i, pp = itos_or_0(p);
      99       75337 :   for(i=2; i<lg(d); i++)
     100             :   {
     101             :     GEN B;
     102       23373 :     if (pp && d[i]%pp==0 && (
     103       21742 :        ((pp&3)==1 && (d[i]&1)) ||
     104       10962 :        ((pp&3)==3 && (d[i]&3)==2) ||
     105       10766 :        (pp==2 && (d[i]&7)==4)))
     106        5180 :     {
     107        5180 :       GEN f=factor_Aurifeuille_prime(p,d[i]);
     108        5180 :       B = factorlim(f, lim);
     109        5180 :       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     110        5180 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     111             :     }
     112             :     else
     113        7322 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     114       12502 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     115             :   }
     116       62835 :   return gerepilecopy(av, A);
     117             : }
     118             : GEN
     119       62835 : factor_pn_1(GEN p, ulong n)
     120       62835 : { return factor_pn_1_limit(p, n, 0); }
     121             : 
     122             : #if 0
     123             : static GEN
     124             : to_mat(GEN p, long e) {
     125             :   GEN B = cgetg(3, t_MAT);
     126             :   gel(B,1) = mkcol(p);
     127             :   gel(B,2) = mkcol(utoipos(e)); return B;
     128             : }
     129             : /* factor phi(n) */
     130             : GEN
     131             : factor_eulerphi(GEN n)
     132             : {
     133             :   pari_sp av = avma;
     134             :   GEN B = NULL, A, P, E, AP, AE;
     135             :   long i, l, v = vali(n);
     136             : 
     137             :   l = lg(n);
     138             :   /* result requires less than this: at most expi(n) primes */
     139             :   (void)new_chunk(bit_accuracy(l) * (l /*p*/ + 3 /*e*/ + 2 /*vectors*/) + 3+2);
     140             :   if (v) { n = shifti(n, -v); v--; }
     141             :   A = Z_factor(n); P = gel(A,1); E = gel(A,2); l = lg(P);
     142             :   for(i = 1; i < l; i++)
     143             :   {
     144             :     GEN p = gel(P,i), q = subis(p,1), fa;
     145             :     long e = itos(gel(E,i)), w;
     146             : 
     147             :     w = vali(q); v += w; q = shifti(q,-w);
     148             :     if (! is_pm1(q))
     149             :     {
     150             :       fa = Z_factor(q);
     151             :       B = B? merge_factor(B, fa, (void*)&cmpii, cmp_nodata): fa;
     152             :     }
     153             :     if (e > 1) {
     154             :       if (B) {
     155             :         gel(B,1) = shallowconcat(gel(B,1), p);
     156             :         gel(B,2) = shallowconcat(gel(B,2), utoipos(e-1));
     157             :       } else
     158             :         B = to_mat(p, e-1);
     159             :     }
     160             :   }
     161             :   avma = av;
     162             :   if (!B) return v? to_mat(gen_2, v): trivial_fact();
     163             :   A = cgetg(3, t_MAT);
     164             :   P = gel(B,1); E = gel(B,2); l = lg(P);
     165             :   AP = cgetg(l+1, t_COL); gel(A,1) = AP; AP++;
     166             :   AE = cgetg(l+1, t_COL); gel(A,2) = AE; AE++;
     167             :   /* prepend "2^v" */
     168             :   gel(AP,0) = gen_2;
     169             :   gel(AE,0) = utoipos(v);
     170             :   for (i = 1; i < l; i++)
     171             :   {
     172             :     gel(AP,i) = icopy(gel(P,i));
     173             :     gel(AE,i) = icopy(gel(E,i));
     174             :   }
     175             :   return A;
     176             : }
     177             : #endif
     178             : 
     179             : /***********************************************************************/
     180             : /**                                                                   **/
     181             : /**         CHECK FACTORIZATION FOR ARITHMETIC FUNCTIONS              **/
     182             : /**                                                                   **/
     183             : /***********************************************************************/
     184             : static int
     185     2665473 : RgV_is_ZVpos(GEN v)
     186             : {
     187     2665473 :   long i, l = lg(v);
     188     8779170 :   for (i = 1; i < l; i++)
     189             :   {
     190     6113704 :     GEN c = gel(v,i);
     191     6113704 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     192             :   }
     193     2665466 :   return 1;
     194             : }
     195             : /* check whether v is a ZV with non-0 entries */
     196             : static int
     197        5166 : RgV_is_ZVnon0(GEN v)
     198             : {
     199        5166 :   long i, l = lg(v);
     200       15785 :   for (i = 1; i < l; i++)
     201             :   {
     202       10668 :     GEN c = gel(v,i);
     203       10668 :     if (typ(c) != t_INT || !signe(c)) return 0;
     204             :   }
     205        5117 :   return 1;
     206             : }
     207             : /* check whether v is a ZV with non-zero entries OR exactly [0] */
     208             : static int
     209        3178 : RgV_is_ZV0(GEN v)
     210             : {
     211        3178 :   long i, l = lg(v);
     212        9422 :   for (i = 1; i < l; i++)
     213             :   {
     214        6335 :     GEN c = gel(v,i);
     215             :     long s;
     216        6335 :     if (typ(c) != t_INT) return 0;
     217        6335 :     s = signe(c);
     218        6335 :     if (!s) return (l == 2);
     219             :   }
     220        3087 :   return 1;
     221             : }
     222             : 
     223             : static int
     224     1336919 : is_Z_factor_i(GEN f)
     225     1336919 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
     226             : int
     227     1328575 : is_Z_factorpos(GEN f)
     228     1328575 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     229             : int
     230        3178 : is_Z_factor(GEN f)
     231        3178 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     232             : /* as is_Z_factorpos, also allow factor(0) */
     233             : int
     234        5166 : is_Z_factornon0(GEN f)
     235        5166 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     236             : GEN
     237        4501 : clean_Z_factor(GEN f)
     238             : {
     239        4501 :   GEN P = gel(f,1);
     240        4501 :   long n = lg(P)-1;
     241        4501 :   if (n && equalim1(gel(P,1)))
     242        2156 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     243        2345 :   return f;
     244             : }
     245             : GEN
     246           0 : fuse_Z_factor(GEN f, GEN B)
     247             : {
     248           0 :   GEN P = gel(f,1), E = gel(f,2), P2,E2;
     249           0 :   long i, l = lg(P);
     250           0 :   if (l == 1) return f;
     251           0 :   for (i = 1; i < l; i++)
     252           0 :     if (absi_cmp(gel(P,i), B) > 0) break;
     253           0 :   if (i == l) return f;
     254             :   /* tail / initial segment */
     255           0 :   P2 = vecslice(P, i, l-1); P = vecslice(P, 1, i-1);
     256           0 :   E2 = vecslice(E, i, l-1); E = vecslice(E, 1, i-1);
     257           0 :   P = shallowconcat(P, mkvec(factorback2(P2,E2)));
     258           0 :   E = shallowconcat(E, mkvec(gen_1));
     259           0 :   return mkmat2(P, E);
     260             : }
     261             : 
     262             : /* n attached to a factorization of a positive integer: either N (t_INT)
     263             :  * a factorization matrix faN, or a t_VEC: [N, faN] */
     264             : GEN
     265         189 : check_arith_pos(GEN n, const char *f) {
     266         189 :   switch(typ(n))
     267             :   {
     268             :     case t_INT:
     269         189 :       if (signe(n) <= 0 ) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
     270         189 :       return NULL;
     271             :     case t_VEC:
     272           0 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     273           0 :       n = gel(n,2); /* fall through */
     274             :     case t_MAT:
     275           0 :       if (!is_Z_factorpos(n)) break;
     276           0 :       return n;
     277             :   }
     278           0 :   pari_err_TYPE(f,n);
     279           0 :   return NULL;
     280             : }
     281             : /* n attached to a factorization of a non-0 integer */
     282             : GEN
     283       14402 : check_arith_non0(GEN n, const char *f) {
     284       14402 :   switch(typ(n))
     285             :   {
     286             :     case t_INT:
     287        9285 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     288        9229 :       return NULL;
     289             :     case t_VEC:
     290           7 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     291           7 :       n = gel(n,2); /* fall through */
     292             :     case t_MAT:
     293        5117 :       if (!is_Z_factornon0(n)) break;
     294        5068 :       return n;
     295             :   }
     296          49 :   pari_err_TYPE(f,n);
     297           0 :   return NULL;
     298             : }
     299             : /* n attached to a factorization of an integer */
     300             : GEN
     301      187551 : check_arith_all(GEN n, const char *f) {
     302      187551 :   switch(typ(n))
     303             :   {
     304             :     case t_INT:
     305      184373 :       return NULL;
     306             :     case t_VEC:
     307         126 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     308         126 :       n = gel(n,2); /* fall through */
     309             :     case t_MAT:
     310        3178 :       if (!is_Z_factor(n)) break;
     311        3178 :       return n;
     312             :   }
     313           0 :   pari_err_TYPE(f,n);
     314           0 :   return NULL;
     315             : }
     316             : 
     317             : /***********************************************************************/
     318             : /**                                                                   **/
     319             : /**                MISCELLANEOUS ARITHMETIC FUNCTIONS                 **/
     320             : /**                (ultimately depend on Z_factor())                  **/
     321             : /**                                                                   **/
     322             : /***********************************************************************/
     323             : /* set P,E from F. Check whether F is an integer and kill "factor" -1 */
     324             : static void
     325      350749 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     326             : {
     327             :   GEN E, P;
     328      350749 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     329      350749 :   P = gel(F,1);
     330      350749 :   E = gel(F,2);
     331      350749 :   RgV_check_ZV(E, "divisors");
     332      350749 :   *isint = RgV_is_ZV(P);
     333      350749 :   if (*isint)
     334             :   {
     335      350735 :     long i, l = lg(P);
     336             :     /* skip -1 */
     337      350735 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     338             :     /* test for 0 */
     339     1143135 :     for (i = 1; i < l; i++)
     340      792407 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     341           7 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     342             :   }
     343      350742 :   *pP = P;
     344      350742 :   *pE = E;
     345      350742 : }
     346             : static void
     347       10640 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     348             : 
     349             : int
     350      361396 : divisors_init(GEN n, GEN *pP, GEN *pE)
     351             : {
     352             :   long i,l;
     353             :   GEN E, P, e;
     354             :   int isint;
     355             : 
     356      361396 :   switch(typ(n))
     357             :   {
     358             :     case t_INT:
     359       10619 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     360       10619 :       set_fact(absi_factor(n), &P,&E);
     361       10619 :       isint = 1; break;
     362             :     case t_VEC:
     363          14 :       if (lg(n) != 3 || typ(gel(n,2)) !=t_MAT) pari_err_TYPE("divisors",n);
     364           7 :       set_fact_check(gel(n,2), &P,&E, &isint);
     365           7 :       break;
     366             :     case t_MAT:
     367      350742 :       set_fact_check(n, &P,&E, &isint);
     368      350735 :       break;
     369             :     default:
     370          21 :       set_fact(factor(n), &P,&E);
     371          21 :       isint = 0; break;
     372             :   }
     373      361382 :   l = lg(P);
     374      361382 :   e = cgetg(l, t_VECSMALL);
     375     1184323 :   for (i=1; i<l; i++)
     376             :   {
     377      822948 :     e[i] = itos(gel(E,i));
     378      822948 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     379             :   }
     380      361375 :   *pP = P; *pE = e; return isint;
     381             : }
     382             : 
     383             : GEN
     384      361354 : divisors(GEN n)
     385             : {
     386      361354 :   pari_sp av = avma;
     387             :   long i, j, l;
     388             :   ulong ndiv;
     389             :   GEN *d, *t, *t1, *t2, *t3, P, E, e;
     390      361354 :   int isint = divisors_init(n, &P, &E);
     391             : 
     392      361333 :   l = lg(E); e = cgetg(l, t_VECSMALL);
     393      361333 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     394      361333 :   ndiv = itou_or_0( zv_prod_Z(e) );
     395      361333 :   if (!ndiv || ndiv & ~LGBITS) pari_err_OVERFLOW("divisors");
     396      361333 :   d = t = (GEN*) cgetg(ndiv+1,t_VEC);
     397      361333 :   *++d = gen_1;
     398      361333 :   if (isint)
     399             :   {
     400     1184064 :     for (i=1; i<l; i++)
     401     2210012 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     402     1387267 :         for (t2=d, t3=t1; t3<t2; ) *++d = mulii(*++t3, gel(P,i));
     403      361319 :     e = ZV_sort((GEN)t);
     404             :   } else {
     405          56 :     for (i=1; i<l; i++)
     406         182 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     407         140 :         for (t2=d, t3=t1; t3<t2; ) *++d = gmul(*++t3, gel(P,i));
     408          14 :     e = (GEN)t;
     409             :   }
     410      361333 :   return gerepileupto(av, e);
     411             : }
     412             : 
     413             : GEN
     414       91460 : divisorsu_fact(GEN P, GEN E)
     415             : {
     416       91460 :   long i, j, l = lg(P);
     417       91460 :   ulong nbdiv = 1;
     418             :   GEN d, t, t1, t2, t3;
     419       91460 :   for (i=1; i<l; i++) nbdiv *= 1+E[i];
     420       91460 :   d = t = cgetg(nbdiv+1,t_VECSMALL);
     421       91460 :   *++d = 1;
     422      134866 :   for (i=1; i<l; i++)
     423       98856 :     for (t1=t,j=E[i]; j; j--,t1=t2)
     424       55450 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     425       91460 :   vecsmall_sort(t); return t;
     426             : }
     427             : GEN
     428       91460 : divisorsu(ulong n)
     429             : {
     430       91460 :   pari_sp av = avma;
     431       91460 :   GEN fa = factoru(n);
     432       91460 :   return gerepileupto(av, divisorsu_fact(gel(fa,1), gel(fa,2)));
     433             : }
     434             : 
     435             : static GEN
     436           0 : corefa(GEN fa)
     437             : {
     438           0 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1;
     439             :   long i;
     440           0 :   for (i=1; i<lg(P); i++)
     441           0 :     if (mod2(gel(E,i))) c = mulii(c,gel(P,i));
     442           0 :   return c;
     443             : }
     444             : static GEN
     445         728 : core2fa(GEN fa)
     446             : {
     447         728 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
     448             :   long i;
     449        1687 :   for (i=1; i<lg(P); i++)
     450             :   {
     451         959 :     long e = itos(gel(E,i));
     452         959 :     if (e & 1)  c = mulii(c, gel(P,i));
     453         959 :     if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
     454             :   }
     455         728 :   return mkvec2(c,f);
     456             : }
     457             : GEN
     458           0 : corepartial(GEN n, long all)
     459             : {
     460           0 :   pari_sp av = avma;
     461           0 :   if (typ(n) != t_INT) pari_err_TYPE("corepartial",n);
     462           0 :   return gerepileuptoint(av, corefa(Z_factor_limit(n,all)));
     463             : }
     464             : GEN
     465         728 : core2partial(GEN n, long all)
     466             : {
     467         728 :   pari_sp av = avma;
     468         728 :   if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
     469         728 :   return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
     470             : }
     471             : static GEN
     472        1484 : core2_i(GEN n)
     473             : {
     474        1484 :   GEN f = core(n);
     475        1484 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     476        1449 :   switch(typ(n))
     477             :   {
     478           7 :     case t_VEC: n = gel(n,1); break;
     479         721 :     case t_MAT: n = factorback(n); break;
     480             :   }
     481        1449 :   return mkvec2(f, sqrtint(diviiexact(n, f)));
     482             : }
     483             : GEN
     484        1477 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
     485             : 
     486             : GEN
     487        3094 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
     488             : 
     489             : static long
     490          14 : _mod4(GEN c) {
     491          14 :   long r, s = signe(c);
     492          14 :   if (!s) return 0;
     493          14 :   r = mod4(c); if (s < 0) r = 4-r;
     494          14 :   return r;
     495             : }
     496             : 
     497             : long
     498        2794 : corediscs(long D, ulong *f)
     499             : {
     500             :   /* D = f^2 dK */
     501        2794 :   long dK = D>=0 ? (long) coreu(D) : -(long) coreu(-(ulong) D);
     502        2794 :   ulong dKmod4 = ((ulong)dK)&3UL;
     503        2794 :   if (dKmod4 == 2 || dKmod4 == 3)
     504         231 :     dK *= 4;
     505        2794 :   if (f) *f = usqrt((ulong)(D/dK));
     506        2794 :   return D;
     507             : }
     508             : 
     509             : GEN
     510           7 : coredisc(GEN n)
     511             : {
     512           7 :   pari_sp av = avma;
     513           7 :   GEN c = core(n);
     514           7 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     515           7 :   return gerepileuptoint(av, shifti(c,2));
     516             : }
     517             : 
     518             : GEN
     519           7 : coredisc2(GEN n)
     520             : {
     521           7 :   pari_sp av = avma;
     522           7 :   GEN y = core2_i(n);
     523           7 :   GEN c = gel(y,1), f = gel(y,2);
     524           7 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
     525           7 :   y = cgetg(3,t_VEC);
     526           7 :   gel(y,1) = shifti(c,2);
     527           7 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
     528             : }
     529             : 
     530             : GEN
     531          14 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
     532             : 
     533             : long
     534         815 : omegau(ulong n)
     535             : {
     536             :   pari_sp av;
     537             :   GEN F;
     538         815 :   if (n == 1UL) return 0;
     539         801 :   av = avma; F = factoru(n);
     540         801 :   avma = av; return lg(gel(F,1))-1;
     541             : }
     542             : long
     543        1589 : omega(GEN n)
     544             : {
     545             :   pari_sp av;
     546             :   GEN F, P;
     547        1589 :   if ((F = check_arith_non0(n,"omega"))) {
     548             :     long n;
     549         721 :     P = gel(F,1); n = lg(P)-1;
     550         721 :     if (n && equalim1(gel(P,1))) n--;
     551         721 :     return n;
     552             :   }
     553         854 :   if (lgefint(n) == 3) return omegau(n[2]);
     554          39 :   av = avma;
     555          39 :   F = absi_factor(n);
     556          39 :   P = gel(F,1); avma = av; return lg(P)-1;
     557             : }
     558             : 
     559             : long
     560         821 : bigomegau(ulong n)
     561             : {
     562             :   pari_sp av;
     563             :   GEN F;
     564         821 :   if (n == 1) return 0;
     565         807 :   av = avma; F = factoru(n);
     566         807 :   avma = av; return zv_sum(gel(F,2));
     567             : }
     568             : long
     569        1589 : bigomega(GEN n)
     570             : {
     571        1589 :   pari_sp av = avma;
     572             :   GEN F, E;
     573        1589 :   if ((F = check_arith_non0(n,"bigomega")))
     574             :   {
     575         721 :     GEN P = gel(F,1);
     576         721 :     long n = lg(P)-1;
     577         721 :     E = gel(F,2);
     578         721 :     if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
     579             :   }
     580         854 :   else if (lgefint(n) == 3)
     581         821 :     return bigomegau(n[2]);
     582             :   else
     583          33 :     E = gel(absi_factor(n), 2);
     584         754 :   E = ZV_to_zv(E);
     585         754 :   avma = av; return zv_sum(E);
     586             : }
     587             : 
     588             : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
     589             : ulong
     590      151350 : eulerphiu_fact(GEN f)
     591             : {
     592      151350 :   GEN P = gel(f,1), E = gel(f,2);
     593      151350 :   long i, m = 1, l = lg(P);
     594      595966 :   for (i = 1; i < l; i++)
     595             :   {
     596      444616 :     ulong p = P[i], e = E[i];
     597      444616 :     if (!e) continue;
     598      444616 :     if (p == 2)
     599      113991 :     { if (e > 1) m <<= e-1; }
     600             :     else
     601             :     {
     602      330625 :       m *= (p-1);
     603      330625 :       if (e > 1) m *= upowuu(p, e-1);
     604             :     }
     605             :   }
     606      151350 :   return m;
     607             : }
     608             : ulong
     609      150090 : eulerphiu(ulong n)
     610             : {
     611      150090 :   pari_sp av = avma;
     612             :   GEN F;
     613      150090 :   if (!n) return 2;
     614      150090 :   F = factoru(n);
     615      150090 :   avma = av; return eulerphiu_fact(F);
     616             : }
     617             : GEN
     618      145607 : eulerphi(GEN n)
     619             : {
     620      145607 :   pari_sp av = avma;
     621             :   GEN Q, F, P, E;
     622             :   long i, l;
     623             : 
     624      145607 :   if ((F = check_arith_all(n,"eulerphi")))
     625             :   {
     626         735 :     F = clean_Z_factor(F);
     627         735 :     n = (typ(n) == t_VEC)? gel(n,1): factorback(F);
     628         735 :     if (lgefint(n) == 3)
     629             :     {
     630             :       ulong e;
     631         721 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
     632         721 :       e = eulerphiu_fact(F);
     633         721 :       avma = av; return utoipos(e);
     634             :     }
     635             :   }
     636      144872 :   else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
     637             :   else
     638          53 :     F = absi_factor(n);
     639          67 :   if (!signe(n)) return gen_2;
     640          39 :   P = gel(F,1);
     641          39 :   E = gel(F,2); l = lg(P);
     642          39 :   Q = cgetg(l, t_VEC);
     643         252 :   for (i = 1; i < l; i++)
     644             :   {
     645         213 :     GEN p = gel(P,i), q;
     646         213 :     ulong v = itou(gel(E,i));
     647         213 :     q = subiu(p,1);
     648         213 :     if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
     649         213 :     gel(Q,i) = q;
     650             :   }
     651          39 :   return gerepileuptoint(av, ZV_prod(Q));
     652             : }
     653             : 
     654             : static GEN
     655         760 : numdiv_aux(GEN F)
     656             : {
     657         760 :   GEN x, E = gel(F,2);
     658         760 :   long i, l = lg(E);
     659         760 :   x = cgetg(l, t_VECSMALL);
     660         760 :   for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
     661         760 :   return x;
     662             : }
     663             : GEN
     664        1897 : numdiv(GEN n)
     665             : {
     666        1897 :   pari_sp av = avma;
     667             :   GEN F, E;
     668             :   long i, l;
     669        1897 :   if ((F = check_arith_non0(n,"numdiv")))
     670             :   {
     671         721 :     F = clean_Z_factor(F);
     672         721 :     E = numdiv_aux(F);
     673             :   }
     674        1162 :   else if (lgefint(n) == 3)
     675             :   {
     676        1123 :     if (n[2] == 1) return gen_1;
     677        1095 :     F = factoru(n[2]);
     678        1095 :     E = gel(F,2); l = lg(E);
     679        1095 :     for (i=1; i<l; i++) E[i]++;
     680             :   }
     681             :   else
     682          39 :     E = numdiv_aux(absi_factor(n));
     683        1855 :   return gerepileuptoint(av, zv_prod_Z(E));
     684             : }
     685             : 
     686             : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
     687             : static GEN
     688        1635 : u_euler_sumdiv(ulong p, long v)
     689             : {
     690        1635 :   GEN u = utoipos(1 + p); /* can't overflow */
     691        1635 :   for (; v > 1; v--) u = addsi(1, mului(p, u));
     692        1635 :   return u;
     693             : }
     694             : /* 1 + q + ... + q^v */
     695             : static GEN
     696        9677 : euler_sumdiv(GEN q, long v)
     697             : {
     698        9677 :   GEN u = addui(1, q);
     699        9677 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
     700        9677 :   return u;
     701             : }
     702             : 
     703             : static GEN
     704         753 : sumdiv_aux(GEN F)
     705             : {
     706         753 :   GEN x, P = gel(F,1), E = gel(F,2);
     707         753 :   long i, l = lg(P);
     708         753 :   x = cgetg(l, t_VEC);
     709         753 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
     710         753 :   return ZV_prod(x);
     711             : }
     712             : GEN
     713        1589 : sumdiv(GEN n)
     714             : {
     715        1589 :   pari_sp av = avma;
     716             :   GEN F, v;
     717             : 
     718        1589 :   if ((F = check_arith_non0(n,"sumdiv")))
     719             :   {
     720         721 :     F = clean_Z_factor(F);
     721         721 :     v = sumdiv_aux(F);
     722             :   }
     723         854 :   else if (lgefint(n) == 3)
     724             :   {
     725         822 :     if (n[2] == 1) return gen_1;
     726         808 :     F = factoru(n[2]);
     727         808 :     v = usumdiv_fact(F);
     728             :   }
     729             :   else
     730          32 :     v = sumdiv_aux(absi_factor(n));
     731        1561 :   return gerepileuptoint(av, v);
     732             : }
     733             : 
     734             : static GEN
     735        1506 : sumdivk_aux(GEN F, long k)
     736             : {
     737        1506 :   GEN x, P = gel(F,1), E = gel(F,2);
     738        1506 :   long i, l = lg(P);
     739        1506 :   x = cgetg(l, t_VEC);
     740        1506 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(powiu(gel(P,i),k), gel(E,i)[2]);
     741        1506 :   return ZV_prod(x);
     742             : }
     743             : GEN
     744        6545 : sumdivk(GEN n, long k)
     745             : {
     746        6545 :   pari_sp av = avma;
     747             :   GEN F, v;
     748             :   long k1;
     749             : 
     750        6545 :   if (!k) return numdiv(n);
     751        6545 :   if (k == 1) return sumdiv(n);
     752        4956 :   if (k ==-1) return gerepileupto(av, gdiv(sumdiv(n), n));
     753        4956 :   k1 = k;
     754        4956 :   if (k < 0)  k = -k;
     755        4956 :   if ((F = check_arith_non0(n,"sumdivk")))
     756             :   {
     757        1442 :     F = clean_Z_factor(F);
     758        1442 :     v = sumdivk_aux(F,k);
     759             :   }
     760        3486 :   else if (lgefint(n) == 3)
     761             :   {
     762        3422 :     if (n[2] == 1) return gen_1;
     763        3324 :     F = factoru(n[2]);
     764        3324 :     v = usumdivk_fact(F,k);
     765             :   }
     766             :   else
     767          64 :     v = sumdivk_aux(absi_factor(n), k);
     768        4830 :   if (k1 > 0) return gerepileuptoint(av, v);
     769           7 :   return gerepileupto(av, gdiv(v, powiu(n,k)));
     770             : }
     771             : 
     772             : GEN
     773         808 : usumdiv_fact(GEN f)
     774             : {
     775         808 :   GEN P = gel(f,1), E = gel(f,2);
     776         808 :   long i, l = lg(P);
     777         808 :   GEN v = cgetg(l, t_VEC);
     778         808 :   for (i=1; i<l; i++) gel(v,i) = u_euler_sumdiv(P[i],E[i]);
     779         808 :   return ZV_prod(v);
     780             : }
     781             : GEN
     782        3380 : usumdivk_fact(GEN f, ulong k)
     783             : {
     784        3380 :   GEN P = gel(f,1), E = gel(f,2);
     785        3380 :   long i, l = lg(P);
     786        3380 :   GEN v = cgetg(l, t_VEC);
     787        3380 :   for (i=1; i<l; i++) gel(v,i) = euler_sumdiv(powuu(P[i],k),E[i]);
     788        3380 :   return ZV_prod(v);
     789             : }
     790             : 
     791             : long
     792         280 : uissquarefree_fact(GEN f)
     793             : {
     794         280 :   GEN E = gel(f,2);
     795         280 :   long i, l = lg(E);
     796         574 :   for (i = 1; i < l; i++)
     797         434 :     if (E[i] > 1) return 0;
     798         140 :   return 1;
     799             : }
     800             : long
     801       11485 : uissquarefree(ulong n)
     802             : {
     803       11485 :   if (!n) return 0;
     804       11485 :   return moebiusu(n)? 1: 0;
     805             : }
     806             : long
     807         772 : Z_issquarefree(GEN n)
     808             : {
     809         772 :   switch(lgefint(n))
     810             :   {
     811           7 :     case 2: return 0;
     812           6 :     case 3: return uissquarefree(n[2]);
     813             :   }
     814         759 :   return moebius(n)? 1: 0;
     815             : }
     816             : long
     817       46970 : issquarefree(GEN x)
     818             : {
     819             :   pari_sp av;
     820             :   GEN d;
     821       46970 :   switch(typ(x))
     822             :   {
     823          14 :     case t_INT: return Z_issquarefree(x);
     824             :     case t_POL:
     825       46956 :       if (!signe(x)) return 0;
     826       46956 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
     827       46956 :       avma = av; return (lg(d) == 3);
     828           0 :     default: pari_err_TYPE("issquarefree",x);
     829           0 :       return 0; /* not reached */
     830             :   }
     831             : }
     832             : 
     833             : /*********************************************************************/
     834             : /**                                                                 **/
     835             : /**                    DIGITS / SUM OF DIGITS                       **/
     836             : /**                                                                 **/
     837             : /*********************************************************************/
     838             : 
     839             : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
     840             : static void
     841     1820756 : set_vexp(GEN v, long l)
     842             : {
     843             :   long m;
     844     3641512 :   if (v[l]) return;
     845      610574 :   v[l] = 1; m = l>>1;
     846      610574 :   set_vexp(v, m);
     847      610574 :   set_vexp(v, l-m);
     848             : }
     849             : 
     850             : /* return all needed B^i for DAC algorithm, for lz digits */
     851             : static GEN
     852      599608 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
     853             : {
     854      599608 :   GEN vB, vexp = const_vecsmall(lz, 0);
     855      599608 :   long i, l = (lz+1) >> 1;
     856      599608 :   vexp[1] = 1;
     857      599608 :   vexp[2] = 1;
     858      599608 :   set_vexp(vexp, lz);
     859      599608 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
     860      599608 :   gel(vB, 1) = T;
     861     1261161 :   for (i = 2; i <= l; i++)
     862      661553 :     if (vexp[i])
     863             :     {
     864      610574 :       long j = i >> 1;
     865      610574 :       GEN B2j = r->sqr(E, gel(vB,j));
     866      610574 :       gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
     867             :     }
     868      599608 :   return vB;
     869             : }
     870             : 
     871             : static void
     872      166283 : gen_digits_dac(GEN x, GEN vB, long l, GEN *z,
     873             :                void *E, GEN div(void *E, GEN a, GEN b, GEN *r))
     874             : {
     875             :   GEN q, r;
     876      166283 :   long m = l>>1;
     877      332566 :   if (l==1) { *z=x; return; }
     878       80047 :   q = div(E, x, gel(vB,m), &r);
     879       80047 :   gen_digits_dac(r, vB, m, z, E, div);
     880       80047 :   gen_digits_dac(q, vB, l-m, z+m, E, div);
     881             : }
     882             : 
     883             : static GEN
     884       28889 : gen_fromdigits_dac(GEN x, GEN vB, long i, long l, void *E,
     885             :                    GEN add(void *E, GEN a, GEN b),
     886             :                    GEN mul(void *E, GEN a, GEN b))
     887             : {
     888             :   GEN a, b;
     889       28889 :   long m = l>>1;
     890       28889 :   if (l==1) return gel(x,i);
     891       13181 :   a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
     892       13181 :   b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
     893       13181 :   return add(E, a, mul(E, b, gel(vB, m)));
     894             : }
     895             : 
     896             : static GEN
     897        6280 : gen_digits_i(GEN x, GEN B, long n, void *E, struct bb_ring *r,
     898             :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
     899             : {
     900             :   GEN z, vB;
     901        6280 :   if (n==1) retmkvec(gcopy(x));
     902        6189 :   vB = get_vB(B, n, E, r);
     903        6189 :   z = cgetg(n+1, t_VEC);
     904        6189 :   gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
     905        6189 :   return z;
     906             : }
     907             : 
     908             : GEN
     909        6230 : gen_digits(GEN x, GEN B, long n, void *E, struct bb_ring *r,
     910             :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
     911             : {
     912        6230 :   pari_sp av = avma;
     913        6230 :   return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
     914             : }
     915             : 
     916             : GEN
     917        2527 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
     918             : {
     919        2527 :   pari_sp av = avma;
     920        2527 :   long n = lg(x)-1;
     921        2527 :   GEN vB = get_vB(B, n, E, r);
     922        2527 :   GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
     923        2527 :   return gerepilecopy(av, z);
     924             : }
     925             : 
     926             : static GEN
     927        1771 : _addii(void *data /* ignored */, GEN x, GEN y)
     928        1771 : { (void)data; return addii(x,y); }
     929             : static GEN
     930      584527 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
     931             : static GEN
     932      150874 : _mulii(void *data /* ignored */, GEN x, GEN y)
     933      150874 :  { (void)data; return mulii(x,y); }
     934             : static GEN
     935         436 : _dvmdii(void *data /* ignored */, GEN x, GEN y, GEN *r)
     936         436 :  { (void)data; return dvmdii(x,y,r); }
     937             : 
     938             : static struct bb_ring Z_ring = { _addii, _mulii, _sqri };
     939             : 
     940             : static GEN
     941         182 : check_basis(GEN B)
     942             : {
     943         182 :   if (!B) return utoipos(10);
     944         161 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
     945         161 :   if (cmpiu(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
     946         161 :   return B;
     947             : }
     948             : 
     949             : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
     950             : static void
     951        3075 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
     952             : {
     953        3075 :   pari_sp av = avma;
     954             :   GEN q,r;
     955             :   long m;
     956        6150 :   if (l==1) { *z=itou(x); return; }
     957        1524 :   m=l>>1;
     958        1524 :   q = dvmdii(x, gel(vB,m), &r);
     959        1524 :   digits_dacsmall(q,vB,l-m,z);
     960        1524 :   digits_dacsmall(r,vB,m,z+l-m);
     961        1524 :   avma = av;
     962             : }
     963             : 
     964             : GEN
     965          56 : digits(GEN x, GEN B)
     966             : {
     967          56 :   pari_sp av=avma;
     968             :   long lz;
     969             :   GEN z, vB;
     970          56 :   if (typ(x)!=t_INT) pari_err_TYPE("digits",x);
     971          56 :   B = check_basis(B);
     972          56 :   if (signe(B)<0) pari_err_DOMAIN("digits","B","<",gen_0,B);
     973          56 :   if (!signe(x))       {avma = av; return cgetg(1,t_VEC); }
     974          49 :   if (absi_cmp(x,B)<0) {avma = av; retmkvec(absi(x)); }
     975          49 :   if (Z_ispow2(B))
     976             :   {
     977          14 :     long k = expi(B);
     978          14 :     if (k == 1) return binaire(x);
     979           7 :     if (k < BITS_IN_LONG)
     980             :     {
     981           0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
     982           0 :       z = binary_2k_nv(x, k);
     983           0 :       avma = av; return Flv_to_ZV(z);
     984             :     }
     985             :     else
     986             :     {
     987           7 :       avma = av; return binary_2k(x, k);
     988             :     }
     989             :   }
     990          35 :   if (signe(x) < 0) x = absi(x);
     991          35 :   lz = logint(x,B,NULL);
     992          35 :   if (lgefint(B)>3)
     993             :   {
     994           8 :     z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
     995           8 :     vecreverse_inplace(z);
     996           8 :     return z;
     997             :   }
     998             :   else
     999             :   {
    1000          27 :     vB = get_vB(B, lz, NULL, &Z_ring);
    1001          27 :     (void)new_chunk(3*lz); /* HACK */
    1002          27 :     z = zero_zv(lz);
    1003          27 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1004          27 :     avma = av; return Flv_to_ZV(z);
    1005             :   }
    1006             : }
    1007             : 
    1008             : static GEN
    1009     1911499 : fromdigitsu_dac(GEN x, GEN vB, long i, long l)
    1010             : {
    1011             :   GEN a, b;
    1012     1911499 :   long m = l>>1;
    1013     1911499 :   if (l==1) return utoi(uel(x,i));
    1014     1577028 :   if (l==2) return addumului(uel(x,i), uel(x,i+1), gel(vB, m));
    1015      660317 :   a = fromdigitsu_dac(x, vB, i, m);
    1016      660317 :   b = fromdigitsu_dac(x, vB, i+m, l-m);
    1017      660317 :   return addii(a, mulii(b, gel(vB, m)));
    1018             : }
    1019             : 
    1020             : GEN
    1021      590865 : fromdigitsu(GEN x, GEN B)
    1022             : {
    1023      590865 :   pari_sp av = avma;
    1024      590865 :   long n = lg(x)-1;
    1025             :   GEN vB, z;
    1026      590865 :   if (n==0) return gen_0;
    1027      590865 :   vB = get_vB(B, n, NULL, &Z_ring);
    1028      590865 :   z = fromdigitsu_dac(x, vB, 1, n);
    1029      590865 :   return gerepileuptoint(av, z);
    1030             : }
    1031             : 
    1032             : static int
    1033          14 : ZV_in_range(GEN v, GEN B)
    1034             : {
    1035          14 :   long i, l = lg(v);
    1036        1659 :   for(i=1; i < l; i++)
    1037             :   {
    1038        1652 :     GEN vi = gel(v, i);
    1039        1652 :     if (signe(vi) < 0 || cmpii(vi, B) >= 0)
    1040           7 :       return 0;
    1041             :   }
    1042           7 :   return 1;
    1043             : }
    1044             : 
    1045             : GEN
    1046          56 : fromdigits(GEN x, GEN B)
    1047             : {
    1048          56 :   pari_sp av = avma;
    1049          56 :   if (typ(x)!=t_VEC || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
    1050          56 :   if (lg(x)==1) return gen_0;
    1051          49 :   B = check_basis(B);
    1052          49 :   if (Z_ispow2(B) && ZV_in_range(x, B))
    1053           7 :     return fromdigits_2k(x, expi(B));
    1054          42 :   x = vecreverse(x);
    1055          42 :   return gerepileuptoint(av, gen_fromdigits(x, B, NULL, &Z_ring));
    1056             : }
    1057             : 
    1058             : static const ulong digsum[] ={
    1059             :   0,1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9,10,2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,
    1060             :   9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,
    1061             :   12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,
    1062             :   12,13,14,15,16,17,18,1,2,3,4,5,6,7,8,9,10,2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,
    1063             :   9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,
    1064             :   12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,
    1065             :   12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,2,3,4,5,6,7,8,9,10,11,3,
    1066             :   4,5,6,7,8,9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,
    1067             :   9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,
    1068             :   9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,
    1069             :   16,17,18,19,20,3,4,5,6,7,8,9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,
    1070             :   11,12,13,14,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,
    1071             :   12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,
    1072             :   19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,21,4,5,6,7,8,9,
    1073             :   10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,
    1074             :   12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,
    1075             :   11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,
    1076             :   18,19,20,21,13,14,15,16,17,18,19,20,21,22,5,6,7,8,9,10,11,12,13,14,6,7,8,9,
    1077             :   10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,
    1078             :   10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,
    1079             :   17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,20,21,22,14,
    1080             :   15,16,17,18,19,20,21,22,23,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,
    1081             :   15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,12,13,
    1082             :   14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,
    1083             :   21,13,14,15,16,17,18,19,20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,17,18,
    1084             :   19,20,21,22,23,24,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,
    1085             :   10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,
    1086             :   17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,20,21,22,14,
    1087             :   15,16,17,18,19,20,21,22,23,15,16,17,18,19,20,21,22,23,24,16,17,18,19,20,21,
    1088             :   22,23,24,25,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,
    1089             :   12,13,14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,
    1090             :   19,20,21,13,14,15,16,17,18,19,20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,
    1091             :   17,18,19,20,21,22,23,24,16,17,18,19,20,21,22,23,24,25,17,18,19,20,21,22,23,
    1092             :   24,25,26,9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,
    1093             :   13,14,15,16,17,18,19,20,12,13,14,15,16,17,18,19,20,21,13,14,15,16,17,18,19,
    1094             :   20,21,22,14,15,16,17,18,19,20,21,22,23,15,16,17,18,19,20,21,22,23,24,16,17,
    1095             :   18,19,20,21,22,23,24,25,17,18,19,20,21,22,23,24,25,26,18,19,20,21,22,23,24,
    1096             :   25,26,27
    1097             : };
    1098             : 
    1099             : ulong
    1100      355152 : sumdigitsu(ulong n)
    1101             : {
    1102      355152 :   ulong s = 0;
    1103      355152 :   while (n) { s += digsum[n % 1000]; n /= 1000; }
    1104      355152 :   return s;
    1105             : }
    1106             : 
    1107             : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1108             : static ulong
    1109          14 : sumdigits_block(ulong *res, long l)
    1110             : {
    1111          14 :   long s = sumdigitsu(*--res);
    1112          14 :   while (--l > 0) s += sumdigitsu(*--res);
    1113          14 :   return s;
    1114             : }
    1115             : 
    1116             : GEN
    1117          35 : sumdigits(GEN n)
    1118             : {
    1119          35 :   pari_sp av = avma;
    1120             :   ulong s, *res;
    1121             :   long l;
    1122             : 
    1123          35 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1124          35 :   l = lgefint(n);
    1125          35 :   switch(l)
    1126             :   {
    1127           7 :     case 2: return gen_0;
    1128          14 :     case 3: return utoipos(sumdigitsu(n[2]));
    1129             :   }
    1130          14 :   res = convi(n, &l);
    1131          14 :   if ((ulong)l < ULONG_MAX / 81)
    1132             :   {
    1133          14 :     s = sumdigits_block(res, l);
    1134          14 :     avma = av; return utoipos(s);
    1135             :   }
    1136             :   else /* Huge. Overflows ulong */
    1137             :   {
    1138           0 :     const long L = (long)(ULONG_MAX / 81);
    1139           0 :     GEN S = gen_0;
    1140           0 :     while (l > L)
    1141             :     {
    1142           0 :       S = addiu(S, sumdigits_block(res, L));
    1143           0 :       res += L; l -= L;
    1144             :     }
    1145           0 :     if (l)
    1146           0 :       S = addiu(S, sumdigits_block(res, l));
    1147           0 :     return gerepileuptoint(av, S);
    1148             :   }
    1149             : }
    1150             : 
    1151             : GEN
    1152         105 : sumdigits0(GEN x, GEN B)
    1153             : {
    1154         105 :   pari_sp av = avma;
    1155             :   GEN z;
    1156             :   long lz;
    1157             : 
    1158         105 :   if (!B) return sumdigits(x);
    1159          77 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1160          77 :   B = check_basis(B);
    1161          77 :   if (Z_ispow2(B))
    1162             :   {
    1163          28 :     long k = expi(B);
    1164          28 :     if (k == 1) { avma = av; return utoi(hammingweight(x)); }
    1165          21 :     if (k < BITS_IN_LONG)
    1166             :     {
    1167          14 :       GEN z = binary_2k_nv(x, k);
    1168          14 :       if (lg(z)-1 > 1L<<(BITS_IN_LONG-k)) /* may overflow */
    1169           0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1170          14 :       avma = av; return utoi(zv_sum(z));
    1171             :     }
    1172           7 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1173             :   }
    1174          49 :   if (!signe(x))       { avma = av; return gen_0; }
    1175          49 :   if (absi_cmp(x,B)<0) { avma = av; return absi(x); }
    1176          49 :   if (equaliu(B,10))   { avma = av; return sumdigits(x); }
    1177          42 :   lz = logint(x,B,NULL);
    1178          42 :   z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
    1179          42 :   return gerepileuptoint(av, ZV_sum(z));
    1180             : }

Generated by: LCOV version 1.11