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 19043-64ae2d8) Lines: 541 581 93.1 %
Date: 2016-06-25 Functions: 75 78 96.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 316 399 79.2 %

           Branch data     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                 :         35 :   return NULL; /* not reached */
      39                 :            : }
      40                 :            : 
      41                 :            : /* NOT memory clean */
      42                 :            : GEN
      43                 :      12481 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
      44                 :            : {
      45                 :      12481 :   long i, j, l = lg(L);
      46                 :      12481 :   GEN e = new_chunk(l), P = new_chunk(l);
      47         [ +  + ]:      87836 :   for (i = j = 1; i < l; i++)
      48                 :            :   {
      49                 :      84140 :     ulong p = uel(L,i);
      50                 :      84140 :     long v = Z_lvalrem(N, p, &N);
      51 [ +  + ][ +  + ]:      84140 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
      52                 :            :   }
      53                 :      12481 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
      54                 :      12481 :   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                 :      28292 : factoru_pow(ulong n)
      65                 :            : {
      66                 :      28292 :   GEN f = cgetg(4,t_VEC);
      67                 :      28291 :   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                 :      28291 :   (void)new_chunk((15 + 1)*3);
      72                 :      28293 :   F = factoru(n);
      73                 :      28293 :   P = gel(F,1);
      74                 :      28293 :   E = gel(F,2); l = lg(P);
      75                 :      28293 :   avma = av;
      76                 :      28293 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
      77                 :      28292 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
      78                 :      28293 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
      79         [ +  + ]:      89613 :   for(i = 1; i < l; i++)
      80                 :            :   {
      81                 :      61320 :     p[i] = P[i];
      82                 :      61320 :     e[i] = E[i];
      83                 :      61320 :     c[i] = upowuu(p[i], e[i]);
      84                 :            :   }
      85                 :      28293 :   return f;
      86                 :            : }
      87                 :            : 
      88                 :            : static GEN
      89                 :      34252 : factorlim(GEN n, ulong lim)
      90         [ -  + ]:      34252 : { 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                 :      16570 : factor_pn_1_limit(GEN p, long n, ulong lim)
      95                 :            : {
      96                 :      16570 :   pari_sp av = avma;
      97                 :      16570 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
      98                 :      16570 :   long i, pp = itos_or_0(p);
      99         [ +  + ]:      29072 :   for(i=2; i<lg(d); i++)
     100                 :            :   {
     101                 :            :     GEN B;
     102 [ +  + ][ +  + ]:      12502 :     if (pp && d[i]%pp==0 && (
                 [ +  + ]
     103 [ -  + ][ +  + ]:      10871 :        ((pp&3)==1 && (d[i]&1)) ||
     104 [ +  + ][ +  + ]:      10864 :        ((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                 :      16570 :   return gerepilecopy(av, A);
     117                 :            : }
     118                 :            : GEN
     119                 :      16570 : factor_pn_1(GEN p, ulong n)
     120                 :      16570 : { 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                 :     790705 : RgV_is_ZVpos(GEN v)
     186                 :            : {
     187                 :     790705 :   long i, l = lg(v);
     188         [ +  + ]:    2341238 :   for (i = 1; i < l; i++)
     189                 :            :   {
     190                 :    1550540 :     GEN c = gel(v,i);
     191 [ +  - ][ +  + ]:    1550540 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     192                 :            :   }
     193                 :     790705 :   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                 :       5166 :   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                 :       3178 :   return 1;
     221                 :            : }
     222                 :            : 
     223                 :            : static int
     224                 :     399535 : is_Z_factor_i(GEN f)
     225 [ +  - ][ +  + ]:     399535 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
                 [ +  + ]
     226                 :            : int
     227                 :     391191 : is_Z_factorpos(GEN f)
     228 [ +  + ][ +  - ]:     391191 : { 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                 :       4501 :   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 associated 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                 :        189 :   return NULL;
     280                 :            : }
     281                 :            : /* n associated 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                 :      14297 :   return NULL;
     298                 :            : }
     299                 :            : /* n associated to a factorization of an integer */
     300                 :            : GEN
     301                 :     187481 : check_arith_all(GEN n, const char *f) {
     302   [ +  +  +  - ]:     187481 :   switch(typ(n))
     303                 :            :   {
     304                 :            :     case t_INT:
     305                 :     184303 :       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                 :     187481 :   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                 :     350756 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     326                 :            : {
     327                 :            :   GEN E, P;
     328         [ -  + ]:     350756 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     329                 :     350756 :   P = gel(F,1);
     330                 :     350756 :   E = gel(F,2);
     331                 :     350756 :   RgV_check_ZV(E, "divisors");
     332                 :     350756 :   *isint = RgV_is_ZV(P);
     333         [ +  + ]:     350756 :   if (*isint)
     334                 :            :   {
     335                 :     350742 :     long i, l = lg(P);
     336                 :            :     /* skip -1 */
     337 [ +  + ][ -  + ]:     350742 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     338                 :            :     /* test for 0 */
     339         [ +  + ]:    1143149 :     for (i = 1; i < l; i++)
     340 [ +  + ][ +  - ]:     792414 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     341                 :          7 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     342                 :            :   }
     343                 :     350749 :   *pP = P;
     344                 :     350749 :   *pE = E;
     345                 :     350749 : }
     346                 :            : static void
     347                 :      10633 : 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         [ -  + ]:      10612 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     360                 :      10612 :       set_fact(absi_factor(n), &P,&E);
     361                 :      10612 :       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                 :     350749 :       set_fact_check(n, &P,&E, &isint);
     368                 :     350742 :       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         [ +  + ]:    1184316 :   for (i=1; i<l; i++)
     376                 :            :   {
     377                 :     822941 :     e[i] = itos(gel(E,i));
     378         [ +  + ]:     822941 :     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         [ +  + ]:    1184113 :   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         [ +  + ]:    1184057 :     for (i=1; i<l; i++)
     401         [ +  + ]:    2210005 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     402         [ +  + ]:    5355693 :         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         [ +  + ]:       1134 :         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                 :     110967 : divisorsu_fact(GEN P, GEN E)
     415                 :            : {
     416                 :     110967 :   long i, j, l = lg(P);
     417                 :     110967 :   ulong nbdiv = 1;
     418                 :            :   GEN d, t, t1, t2, t3;
     419         [ +  + ]:     259184 :   for (i=1; i<l; i++) nbdiv *= 1+E[i];
     420                 :     110967 :   d = t = cgetg(nbdiv+1,t_VECSMALL);
     421                 :     110967 :   *++d = 1;
     422         [ +  + ]:     259184 :   for (i=1; i<l; i++)
     423         [ +  + ]:     337836 :     for (t1=t,j=E[i]; j; j--,t1=t2)
     424         [ +  + ]:     496657 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     425                 :     110967 :   vecsmall_sort(t); return t;
     426                 :            : }
     427                 :            : GEN
     428                 :     110967 : divisorsu(ulong n)
     429                 :            : {
     430                 :     110967 :   pari_sp av = avma;
     431                 :     110967 :   GEN fa = factoru(n);
     432                 :     110967 :   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                 :       1484 :   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                 :        934 : corediscs(long D, ulong *f)
     499                 :            : {
     500                 :            :   /* D = f^2 dK */
     501         [ -  + ]:        934 :   long dK = D>=0 ? (long) coreu(D) : -(long) coreu(-(ulong) D);
     502                 :        934 :   ulong dKmod4 = ((ulong)dK)&3UL;
     503 [ +  + ][ +  + ]:        934 :   if (dKmod4 == 2 || dKmod4 == 3)
     504                 :        105 :     dK *= 4;
     505         [ +  - ]:        934 :   if (f) *f = usqrt((ulong)(D/dK));
     506                 :        934 :   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                 :        815 :   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                 :       1575 :   P = gel(F,1); avma = av; return lg(P)-1;
     557                 :            : }
     558                 :            : 
     559                 :            : long
     560                 :     284244 : bigomegau(ulong n)
     561                 :            : {
     562                 :            :   pari_sp av;
     563                 :            :   GEN F;
     564         [ +  + ]:     284244 :   if (n == 1) return 0;
     565                 :     218458 :   av = avma; F = factoru(n);
     566                 :     284244 :   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                 :       1575 :   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                 :     145607 :   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         [ +  + ]:       2079 :   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         [ +  + ]:       3129 :     for (i=1; i<l; i++) E[i]++;
     680                 :            :   }
     681                 :            :   else
     682                 :         39 :     E = numdiv_aux(absi_factor(n));
     683                 :       1883 :   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         [ +  + ]:       2083 :   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         [ +  + ]:      13093 :   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         [ +  + ]:       2058 :   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                 :       1575 :   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         [ +  + ]:       4116 :   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                 :       6503 :   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         [ +  + ]:       2443 :   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         [ +  + ]:       9142 :   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                 :        280 :   return 1;
     799                 :            : }
     800                 :            : long
     801                 :       9744 : uissquarefree(ulong n)
     802                 :            : {
     803         [ -  + ]:       9744 :   if (!n) return 0;
     804                 :       9744 :   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                 :        772 :   return moebius(n)? 1: 0;
     815                 :            : }
     816                 :            : long
     817                 :      46963 : issquarefree(GEN x)
     818                 :            : {
     819                 :            :   pari_sp av;
     820                 :            :   GEN d;
     821      [ +  +  - ]:      46963 :   switch(typ(x))
     822                 :            :   {
     823                 :         14 :     case t_INT: return Z_issquarefree(x);
     824                 :            :     case t_POL:
     825         [ -  + ]:      46949 :       if (!signe(x)) return 0;
     826                 :      46949 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
     827                 :      46949 :       avma = av; return (lg(d) == 3);
     828                 :          0 :     default: pari_err_TYPE("issquarefree",x);
     829                 :      46963 :       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                 :    1819212 : set_vexp(GEN v, long l)
     842                 :            : {
     843                 :            :   long m;
     844         [ +  + ]:    2429511 :   if (v[l]) return;
     845                 :     610299 :   v[l] = 1; m = l>>1;
     846                 :     610299 :   set_vexp(v, m);
     847                 :     610299 :   set_vexp(v, l-m);
     848                 :            : }
     849                 :            : 
     850                 :            : /* return all needed B^i for DAC algorithm, for lz digits */
     851                 :            : static GEN
     852                 :     598614 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
     853                 :            : {
     854                 :     598614 :   GEN vB, vexp = const_vecsmall(lz, 0);
     855                 :     598614 :   long i, l = (lz+1) >> 1;
     856                 :     598614 :   vexp[1] = 1;
     857                 :     598614 :   vexp[2] = 1;
     858                 :     598614 :   set_vexp(vexp, lz);
     859                 :     598614 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
     860                 :     598614 :   gel(vB, 1) = T;
     861         [ +  + ]:    1259301 :   for (i = 2; i <= l; i++)
     862         [ +  + ]:     660687 :     if (vexp[i])
     863                 :            :     {
     864                 :     610299 :       long j = i >> 1;
     865                 :     610299 :       GEN B2j = r->sqr(E, gel(vB,j));
     866         [ +  + ]:     610299 :       gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
     867                 :            :     }
     868                 :     598614 :   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         [ +  + ]:     246330 :   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                 :      28805 : 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                 :      28805 :   long m = l>>1;
     890         [ +  + ]:      28805 :   if (l==1) return gel(x,i);
     891                 :      13146 :   a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
     892                 :      13146 :   b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
     893                 :      28805 :   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                 :       6280 :   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                 :       2513 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
     918                 :            : {
     919                 :       2513 :   pari_sp av = avma;
     920                 :       2513 :   long n = lg(x)-1;
     921                 :       2513 :   GEN vB = get_vB(B, n, E, r);
     922                 :       2513 :   GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
     923                 :       2513 :   return gerepilecopy(av, z);
     924                 :            : }
     925                 :            : 
     926                 :            : static GEN
     927                 :       1736 : _addii(void *data /* ignored */, GEN x, GEN y)
     928                 :       1736 : { (void)data; return addii(x,y); }
     929                 :            : static GEN
     930                 :     584252 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
     931                 :            : static GEN
     932                 :     150734 : _mulii(void *data /* ignored */, GEN x, GEN y)
     933                 :     150734 :  { (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                 :        168 : check_basis(GEN B)
     942                 :            : {
     943         [ +  + ]:        168 :   if (!B) return utoipos(10);
     944         [ -  + ]:        147 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
     945         [ -  + ]:        147 :   if (cmpis(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
     946                 :        168 :   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         [ +  + ]:       4599 :   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(x))       {avma = av; return cgetg(1,t_VEC); }
     973         [ -  + ]:         49 :   if (absi_cmp(x,B)<0) {avma = av; retmkvec(absi(x)); }
     974         [ +  + ]:         49 :   if (Z_ispow2(B))
     975                 :            :   {
     976                 :         14 :     long k = expi(B);
     977         [ +  + ]:         14 :     if (k == 1) return binaire(x);
     978         [ -  + ]:          7 :     if (k < BITS_IN_LONG)
     979                 :            :     {
     980                 :          0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
     981                 :          0 :       z = binary_2k_nv(x, k);
     982                 :          0 :       avma = av; return Flv_to_ZV(z);
     983                 :            :     }
     984                 :            :     else
     985                 :            :     {
     986                 :          7 :       avma = av; return binary_2k(x, k);
     987                 :            :     }
     988                 :            :   }
     989         [ -  + ]:         35 :   if (signe(x) < 0) x = absi(x);
     990                 :         35 :   lz = logint(x,B,NULL);
     991         [ +  + ]:         35 :   if (lgefint(B)>3)
     992                 :            :   {
     993                 :          8 :     z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
     994                 :          8 :     vecreverse_inplace(z);
     995                 :          8 :     return z;
     996                 :            :   }
     997                 :            :   else
     998                 :            :   {
     999                 :         27 :     vB = get_vB(B, lz, NULL, &Z_ring);
    1000                 :         27 :     (void)new_chunk(3*lz); /* HACK */
    1001                 :         27 :     z = zero_zv(lz);
    1002                 :         27 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1003                 :         56 :     avma = av; return Flv_to_ZV(z);
    1004                 :            :   }
    1005                 :            : }
    1006                 :            : 
    1007                 :            : static GEN
    1008                 :    1908743 : fromdigitsu_dac(GEN x, GEN vB, long i, long l)
    1009                 :            : {
    1010                 :            :   GEN a, b;
    1011                 :    1908743 :   long m = l>>1;
    1012         [ +  + ]:    1908743 :   if (l==1) return utoi(uel(x,i));
    1013         [ +  + ]:    1574687 :   if (l==2) return addumului(uel(x,i), uel(x,i+1), gel(vB, m));
    1014                 :     659429 :   a = fromdigitsu_dac(x, vB, i, m);
    1015                 :     659429 :   b = fromdigitsu_dac(x, vB, i+m, l-m);
    1016                 :    1908743 :   return addii(a, mulii(b, gel(vB, m)));
    1017                 :            : }
    1018                 :            : 
    1019                 :            : GEN
    1020                 :     589885 : fromdigitsu(GEN x, GEN B)
    1021                 :            : {
    1022                 :     589885 :   pari_sp av = avma;
    1023                 :     589885 :   long n = lg(x)-1;
    1024                 :            :   GEN vB, z;
    1025         [ -  + ]:     589885 :   if (n==0) return gen_0;
    1026                 :     589885 :   vB = get_vB(B, n, NULL, &Z_ring);
    1027                 :     589885 :   z = fromdigitsu_dac(x, vB, 1, n);
    1028                 :     589885 :   return gerepileuptoint(av, z);
    1029                 :            : }
    1030                 :            : 
    1031                 :            : GEN
    1032                 :         42 : fromdigits(GEN x, GEN B)
    1033                 :            : {
    1034                 :         42 :   pari_sp av = avma;
    1035 [ +  - ][ -  + ]:         42 :   if (typ(x)!=t_VEC || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
    1036         [ +  + ]:         42 :   if (lg(x)==1) return gen_0;
    1037                 :         35 :   B = check_basis(B);
    1038         [ +  + ]:         35 :   if (Z_ispow2(B))
    1039                 :          7 :     return fromdigits_2k(x, expi(B));
    1040                 :         28 :   x = vecreverse(x);
    1041                 :         42 :   return gerepileuptoint(av, gen_fromdigits(x, B, NULL, &Z_ring));
    1042                 :            : }
    1043                 :            : 
    1044                 :            : static const ulong digsum[] ={
    1045                 :            :   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,
    1046                 :            :   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,
    1047                 :            :   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,
    1048                 :            :   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,
    1049                 :            :   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,
    1050                 :            :   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,
    1051                 :            :   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,
    1052                 :            :   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,
    1053                 :            :   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,
    1054                 :            :   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,
    1055                 :            :   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,
    1056                 :            :   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,
    1057                 :            :   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,
    1058                 :            :   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,
    1059                 :            :   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,
    1060                 :            :   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,
    1061                 :            :   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,
    1062                 :            :   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,
    1063                 :            :   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,
    1064                 :            :   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,
    1065                 :            :   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,
    1066                 :            :   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,
    1067                 :            :   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,
    1068                 :            :   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,
    1069                 :            :   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,
    1070                 :            :   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,
    1071                 :            :   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,
    1072                 :            :   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,
    1073                 :            :   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,
    1074                 :            :   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,
    1075                 :            :   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,
    1076                 :            :   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,
    1077                 :            :   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,
    1078                 :            :   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,
    1079                 :            :   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,
    1080                 :            :   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,
    1081                 :            :   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,
    1082                 :            :   25,26,27
    1083                 :            : };
    1084                 :            : 
    1085                 :            : ulong
    1086                 :     355152 : sumdigitsu(ulong n)
    1087                 :            : {
    1088                 :     355152 :   ulong s = 0;
    1089         [ +  + ]:    1361843 :   while (n) { s += digsum[n % 1000]; n /= 1000; }
    1090                 :     355152 :   return s;
    1091                 :            : }
    1092                 :            : 
    1093                 :            : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1094                 :            : static ulong
    1095                 :         14 : sumdigits_block(ulong *res, long l)
    1096                 :            : {
    1097                 :         14 :   long s = sumdigitsu(*--res);
    1098         [ +  + ]:     355138 :   while (--l > 0) s += sumdigitsu(*--res);
    1099                 :         14 :   return s;
    1100                 :            : }
    1101                 :            : 
    1102                 :            : GEN
    1103                 :         35 : sumdigits(GEN n)
    1104                 :            : {
    1105                 :         35 :   pari_sp av = avma;
    1106                 :            :   ulong s, *res;
    1107                 :            :   long l;
    1108                 :            : 
    1109         [ -  + ]:         35 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1110                 :         35 :   l = lgefint(n);
    1111      [ +  +  + ]:         35 :   switch(l)
    1112                 :            :   {
    1113                 :          7 :     case 2: return gen_0;
    1114                 :         14 :     case 3: return utoipos(sumdigitsu(n[2]));
    1115                 :            :   }
    1116                 :         14 :   res = convi(n, &l);
    1117         [ +  - ]:         14 :   if ((ulong)l < ULONG_MAX / 81)
    1118                 :            :   {
    1119                 :         14 :     s = sumdigits_block(res, l);
    1120                 :         14 :     avma = av; return utoipos(s);
    1121                 :            :   }
    1122                 :            :   else /* Huge. Overflows ulong */
    1123                 :            :   {
    1124                 :          0 :     const long L = (long)(ULONG_MAX / 81);
    1125                 :          0 :     GEN S = gen_0;
    1126         [ #  # ]:          0 :     while (l > L)
    1127                 :            :     {
    1128                 :          0 :       S = addiu(S, sumdigits_block(res, L));
    1129                 :          0 :       res += L; l -= L;
    1130                 :            :     }
    1131         [ #  # ]:          0 :     if (l)
    1132                 :          0 :       S = addiu(S, sumdigits_block(res, l));
    1133                 :         35 :     return gerepileuptoint(av, S);
    1134                 :            :   }
    1135                 :            : }
    1136                 :            : 
    1137                 :            : GEN
    1138                 :        105 : sumdigits0(GEN x, GEN B)
    1139                 :            : {
    1140                 :        105 :   pari_sp av = avma;
    1141                 :            :   GEN z;
    1142                 :            :   long lz;
    1143                 :            : 
    1144         [ +  + ]:        105 :   if (!B) return sumdigits(x);
    1145         [ -  + ]:         77 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1146                 :         77 :   B = check_basis(B);
    1147         [ +  + ]:         77 :   if (Z_ispow2(B))
    1148                 :            :   {
    1149                 :         28 :     long k = expi(B);
    1150         [ +  + ]:         28 :     if (k == 1) { avma = av; return utoi(hammingweight(x)); }
    1151         [ +  + ]:         21 :     if (k < BITS_IN_LONG)
    1152                 :            :     {
    1153                 :         14 :       GEN z = binary_2k_nv(x, k);
    1154         [ -  + ]:         14 :       if (lg(z)-1 > 1L<<(BITS_IN_LONG-k)) /* may overflow */
    1155                 :          0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1156                 :         14 :       avma = av; return utoi(zv_sum(z));
    1157                 :            :     }
    1158                 :          7 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1159                 :            :   }
    1160         [ -  + ]:         49 :   if (!signe(x))       { avma = av; return gen_0; }
    1161         [ -  + ]:         49 :   if (absi_cmp(x,B)<0) { avma = av; return absi(x); }
    1162         [ +  + ]:         49 :   if (equaliu(B,10))   { avma = av; return sumdigits(x); }
    1163                 :         42 :   lz = logint(x,B,NULL);
    1164                 :         42 :   z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
    1165                 :        105 :   return gerepileuptoint(av, ZV_sum(z));
    1166                 :            : }

Generated by: LCOV version 1.9