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.12.0 lcov report (development 23332-367b47754) Lines: 590 671 87.9 %
Date: 2018-12-10 05:41:52 Functions: 87 95 91.6 %
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          42 : boundfact(GEN n, ulong lim)
      25             : {
      26          42 :   switch(typ(n))
      27             :   {
      28          28 :     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(a,b,(void*)&cmpii,cmp_nodata));
      35             :     }
      36             :   }
      37           0 :   pari_err_TYPE("boundfact",n);
      38             :   return NULL; /* LCOV_EXCL_LINE */
      39             : }
      40             : 
      41             : /* NOT memory clean */
      42             : GEN
      43       12183 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
      44             : {
      45       12183 :   long i, j, l = lg(L);
      46       12183 :   GEN e = new_chunk(l), P = new_chunk(l);
      47      112719 :   for (i = j = 1; i < l; i++)
      48             :   {
      49      107498 :     ulong p = uel(L,i);
      50      107498 :     long v = Z_lvalrem(N, p, &N);
      51      107498 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
      52             :   }
      53       12183 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
      54       12183 :   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       31543 : factoru_pow(ulong n)
      65             : {
      66       31543 :   GEN f = cgetg(4,t_VEC);
      67       31538 :   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       31538 :   (void)new_chunk((15 + 1)*3);
      72       31538 :   F = factoru(n);
      73       31545 :   P = gel(F,1);
      74       31545 :   E = gel(F,2); l = lg(P);
      75       31545 :   set_avma(av);
      76       31546 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
      77       31547 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
      78       31542 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
      79       96573 :   for(i = 1; i < l; i++)
      80             :   {
      81       65034 :     p[i] = P[i];
      82       65034 :     e[i] = E[i];
      83       65034 :     c[i] = upowuu(p[i], e[i]);
      84             :   }
      85       31539 :   return f;
      86             : }
      87             : 
      88             : static GEN
      89       95070 : factorlim(GEN n, ulong lim)
      90       95070 : { 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       72768 : factor_pn_1_limit(GEN p, long n, ulong lim)
      95             : {
      96       72768 :   pari_sp av = avma;
      97       72768 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
      98       72768 :   long i, pp = itos_or_0(p);
      99       89876 :   for(i=2; i<lg(d); i++)
     100             :   {
     101             :     GEN B;
     102       30226 :     if (pp && d[i]%pp==0 && (
     103       26236 :        ((pp&3)==1 && (d[i]&1)) ||
     104       13237 :        ((pp&3)==3 && (d[i]&3)==2) ||
     105       12985 :        (pp==2 && (d[i]&7)==4)))
     106        5194 :     {
     107        5194 :       GEN f=factor_Aurifeuille_prime(p,d[i]);
     108        5194 :       B = factorlim(f, lim);
     109        5194 :       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     110        5194 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     111             :     }
     112             :     else
     113       11914 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     114       17108 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     115             :   }
     116       72768 :   return gerepilecopy(av, A);
     117             : }
     118             : GEN
     119       72768 : factor_pn_1(GEN p, ulong n)
     120       72768 : { 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 = subiu(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             :   set_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             : int
     185     6144413 : RgV_is_ZVpos(GEN v)
     186             : {
     187     6144413 :   long i, l = lg(v);
     188    18606307 :   for (i = 1; i < l; i++)
     189             :   {
     190    12467305 :     GEN c = gel(v,i);
     191    12467305 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     192             :   }
     193     6139002 :   return 1;
     194             : }
     195             : /* check whether v is a ZV with non-0 entries */
     196             : int
     197        6811 : RgV_is_ZVnon0(GEN v)
     198             : {
     199        6811 :   long i, l = lg(v);
     200       20384 :   for (i = 1; i < l; i++)
     201             :   {
     202       13629 :     GEN c = gel(v,i);
     203       13629 :     if (typ(c) != t_INT || !signe(c)) return 0;
     204             :   }
     205        6755 :   return 1;
     206             : }
     207             : /* check whether v is a ZV with non-zero entries OR exactly [0] */
     208             : static int
     209       17052 : RgV_is_ZV0(GEN v)
     210             : {
     211       17052 :   long i, l = lg(v);
     212       43337 :   for (i = 1; i < l; i++)
     213             :   {
     214       26397 :     GEN c = gel(v,i);
     215             :     long s;
     216       26397 :     if (typ(c) != t_INT) return 0;
     217       26397 :     s = signe(c);
     218       26397 :     if (!s) return (l == 2);
     219             :   }
     220       16940 :   return 1;
     221             : }
     222             : 
     223             : static int
     224       40404 : RgV_is_prV(GEN v)
     225             : {
     226       40404 :   long l = lg(v), i;
     227       41643 :   for (i = 1; i < l; i++)
     228       40733 :     if (!checkprid_i(gel(v,i))) return 0;
     229         910 :   return 1;
     230             : }
     231             : int
     232       44135 : is_nf_factor(GEN F)
     233             : {
     234       87276 :   return typ(F) == t_MAT && lg(F) == 3
     235       84532 :     && RgV_is_prV(gel(F,1)) && RgV_is_ZVpos(gel(F,2));
     236             : }
     237             : int
     238          14 : is_nf_extfactor(GEN F)
     239             : {
     240          28 :   return typ(F) == t_MAT && lg(F) == 3
     241          21 :     && RgV_is_prV(gel(F,1)) && RgV_is_ZV(gel(F,2));
     242             : }
     243             : 
     244             : static int
     245     3080904 : is_Z_factor_i(GEN f)
     246     3080904 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
     247             : int
     248     3057034 : is_Z_factorpos(GEN f)
     249     3057034 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     250             : int
     251       17052 : is_Z_factor(GEN f)
     252       17052 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     253             : /* as is_Z_factorpos, also allow factor(0) */
     254             : int
     255        6818 : is_Z_factornon0(GEN f)
     256        6818 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     257             : GEN
     258       15589 : clean_Z_factor(GEN f)
     259             : {
     260       15589 :   GEN P = gel(f,1);
     261       15589 :   long n = lg(P)-1;
     262       15589 :   if (n && equalim1(gel(P,1)))
     263        2513 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     264       13076 :   return f;
     265             : }
     266             : GEN
     267           0 : fuse_Z_factor(GEN f, GEN B)
     268             : {
     269           0 :   GEN P = gel(f,1), E = gel(f,2), P2,E2;
     270           0 :   long i, l = lg(P);
     271           0 :   if (l == 1) return f;
     272           0 :   for (i = 1; i < l; i++)
     273           0 :     if (abscmpii(gel(P,i), B) > 0) break;
     274           0 :   if (i == l) return f;
     275             :   /* tail / initial segment */
     276           0 :   P2 = vecslice(P, i, l-1); P = vecslice(P, 1, i-1);
     277           0 :   E2 = vecslice(E, i, l-1); E = vecslice(E, 1, i-1);
     278           0 :   P = shallowconcat(P, mkvec(factorback2(P2,E2)));
     279           0 :   E = shallowconcat(E, mkvec(gen_1));
     280           0 :   return mkmat2(P, E);
     281             : }
     282             : 
     283             : /* n attached to a factorization of a positive integer: either N (t_INT)
     284             :  * a factorization matrix faN, or a t_VEC: [N, faN] */
     285             : GEN
     286         210 : check_arith_pos(GEN n, const char *f) {
     287         210 :   switch(typ(n))
     288             :   {
     289             :     case t_INT:
     290         203 :       if (signe(n) <= 0) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
     291         203 :       return NULL;
     292             :     case t_VEC:
     293           0 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT || signe(gel(n,1)) <= 0) break;
     294           0 :       n = gel(n,2); /* fall through */
     295             :     case t_MAT:
     296           7 :       if (!is_Z_factorpos(n)) break;
     297           7 :       return n;
     298             :   }
     299           0 :   pari_err_TYPE(f,n);
     300           0 :   return NULL;
     301             : }
     302             : /* n attached to a factorization of a non-0 integer */
     303             : GEN
     304     2782133 : check_arith_non0(GEN n, const char *f) {
     305     2782133 :   switch(typ(n))
     306             :   {
     307             :     case t_INT:
     308     2775555 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     309     2775495 :       return NULL;
     310             :     case t_VEC:
     311          14 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT || !signe(gel(n,1))) break;
     312          14 :       n = gel(n,2); /* fall through */
     313             :     case t_MAT:
     314        6580 :       if (!is_Z_factornon0(n)) break;
     315        6524 :       return n;
     316             :   }
     317          54 :   pari_err_TYPE(f,n);
     318           0 :   return NULL;
     319             : }
     320             : /* n attached to a factorization of an integer */
     321             : GEN
     322      262164 : check_arith_all(GEN n, const char *f) {
     323      262164 :   switch(typ(n))
     324             :   {
     325             :     case t_INT:
     326      245112 :       return NULL;
     327             :     case t_VEC:
     328       11179 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     329       11179 :       n = gel(n,2); /* fall through */
     330             :     case t_MAT:
     331       17052 :       if (!is_Z_factor(n)) break;
     332       17052 :       return n;
     333             :   }
     334           0 :   pari_err_TYPE(f,n);
     335           0 :   return NULL;
     336             : }
     337             : 
     338             : /***********************************************************************/
     339             : /**                                                                   **/
     340             : /**                MISCELLANEOUS ARITHMETIC FUNCTIONS                 **/
     341             : /**                (ultimately depend on Z_factor())                  **/
     342             : /**                                                                   **/
     343             : /***********************************************************************/
     344             : /* set P,E from F. Check whether F is an integer and kill "factor" -1 */
     345             : static void
     346      350749 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     347             : {
     348             :   GEN E, P;
     349      350749 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     350      350749 :   P = gel(F,1);
     351      350749 :   E = gel(F,2);
     352      350749 :   RgV_check_ZV(E, "divisors");
     353      350749 :   *isint = RgV_is_ZV(P);
     354      350749 :   if (*isint)
     355             :   {
     356      350735 :     long i, l = lg(P);
     357             :     /* skip -1 */
     358      350735 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     359             :     /* test for 0 */
     360     1143135 :     for (i = 1; i < l; i++)
     361      792407 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     362           7 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     363             :   }
     364      350742 :   *pP = P;
     365      350742 :   *pE = E;
     366      350742 : }
     367             : static void
     368       11515 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     369             : 
     370             : int
     371      362271 : divisors_init(GEN n, GEN *pP, GEN *pE)
     372             : {
     373             :   long i,l;
     374             :   GEN E, P, e;
     375             :   int isint;
     376             : 
     377      362271 :   switch(typ(n))
     378             :   {
     379             :     case t_INT:
     380       11494 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     381       11494 :       set_fact(absZ_factor(n), &P,&E);
     382       11494 :       isint = 1; break;
     383             :     case t_VEC:
     384          14 :       if (lg(n) != 3 || typ(gel(n,2)) !=t_MAT) pari_err_TYPE("divisors",n);
     385           7 :       set_fact_check(gel(n,2), &P,&E, &isint);
     386           7 :       break;
     387             :     case t_MAT:
     388      350742 :       set_fact_check(n, &P,&E, &isint);
     389      350735 :       break;
     390             :     default:
     391          21 :       set_fact(factor(n), &P,&E);
     392          21 :       isint = 0; break;
     393             :   }
     394      362257 :   l = lg(P);
     395      362257 :   e = cgetg(l, t_VECSMALL);
     396     1186598 :   for (i=1; i<l; i++)
     397             :   {
     398      824348 :     e[i] = itos(gel(E,i));
     399      824348 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     400             :   }
     401      362250 :   *pP = P; *pE = e; return isint;
     402             : }
     403             : 
     404             : static long
     405      362208 : ndiv(GEN E)
     406             : {
     407             :   long i, l;
     408      362208 :   GEN e = cgetg_copy(E, &l); /* left on stack */
     409             :   ulong n;
     410      362208 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     411      362208 :   n = itou_or_0( zv_prod_Z(e) );
     412      362208 :   if (!n || n & ~LGBITS) pari_err_OVERFLOW("divisors");
     413      362208 :   return n;
     414             : }
     415             : static int
     416        4354 : cmpi1(void *E, GEN a, GEN b) { (void)E; return cmpii(gel(a,1), gel(b,1)); }
     417             : /* P a t_COL of objects, E a t_VECSMALL of exponents, return cleaned-up
     418             :  * factorization (removing 0 exponents) as a t_MAT with 2 cols. */
     419             : static GEN
     420        3486 : fa_clean(GEN P, GEN E)
     421             : {
     422        3486 :   long i, j, l = lg(E);
     423        3486 :   GEN Q = cgetg(l, t_COL);
     424       10311 :   for (i = j = 1; i < l; i++)
     425        6825 :     if (E[i]) { gel(Q,j) = gel(P,i); E[j] = E[i]; j++; }
     426        3486 :   setlg(Q,j);
     427        3486 :   setlg(E,j); return mkmat2(Q,Flc_to_ZC(E));
     428             : }
     429             : GEN
     430         721 : divisors_factored(GEN N)
     431             : {
     432         721 :   pari_sp av = avma;
     433             :   GEN *d, *t1, *t2, *t3, D, P, E;
     434         721 :   int isint = divisors_init(N, &P, &E);
     435         721 :   GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
     436         721 :   long i, j, l, n = ndiv(E);
     437             : 
     438         721 :   D = cgetg(n+1,t_VEC); d = (GEN*)D;
     439         721 :   l = lg(E);
     440         721 :   *++d = mkvec2(gen_1, const_vecsmall(l-1,0));
     441        1946 :   for (i=1; i<l; i++)
     442        2947 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     443        6209 :       for (t2=d, t3=t1; t3<t2; )
     444             :       {
     445             :         GEN a, b;
     446        2765 :         a = mul(gel(*++t3,1), gel(P,i));
     447        2765 :         b = leafcopy(gel(*t3,2)); b[i]++;
     448        2765 :         *++d = mkvec2(a,b);
     449             :       }
     450         721 :   if (isint) gen_sort_inplace(D,NULL,&cmpi1,NULL);
     451         721 :   for (i = 1; i <= n; i++) gmael(D,i,2) = fa_clean(P, gmael(D,i,2));
     452         721 :   return gerepilecopy(av, D);
     453             : }
     454             : static int
     455           0 : cmpu1(void *E, GEN va, GEN vb)
     456           0 : { long a = va[1], b = vb[1]; (void)E; return a>b? 1: (a<b? -1: 0); }
     457             : static GEN
     458           0 : fa_clean_u(GEN P, GEN E)
     459             : {
     460           0 :   long i, j, l = lg(E);
     461           0 :   GEN Q = cgetg(l, t_VECSMALL);
     462           0 :   for (i = j = 1; i < l; i++)
     463           0 :     if (E[i]) { Q[j] = P[i]; E[j] = E[i]; j++; }
     464           0 :   setlg(Q,j);
     465           0 :   setlg(E,j); return mkmat2(Q,E);
     466             : }
     467             : GEN
     468           0 : divisorsu_fact_factored(GEN fa)
     469             : {
     470           0 :   pari_sp av = avma;
     471           0 :   GEN *d, *t1, *t2, *t3, vD, D, P = gel(fa,1), E = gel(fa,2);
     472           0 :   long i, j, l, n = ndiv(E);
     473             : 
     474           0 :   D = cgetg(n+1,t_VEC); d = (GEN*)D;
     475           0 :   l = lg(E);
     476           0 :   *++d = mkvec2((GEN)1, const_vecsmall(l-1,0));
     477           0 :   for (i=1; i<l; i++)
     478           0 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     479           0 :       for (t2=d, t3=t1; t3<t2; )
     480             :       {
     481             :         ulong a;
     482             :         GEN b;
     483           0 :         a = (*++t3)[1] * P[i];
     484           0 :         b = leafcopy(gel(*t3,2)); b[i]++;
     485           0 :         *++d = mkvec2((GEN)a,b);
     486             :       }
     487           0 :   gen_sort_inplace(D,NULL,&cmpu1,NULL);
     488           0 :   vD = cgetg(n+1, t_VECSMALL);
     489           0 :   for (i = 1; i <= n; i++)
     490             :   {
     491           0 :     vD[i] = umael(D,i,1);
     492           0 :     gel(D,i) = fa_clean_u(P, gmael(D,i,2));
     493             :   }
     494           0 :   return gerepilecopy(av, mkvec2(vD,D));
     495             : }
     496             : GEN
     497      361508 : divisors(GEN N)
     498             : {
     499             :   long i, j, l;
     500             :   GEN *d, *t1, *t2, *t3, D, P, E;
     501      361508 :   int isint = divisors_init(N, &P, &E);
     502      361487 :   GEN (*mul)(GEN,GEN) = isint? mulii: gmul;
     503             : 
     504      361487 :   D = cgetg(ndiv(E)+1,t_VEC); d = (GEN*)D;
     505      361487 :   l = lg(E);
     506      361487 :   *++d = gen_1;
     507     1184449 :   for (i=1; i<l; i++)
     508     2210866 :     for (t1=(GEN*)D,j=E[i]; j; j--,t1=t2)
     509     1387904 :       for (t2=d, t3=t1; t3<t2; ) *++d = mul(*++t3, gel(P,i));
     510      361487 :   if (isint) ZV_sort_inplace(D);
     511      361487 :   return D;
     512             : }
     513             : GEN
     514         784 : divisors0(GEN N, long flag)
     515             : {
     516         784 :   if (flag && flag != 1) pari_err_FLAG("divisors");
     517         784 :   return flag? divisors_factored(N): divisors(N);
     518             : }
     519             : 
     520             : GEN
     521        8253 : divisorsu_moebius(GEN P)
     522             : {
     523             :   GEN d, t, t2, t3;
     524        8253 :   long i, l = lg(P);
     525        8253 :   d = t = cgetg((1 << (l-1)) + 1, t_VECSMALL);
     526        8253 :   *++d = 1;
     527       19600 :   for (i=1; i<l; i++)
     528       11347 :     for (t2=d, t3=t; t3<t2; ) *(++d) = *(++t3) * -P[i];
     529        8253 :   return t;
     530             : }
     531             : GEN
     532    12108101 : divisorsu_fact(GEN fa)
     533             : {
     534    12108101 :   GEN d, t, t1, t2, t3, P = gel(fa,1), E = gel(fa,2);
     535    12108101 :   long i, j, l = lg(P);
     536    12108101 :   d = t = cgetg(numdivu_fact(fa) + 1,t_VECSMALL);
     537    12108101 :   *++d = 1;
     538    40973005 :   for (i=1; i<l; i++)
     539    62821785 :     for (t1=t,j=E[i]; j; j--,t1=t2)
     540    33956881 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     541    12108101 :   vecsmall_sort(t); return t;
     542             : }
     543             : GEN
     544      137274 : divisorsu(ulong N)
     545             : {
     546      137274 :   pari_sp av = avma;
     547      137274 :   return gerepileupto(av, divisorsu_fact(factoru(N)));
     548             : }
     549             : 
     550             : static GEN
     551           0 : corefa(GEN fa)
     552             : {
     553           0 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1;
     554             :   long i;
     555           0 :   for (i=1; i<lg(P); i++)
     556           0 :     if (mod2(gel(E,i))) c = mulii(c,gel(P,i));
     557           0 :   return c;
     558             : }
     559             : static GEN
     560           0 : core2fa(GEN fa)
     561             : {
     562           0 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
     563             :   long i;
     564           0 :   for (i=1; i<lg(P); i++)
     565             :   {
     566           0 :     long e = itos(gel(E,i));
     567           0 :     if (e & 1)  c = mulii(c, gel(P,i));
     568           0 :     if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
     569             :   }
     570           0 :   return mkvec2(c,f);
     571             : }
     572             : GEN
     573           0 : corepartial(GEN n, long all)
     574             : {
     575           0 :   pari_sp av = avma;
     576           0 :   if (typ(n) != t_INT) pari_err_TYPE("corepartial",n);
     577           0 :   return gerepileuptoint(av, corefa(Z_factor_limit(n,all)));
     578             : }
     579             : GEN
     580           0 : core2partial(GEN n, long all)
     581             : {
     582           0 :   pari_sp av = avma;
     583           0 :   if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
     584           0 :   return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
     585             : }
     586             : /* given an arithmetic function argument, return the underlying integer */
     587             : static GEN
     588        3500 : arith_n(GEN A)
     589             : {
     590        3500 :   switch(typ(A))
     591             :   {
     592         728 :     case t_INT: return A;
     593         595 :     case t_VEC: return gel(A,1);
     594        2177 :     default: return factorback(A);
     595             :   }
     596             : }
     597             : static GEN
     598        1484 : core2_i(GEN n)
     599             : {
     600        1484 :   GEN f = core(n);
     601        1484 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     602        1449 :   return mkvec2(f, sqrtint(diviiexact(arith_n(n), f)));
     603             : }
     604             : GEN
     605        1477 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
     606             : 
     607             : GEN
     608        3094 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
     609             : 
     610             : static long
     611         329 : _mod4(GEN c) {
     612         329 :   long r, s = signe(c);
     613         329 :   if (!s) return 0;
     614         329 :   r = mod4(c); if (s < 0) r = 4-r;
     615         329 :   return r;
     616             : }
     617             : 
     618             : long
     619        3816 : corediscs(long D, ulong *f)
     620             : {
     621             :   /* D = f^2 dK */
     622        3816 :   long dK = D>=0 ? (long) coreu(D) : -(long) coreu(-(ulong) D);
     623        3816 :   ulong dKmod4 = ((ulong)dK)&3UL;
     624        3816 :   if (dKmod4 == 2 || dKmod4 == 3)
     625         441 :     dK *= 4;
     626        3816 :   if (f) *f = usqrt((ulong)(D/dK));
     627        3816 :   return D;
     628             : }
     629             : 
     630             : GEN
     631         322 : coredisc(GEN n)
     632             : {
     633         322 :   pari_sp av = avma;
     634         322 :   GEN c = core(n);
     635         322 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     636          70 :   return gerepileuptoint(av, shifti(c,2));
     637             : }
     638             : 
     639             : GEN
     640           7 : coredisc2(GEN n)
     641             : {
     642           7 :   pari_sp av = avma;
     643           7 :   GEN y = core2_i(n);
     644           7 :   GEN c = gel(y,1), f = gel(y,2);
     645           7 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
     646           7 :   y = cgetg(3,t_VEC);
     647           7 :   gel(y,1) = shifti(c,2);
     648           7 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
     649             : }
     650             : 
     651             : GEN
     652          14 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
     653             : 
     654             : long
     655         815 : omegau(ulong n)
     656             : {
     657             :   pari_sp av;
     658         815 :   if (n == 1UL) return 0;
     659         801 :   av = avma; return gc_long(av, nbrows(factoru(n)));
     660             : }
     661             : long
     662        1589 : omega(GEN n)
     663             : {
     664             :   pari_sp av;
     665             :   GEN F, P;
     666        1589 :   if ((F = check_arith_non0(n,"omega"))) {
     667             :     long n;
     668         721 :     P = gel(F,1); n = lg(P)-1;
     669         721 :     if (n && equalim1(gel(P,1))) n--;
     670         721 :     return n;
     671             :   }
     672         854 :   if (lgefint(n) == 3) return omegau(n[2]);
     673          39 :   av = avma;
     674          39 :   F = absZ_factor(n);
     675          39 :   return gc_long(av, nbrows(F));
     676             : }
     677             : 
     678             : long
     679         821 : bigomegau(ulong n)
     680             : {
     681             :   pari_sp av;
     682         821 :   if (n == 1) return 0;
     683         807 :   av = avma; return gc_long(av, zv_sum(gel(factoru(n),2)));
     684             : }
     685             : long
     686        1589 : bigomega(GEN n)
     687             : {
     688        1589 :   pari_sp av = avma;
     689             :   GEN F, E;
     690        1589 :   if ((F = check_arith_non0(n,"bigomega")))
     691             :   {
     692         721 :     GEN P = gel(F,1);
     693         721 :     long n = lg(P)-1;
     694         721 :     E = gel(F,2);
     695         721 :     if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
     696             :   }
     697         854 :   else if (lgefint(n) == 3)
     698         821 :     return bigomegau(n[2]);
     699             :   else
     700          33 :     E = gel(absZ_factor(n), 2);
     701         754 :   E = ZV_to_zv(E);
     702         754 :   return gc_long(av, zv_sum(E));
     703             : }
     704             : 
     705             : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
     706             : ulong
     707      280565 : eulerphiu_fact(GEN f)
     708             : {
     709      280565 :   GEN P = gel(f,1), E = gel(f,2);
     710      280565 :   long i, m = 1, l = lg(P);
     711      882593 :   for (i = 1; i < l; i++)
     712             :   {
     713      602028 :     ulong p = P[i], e = E[i];
     714      602028 :     if (!e) continue;
     715      602028 :     if (p == 2)
     716      229654 :     { if (e > 1) m <<= e-1; }
     717             :     else
     718             :     {
     719      372374 :       m *= (p-1);
     720      372374 :       if (e > 1) m *= upowuu(p, e-1);
     721             :     }
     722             :   }
     723      280565 :   return m;
     724             : }
     725             : ulong
     726      150517 : eulerphiu(ulong n)
     727             : {
     728             :   pari_sp av;
     729      150517 :   if (!n) return 2;
     730      150517 :   av = avma; return gc_long(av, eulerphiu_fact(factoru(n)));
     731             : }
     732             : GEN
     733      146195 : eulerphi(GEN n)
     734             : {
     735      146195 :   pari_sp av = avma;
     736             :   GEN Q, F, P, E;
     737             :   long i, l;
     738             : 
     739      146195 :   if ((F = check_arith_all(n,"eulerphi")))
     740             :   {
     741        1323 :     F = clean_Z_factor(F);
     742        1323 :     n = arith_n(n);
     743        1323 :     if (lgefint(n) == 3)
     744             :     {
     745             :       ulong e;
     746        1309 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
     747        1309 :       e = eulerphiu_fact(F);
     748        1309 :       set_avma(av); return utoipos(e);
     749             :     }
     750             :   }
     751      144872 :   else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
     752             :   else
     753          53 :     F = absZ_factor(n);
     754          67 :   if (!signe(n)) return gen_2;
     755          39 :   P = gel(F,1);
     756          39 :   E = gel(F,2); l = lg(P);
     757          39 :   Q = cgetg(l, t_VEC);
     758         252 :   for (i = 1; i < l; i++)
     759             :   {
     760         213 :     GEN p = gel(P,i), q;
     761         213 :     ulong v = itou(gel(E,i));
     762         213 :     q = subiu(p,1);
     763         213 :     if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
     764         213 :     gel(Q,i) = q;
     765             :   }
     766          39 :   return gerepileuptoint(av, ZV_prod(Q));
     767             : }
     768             : 
     769             : long
     770    12111583 : numdivu_fact(GEN fa)
     771             : {
     772    12111583 :   GEN E = gel(fa,2);
     773    12111583 :   long n = 1, i, l = lg(E);
     774    12111583 :   for (i = 1; i < l; i++) n *= E[i]+1;
     775    12111583 :   return n;
     776             : }
     777             : long
     778         815 : numdivu(long N)
     779             : {
     780             :   pari_sp av;
     781         815 :   if (N == 1) return 1;
     782         801 :   av = avma; return gc_long(av, numdivu_fact(factoru(N)));
     783             : }
     784             : static GEN
     785         760 : numdiv_aux(GEN F)
     786             : {
     787         760 :   GEN x, E = gel(F,2);
     788         760 :   long i, l = lg(E);
     789         760 :   x = cgetg(l, t_VECSMALL);
     790         760 :   for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
     791         760 :   return x;
     792             : }
     793             : GEN
     794        1589 : numdiv(GEN n)
     795             : {
     796        1589 :   pari_sp av = avma;
     797             :   GEN F, E;
     798        1589 :   if ((F = check_arith_non0(n,"numdiv")))
     799             :   {
     800         721 :     F = clean_Z_factor(F);
     801         721 :     E = numdiv_aux(F);
     802             :   }
     803         854 :   else if (lgefint(n) == 3)
     804         815 :     return utoipos(numdivu(n[2]));
     805             :   else
     806          39 :     E = numdiv_aux(absZ_factor(n));
     807         760 :   return gerepileuptoint(av, zv_prod_Z(E));
     808             : }
     809             : 
     810             : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
     811             : static GEN
     812       15950 : u_euler_sumdiv(ulong p, long v)
     813             : {
     814       15950 :   GEN u = utoipos(1 + p); /* can't overflow */
     815       15950 :   for (; v > 1; v--) u = addui(1, mului(p, u));
     816       15950 :   return u;
     817             : }
     818             : /* 1 + q + ... + q^v */
     819             : static GEN
     820       23138 : euler_sumdiv(GEN q, long v)
     821             : {
     822       23138 :   GEN u = addui(1, q);
     823       23138 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
     824       23138 :   return u;
     825             : }
     826             : 
     827             : static GEN
     828        1474 : sumdiv_aux(GEN F)
     829             : {
     830        1474 :   GEN x, P = gel(F,1), E = gel(F,2);
     831        1474 :   long i, l = lg(P);
     832        1474 :   x = cgetg(l, t_VEC);
     833        1474 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
     834        1474 :   return ZV_prod(x);
     835             : }
     836             : GEN
     837        3017 : sumdiv(GEN n)
     838             : {
     839        3017 :   pari_sp av = avma;
     840             :   GEN F, v;
     841             : 
     842        3017 :   if ((F = check_arith_non0(n,"sumdiv")))
     843             :   {
     844        1442 :     F = clean_Z_factor(F);
     845        1442 :     v = sumdiv_aux(F);
     846             :   }
     847        1561 :   else if (lgefint(n) == 3)
     848             :   {
     849        1529 :     if (n[2] == 1) return gen_1;
     850        1501 :     F = factoru(n[2]);
     851        1501 :     v = usumdiv_fact(F);
     852             :   }
     853             :   else
     854          32 :     v = sumdiv_aux(absZ_factor(n));
     855        2975 :   return gerepileuptoint(av, v);
     856             : }
     857             : 
     858             : static GEN
     859        1506 : sumdivk_aux(GEN F, long k)
     860             : {
     861        1506 :   GEN x, P = gel(F,1), E = gel(F,2);
     862        1506 :   long i, l = lg(P);
     863        1506 :   x = cgetg(l, t_VEC);
     864        1506 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(powiu(gel(P,i),k), gel(E,i)[2]);
     865        1506 :   return ZV_prod(x);
     866             : }
     867             : GEN
     868       10423 : sumdivk(GEN n, long k)
     869             : {
     870       10423 :   pari_sp av = avma;
     871             :   GEN F, v;
     872             :   long k1;
     873             : 
     874       10423 :   if (!k) return numdiv(n);
     875       10423 :   if (k == 1) return sumdiv(n);
     876        8834 :   if ((F = check_arith_non0(n,"sumdivk"))) F = clean_Z_factor(F);
     877        8792 :   k1 = k; if (k < 0)  k = -k;
     878        8792 :   if (k == 1)
     879        1428 :     v = sumdiv(F? F: n);
     880             :   else
     881             :   {
     882        7364 :     if (F)
     883        1442 :       v = sumdivk_aux(F,k);
     884        5922 :     else if (lgefint(n) == 3)
     885             :     {
     886        5858 :       if (n[2] == 1) return gen_1;
     887        5697 :       F = factoru(n[2]);
     888        5697 :       v = usumdivk_fact(F,k);
     889             :     }
     890             :     else
     891          64 :       v = sumdivk_aux(absZ_factor(n), k);
     892        7203 :     if (k1 > 0) return gerepileuptoint(av, v);
     893             :   }
     894             : 
     895        1435 :   if (F) n = arith_n(n);
     896        1435 :   if (k != 1) n = powiu(n,k);
     897        1435 :   return gerepileupto(av, gdiv(v, n));
     898             : }
     899             : 
     900             : GEN
     901       11518 : usumdiv_fact(GEN f)
     902             : {
     903       11518 :   GEN P = gel(f,1), E = gel(f,2);
     904       11518 :   long i, l = lg(P);
     905       11518 :   GEN v = cgetg(l, t_VEC);
     906       11518 :   for (i=1; i<l; i++) gel(v,i) = u_euler_sumdiv(P[i],E[i]);
     907       11518 :   return ZV_prod(v);
     908             : }
     909             : GEN
     910       12732 : usumdivk_fact(GEN f, ulong k)
     911             : {
     912       12732 :   GEN P = gel(f,1), E = gel(f,2);
     913       12732 :   long i, l = lg(P);
     914       12732 :   GEN v = cgetg(l, t_VEC);
     915       12732 :   for (i=1; i<l; i++) gel(v,i) = euler_sumdiv(powuu(P[i],k),E[i]);
     916       12732 :   return ZV_prod(v);
     917             : }
     918             : 
     919             : long
     920        1246 : uissquarefree_fact(GEN f)
     921             : {
     922        1246 :   GEN E = gel(f,2);
     923        1246 :   long i, l = lg(E);
     924        2506 :   for (i = 1; i < l; i++)
     925        1883 :     if (E[i] > 1) return 0;
     926         623 :   return 1;
     927             : }
     928             : long
     929       21768 : uissquarefree(ulong n)
     930             : {
     931       21768 :   if (!n) return 0;
     932       21768 :   return moebiusu(n)? 1: 0;
     933             : }
     934             : long
     935        2179 : Z_issquarefree(GEN n)
     936             : {
     937        2179 :   switch(lgefint(n))
     938             :   {
     939          14 :     case 2: return 0;
     940        1406 :     case 3: return uissquarefree(n[2]);
     941             :   }
     942         759 :   return moebius(n)? 1: 0;
     943             : }
     944             : 
     945             : static int
     946        1407 : fa_issquarefree(GEN F)
     947             : {
     948        1407 :   GEN P = gel(F,1), E = gel(F,2);
     949        1407 :   long i, s, l = lg(P);
     950        1407 :   if (l == 1) return 1;
     951        1400 :   s = signe(gel(P,1)); /* = signe(x) */
     952        1400 :   if (!s) return 0;
     953        3577 :   for(i = 1; i < l; i++)
     954        2730 :     if (!equali1(gel(E,i))) return 0;
     955         847 :   return 1;
     956             : }
     957             : long
     958        4046 : issquarefree(GEN x)
     959             : {
     960             :   pari_sp av;
     961             :   GEN d;
     962        4046 :   switch(typ(x))
     963             :   {
     964        1421 :     case t_INT: return Z_issquarefree(x);
     965             :     case t_POL:
     966        1218 :       if (!signe(x)) return 0;
     967        1218 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
     968        1218 :       return gc_long(av, lg(d)==3);
     969             :     case t_VEC:
     970        1407 :     case t_MAT: return fa_issquarefree(check_arith_all(x,"issquarefree"));
     971           0 :     default: pari_err_TYPE("issquarefree",x);
     972             :       return 0; /* LCOV_EXCL_LINE */
     973             :   }
     974             : }
     975             : 
     976             : /*********************************************************************/
     977             : /**                                                                 **/
     978             : /**                    DIGITS / SUM OF DIGITS                       **/
     979             : /**                                                                 **/
     980             : /*********************************************************************/
     981             : 
     982             : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
     983             : static void
     984     2920618 : set_vexp(GEN v, long l)
     985             : {
     986             :   long m;
     987     2920618 :   if (v[l]) return;
     988      992336 :   v[l] = 1; m = l>>1;
     989      992336 :   set_vexp(v, m);
     990      992336 :   set_vexp(v, l-m);
     991             : }
     992             : 
     993             : /* return all needed B^i for DAC algorithm, for lz digits */
     994             : static GEN
     995      935946 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
     996             : {
     997      935946 :   GEN vB, vexp = const_vecsmall(lz, 0);
     998      935946 :   long i, l = (lz+1) >> 1;
     999      935946 :   vexp[1] = 1;
    1000      935946 :   vexp[2] = 1;
    1001      935946 :   set_vexp(vexp, lz);
    1002      935946 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
    1003      935946 :   gel(vB, 1) = T;
    1004     2000494 :   for (i = 2; i <= l; i++)
    1005     1064548 :     if (vexp[i])
    1006             :     {
    1007      992336 :       long j = i >> 1;
    1008      992336 :       GEN B2j = r->sqr(E, gel(vB,j));
    1009      992336 :       gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
    1010             :     }
    1011      935946 :   return vB;
    1012             : }
    1013             : 
    1014             : static void
    1015      206435 : gen_digits_dac(GEN x, GEN vB, long l, GEN *z,
    1016             :                void *E, GEN div(void *E, GEN a, GEN b, GEN *r))
    1017             : {
    1018             :   GEN q, r;
    1019      206435 :   long m = l>>1;
    1020      206435 :   if (l==1) { *z=x; return; }
    1021       98072 :   q = div(E, x, gel(vB,m), &r);
    1022       98072 :   gen_digits_dac(r, vB, m, z, E, div);
    1023       98072 :   gen_digits_dac(q, vB, l-m, z+m, E, div);
    1024             : }
    1025             : 
    1026             : static GEN
    1027       28896 : gen_fromdigits_dac(GEN x, GEN vB, long i, long l, void *E,
    1028             :                    GEN add(void *E, GEN a, GEN b),
    1029             :                    GEN mul(void *E, GEN a, GEN b))
    1030             : {
    1031             :   GEN a, b;
    1032       28896 :   long m = l>>1;
    1033       28896 :   if (l==1) return gel(x,i);
    1034       13181 :   a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
    1035       13181 :   b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
    1036       13181 :   return add(E, a, mul(E, b, gel(vB, m)));
    1037             : }
    1038             : 
    1039             : static GEN
    1040       10795 : gen_digits_i(GEN x, GEN B, long n, void *E, struct bb_ring *r,
    1041             :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
    1042             : {
    1043             :   GEN z, vB;
    1044       10795 :   if (n==1) retmkvec(gcopy(x));
    1045       10291 :   vB = get_vB(B, n, E, r);
    1046       10291 :   z = cgetg(n+1, t_VEC);
    1047       10291 :   gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
    1048       10291 :   return z;
    1049             : }
    1050             : 
    1051             : GEN
    1052       10745 : gen_digits(GEN x, GEN B, long n, void *E, struct bb_ring *r,
    1053             :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
    1054             : {
    1055       10745 :   pari_sp av = avma;
    1056       10745 :   return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
    1057             : }
    1058             : 
    1059             : GEN
    1060        2534 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
    1061             : {
    1062        2534 :   pari_sp av = avma;
    1063        2534 :   long n = lg(x)-1;
    1064        2534 :   GEN vB = get_vB(B, n, E, r);
    1065        2534 :   GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
    1066        2534 :   return gerepilecopy(av, z);
    1067             : }
    1068             : 
    1069             : static GEN
    1070        1771 : _addii(void *data /* ignored */, GEN x, GEN y)
    1071        1771 : { (void)data; return addii(x,y); }
    1072             : static GEN
    1073      958862 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
    1074             : static GEN
    1075      249302 : _mulii(void *data /* ignored */, GEN x, GEN y)
    1076      249302 :  { (void)data; return mulii(x,y); }
    1077             : static GEN
    1078         436 : _dvmdii(void *data /* ignored */, GEN x, GEN y, GEN *r)
    1079         436 :  { (void)data; return dvmdii(x,y,r); }
    1080             : 
    1081             : static struct bb_ring Z_ring = { _addii, _mulii, _sqri };
    1082             : 
    1083             : static GEN
    1084         196 : check_basis(GEN B)
    1085             : {
    1086         196 :   if (!B) return utoipos(10);
    1087         175 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
    1088         175 :   if (abscmpiu(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
    1089         175 :   return B;
    1090             : }
    1091             : 
    1092             : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
    1093             : static void
    1094        3117 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
    1095             : {
    1096        3117 :   pari_sp av = avma;
    1097             :   GEN q,r;
    1098             :   long m;
    1099        3117 :   if (l==1) { *z=itou(x); return; }
    1100        1538 :   m=l>>1;
    1101        1538 :   q = dvmdii(x, gel(vB,m), &r);
    1102        1538 :   digits_dacsmall(q,vB,l-m,z);
    1103        1538 :   digits_dacsmall(r,vB,m,z+l-m);
    1104        1538 :   set_avma(av);
    1105             : }
    1106             : 
    1107             : GEN
    1108          70 : digits(GEN x, GEN B)
    1109             : {
    1110          70 :   pari_sp av=avma;
    1111             :   long lz;
    1112             :   GEN z, vB;
    1113          70 :   if (typ(x)!=t_INT) pari_err_TYPE("digits",x);
    1114          70 :   B = check_basis(B);
    1115          70 :   if (signe(B)<0) pari_err_DOMAIN("digits","B","<",gen_0,B);
    1116          70 :   if (!signe(x))       {set_avma(av); return cgetg(1,t_VEC); }
    1117          63 :   if (abscmpii(x,B)<0) {set_avma(av); retmkvec(absi(x)); }
    1118          63 :   if (Z_ispow2(B))
    1119             :   {
    1120          14 :     long k = expi(B);
    1121          14 :     if (k == 1) return binaire(x);
    1122           7 :     if (k < BITS_IN_LONG)
    1123             :     {
    1124           0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
    1125           0 :       z = binary_2k_nv(x, k);
    1126           0 :       set_avma(av); return Flv_to_ZV(z);
    1127             :     }
    1128             :     else
    1129             :     {
    1130           7 :       set_avma(av); return binary_2k(x, k);
    1131             :     }
    1132             :   }
    1133          49 :   x = absi_shallow(x);
    1134          49 :   lz = logint(x,B) + 1;
    1135          49 :   if (lgefint(B)>3)
    1136             :   {
    1137           8 :     z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
    1138           8 :     vecreverse_inplace(z); return z;
    1139             :   }
    1140             :   else
    1141             :   {
    1142          41 :     vB = get_vB(B, lz, NULL, &Z_ring);
    1143          41 :     (void)new_chunk(3*lz); /* HACK */
    1144          41 :     z = zero_zv(lz);
    1145          41 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1146          41 :     set_avma(av); return Flv_to_ZV(z);
    1147             :   }
    1148             : }
    1149             : 
    1150             : static GEN
    1151     3091822 : fromdigitsu_dac(GEN x, GEN vB, long i, long l)
    1152             : {
    1153             :   GEN a, b;
    1154     3091822 :   long m = l>>1;
    1155     3091822 :   if (l==1) return utoi(uel(x,i));
    1156     2551051 :   if (l==2) return addumului(uel(x,i), uel(x,i+1), gel(vB, m));
    1157     1084371 :   a = fromdigitsu_dac(x, vB, i, m);
    1158     1084371 :   b = fromdigitsu_dac(x, vB, i+m, l-m);
    1159     1084371 :   return addii(a, mulii(b, gel(vB, m)));
    1160             : }
    1161             : 
    1162             : GEN
    1163      923080 : fromdigitsu(GEN x, GEN B)
    1164             : {
    1165      923080 :   pari_sp av = avma;
    1166      923080 :   long n = lg(x)-1;
    1167             :   GEN vB, z;
    1168      923080 :   if (n==0) return gen_0;
    1169      923080 :   vB = get_vB(B, n, NULL, &Z_ring);
    1170      923080 :   z = fromdigitsu_dac(x, vB, 1, n);
    1171      923080 :   return gerepileuptoint(av, z);
    1172             : }
    1173             : 
    1174             : static int
    1175          14 : ZV_in_range(GEN v, GEN B)
    1176             : {
    1177          14 :   long i, l = lg(v);
    1178        1659 :   for(i=1; i < l; i++)
    1179             :   {
    1180        1652 :     GEN vi = gel(v, i);
    1181        1652 :     if (signe(vi) < 0 || cmpii(vi, B) >= 0)
    1182           7 :       return 0;
    1183             :   }
    1184           7 :   return 1;
    1185             : }
    1186             : 
    1187             : GEN
    1188          56 : fromdigits(GEN x, GEN B)
    1189             : {
    1190          56 :   pari_sp av = avma;
    1191          56 :   if (typ(x)!=t_VEC || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
    1192          56 :   if (lg(x)==1) return gen_0;
    1193          49 :   B = check_basis(B);
    1194          49 :   if (Z_ispow2(B) && ZV_in_range(x, B))
    1195           7 :     return fromdigits_2k(x, expi(B));
    1196          42 :   x = vecreverse(x);
    1197          42 :   return gerepileuptoint(av, gen_fromdigits(x, B, NULL, &Z_ring));
    1198             : }
    1199             : 
    1200             : static const ulong digsum[] ={
    1201             :   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,
    1202             :   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,
    1203             :   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,
    1204             :   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,
    1205             :   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,
    1206             :   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,
    1207             :   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,
    1208             :   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,
    1209             :   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,
    1210             :   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,
    1211             :   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,
    1212             :   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,
    1213             :   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,
    1214             :   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,
    1215             :   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,
    1216             :   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,
    1217             :   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,
    1218             :   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,
    1219             :   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,
    1220             :   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,
    1221             :   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,
    1222             :   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,
    1223             :   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,
    1224             :   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,
    1225             :   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,
    1226             :   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,
    1227             :   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,
    1228             :   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,
    1229             :   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,
    1230             :   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,
    1231             :   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,
    1232             :   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,
    1233             :   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,
    1234             :   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,
    1235             :   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,
    1236             :   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,
    1237             :   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,
    1238             :   25,26,27
    1239             : };
    1240             : 
    1241             : ulong
    1242      355152 : sumdigitsu(ulong n)
    1243             : {
    1244      355152 :   ulong s = 0;
    1245      355152 :   while (n) { s += digsum[n % 1000]; n /= 1000; }
    1246      355152 :   return s;
    1247             : }
    1248             : 
    1249             : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1250             : static ulong
    1251          14 : sumdigits_block(ulong *res, long l)
    1252             : {
    1253          14 :   long s = sumdigitsu(*--res);
    1254          14 :   while (--l > 0) s += sumdigitsu(*--res);
    1255          14 :   return s;
    1256             : }
    1257             : 
    1258             : GEN
    1259          35 : sumdigits(GEN n)
    1260             : {
    1261          35 :   pari_sp av = avma;
    1262             :   ulong s, *res;
    1263             :   long l;
    1264             : 
    1265          35 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1266          35 :   l = lgefint(n);
    1267          35 :   switch(l)
    1268             :   {
    1269           7 :     case 2: return gen_0;
    1270          14 :     case 3: return utoipos(sumdigitsu(n[2]));
    1271             :   }
    1272          14 :   res = convi(n, &l);
    1273          14 :   if ((ulong)l < ULONG_MAX / 81)
    1274             :   {
    1275          14 :     s = sumdigits_block(res, l);
    1276          14 :     set_avma(av); return utoipos(s);
    1277             :   }
    1278             :   else /* Huge. Overflows ulong */
    1279             :   {
    1280           0 :     const long L = (long)(ULONG_MAX / 81);
    1281           0 :     GEN S = gen_0;
    1282           0 :     while (l > L)
    1283             :     {
    1284           0 :       S = addiu(S, sumdigits_block(res, L));
    1285           0 :       res += L; l -= L;
    1286             :     }
    1287           0 :     if (l)
    1288           0 :       S = addiu(S, sumdigits_block(res, l));
    1289           0 :     return gerepileuptoint(av, S);
    1290             :   }
    1291             : }
    1292             : 
    1293             : GEN
    1294         105 : sumdigits0(GEN x, GEN B)
    1295             : {
    1296         105 :   pari_sp av = avma;
    1297             :   GEN z;
    1298             :   long lz;
    1299             : 
    1300         105 :   if (!B) return sumdigits(x);
    1301          77 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1302          77 :   B = check_basis(B);
    1303          77 :   if (Z_ispow2(B))
    1304             :   {
    1305          28 :     long k = expi(B);
    1306          28 :     if (k == 1) { set_avma(av); return utoi(hammingweight(x)); }
    1307          21 :     if (k < BITS_IN_LONG)
    1308             :     {
    1309          14 :       GEN z = binary_2k_nv(x, k);
    1310          14 :       if (lg(z)-1 > 1L<<(BITS_IN_LONG-k)) /* may overflow */
    1311           0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1312          14 :       set_avma(av); return utoi(zv_sum(z));
    1313             :     }
    1314           7 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1315             :   }
    1316          49 :   if (!signe(x))       { set_avma(av); return gen_0; }
    1317          49 :   if (abscmpii(x,B)<0) { set_avma(av); return absi(x); }
    1318          49 :   if (absequaliu(B,10))   { set_avma(av); return sumdigits(x); }
    1319          42 :   lz = logint(x,B) + 1;
    1320          42 :   z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
    1321          42 :   return gerepileuptoint(av, ZV_sum(z));
    1322             : }

Generated by: LCOV version 1.13