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-bordeaux1.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 16937-4bd9b4e) Lines: 608 720 84.4 %
Date: 2014-10-24 Functions: 77 86 89.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 332 479 69.3 %

           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                 :            : static ulong _maxprime = 0;
      24                 :            : static ulong diffptrlen;
      25                 :            : 
      26                 :            : /* Building/Rebuilding the diffptr table. The actual work is done by the
      27                 :            :  * following two subroutines;  the user entry point is the function
      28                 :            :  * initprimes() below.  initprimes1() is the old algorithm, called when
      29                 :            :  * maxnum (size) is moderate. Must be called after pari_init_stack() )*/
      30                 :            : static void
      31                 :       1201 : initprimes1(ulong size, long *lenp, ulong *lastp, byteptr p1)
      32                 :            : {
      33                 :       1201 :   pari_sp av = avma;
      34                 :            :   long k;
      35                 :       1201 :   byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
      36                 :            : 
      37         [ +  + ]:      12010 :   for (r=q=p,k=1; r<=fin; )
      38                 :            :   {
      39         [ +  + ]:      16814 :     do { r+=k; k+=2; r+=k; } while (*++q);
      40         [ +  + ]:     381918 :     for (s=r; s<=fin; s+=k) *s = 1;
      41                 :            :   }
      42                 :       1201 :   r = p1; *r++ = 2; *r++ = 1; /* 2 and 3 */
      43                 :       1201 :   for (s=q=p+1; ; s=q)
      44                 :            :   {
      45         [ +  + ]:     423953 :     do q++; while (*q);
      46         [ +  + ]:     150125 :     if (q > fin) break;
      47                 :     148924 :     *r++ = (unsigned char) ((q-s) << 1);
      48                 :     148924 :   }
      49                 :       1201 :   *r++ = 0;
      50                 :       1201 :   *lenp = r - p1;
      51                 :       1201 :   *lastp = ((s - p) << 1) + 1;
      52                 :       1201 :   avma = av;
      53                 :       1201 : }
      54                 :            : 
      55                 :            : /*  Timing in ms (Athlon/850; reports 512K of secondary cache; looks
      56                 :            :     like there is 64K of quickier cache too).
      57                 :            : 
      58                 :            :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
      59                 :            :       =================================================
      60                 :            :       16K       1.1053  1.1407  1.2589  1.4368   1.6086
      61                 :            :       24K       1.0000  1.0625  1.1320  1.2443   1.3095
      62                 :            :       32K       1.0000  1.0469  1.0761  1.1336   1.1776
      63                 :            :       48K       1.0000  1.0000  1.0254  1.0445   1.0546
      64                 :            :       50K       1.0000  1.0000  1.0152  1.0345   1.0464
      65                 :            :       52K       1.0000  1.0000  1.0203  1.0273   1.0362
      66                 :            :       54K       1.0000  1.0000  1.0812  1.0216   1.0281
      67                 :            :       56K       1.0526  1.0000  1.0051  1.0144   1.0205
      68                 :            :       58K       1.0000  1.0000  1.0000  1.0086   1.0123
      69                 :            :       60K       0.9473  0.9844  1.0051  1.0014   1.0055
      70                 :            :       62K       1.0000  0.9844  0.9949  0.9971   0.9993
      71                 :            :       64K       1.0000  1.0000  1.0000  1.0000   1.0000
      72                 :            :       66K       1.2632  1.2187  1.2183  1.2055   1.1953
      73                 :            :       68K       1.4211  1.4844  1.4721  1.4425   1.4188
      74                 :            :       70K       1.7368  1.7188  1.7107  1.6767   1.6421
      75                 :            :       72K       1.9474  1.9531  1.9594  1.9023   1.8573
      76                 :            :       74K       2.2105  2.1875  2.1827  2.1207   2.0650
      77                 :            :       76K       2.4211  2.4219  2.4010  2.3305   2.2644
      78                 :            :       78K       2.5789  2.6250  2.6091  2.5330   2.4571
      79                 :            :       80K       2.8421  2.8125  2.8223  2.7213   2.6380
      80                 :            :       84K       3.1053  3.1875  3.1776  3.0819   2.9802
      81                 :            :       88K       3.5263  3.5312  3.5228  3.4124   3.2992
      82                 :            :       92K       3.7895  3.8438  3.8375  3.7213   3.5971
      83                 :            :       96K       4.0000  4.1093  4.1218  3.9986   3.9659
      84                 :            :       112K      4.3684  4.5781  4.5787  4.4583   4.6115
      85                 :            :       128K      4.7368  4.8750  4.9188  4.8075   4.8997
      86                 :            :       192K      5.5263  5.7188  5.8020  5.6911   5.7064
      87                 :            :       256K      6.0000  6.2187  6.3045  6.1954   6.1033
      88                 :            :       384K      6.7368  6.9531  7.0405  6.9181   6.7912
      89                 :            :       512K      7.3158  7.5156  7.6294  7.5000   7.4654
      90                 :            :       768K      9.1579  9.4531  9.6395  9.5014   9.1075
      91                 :            :       1024K    10.368  10.7497 10.9999 10.878   10.8201
      92                 :            :       1536K    12.579  13.3124 13.7660 13.747   13.4739
      93                 :            :       2048K    13.737  14.4839 15.0509 15.151   15.1282
      94                 :            :       3076K    14.789  15.5780 16.2993 16.513   16.3365
      95                 :            : 
      96                 :            :     Now the same number relative to the model
      97                 :            : 
      98                 :            :     (1 + 0.36*sqrt(primelimit)/arena) * (arena <= 64 ? 1.05 : (arena-64)**0.38)
      99                 :            : 
     100                 :            :      [SLOW2_IN_ROOTS = 0.36, ALPHA = 0.38]
     101                 :            : 
     102                 :            :       arena|    30m     100m    300m    1000m    2000m  <-- primelimit
     103                 :            :       =================================================
     104                 :            :         16K    1.014    0.9835  0.9942  0.9889  1.004
     105                 :            :         24K    0.9526   0.9758  0.9861  0.9942  0.981
     106                 :            :         32K    0.971    0.9939  0.9884  0.9849  0.9806
     107                 :            :         48K    0.9902   0.9825  0.996   0.9945  0.9885
     108                 :            :         50K    0.9917   0.9853  0.9906  0.9926  0.9907
     109                 :            :         52K    0.9932   0.9878  0.9999  0.9928  0.9903
     110                 :            :         54K    0.9945   0.9902  1.064   0.9939  0.9913
     111                 :            :         56K    1.048    0.9924  0.9925  0.993   0.9921
     112                 :            :         58K    0.9969   0.9945  0.9909  0.9932  0.9918
     113                 :            :         60K    0.9455   0.9809  0.9992  0.9915  0.9923
     114                 :            :         62K    0.9991   0.9827  0.9921  0.9924  0.9929
     115                 :            :         64K    1        1       1       1       1
     116                 :            :         66K    1.02     0.9849  0.9857  0.9772  0.9704
     117                 :            :         68K    0.8827   0.9232  0.9176  0.9025  0.8903
     118                 :            :         70K    0.9255   0.9177  0.9162  0.9029  0.8881
     119                 :            :         72K    0.9309   0.936   0.9429  0.9219  0.9052
     120                 :            :         74K    0.9715   0.9644  0.967   0.9477  0.9292
     121                 :            :         76K    0.9935   0.9975  0.9946  0.9751  0.9552
     122                 :            :         78K    0.9987   1.021   1.021   1.003   0.9819
     123                 :            :         80K    1.047    1.041   1.052   1.027   1.006
     124                 :            :         84K    1.052    1.086   1.092   1.075   1.053
     125                 :            :         88K    1.116    1.125   1.133   1.117   1.096
     126                 :            :         92K    1.132    1.156   1.167   1.155   1.134
     127                 :            :         96K    1.137    1.177   1.195   1.185   1.196
     128                 :            :        112K    1.067    1.13    1.148   1.15    1.217
     129                 :            :        128K    1.04     1.083   1.113   1.124   1.178
     130                 :            :        192K    0.9368   0.985   1.025   1.051   1.095
     131                 :            :        256K    0.8741   0.9224  0.9619  0.995   1.024
     132                 :            :        384K    0.8103   0.8533  0.8917  0.9282  0.9568
     133                 :            :        512K    0.7753   0.8135  0.8537  0.892   0.935
     134                 :            :        768K    0.8184   0.8638  0.9121  0.9586  0.9705
     135                 :            :       1024K    0.8241   0.8741  0.927   0.979   1.03
     136                 :            :       1536K    0.8505   0.9212  0.9882  1.056   1.096
     137                 :            :       2048K    0.8294   0.8954  0.9655  1.041   1.102
     138                 :            : 
     139                 :            : */
     140                 :            : 
     141                 :            : #ifndef SLOW2_IN_ROOTS
     142                 :            :   /* SLOW2_IN_ROOTS below 3: some slowdown starts to be noticable
     143                 :            :    * when things fit into the cache on Sparc.
     144                 :            :    * The choice of 2.6 gives a slowdown of 1-2% on UltraSparcII,
     145                 :            :    * but makes calculations for "maximum" of 436273009
     146                 :            :    * fit into 256K cache (still common for some architectures).
     147                 :            :    *
     148                 :            :    * One may change it when small caches become uncommon, but the gain
     149                 :            :    * is not going to be very noticable... */
     150                 :            : #  ifdef i386           /* gcc defines this? */
     151                 :            : #    define SLOW2_IN_ROOTS      0.36
     152                 :            : #  else
     153                 :            : #    define SLOW2_IN_ROOTS      2.6
     154                 :            : #  endif
     155                 :            : #endif
     156                 :            : #ifndef CACHE_ARENA
     157                 :            : #  ifdef i386           /* gcc defines this? */
     158                 :            :    /* Due to smaller SLOW2_IN_ROOTS, smaller arena is OK; fit L1 cache */
     159                 :            : #    define CACHE_ARENA (63 * 1024UL) /* No slowdown even with 64K L1 cache */
     160                 :            : #  else
     161                 :            : #    define CACHE_ARENA (200 * 1024UL) /* No slowdown even with 256K L2 cache */
     162                 :            : #  endif
     163                 :            : #endif
     164                 :            : 
     165                 :            : #define CACHE_ALPHA     (0.38)          /* Cache performance model parameter */
     166                 :            : #define CACHE_CUTOFF    (0.018)         /* Cache performance not smooth here */
     167                 :            : 
     168                 :            : static double slow2_in_roots = SLOW2_IN_ROOTS;
     169                 :            : 
     170                 :            : typedef struct {
     171                 :            :     ulong arena;
     172                 :            :     double power;
     173                 :            :     double cutoff;
     174                 :            : } cache_model_t;
     175                 :            : 
     176                 :            : static cache_model_t cache_model = { CACHE_ARENA, CACHE_ALPHA, CACHE_CUTOFF };
     177                 :            : 
     178                 :            : /* Assume that some calculation requires a chunk of memory to be
     179                 :            :    accessed often in more or less random fashion (as in sieving).
     180                 :            :    Assume that the calculation can be done in steps by subdividing the
     181                 :            :    chunk into smaller subchunks (arenas) and treating them
     182                 :            :    separately.  Assume that the overhead of subdivision is equivalent
     183                 :            :    to the number of arenas.
     184                 :            : 
     185                 :            :    Find an optimal size of the arena taking into account the overhead
     186                 :            :    of subdivision, and the overhead of arena not fitting into the
     187                 :            :    cache.  Assume that arenas of size slow2_in_roots slows down the
     188                 :            :    calculation 2x (comparing to very big arenas; when cache hits do
     189                 :            :    not matter).  Since cache performance varies wildly with
     190                 :            :    architecture, load, and wheather (especially with cache coloring
     191                 :            :    enabled), use an idealized cache model based on benchmarks above.
     192                 :            : 
     193                 :            :    Assume that an independent region of FIXED_TO_CACHE bytes is accessed
     194                 :            :    very often concurrently with the arena access.
     195                 :            :  */
     196                 :            : static ulong
     197                 :       1201 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
     198                 :            :                 cache_model_t *cache_model)
     199                 :            : {
     200                 :       1201 :   ulong asize, cache_arena = cache_model->arena;
     201                 :            :   double Xmin, Xmax, A, B, C1, C2, D, V;
     202                 :       1201 :   double alpha = cache_model->power, cut_off = cache_model->cutoff;
     203                 :            : 
     204                 :            :   /* Estimated relative slowdown,
     205                 :            :      with overhead = max((fixed_to_cache+arena)/cache_arena - 1, 0):
     206                 :            : 
     207                 :            :      1 + slow2_size/arena due to initialization overhead;
     208                 :            : 
     209                 :            :      max(1, 4.63 * overhead^0.38 ) due to footprint > cache size.
     210                 :            : 
     211                 :            :      [The latter is hard to substantiate theoretically, but this
     212                 :            :      function describes benchmarks pretty close; it does not hurt that
     213                 :            :      one can minimize it explicitly too ;-).  The switch between
     214                 :            :      different choices of max() happens when overhead=0.018.]
     215                 :            : 
     216                 :            :      Thus the problem is minimizing (1 + slow2_size/arena)*overhead**0.29.
     217                 :            :      This boils down to F=((X+A)/(X+B))X^alpha, X=overhead,
     218                 :            :      B = (1 - fixed_to_cache/cache_arena), A = B +
     219                 :            :      slow2_size/cache_arena, alpha = 0.38, and X>=0.018, X>-B.
     220                 :            : 
     221                 :            :      We need to find the rightmost root of (X+A)*(X+B) - alpha(A-B)X to the
     222                 :            :      right of 0.018 (if such exists and is below Xmax).  Then we manually
     223                 :            :      check the remaining region [0, 0.018].
     224                 :            : 
     225                 :            :      Since we cannot trust the purely-experimental cache-hit slowdown
     226                 :            :      function, as a sanity check always prefer fitting into the
     227                 :            :      cache (or "almost fitting") if F-law predicts that the larger
     228                 :            :      value of the arena provides less than 10% speedup.
     229                 :            : 
     230                 :            :    */
     231                 :            : 
     232                 :            :   /* The simplest case: we fit into cache */
     233         [ -  + ]:       1201 :   if (total + fixed_to_cache <= cache_arena)
     234                 :          0 :       return total;
     235                 :            :   /* The simple case: fitting into cache doesn't slow us down more than 10% */
     236         [ +  - ]:       1201 :   if (cache_arena - fixed_to_cache > 10 * slow2_size) {
     237                 :       1201 :       asize = cache_arena - fixed_to_cache;
     238         [ -  + ]:       1201 :       if (asize > total) asize = total; /* Automatically false... */
     239                 :       1201 :       return asize;
     240                 :            :   }
     241                 :            :   /* Slowdown of not fitting into cache is significant.  Try to optimize.
     242                 :            :      Do not be afraid to spend some time on optimization - in trivial
     243                 :            :      cases we do not reach this point; any gain we get should
     244                 :            :      compensate the time spent on optimization.  */
     245                 :            : 
     246                 :          0 :   B = (1 - ((double)fixed_to_cache)/cache_arena);
     247                 :          0 :   A = B + ((double)slow2_size)/cache_arena;
     248                 :          0 :   C2 = A*B;
     249                 :          0 :   C1 = (A + B - 1/alpha*(A - B))/2;
     250                 :          0 :   D = C1*C1 - C2;
     251         [ #  # ]:          0 :   if (D > 0)
     252                 :          0 :       V = cut_off*cut_off + 2*C1*cut_off + C2; /* Value at CUT_OFF */
     253                 :            :   else
     254                 :          0 :       V = 0;                            /* Peacify the warning */
     255                 :          0 :   Xmin = cut_off;
     256                 :          0 :   Xmax = ((double)total - fixed_to_cache)/cache_arena; /* Two candidates */
     257                 :            : 
     258 [ #  # ][ #  # ]:          0 :   if ( D <= 0 || (V >= 0 && C1 + cut_off >= 0) ) /* slowdown increasing */
                 [ #  # ]
     259                 :          0 :       Xmax = cut_off;                   /* Only one candidate */
     260 [ #  # ][ #  # ]:          0 :   else if (V >= 0 &&                    /* slowdown concave down */
     261         [ #  # ]:          0 :            ((Xmax + C1) <= 0 || (Xmax*Xmax + 2*C1*Xmax + C2) <= 0))
     262                 :            :       /* DO NOTHING */;                 /* Keep both candidates */
     263 [ #  # ][ #  # ]:          0 :   else if (V <= 0 && (Xmax*Xmax + 2*C1*Xmax + C2) <= 0) /* slowdown decreasing */
     264                 :          0 :       Xmin = cut_off;                   /* Only one candidate */
     265                 :            :   else /* Now we know: 2 roots, the largest is in CUT_OFF..Xmax */
     266                 :          0 :       Xmax = sqrt(D) - C1;
     267         [ #  # ]:          0 :   if (Xmax != Xmin) {   /* Xmin == CUT_OFF; Check which one is better */
     268                 :          0 :       double v1 = (cut_off + A)/(cut_off + B);
     269                 :          0 :       double v2 = 2.33 * (Xmax + A)/(Xmax + B) * pow(Xmax, alpha);
     270                 :            : 
     271         [ #  # ]:          0 :       if (1.1 * v2 >= v1) /* Prefer fitting into the cache if slowdown < 10% */
     272                 :          0 :           V = v1;
     273                 :            :       else {
     274                 :          0 :           Xmin = Xmax;
     275                 :          0 :           V = v2;
     276                 :            :       }
     277         [ #  # ]:          0 :   } else if (B > 0)                     /* We need V */
     278                 :          0 :       V = 2.33 * (Xmin + A)/(Xmin + B) * pow(Xmin, alpha);
     279 [ #  # ][ #  # ]:          0 :   if (B > 0 && 1.1 * V > A/B)  /* Now Xmin is the minumum.  Compare with 0 */
     280                 :          0 :       Xmin = 0;
     281                 :            : 
     282                 :          0 :   asize = (ulong)((1 + Xmin)*cache_arena - fixed_to_cache);
     283         [ #  # ]:          0 :   if (asize > total) asize = total;     /* May happen due to approximations */
     284                 :       1201 :   return asize;
     285                 :            : }
     286                 :            : 
     287                 :            : /* Use as in
     288                 :            :     install(set_optimize,lLDG)          \\ Through some M too?
     289                 :            :     set_optimize(2,1) \\ disable dependence on limit
     290                 :            :     \\ 1: how much cache usable, 2: slowdown of setup, 3: alpha, 4: cutoff
     291                 :            :     \\ 2,3,4 are in units of 0.001
     292                 :            : 
     293                 :            :     { time_primes_arena(ar,limit) =     \\ ar = arena size in K
     294                 :            :         set_optimize(1,floor(ar*1024));
     295                 :            :         default(primelimit, 200 000);   \\ 100000 results in *larger* malloc()!
     296                 :            :         gettime;
     297                 :            :         default(primelimit, floor(limit));
     298                 :            :         if(ar >= 1, ar=floor(ar));
     299                 :            :         print("arena "ar"K => "gettime"ms");
     300                 :            :     }
     301                 :            : */
     302                 :            : long
     303                 :          0 : set_optimize(long what, GEN g)
     304                 :            : {
     305                 :          0 :   long ret = 0;
     306                 :            : 
     307   [ #  #  #  #  :          0 :   switch (what) {
                      # ]
     308                 :            :   case 1:
     309                 :          0 :     ret = (long)cache_model.arena;
     310                 :          0 :     break;
     311                 :            :   case 2:
     312                 :          0 :     ret = (long)(slow2_in_roots * 1000);
     313                 :          0 :     break;
     314                 :            :   case 3:
     315                 :          0 :     ret = (long)(cache_model.power * 1000);
     316                 :          0 :     break;
     317                 :            :   case 4:
     318                 :          0 :     ret = (long)(cache_model.cutoff * 1000);
     319                 :          0 :     break;
     320                 :            :   default:
     321                 :          0 :     pari_err_BUG("set_optimize");
     322                 :          0 :     break;
     323                 :            :   }
     324         [ #  # ]:          0 :   if (g != NULL) {
     325                 :          0 :     ulong val = itou(g);
     326                 :            : 
     327   [ #  #  #  #  :          0 :     switch (what) {
                      # ]
     328                 :          0 :     case 1: cache_model.arena = val; break;
     329                 :          0 :     case 2: slow2_in_roots     = (double)val / 1000.; break;
     330                 :          0 :     case 3: cache_model.power  = (double)val / 1000.; break;
     331                 :          0 :     case 4: cache_model.cutoff = (double)val / 1000.; break;
     332                 :            :     }
     333                 :            :   }
     334                 :          0 :   return ret;
     335                 :            : }
     336                 :            : 
     337                 :            : /* s is odd; prime differences (starting from 5-3=2) start at known_primes[2],
     338                 :            :   terminated by a 0 byte. Checks n odd numbers starting at 'start', setting
     339                 :            :   bytes starting at data to 0 (composite) or 1 (prime) */
     340                 :            : static void
     341                 :       2736 : sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulong n)
     342                 :            : {
     343                 :       2736 :   ulong p, cnt = n-1, start = s, delta = 1;
     344                 :            :   byteptr q;
     345                 :            : 
     346                 :       2736 :   memset(data, 0, n);
     347                 :       2736 :   start >>= 1;  /* (start - 1)/2 */
     348                 :       2736 :   start += n; /* Corresponds to the end */
     349                 :            :   /* data corresponds to start, q runs over primediffs */
     350         [ +  + ]:     320066 :   for (q = known_primes + 1, p = 3; delta; delta = *++q, p += delta)
     351                 :            :   { /* first odd number >= start > p and divisible by p
     352                 :            :        = last odd number <= start + 2p - 2 and 0 (mod p)
     353                 :            :        = p + last number <= start + p - 2 and 0 (mod 2p)
     354                 :            :        = p + start+p-2 - (start+p-2) % 2p
     355                 :            :        = start + 2(p - 1 - ((start-1)/2 + (p-1)/2) % p). */
     356                 :     317330 :     long off = cnt - ((start+(p>>1)) % p);
     357         [ +  + ]:  489638736 :     while (off >= 0) { data[off] = 1; off -= p; }
     358                 :            :   }
     359                 :       2736 : }
     360                 :            : 
     361                 :            : /* assume maxnum <= 436273289 < 2^29 */
     362                 :            : static void
     363                 :       1201 : initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
     364                 :            : {
     365                 :       1201 :   pari_sp av = avma, bot = pari_mainstack->bot;
     366                 :            :   long alloced, psize;
     367                 :            :   byteptr q, end, p, end1, plast, curdiff;
     368                 :            :   ulong last, remains, curlow, rootnum, asize;
     369                 :            :   ulong prime_above;
     370                 :            :   byteptr p_prime_above;
     371                 :            : 
     372                 :       1201 :   maxnum |= 1; /* make it odd. */
     373                 :            :   /* base case */
     374         [ -  + ]:       2402 :   if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, lastp, p1); return; }
     375                 :            : 
     376                 :            :   /* Checked to be enough up to 40e6, attained at 155893 */
     377                 :       1201 :   rootnum = usqrt(maxnum) | 1;
     378                 :       1201 :   initprimes1(rootnum>>1, &psize, &last, p1);
     379                 :       1201 :   end1 = p1 + psize - 1;
     380                 :       1201 :   remains = (maxnum - last) >> 1; /* number of odd numbers to check */
     381                 :            : 
     382                 :            :   /* we access primes array of psize too; but we access it consecutively,
     383                 :            :    * thus we do not include it in fixed_to_cache */
     384                 :       1201 :   asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
     385                 :            :                           &cache_model) - 1;
     386                 :            :   /* enough room on the stack ? */
     387                 :       1201 :   alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
     388         [ -  + ]:       1201 :   if (alloced)
     389                 :          0 :     p = (byteptr)pari_malloc(asize+1);
     390                 :            :   else
     391                 :       1201 :     p = (byteptr)stack_malloc(asize+1);
     392                 :       1201 :   end = p + asize; /* the 0 sentinel goes at end. */
     393                 :       1201 :   curlow = last + 2; /* First candidate: know primes up to last (odd). */
     394                 :       1201 :   curdiff = end1;
     395                 :            : 
     396                 :            :   /* During each iteration p..end-1 represents a range of odd
     397                 :            :      numbers.  plast is a pointer which represents the last prime seen,
     398                 :            :      it may point before p..end-1. */
     399                 :       1201 :   plast = p - 1;
     400                 :       1201 :   p_prime_above = p1 + 2;
     401                 :       1201 :   prime_above = 3;
     402         [ +  + ]:       3937 :   while (remains)
     403                 :            :   { /* cycle over arenas; performance not crucial */
     404                 :            :     unsigned char was_delta;
     405         [ +  + ]:       2736 :     if (asize > remains) { asize = remains; end = p + asize; }
     406                 :            :     /* Fake the upper limit appropriate for the given arena */
     407 [ +  + ][ +  - ]:     152861 :     while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
     408                 :     150125 :       prime_above += *p_prime_above++;
     409                 :       2736 :     was_delta = *p_prime_above;
     410                 :       2736 :     *p_prime_above = 0; /* sentinel for sieve_chunk */
     411                 :       2736 :     sieve_chunk(p1, curlow, p, asize);
     412                 :       2736 :     *p_prime_above = was_delta; /* restore */
     413                 :            : 
     414                 :       2736 :     p[asize] = 0; /* sentinel */
     415                 :       2736 :     for (q = p; ; plast = q++)
     416                 :            :     { /* q runs over addresses corresponding to primes */
     417         [ +  + ]:  299832386 :       while (*q) q++; /* use sentinel at end */
     418         [ +  + ]:   49738548 :       if (q >= end) break;
     419                 :   49735812 :       *curdiff++ = (unsigned char)(q-plast) << 1; /* < 255 for q < 436273291 */
     420                 :   49735812 :     }
     421                 :       2736 :     plast -= asize;
     422                 :       2736 :     remains -= asize;
     423                 :       2736 :     curlow += (asize<<1);
     424                 :            :   }
     425                 :       1201 :   last = curlow - ((p - plast) << 1);
     426                 :       1201 :   *curdiff++ = 0; /* sentinel */
     427                 :       1201 :   *lenp = curdiff - p1;
     428                 :       1201 :   *lastp = last;
     429         [ -  + ]:       1201 :   if (alloced) pari_free(p); else avma = av;
     430                 :            : }
     431                 :            : 
     432                 :            : ulong
     433         [ +  - ]:   13094231 : maxprime(void) { return diffptr ? _maxprime : 0; }
     434                 :            : 
     435                 :            : void
     436         [ #  # ]:          0 : maxprime_check(ulong c) { if (_maxprime < c) pari_err_MAXPRIME(c); }
     437                 :            : 
     438                 :            : /* We ensure 65302 <= maxnum <= 436273289: the LHS ensures modular function
     439                 :            :  * have enough fast primes to work, the RHS ensures that p_{n+1} - p_n < 255
     440                 :            :  * (N.B. RHS would be incorrect since initprimes0 would make it odd, thereby
     441                 :            :  * increasing it by 1) */
     442                 :            : byteptr
     443                 :       1201 : initprimes(ulong maxnum, long *lenp, ulong *lastp)
     444                 :            : {
     445                 :            :   byteptr t;
     446         [ -  + ]:       1201 :   if (maxnum < 65537)
     447                 :          0 :     maxnum = 65537;
     448         [ -  + ]:       1201 :   else if (maxnum > 436273289)
     449                 :          0 :     maxnum = 436273289;
     450                 :       1201 :   t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
     451                 :       1201 :   initprimes0(maxnum, lenp, lastp, t);
     452                 :       1201 :   return (byteptr)pari_realloc(t, *lenp);
     453                 :            : }
     454                 :            : 
     455                 :            : void
     456                 :       1201 : initprimetable(ulong maxnum)
     457                 :            : {
     458                 :            :   long len;
     459                 :            :   ulong last;
     460                 :       1201 :   byteptr p = initprimes(maxnum, &len, &last), old = diffptr;
     461                 :       1201 :   diffptrlen = minss(diffptrlen, len);
     462                 :       1201 :   _maxprime  = minss(_maxprime,last); /*Protect against ^C*/
     463                 :       1201 :   diffptr = p; diffptrlen = len; _maxprime = last;
     464         [ -  + ]:       1201 :   if (old) free(old);
     465                 :       1201 : }
     466                 :            : 
     467                 :            : /* all init_primepointer_xx routines set *ptr to the corresponding place
     468                 :            :  * in prime table */
     469                 :            : /* smallest p >= a */
     470                 :            : ulong
     471                 :          0 : init_primepointer_geq(ulong a, byteptr *pd)
     472                 :            : {
     473                 :            :   ulong n, p;
     474                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     475                 :          0 :   return p;
     476                 :            : }
     477                 :            : /* largest p < a */
     478                 :            : ulong
     479                 :    6080756 : init_primepointer_lt(ulong a, byteptr *pd)
     480                 :            : {
     481                 :            :   ulong n, p;
     482                 :    6080756 :   prime_table_next_p(a, pd, &p, &n);
     483                 :    6080756 :   PREC_PRIME_VIADIFF(p, *pd);
     484                 :    6080756 :   return p;
     485                 :            : }
     486                 :            : /* largest p <= a */
     487                 :            : ulong
     488                 :          0 : init_primepointer_leq(ulong a, byteptr *pd)
     489                 :            : {
     490                 :            :   ulong n, p;
     491                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     492         [ #  # ]:          0 :   if (p != a) PREC_PRIME_VIADIFF(p, *pd);
     493                 :          0 :   return p;
     494                 :            : }
     495                 :            : /* smallest p > a */
     496                 :            : ulong
     497                 :          0 : init_primepointer_gt(ulong a, byteptr *pd)
     498                 :            : {
     499                 :            :   ulong n, p;
     500                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     501         [ #  # ]:          0 :   if (p == a) NEXT_PRIME_VIADIFF(p, *pd);
     502                 :          0 :   return p;
     503                 :            : }
     504                 :            : 
     505                 :            : GEN
     506                 :         35 : boundfact(GEN n, ulong lim)
     507                 :            : {
     508      [ +  +  - ]:         35 :   switch(typ(n))
     509                 :            :   {
     510                 :         21 :     case t_INT: return Z_factor_limit(n,lim);
     511                 :            :     case t_FRAC: {
     512                 :         14 :       pari_sp av = avma;
     513                 :         14 :       GEN a = Z_factor_limit(gel(n,1),lim);
     514                 :         14 :       GEN b = Z_factor_limit(gel(n,2),lim);
     515                 :         14 :       gel(b,2) = ZC_neg(gel(b,2));
     516                 :         14 :       return gerepilecopy(av, merge_factor_i(a,b));
     517                 :            :     }
     518                 :            :   }
     519                 :          0 :   pari_err_TYPE("boundfact",n);
     520                 :         35 :   return NULL; /* not reached */
     521                 :            : }
     522                 :            : 
     523                 :            : /* NOT memory clean */
     524                 :            : GEN
     525                 :       9541 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
     526                 :            : {
     527                 :       9541 :   long i, j, l = lg(L);
     528                 :       9541 :   GEN e = new_chunk(l), P = new_chunk(l);
     529         [ +  + ]:      58618 :   for (i = j = 1; i < l; i++)
     530                 :            :   {
     531                 :      56210 :     ulong p = uel(L,i);
     532                 :      56210 :     long v = Z_lvalrem(N, p, &N);
     533 [ +  + ][ +  + ]:      56210 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
     534                 :            :   }
     535                 :       9541 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
     536                 :       9541 :   e[0] = evaltyp(t_VECSMALL) | evallg(j); *pe = e; return N;
     537                 :            : }
     538                 :            : /***********************************************************************/
     539                 :            : /**                                                                   **/
     540                 :            : /**                    SIMPLE FACTORISATIONS                          **/
     541                 :            : /**                                                                   **/
     542                 :            : /***********************************************************************/
     543                 :            : /* Factor n and output [p,e,c] where
     544                 :            :  * p, e and c are vecsmall with n = prod{p[i]^e[i]} and c[i] = p[i]^e[i] */
     545                 :            : GEN
     546                 :      22206 : factoru_pow(ulong n)
     547                 :            : {
     548                 :      22206 :   GEN f = cgetg(4,t_VEC);
     549                 :      22206 :   pari_sp av = avma;
     550                 :            :   GEN F, P, E, p, e, c;
     551                 :            :   long i, l;
     552                 :            :   /* enough room to store <= 15 * [p,e,p^e] (OK if n < 2^64) */
     553                 :      22206 :   (void)new_chunk((15 + 1)*3);
     554                 :      22206 :   F = factoru(n);
     555                 :      22206 :   P = gel(F,1);
     556                 :      22206 :   E = gel(F,2); l = lg(P);
     557                 :      22206 :   avma = av;
     558                 :      22206 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
     559                 :      22206 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
     560                 :      22206 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
     561         [ +  + ]:      71565 :   for(i = 1; i < l; i++)
     562                 :            :   {
     563                 :      49359 :     p[i] = P[i];
     564                 :      49359 :     e[i] = E[i];
     565                 :      49359 :     c[i] = upowuu(p[i], e[i]);
     566                 :            :   }
     567                 :      22206 :   return f;
     568                 :            : }
     569                 :            : 
     570                 :            : static GEN
     571                 :      23730 : factorlim(GEN n, ulong lim)
     572         [ -  + ]:      23730 : { return lim? Z_factor_limit(n, lim): Z_factor(n); }
     573                 :            : /* factor p^n - 1, assuming p prime. If lim != 0, limit factorization to
     574                 :            :  * primes <= lim */
     575                 :            : GEN
     576                 :       6622 : factor_pn_1_limit(GEN p, long n, ulong lim)
     577                 :            : {
     578                 :       6622 :   pari_sp av = avma;
     579                 :       6622 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
     580                 :       6622 :   long i, pp = itos_or_0(p);
     581         [ +  + ]:      18564 :   for(i=2; i<lg(d); i++)
     582                 :            :   {
     583                 :            :     GEN B;
     584 [ +  + ][ +  + ]:      11942 :     if (pp && d[i]%pp==0 && (
                 [ -  + ]
     585 [ #  # ][ +  + ]:      10717 :        ((pp&3)==1 && (d[i]&1)) ||
     586 [ +  + ][ +  + ]:      10717 :        ((pp&3)==3 && (d[i]&3)==2) ||
     587         [ +  + ]:      10640 :        (pp==2 && (d[i]&7)==4)))
     588                 :       5166 :     {
     589                 :       5166 :       GEN f=factor_Aurifeuille_prime(p,d[i]);
     590                 :       5166 :       B = factorlim(f, lim);
     591                 :       5166 :       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     592                 :       5166 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     593                 :            :     }
     594                 :            :     else
     595                 :       6776 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     596                 :      11942 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     597                 :            :   }
     598                 :       6622 :   return gerepilecopy(av, A);
     599                 :            : }
     600                 :            : GEN
     601                 :       6622 : factor_pn_1(GEN p, ulong n)
     602                 :       6622 : { return factor_pn_1_limit(p, n, 0); }
     603                 :            : 
     604                 :            : #if 0
     605                 :            : static GEN
     606                 :            : to_mat(GEN p, long e) {
     607                 :            :   GEN B = cgetg(3, t_MAT);
     608                 :            :   gel(B,1) = mkcol(p);
     609                 :            :   gel(B,2) = mkcol(utoipos(e)); return B;
     610                 :            : }
     611                 :            : /* factor phi(n) */
     612                 :            : GEN
     613                 :            : factor_eulerphi(GEN n)
     614                 :            : {
     615                 :            :   pari_sp av = avma;
     616                 :            :   GEN B = NULL, A, P, E, AP, AE;
     617                 :            :   long i, l, v = vali(n);
     618                 :            : 
     619                 :            :   l = lg(n);
     620                 :            :   /* result requires less than this: at most expi(n) primes */
     621                 :            :   (void)new_chunk(bit_accuracy(l) * (l /*p*/ + 3 /*e*/ + 2 /*vectors*/) + 3+2);
     622                 :            :   if (v) { n = shifti(n, -v); v--; }
     623                 :            :   A = Z_factor(n); P = gel(A,1); E = gel(A,2); l = lg(P);
     624                 :            :   for(i = 1; i < l; i++)
     625                 :            :   {
     626                 :            :     GEN p = gel(P,i), q = subis(p,1), fa;
     627                 :            :     long e = itos(gel(E,i)), w;
     628                 :            : 
     629                 :            :     w = vali(q); v += w; q = shifti(q,-w);
     630                 :            :     if (! is_pm1(q))
     631                 :            :     {
     632                 :            :       fa = Z_factor(q);
     633                 :            :       B = B? merge_factor(B, fa, (void*)&cmpii, cmp_nodata): fa;
     634                 :            :     }
     635                 :            :     if (e > 1) {
     636                 :            :       if (B) {
     637                 :            :         gel(B,1) = shallowconcat(gel(B,1), p);
     638                 :            :         gel(B,2) = shallowconcat(gel(B,2), utoipos(e-1));
     639                 :            :       } else
     640                 :            :         B = to_mat(p, e-1);
     641                 :            :     }
     642                 :            :   }
     643                 :            :   avma = av;
     644                 :            :   if (!B) return v? to_mat(gen_2, v): trivial_fact();
     645                 :            :   A = cgetg(3, t_MAT);
     646                 :            :   P = gel(B,1); E = gel(B,2); l = lg(P);
     647                 :            :   AP = cgetg(l+1, t_COL); gel(A,1) = AP; AP++;
     648                 :            :   AE = cgetg(l+1, t_COL); gel(A,2) = AE; AE++;
     649                 :            :   /* prepend "2^v" */
     650                 :            :   gel(AP,0) = gen_2;
     651                 :            :   gel(AE,0) = utoipos(v);
     652                 :            :   for (i = 1; i < l; i++)
     653                 :            :   {
     654                 :            :     gel(AP,i) = icopy(gel(P,i));
     655                 :            :     gel(AE,i) = icopy(gel(E,i));
     656                 :            :   }
     657                 :            :   return A;
     658                 :            : }
     659                 :            : #endif
     660                 :            : 
     661                 :            : /***********************************************************************/
     662                 :            : /**                                                                   **/
     663                 :            : /**         CHECK FACTORIZATION FOR ARITHMETIC FUNCTIONS              **/
     664                 :            : /**                                                                   **/
     665                 :            : /***********************************************************************/
     666                 :            : static int
     667                 :     271232 : RgV_is_ZVpos(GEN v)
     668                 :            : {
     669                 :     271232 :   long i, l = lg(v);
     670         [ +  + ]:     709495 :   for (i = 1; i < l; i++)
     671                 :            :   {
     672                 :     438270 :     GEN c = gel(v,i);
     673 [ +  - ][ +  + ]:     438270 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     674                 :            :   }
     675                 :     271232 :   return 1;
     676                 :            : }
     677                 :            : /* check whether v is a ZV with non-0 entries */
     678                 :            : static int
     679                 :        259 : RgV_is_ZVnon0(GEN v)
     680                 :            : {
     681                 :        259 :   long i, l = lg(v);
     682         [ +  + ]:        868 :   for (i = 1; i < l; i++)
     683                 :            :   {
     684                 :        658 :     GEN c = gel(v,i);
     685 [ +  - ][ +  + ]:        658 :     if (typ(c) != t_INT || !signe(c)) return 0;
     686                 :            :   }
     687                 :        259 :   return 1;
     688                 :            : }
     689                 :            : /* check whether v is a ZV with non-zero entries OR exactly [0] */
     690                 :            : static int
     691                 :        210 : RgV_is_ZV0(GEN v)
     692                 :            : {
     693                 :        210 :   long i, l = lg(v);
     694         [ +  + ]:        581 :   for (i = 1; i < l; i++)
     695                 :            :   {
     696                 :        427 :     GEN c = gel(v,i);
     697                 :            :     long s;
     698         [ -  + ]:        427 :     if (typ(c) != t_INT) return 0;
     699                 :        427 :     s = signe(c);
     700         [ +  + ]:        427 :     if (!s) return (l == 2);
     701                 :            :   }
     702                 :        210 :   return 1;
     703                 :            : }
     704                 :            : 
     705                 :            : static int
     706                 :     135861 : is_Z_factor_i(GEN f)
     707 [ +  - ][ +  + ]:     135861 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
                 [ +  + ]
     708                 :            : int
     709                 :     135392 : is_Z_factorpos(GEN f)
     710 [ +  + ][ +  - ]:     135392 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     711                 :            : int
     712                 :        210 : is_Z_factor(GEN f)
     713 [ +  - ][ +  - ]:        210 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     714                 :            : /* as is_Z_factorpos, also allow factor(0) */
     715                 :            : int
     716                 :        259 : is_Z_factornon0(GEN f)
     717 [ +  - ][ +  + ]:        259 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     718                 :            : GEN
     719                 :        168 : clean_Z_factor(GEN f)
     720                 :            : {
     721                 :        168 :   GEN P = gel(f,1);
     722                 :        168 :   long n = lg(P)-1;
     723 [ +  + ][ +  + ]:        168 :   if (n && equalim1(gel(P,1)))
     724                 :         56 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     725                 :        168 :   return f;
     726                 :            : }
     727                 :            : GEN
     728                 :          0 : fuse_Z_factor(GEN f, GEN B)
     729                 :            : {
     730                 :          0 :   GEN P = gel(f,1), E = gel(f,2), P2,E2;
     731                 :          0 :   long i, l = lg(P);
     732         [ #  # ]:          0 :   if (l == 1) return f;
     733         [ #  # ]:          0 :   for (i = 1; i < l; i++)
     734         [ #  # ]:          0 :     if (absi_cmp(gel(P,i), B) > 0) break;
     735         [ #  # ]:          0 :   if (i == l) return f;
     736                 :            :   /* tail / initial segment */
     737                 :          0 :   P2 = vecslice(P, i, l-1); P = vecslice(P, 1, i-1);
     738                 :          0 :   E2 = vecslice(E, i, l-1); E = vecslice(E, 1, i-1);
     739                 :          0 :   P = shallowconcat(P, mkvec(factorback2(P2,E2)));
     740                 :          0 :   E = shallowconcat(E, mkvec(gen_1));
     741                 :          0 :   return mkmat2(P, E);
     742                 :            : }
     743                 :            : 
     744                 :            : /* n associated to a factorization of a positive integer: either N (t_INT)
     745                 :            :  * a factorization matrix faN, or a t_VEC: [N, faN] */
     746                 :            : GEN
     747                 :          0 : check_arith_pos(GEN n, const char *f) {
     748   [ #  #  #  # ]:          0 :   switch(typ(n))
     749                 :            :   {
     750                 :            :     case t_INT:
     751         [ #  # ]:          0 :       if (signe(n) <= 0 ) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
     752                 :          0 :       return NULL;
     753                 :            :     case t_VEC:
     754 [ #  # ][ #  # ]:          0 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     755                 :          0 :       n = gel(n,2); /* fall through */
     756                 :            :     case t_MAT:
     757         [ #  # ]:          0 :       if (!is_Z_factorpos(n)) break;
     758                 :          0 :       return n;
     759                 :            :   }
     760                 :          0 :   pari_err_TYPE(f,n);
     761                 :          0 :   return NULL;
     762                 :            : }
     763                 :            : /* n associated to a factorization of a non-0 integer */
     764                 :            : GEN
     765                 :       2494 : check_arith_non0(GEN n, const char *f) {
     766   [ +  +  +  - ]:       2494 :   switch(typ(n))
     767                 :            :   {
     768                 :            :     case t_INT:
     769         [ +  + ]:       2284 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     770                 :       2228 :       return NULL;
     771                 :            :     case t_VEC:
     772 [ +  - ][ +  - ]:          7 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     773                 :          7 :       n = gel(n,2); /* fall through */
     774                 :            :     case t_MAT:
     775         [ +  + ]:        210 :       if (!is_Z_factornon0(n)) break;
     776                 :        161 :       return n;
     777                 :            :   }
     778                 :         49 :   pari_err_TYPE(f,n);
     779                 :       2389 :   return NULL;
     780                 :            : }
     781                 :            : /* n associated to a factorization of an integer */
     782                 :            : GEN
     783                 :     169414 : check_arith_all(GEN n, const char *f) {
     784   [ +  +  +  - ]:     169414 :   switch(typ(n))
     785                 :            :   {
     786                 :            :     case t_INT:
     787                 :     169204 :       return NULL;
     788                 :            :     case t_VEC:
     789 [ +  - ][ +  - ]:          7 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     790                 :          7 :       n = gel(n,2); /* fall through */
     791                 :            :     case t_MAT:
     792         [ -  + ]:        210 :       if (!is_Z_factor(n)) break;
     793                 :        210 :       return n;
     794                 :            :   }
     795                 :          0 :   pari_err_TYPE(f,n);
     796                 :     169414 :   return NULL;
     797                 :            : }
     798                 :            : 
     799                 :            : /***********************************************************************/
     800                 :            : /**                                                                   **/
     801                 :            : /**                MISCELLANEOUS ARITHMETIC FUNCTIONS                 **/
     802                 :            : /**                (ultimately depend on Z_factor())                  **/
     803                 :            : /**                                                                   **/
     804                 :            : /***********************************************************************/
     805                 :            : /* set P,E from F. Check whether F is an integer and kill "factor" -1 */
     806                 :            : static void
     807                 :     350749 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     808                 :            : {
     809                 :            :   GEN E, P;
     810         [ -  + ]:     350749 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     811                 :     350749 :   P = gel(F,1);
     812                 :     350749 :   E = gel(F,2);
     813                 :     350749 :   RgV_check_ZV(E, "divisors");
     814                 :     350749 :   *isint = RgV_is_ZV(P);
     815         [ +  + ]:     350749 :   if (*isint)
     816                 :            :   {
     817                 :     350735 :     long i, l = lg(P);
     818                 :            :     /* skip -1 */
     819 [ +  + ][ -  + ]:     350735 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     820                 :            :     /* test for 0 */
     821         [ +  + ]:    1143135 :     for (i = 1; i < l; i++)
     822 [ +  + ][ +  - ]:     792407 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     823                 :          7 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     824                 :            :   }
     825                 :     350742 :   *pP = P;
     826                 :     350742 :   *pE = E;
     827                 :     350742 : }
     828                 :            : static void
     829                 :      10612 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     830                 :            : 
     831                 :            : int
     832                 :     361361 : divisors_init(GEN n, GEN *pP, GEN *pE)
     833                 :            : {
     834                 :            :   long i,l;
     835                 :            :   GEN E, P, e;
     836                 :            :   int isint;
     837                 :            : 
     838   [ +  +  +  + ]:     361361 :   switch(typ(n))
     839                 :            :   {
     840                 :            :     case t_INT:
     841         [ -  + ]:      10591 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     842                 :      10591 :       set_fact(absi_factor(n), &P,&E);
     843                 :      10591 :       isint = 1; break;
     844                 :            :     case t_VEC:
     845         [ -  + ]:          7 :       if (lg(n) != 3) pari_err_TYPE("divisors",n);
     846                 :          7 :       set_fact_check(gel(n,2), &P,&E, &isint);
     847                 :          7 :       break;
     848                 :            :     case t_MAT:
     849                 :     350742 :       set_fact_check(n, &P,&E, &isint);
     850                 :     350735 :       break;
     851                 :            :     default:
     852                 :         21 :       set_fact(factor(n), &P,&E);
     853                 :         21 :       isint = 0; break;
     854                 :            :   }
     855                 :     361354 :   l = lg(P);
     856                 :     361354 :   e = cgetg(l, t_VECSMALL);
     857         [ +  + ]:    1184253 :   for (i=1; i<l; i++)
     858                 :            :   {
     859                 :     822906 :     e[i] = itos(gel(E,i));
     860         [ +  + ]:     822906 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     861                 :            :   }
     862                 :     361347 :   *pP = P; *pE = e; return isint;
     863                 :            : }
     864                 :            : 
     865                 :            : GEN
     866                 :     361319 : divisors(GEN n)
     867                 :            : {
     868                 :     361319 :   pari_sp av = avma;
     869                 :            :   long i, j, l;
     870                 :            :   ulong ndiv;
     871                 :            :   GEN *d, *t, *t1, *t2, *t3, P, E, e;
     872                 :     361319 :   int isint = divisors_init(n, &P, &E);
     873                 :            : 
     874                 :     361305 :   l = lg(E); e = cgetg(l, t_VECSMALL);
     875         [ +  + ]:    1184050 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     876                 :     361305 :   ndiv = itou_or_0( zv_prod_Z(e) );
     877 [ +  - ][ -  + ]:     361305 :   if (!ndiv || ndiv & ~LGBITS) pari_err_OVERFLOW("divisors");
     878                 :     361305 :   d = t = (GEN*) cgetg(ndiv+1,t_VEC);
     879                 :     361305 :   *++d = gen_1;
     880         [ +  + ]:     361305 :   if (isint)
     881                 :            :   {
     882         [ +  + ]:    1183994 :     for (i=1; i<l; i++)
     883         [ +  + ]:    2209928 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     884         [ +  + ]:    5355595 :         for (t2=d, t3=t1; t3<t2; ) *++d = mulii(*++t3, gel(P,i));
     885                 :     361291 :     e = ZV_sort((GEN)t);
     886                 :            :   } else {
     887         [ +  + ]:         56 :     for (i=1; i<l; i++)
     888         [ +  + ]:        182 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     889         [ +  + ]:       1134 :         for (t2=d, t3=t1; t3<t2; ) *++d = gmul(*++t3, gel(P,i));
     890                 :         14 :     e = (GEN)t;
     891                 :            :   }
     892                 :     361305 :   return gerepileupto(av, e);
     893                 :            : }
     894                 :            : 
     895                 :            : GEN
     896                 :      30115 : divisorsu(ulong n)
     897                 :            : {
     898                 :      30115 :   pari_sp av = avma;
     899                 :            :   long i, j, l;
     900                 :            :   ulong nbdiv;
     901                 :      30115 :   GEN d, t, t1, t2, t3, P, e, fa = factoru(n);
     902                 :            : 
     903                 :      30115 :   P = gel(fa,1); l = lg(P);
     904                 :      30115 :   e = gel(fa,2);
     905                 :      30115 :   nbdiv = 1;
     906         [ +  + ]:      66253 :   for (i=1; i<l; i++) nbdiv *= 1+e[i];
     907                 :      30115 :   d = t = cgetg(nbdiv+1,t_VECSMALL);
     908                 :      30115 :   *++d = 1;
     909         [ +  + ]:      66253 :   for (i=1; i<l; i++)
     910         [ +  + ]:      81819 :     for (t1=t,j=e[i]; j; j--,t1=t2)
     911         [ +  + ]:     117039 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     912                 :      30115 :   vecsmall_sort(t);
     913                 :      30115 :   return gerepileupto(av, t);
     914                 :            : }
     915                 :            : 
     916                 :            : static GEN
     917                 :          0 : corefa(GEN fa)
     918                 :            : {
     919                 :          0 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1;
     920                 :            :   long i;
     921         [ #  # ]:          0 :   for (i=1; i<lg(P); i++)
     922         [ #  # ]:          0 :     if (mod2(gel(E,i))) c = mulii(c,gel(P,i));
     923                 :          0 :   return c;
     924                 :            : }
     925                 :            : static GEN
     926                 :        154 : core2fa(GEN fa)
     927                 :            : {
     928                 :        154 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
     929                 :            :   long i;
     930         [ +  + ]:        413 :   for (i=1; i<lg(P); i++)
     931                 :            :   {
     932                 :        259 :     long e = itos(gel(E,i));
     933         [ +  + ]:        259 :     if (e & 1)  c = mulii(c, gel(P,i));
     934         [ +  + ]:        259 :     if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
     935                 :            :   }
     936                 :        154 :   return mkvec2(c,f);
     937                 :            : }
     938                 :            : GEN
     939                 :          0 : corepartial(GEN n, long all)
     940                 :            : {
     941                 :          0 :   pari_sp av = avma;
     942         [ #  # ]:          0 :   if (typ(n) != t_INT) pari_err_TYPE("corepartial",n);
     943                 :          0 :   return gerepileuptoint(av, corefa(Z_factor_limit(n,all)));
     944                 :            : }
     945                 :            : GEN
     946                 :        154 : core2partial(GEN n, long all)
     947                 :            : {
     948                 :        154 :   pari_sp av = avma;
     949         [ -  + ]:        154 :   if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
     950                 :        154 :   return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
     951                 :            : }
     952                 :            : static GEN
     953                 :         70 : core2_i(GEN n)
     954                 :            : {
     955                 :         70 :   GEN f = core(n);
     956         [ +  + ]:         70 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     957      [ +  +  + ]:         49 :   switch(typ(n))
     958                 :            :   {
     959                 :          7 :     case t_VEC: n = gel(n,1); break;
     960                 :         21 :     case t_MAT: n = factorback(n); break;
     961                 :            :   }
     962                 :         70 :   return mkvec2(f, sqrtint(diviiexact(n, f)));
     963                 :            : }
     964                 :            : GEN
     965                 :         63 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
     966                 :            : 
     967                 :            : GEN
     968         [ +  + ]:        266 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
     969                 :            : 
     970                 :            : static long
     971                 :         14 : _mod4(GEN c) {
     972                 :         14 :   long r, s = signe(c);
     973         [ -  + ]:         14 :   if (!s) return 0;
     974         [ -  + ]:         14 :   r = mod4(c); if (s < 0) r = 4-r;
     975                 :         14 :   return r;
     976                 :            : }
     977                 :            : 
     978                 :            : GEN
     979                 :          7 : coredisc(GEN n)
     980                 :            : {
     981                 :          7 :   pari_sp av = avma;
     982                 :          7 :   GEN c = core(n);
     983         [ -  + ]:          7 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     984                 :          7 :   return gerepileuptoint(av, shifti(c,2));
     985                 :            : }
     986                 :            : 
     987                 :            : GEN
     988                 :          7 : coredisc2(GEN n)
     989                 :            : {
     990                 :          7 :   pari_sp av = avma;
     991                 :          7 :   GEN y = core2_i(n);
     992                 :          7 :   GEN c = gel(y,1), f = gel(y,2);
     993         [ -  + ]:          7 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
     994                 :          7 :   y = cgetg(3,t_VEC);
     995                 :          7 :   gel(y,1) = shifti(c,2);
     996                 :          7 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
     997                 :            : }
     998                 :            : 
     999                 :            : GEN
    1000         [ +  + ]:         14 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
    1001                 :            : 
    1002                 :            : long
    1003                 :        189 : omega(GEN n)
    1004                 :            : {
    1005                 :        189 :   pari_sp av = avma;
    1006                 :            :   GEN F, P;
    1007         [ +  + ]:        189 :   if ((F = check_arith_non0(n,"omega"))) {
    1008                 :            :     long n;
    1009                 :         21 :     P = gel(F,1); n = lg(P)-1;
    1010 [ +  + ][ +  + ]:         21 :     if (n && equalim1(gel(P,1))) n--;
    1011                 :         21 :     return n;
    1012                 :            :   }
    1013         [ +  + ]:        154 :   else if (lgefint(n) == 3)
    1014                 :            :   {
    1015         [ -  + ]:        115 :     if (n[2] == 1) return 0;
    1016                 :        115 :     F = factoru(n[2]);
    1017                 :            :   }
    1018                 :            :   else
    1019                 :         39 :     F = absi_factor(n);
    1020                 :        175 :   P = gel(F,1); avma = av; return lg(P)-1;
    1021                 :            : }
    1022                 :            : 
    1023                 :            : long
    1024                 :        189 : bigomega(GEN n)
    1025                 :            : {
    1026                 :        189 :   pari_sp av = avma;
    1027                 :            :   GEN F, E;
    1028         [ +  + ]:        189 :   if ((F = check_arith_non0(n,"bigomega")))
    1029                 :            :   {
    1030                 :         21 :     GEN P = gel(F,1);
    1031                 :         21 :     long n = lg(P)-1;
    1032                 :         21 :     E = gel(F,2);
    1033 [ +  + ][ +  + ]:         21 :     if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
    1034                 :         21 :     E = ZV_to_zv(E);
    1035                 :            :   }
    1036         [ +  + ]:        154 :   else if (lgefint(n) == 3)
    1037                 :            :   {
    1038         [ -  + ]:        121 :     if (n[2] == 1) return 0;
    1039                 :        121 :     F = factoru(n[2]);
    1040                 :        121 :     E = gel(F,2);
    1041                 :            :   }
    1042                 :            :   else
    1043                 :         33 :     E = ZV_to_zv(gel(absi_factor(n), 2));
    1044                 :        175 :   avma = av; return zv_sum(E);
    1045                 :            : }
    1046                 :            : 
    1047                 :            : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
    1048                 :            : ulong
    1049                 :     149187 : eulerphiu_fact(GEN f)
    1050                 :            : {
    1051                 :     149187 :   GEN P = gel(f,1), E = gel(f,2);
    1052                 :     149187 :   long i, m = 1, l = lg(P);
    1053         [ +  + ]:     590723 :   for (i = 1; i < l; i++)
    1054                 :            :   {
    1055                 :     441536 :     ulong p = P[i], e = E[i];
    1056         [ -  + ]:     441536 :     if (!e) continue;
    1057         [ +  + ]:     441536 :     if (p == 2)
    1058         [ +  + ]:     112808 :     { if (e > 1) m <<= e-1; }
    1059                 :            :     else
    1060                 :            :     {
    1061                 :     328728 :       m *= (p-1);
    1062         [ +  + ]:     328728 :       if (e > 1) m *= upowuu(p, e-1);
    1063                 :            :     }
    1064                 :            :   }
    1065                 :     149187 :   return m;
    1066                 :            : }
    1067                 :            : ulong
    1068                 :     148816 : eulerphiu(ulong n)
    1069                 :            : {
    1070                 :     148816 :   pari_sp av = avma;
    1071                 :            :   GEN F;
    1072         [ -  + ]:     148816 :   if (!n) return 2;
    1073                 :     148816 :   F = factoru(n);
    1074                 :     148816 :   avma = av; return eulerphiu_fact(F);
    1075                 :            : }
    1076                 :            : GEN
    1077                 :     144193 : eulerphi(GEN n)
    1078                 :            : {
    1079                 :     144193 :   pari_sp av = avma;
    1080                 :            :   GEN Q, F, P, E;
    1081                 :            :   long i, l;
    1082                 :            : 
    1083         [ +  + ]:     144193 :   if ((F = check_arith_all(n,"eulerphi")))
    1084                 :            :   {
    1085                 :         28 :     F = clean_Z_factor(F);
    1086         [ -  + ]:         28 :     n = (typ(n) == t_VEC)? gel(n,1): factorback(F);
    1087         [ +  + ]:         28 :     if (lgefint(n) == 3)
    1088                 :            :     {
    1089                 :            :       ulong e;
    1090                 :         21 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
    1091                 :         21 :       e = eulerphiu_fact(F);
    1092                 :         21 :       avma = av; return utoipos(e);
    1093                 :            :     }
    1094                 :            :   }
    1095         [ +  + ]:     144165 :   else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
    1096                 :            :   else
    1097                 :         46 :     F = absi_factor(n);
    1098         [ +  + ]:         53 :   if (!signe(n)) return gen_2;
    1099                 :         39 :   P = gel(F,1);
    1100                 :         39 :   E = gel(F,2); l = lg(P);
    1101                 :         39 :   Q = cgetg(l, t_VEC);
    1102         [ +  + ]:        252 :   for (i = 1; i < l; i++)
    1103                 :            :   {
    1104                 :        213 :     GEN p = gel(P,i), q;
    1105                 :        213 :     ulong v = itou(gel(E,i));
    1106                 :        213 :     q = subiu(p,1);
    1107 [ +  + ][ +  + ]:        213 :     if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
    1108                 :        213 :     gel(Q,i) = q;
    1109                 :            :   }
    1110                 :     144193 :   return gerepileuptoint(av, ZV_prod(Q));
    1111                 :            : }
    1112                 :            : 
    1113                 :            : static GEN
    1114                 :         60 : numdiv_aux(GEN F)
    1115                 :            : {
    1116                 :         60 :   GEN x, E = gel(F,2);
    1117                 :         60 :   long i, l = lg(E);
    1118                 :         60 :   x = cgetg(l, t_VECSMALL);
    1119         [ +  + ]:        301 :   for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
    1120                 :         60 :   return x;
    1121                 :            : }
    1122                 :            : GEN
    1123                 :        189 : numdiv(GEN n)
    1124                 :            : {
    1125                 :        189 :   pari_sp av = avma;
    1126                 :            :   GEN F, E;
    1127                 :            :   long i, l;
    1128         [ +  + ]:        189 :   if ((F = check_arith_non0(n,"numdiv")))
    1129                 :            :   {
    1130                 :         21 :     F = clean_Z_factor(F);
    1131                 :         21 :     E = numdiv_aux(F);
    1132                 :            :   }
    1133         [ +  + ]:        154 :   else if (lgefint(n) == 3)
    1134                 :            :   {
    1135         [ -  + ]:        115 :     if (n[2] == 1) return gen_1;
    1136                 :        115 :     F = factoru(n[2]);
    1137                 :        115 :     E = gel(F,2); l = lg(E);
    1138         [ +  + ]:        658 :     for (i=1; i<l; i++) E[i]++;
    1139                 :            :   }
    1140                 :            :   else
    1141                 :         39 :     E = numdiv_aux(absi_factor(n));
    1142                 :        175 :   return gerepileuptoint(av, zv_prod_Z(E));
    1143                 :            : }
    1144                 :            : 
    1145                 :            : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
    1146                 :            : static GEN
    1147                 :        557 : u_euler_sumdiv(ulong p, long v)
    1148                 :            : {
    1149                 :        557 :   GEN u = utoipos(1 + p); /* can't overflow */
    1150         [ +  + ]:        571 :   for (; v > 1; v--) u = addsi(1, mului(p, u));
    1151                 :        557 :   return u;
    1152                 :            : }
    1153                 :            : /* 1 + q + ... + q^v */
    1154                 :            : static GEN
    1155                 :       1970 : euler_sumdiv(GEN q, long v)
    1156                 :            : {
    1157                 :       1970 :   GEN u = addui(1, q);
    1158         [ +  + ]:       2362 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
    1159                 :       1970 :   return u;
    1160                 :            : }
    1161                 :            : static GEN
    1162                 :       1289 : u_euler_sumdivk(ulong p, long v, long k) { return euler_sumdiv(powuu(p,k), v); }
    1163                 :            : static GEN
    1164                 :        454 : euler_sumdivk(GEN p, long v, long k) { return euler_sumdiv(powiu(p,k), v); }
    1165                 :            : 
    1166                 :            : static GEN
    1167                 :         53 : sumdiv_aux(GEN F)
    1168                 :            : {
    1169                 :         53 :   GEN x, P = gel(F,1), E = gel(F,2);
    1170                 :         53 :   long i, l = lg(P);
    1171                 :         53 :   x = cgetg(l, t_VEC);
    1172         [ +  + ]:        280 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
    1173                 :         53 :   return x;
    1174                 :            : }
    1175                 :            : GEN
    1176                 :        189 : sumdiv(GEN n)
    1177                 :            : {
    1178                 :        189 :   pari_sp av = avma;
    1179                 :            :   GEN F, P, E;
    1180                 :            :   long i, l;
    1181                 :            : 
    1182         [ +  + ]:        189 :   if ((F = check_arith_non0(n,"sumdiv")))
    1183                 :            :   {
    1184                 :         21 :     F = clean_Z_factor(F);
    1185                 :         21 :     P = sumdiv_aux(F);
    1186                 :            :   }
    1187         [ +  + ]:        154 :   else if (lgefint(n) == 3)
    1188                 :            :   {
    1189         [ -  + ]:        122 :     if (n[2] == 1) return gen_1;
    1190                 :        122 :     F = factoru(n[2]);
    1191                 :        122 :     P = gel(F,1);
    1192                 :        122 :     E = gel(F,2); l = lg(P);
    1193         [ +  + ]:        679 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdiv(P[i], E[i]);
    1194                 :            :   }
    1195                 :            :   else
    1196                 :         32 :     P = sumdiv_aux(absi_factor(n));
    1197                 :        175 :   return gerepileuptoint(av, ZV_prod(P));
    1198                 :            : }
    1199                 :            : 
    1200                 :            : static GEN
    1201                 :        106 : sumdivk_aux(GEN F, long k)
    1202                 :            : {
    1203                 :        106 :   GEN x, P = gel(F,1), E = gel(F,2);
    1204                 :        106 :   long i, l = lg(P);
    1205                 :        106 :   x = cgetg(l, t_VEC);
    1206         [ +  + ]:        560 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdivk(gel(P,i), gel(E,i)[2], k);
    1207                 :        106 :   return x;
    1208                 :            : }
    1209                 :            : GEN
    1210                 :        672 : sumdivk(GEN n, long k)
    1211                 :            : {
    1212                 :        672 :   pari_sp av = avma;
    1213                 :            :   GEN E, F, P;
    1214                 :            :   long i, l, k1;
    1215                 :            : 
    1216         [ -  + ]:        672 :   if (!k) return numdiv(n);
    1217         [ +  + ]:        672 :   if (k == 1) return sumdiv(n);
    1218         [ -  + ]:        483 :   if (k ==-1) return gerepileupto(av, gdiv(sumdiv(n), n));
    1219                 :        483 :   k1 = k;
    1220         [ +  + ]:        483 :   if (k < 0)  k = -k;
    1221         [ +  + ]:        483 :   if ((F = check_arith_non0(n,"sumdivk")))
    1222                 :            :   {
    1223                 :         42 :     F = clean_Z_factor(F);
    1224                 :         42 :     P = sumdivk_aux(F,k);
    1225                 :            :   }
    1226         [ +  + ]:        413 :   else if (lgefint(n) == 3)
    1227                 :            :   {
    1228         [ +  + ]:        349 :     if (n[2] == 1) return gen_1;
    1229                 :        342 :     F = factoru(n[2]);
    1230                 :        342 :     P = gel(F,1);
    1231                 :        342 :     E = gel(F,2); l = lg(P);
    1232         [ +  + ]:       1589 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdivk(P[i], E[i], k);
    1233                 :            :   }
    1234                 :            :   else
    1235                 :         64 :     P = sumdivk_aux(absi_factor(n), k);
    1236                 :        448 :   P = ZV_prod(P);
    1237         [ +  + ]:        448 :   if (k1 > 0) return gerepileuptoint(av, P);
    1238                 :        630 :   return gerepileupto(av, gdiv(P, powiu(n,k)));
    1239                 :            : }
    1240                 :            : 
    1241                 :            : /* K t_VECSMALL of k >= 0 */
    1242                 :            : GEN
    1243                 :         28 : usumdivkvec(ulong n, GEN K)
    1244                 :            : {
    1245                 :         28 :   pari_sp av = avma;
    1246                 :         28 :   GEN F = factoru(n), P = gel(F,1), E = gel(F,2), Z, S;
    1247                 :         28 :   long i,j, l = lg(P), lK = lg(K);
    1248                 :         28 :   Z = cgetg(l, t_VEC);
    1249                 :         28 :   S = cgetg(lK, t_VEC);
    1250         [ +  + ]:         84 :   for (j=1; j<lK; j++)
    1251                 :            :   {
    1252                 :         56 :     long k = K[j];
    1253         [ +  + ]:         98 :     for (i=1; i<l; i++) gel(Z,i) = u_euler_sumdivk(P[i], E[i], k);
    1254                 :         56 :     gel(S,j) = ZV_prod(Z);
    1255                 :            :   }
    1256                 :         28 :   return gerepilecopy(av, S);
    1257                 :            : }
    1258                 :            : 
    1259                 :            : long
    1260                 :        168 : uissquarefree_fact(GEN f)
    1261                 :            : {
    1262                 :        168 :   GEN E = gel(f,2);
    1263                 :        168 :   long i, l = lg(E);
    1264         [ +  + ]:        308 :   for (i = 1; i < l; i++)
    1265         [ +  + ]:        238 :     if (E[i] > 1) return 0;
    1266                 :        168 :   return 1;
    1267                 :            : }
    1268                 :            : long
    1269                 :       7978 : uissquarefree(ulong n)
    1270                 :            : {
    1271         [ -  + ]:       7978 :   if (!n) return 0;
    1272                 :       7978 :   return moebiusu(n)? 1: 0;
    1273                 :            : }
    1274                 :            : long
    1275                 :        764 : Z_issquarefree(GEN n)
    1276                 :            : {
    1277      [ +  +  + ]:        764 :   switch(lgefint(n))
    1278                 :            :   {
    1279                 :          7 :     case 2: return 0;
    1280                 :          6 :     case 3: return uissquarefree(n[2]);
    1281                 :            :   }
    1282                 :        764 :   return moebius(n)? 1: 0;
    1283                 :            : }
    1284                 :            : long
    1285                 :         49 : issquarefree(GEN x)
    1286                 :            : {
    1287                 :            :   pari_sp av;
    1288                 :            :   GEN d;
    1289      [ +  +  - ]:         49 :   switch(typ(x))
    1290                 :            :   {
    1291                 :         14 :     case t_INT: return Z_issquarefree(x);
    1292                 :            :     case t_POL:
    1293         [ -  + ]:         35 :       if (!signe(x)) return 0;
    1294                 :         35 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
    1295                 :         35 :       avma = av; return (lg(d) == 3);
    1296                 :          0 :     default: pari_err_TYPE("issquarefree",x);
    1297                 :         49 :       return 0; /* not reached */
    1298                 :            :   }
    1299                 :            : }
    1300                 :            : 
    1301                 :            : /*********************************************************************/
    1302                 :            : /**                                                                 **/
    1303                 :            : /**                    DIGITS / SUM OF DIGITS                       **/
    1304                 :            : /**                                                                 **/
    1305                 :            : /*********************************************************************/
    1306                 :            : 
    1307                 :            : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
    1308                 :            : static void
    1309                 :      33509 : set_vexp(GEN v, long l)
    1310                 :            : {
    1311                 :            :   long m;
    1312         [ +  + ]:      47712 :   if (v[l]) return;
    1313                 :      14203 :   v[l] = 1; m = l>>1;
    1314                 :      14203 :   set_vexp(v, m);
    1315                 :      14203 :   set_vexp(v, l-m);
    1316                 :            : }
    1317                 :            : 
    1318                 :            : /* return all needed B^i for DAC algorithm, for lz digits */
    1319                 :            : static GEN
    1320                 :       5103 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
    1321                 :            : {
    1322                 :       5103 :   GEN vB, vexp = const_vecsmall(lz, 0);
    1323                 :       5103 :   long i, l = (lz+1) >> 1;
    1324                 :       5103 :   vexp[1] = 1;
    1325                 :       5103 :   vexp[2] = 1;
    1326                 :       5103 :   set_vexp(vexp, lz);
    1327                 :       5103 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
    1328                 :       5103 :   gel(vB, 1) = T;
    1329         [ +  + ]:      29057 :   for (i = 2; i <= l; i++)
    1330         [ +  + ]:      23954 :     if (vexp[i])
    1331                 :            :     {
    1332                 :      14203 :       long j = i >> 1;
    1333                 :      14203 :       GEN B2j = r->sqr(E, gel(vB,j));
    1334         [ +  + ]:      14203 :       gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
    1335                 :            :     }
    1336                 :       5103 :   return vB;
    1337                 :            : }
    1338                 :            : 
    1339                 :            : static void
    1340                 :      84299 : gen_digits_dac(GEN x, GEN vB, long l, GEN *z,
    1341                 :            :                void *E, GEN div(void *E, GEN a, GEN b, GEN *r))
    1342                 :            : {
    1343                 :            :   GEN q, r;
    1344                 :      84299 :   long m = l>>1;
    1345         [ +  + ]:     124670 :   if (l==1) { *z=x; return; }
    1346                 :      40371 :   q = div(E, x, gel(vB,m), &r);
    1347                 :      40371 :   gen_digits_dac(r, vB, m, z, E, div);
    1348                 :      40371 :   gen_digits_dac(q, vB, l-m, z+m, E, div);
    1349                 :            : }
    1350                 :            : 
    1351                 :            : static GEN
    1352                 :      18830 : gen_fromdigits_dac(GEN x, GEN vB, long i, long l, void *E,
    1353                 :            :                    GEN add(void *E, GEN a, GEN b),
    1354                 :            :                    GEN mul(void *E, GEN a, GEN b))
    1355                 :            : {
    1356                 :            :   GEN a, b;
    1357                 :      18830 :   long m = l>>1;
    1358         [ +  + ]:      18830 :   if (l==1) return gel(x,i);
    1359                 :       8652 :   a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
    1360                 :       8652 :   b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
    1361                 :      18830 :   return add(E, a, mul(E, b, gel(vB, m)));
    1362                 :            : }
    1363                 :            : 
    1364                 :            : static GEN
    1365                 :       3578 : gen_digits_i(GEN x, GEN B, long n, void *E, struct bb_ring *r,
    1366                 :            :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
    1367                 :            : {
    1368                 :            :   GEN z, vB;
    1369         [ +  + ]:       3578 :   if (n==1) retmkvec(gcopy(x));
    1370                 :       3557 :   vB = get_vB(B, n, E, r);
    1371                 :       3557 :   z = cgetg(n+1, t_VEC);
    1372                 :       3557 :   gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
    1373                 :       3578 :   return z;
    1374                 :            : }
    1375                 :            : 
    1376                 :            : GEN
    1377                 :       3528 : gen_digits(GEN x, GEN B, long n, void *E, struct bb_ring *r,
    1378                 :            :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
    1379                 :            : {
    1380                 :       3528 :   pari_sp av = avma;
    1381                 :       3528 :   return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
    1382                 :            : }
    1383                 :            : 
    1384                 :            : GEN
    1385                 :       1526 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
    1386                 :            : {
    1387                 :       1526 :   pari_sp av = avma;
    1388                 :       1526 :   long n = lg(x)-1;
    1389                 :       1526 :   GEN vB = get_vB(B, n, E, r);
    1390                 :       1526 :   GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
    1391                 :       1526 :   return gerepilecopy(av, z);
    1392                 :            : }
    1393                 :            : 
    1394                 :            : static GEN
    1395                 :       3360 : _addii(void *data /* ignored */, GEN x, GEN y)
    1396                 :       3360 : { (void)data; return addii(x,y); }
    1397                 :            : static GEN
    1398                 :        553 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
    1399                 :            : static GEN
    1400                 :       3619 : _mulii(void *data /* ignored */, GEN x, GEN y)
    1401                 :       3619 :  { (void)data; return mulii(x,y); }
    1402                 :            : static GEN
    1403                 :        436 : _dvmdii(void *data /* ignored */, GEN x, GEN y, GEN *r)
    1404                 :        436 :  { (void)data; return dvmdii(x,y,r); }
    1405                 :            : 
    1406                 :            : static struct bb_ring Z_ring = { _addii, _mulii, _sqri };
    1407                 :            : 
    1408                 :            : static GEN
    1409                 :        168 : check_basis(GEN B)
    1410                 :            : {
    1411         [ +  + ]:        168 :   if (!B) return utoipos(10);
    1412         [ -  + ]:        140 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
    1413         [ -  + ]:        140 :   if (cmpis(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
    1414                 :        168 :   return B;
    1415                 :            : }
    1416                 :            : 
    1417                 :            : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
    1418                 :            : static void
    1419                 :       3054 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
    1420                 :            : {
    1421                 :       3054 :   pari_sp av = avma;
    1422                 :            :   GEN q,r;
    1423                 :            :   long m;
    1424         [ +  + ]:       4571 :   if (l==1) { *z=itou(x); return; }
    1425                 :       1517 :   m=l>>1;
    1426                 :       1517 :   q = dvmdii(x, gel(vB,m), &r);
    1427                 :       1517 :   digits_dacsmall(q,vB,l-m,z);
    1428                 :       1517 :   digits_dacsmall(r,vB,m,z+l-m);
    1429                 :       1517 :   avma = av;
    1430                 :            : }
    1431                 :            : 
    1432                 :            : GEN
    1433                 :         49 : digits(GEN x, GEN B)
    1434                 :            : {
    1435                 :         49 :   pari_sp av=avma;
    1436                 :            :   long lz;
    1437                 :            :   GEN z, vB;
    1438         [ -  + ]:         49 :   if (typ(x)!=t_INT) pari_err_TYPE("digits",x);
    1439                 :         49 :   B = check_basis(B);
    1440         [ +  + ]:         49 :   if (!signe(x))       {avma = av; return cgetg(1,t_VEC); }
    1441         [ -  + ]:         42 :   if (absi_cmp(x,B)<0) {avma = av; retmkvec(absi(x)); }
    1442         [ +  + ]:         42 :   if (Z_ispow2(B))
    1443                 :            :   {
    1444                 :         14 :     long k = expi(B);
    1445         [ +  + ]:         14 :     if (k == 1) return binaire(x);
    1446         [ -  + ]:          7 :     if (k < BITS_IN_LONG)
    1447                 :            :     {
    1448                 :          0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
    1449                 :          0 :       z = binary_2k_zv(x, k);
    1450                 :          0 :       avma = av; return Flv_to_ZV(z);
    1451                 :            :     }
    1452                 :            :     else
    1453                 :            :     {
    1454                 :          7 :       avma = av; return binary_2k(x, k);
    1455                 :            :     }
    1456                 :            :   }
    1457         [ -  + ]:         28 :   if (signe(x) < 0) x = absi(x);
    1458                 :         28 :   lz = logint(x,B,NULL);
    1459         [ +  + ]:         28 :   if (lgefint(B)>3)
    1460                 :            :   {
    1461                 :          8 :     z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
    1462                 :          8 :     vecreverse_inplace(z);
    1463                 :          8 :     return z;
    1464                 :            :   }
    1465                 :            :   else
    1466                 :            :   {
    1467                 :         20 :     vB = get_vB(B, lz, NULL, &Z_ring);
    1468                 :         20 :     (void)new_chunk(3*lz); /* HACK */
    1469                 :         20 :     z = zero_zv(lz);
    1470                 :         20 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1471                 :         49 :     avma = av; return vecsmall_to_vec(z);
    1472                 :            :   }
    1473                 :            : }
    1474                 :            : 
    1475                 :            : GEN
    1476                 :         42 : fromdigits(GEN x, GEN B)
    1477                 :            : {
    1478                 :         42 :   pari_sp av = avma;
    1479 [ +  - ][ -  + ]:         42 :   if (typ(x)!=t_VEC || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
    1480                 :         42 :   B = check_basis(B);
    1481         [ +  + ]:         42 :   if (lg(x)==1) { avma = av; return gen_0; }
    1482                 :         35 :   x = vecreverse(x);
    1483                 :         42 :   return gerepileupto(av, gen_fromdigits(x, B, NULL, &Z_ring));
    1484                 :            : }
    1485                 :            : 
    1486                 :            : static ulong DS[] ={
    1487                 :            :   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,
    1488                 :            :   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,
    1489                 :            :   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,
    1490                 :            :   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,
    1491                 :            :   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,
    1492                 :            :   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,
    1493                 :            :   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,
    1494                 :            :   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,
    1495                 :            :   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,
    1496                 :            :   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,
    1497                 :            :   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,
    1498                 :            :   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,
    1499                 :            :   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,
    1500                 :            :   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,
    1501                 :            :   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,
    1502                 :            :   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,
    1503                 :            :   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,
    1504                 :            :   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,
    1505                 :            :   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,
    1506                 :            :   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,
    1507                 :            :   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,
    1508                 :            :   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,
    1509                 :            :   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,
    1510                 :            :   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,
    1511                 :            :   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,
    1512                 :            :   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,
    1513                 :            :   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,
    1514                 :            :   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,
    1515                 :            :   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,
    1516                 :            :   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,
    1517                 :            :   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,
    1518                 :            :   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,
    1519                 :            :   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,
    1520                 :            :   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,
    1521                 :            :   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,
    1522                 :            :   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,
    1523                 :            :   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,
    1524                 :            :   25,26,27
    1525                 :            : };
    1526                 :            : 
    1527                 :            : ulong
    1528                 :     355152 : sumdigitsu(ulong n)
    1529                 :            : {
    1530                 :     355152 :   ulong s = 0;
    1531         [ +  + ]:    1361843 :   while (n) { s += DS[n % 1000]; n /= 1000; }
    1532                 :     355152 :   return s;
    1533                 :            : }
    1534                 :            : 
    1535                 :            : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1536                 :            : static ulong
    1537                 :         14 : sumdigits_block(ulong *res, long l)
    1538                 :            : {
    1539                 :         14 :   long s = sumdigitsu(*--res);
    1540         [ +  + ]:     355138 :   while (--l > 0) s += sumdigitsu(*--res);
    1541                 :         14 :   return s;
    1542                 :            : }
    1543                 :            : 
    1544                 :            : GEN
    1545                 :         35 : sumdigits(GEN n)
    1546                 :            : {
    1547                 :         35 :   pari_sp av = avma;
    1548                 :            :   ulong s, *res;
    1549                 :            :   long l;
    1550                 :            : 
    1551         [ -  + ]:         35 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1552                 :         35 :   l = lgefint(n);
    1553      [ +  +  + ]:         35 :   switch(l)
    1554                 :            :   {
    1555                 :          7 :     case 2: return gen_0;
    1556                 :         14 :     case 3: return utoipos(sumdigitsu(n[2]));
    1557                 :            :   }
    1558                 :         14 :   res = convi(n, &l);
    1559         [ +  - ]:         14 :   if ((ulong)l < ULONG_MAX / 81)
    1560                 :            :   {
    1561                 :         14 :     s = sumdigits_block(res, l);
    1562                 :         14 :     avma = av; return utoipos(s);
    1563                 :            :   }
    1564                 :            :   else /* Huge. Overflows ulong */
    1565                 :            :   {
    1566                 :          0 :     const long L = (long)(ULONG_MAX / 81);
    1567                 :          0 :     GEN S = gen_0;
    1568         [ #  # ]:          0 :     while (l > L)
    1569                 :            :     {
    1570                 :          0 :       S = addiu(S, sumdigits_block(res, L));
    1571                 :          0 :       res += L; l -= L;
    1572                 :            :     }
    1573         [ #  # ]:          0 :     if (l)
    1574                 :          0 :       S = addiu(S, sumdigits_block(res, l));
    1575                 :         35 :     return gerepileuptoint(av, S);
    1576                 :            :   }
    1577                 :            : }
    1578                 :            : 
    1579                 :            : GEN
    1580                 :        105 : sumdigits0(GEN x, GEN B)
    1581                 :            : {
    1582                 :        105 :   pari_sp av = avma;
    1583                 :            :   GEN z;
    1584                 :            :   long lz;
    1585                 :            : 
    1586         [ +  + ]:        105 :   if (!B) return sumdigits(x);
    1587         [ -  + ]:         77 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1588                 :         77 :   B = check_basis(B);
    1589         [ +  + ]:         77 :   if (Z_ispow2(B))
    1590                 :            :   {
    1591                 :         28 :     long k = expi(B);
    1592         [ +  + ]:         28 :     if (k == 1) { avma = av; return utoi(hammingweight(x)); }
    1593         [ +  + ]:         21 :     if (k < BITS_IN_LONG)
    1594                 :            :     {
    1595                 :         14 :       GEN z = binary_2k_zv(x, k);
    1596         [ -  + ]:         14 :       if (lg(z)-1 > 1<<(BITS_IN_LONG-k)) /* may overflow */
    1597                 :          0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1598                 :         14 :       avma = av; return utoi(zv_sum(z));
    1599                 :            :     }
    1600                 :          7 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1601                 :            :   }
    1602         [ -  + ]:         49 :   if (!signe(x))       { avma = av; return gen_0; }
    1603         [ -  + ]:         49 :   if (absi_cmp(x,B)<0) { avma = av; return absi(x); }
    1604         [ +  + ]:         49 :   if (equaliu(B,10))   { avma = av; return sumdigits(x); }
    1605                 :         42 :   lz = logint(x,B,NULL);
    1606                 :         42 :   z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
    1607                 :        105 :   return gerepileuptoint(av, ZV_sum(z));
    1608                 :            : }

Generated by: LCOV version 1.9