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 17110-9967e23) Lines: 615 727 84.6 %
Date: 2014-11-26 Functions: 78 87 89.7 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 338 487 69.4 %

           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                 :        948 : initprimes1(ulong size, long *lenp, ulong *lastp, byteptr p1)
      32                 :            : {
      33                 :        948 :   pari_sp av = avma;
      34                 :            :   long k;
      35                 :        948 :   byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
      36                 :            : 
      37         [ +  + ]:       9480 :   for (r=q=p,k=1; r<=fin; )
      38                 :            :   {
      39         [ +  + ]:      13272 :     do { r+=k; k+=2; r+=k; } while (*++q);
      40         [ +  + ]:     301464 :     for (s=r; s<=fin; s+=k) *s = 1;
      41                 :            :   }
      42                 :        948 :   r = p1; *r++ = 2; *r++ = 1; /* 2 and 3 */
      43                 :        948 :   for (s=q=p+1; ; s=q)
      44                 :            :   {
      45         [ +  + ]:     334644 :     do q++; while (*q);
      46         [ +  + ]:     118500 :     if (q > fin) break;
      47                 :     117552 :     *r++ = (unsigned char) ((q-s) << 1);
      48                 :     117552 :   }
      49                 :        948 :   *r++ = 0;
      50                 :        948 :   *lenp = r - p1;
      51                 :        948 :   *lastp = ((s - p) << 1) + 1;
      52                 :        948 :   avma = av;
      53                 :        948 : }
      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                 :        948 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
     198                 :            :                 cache_model_t *cache_model)
     199                 :            : {
     200                 :        948 :   ulong asize, cache_arena = cache_model->arena;
     201                 :            :   double Xmin, Xmax, A, B, C1, C2, D, V;
     202                 :        948 :   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         [ -  + ]:        948 :   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         [ +  - ]:        948 :   if (cache_arena - fixed_to_cache > 10 * slow2_size) {
     237                 :        948 :       asize = cache_arena - fixed_to_cache;
     238         [ -  + ]:        948 :       if (asize > total) asize = total; /* Automatically false... */
     239                 :        948 :       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                 :        948 :   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                 :       2240 : sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulong n)
     342                 :            : {
     343                 :       2240 :   ulong p, cnt = n-1, start = s, delta = 1;
     344                 :            :   byteptr q;
     345                 :            : 
     346                 :       2240 :   memset(data, 0, n);
     347                 :       2240 :   start >>= 1;  /* (start - 1)/2 */
     348                 :       2240 :   start += n; /* Corresponds to the end */
     349                 :            :   /* data corresponds to start, q runs over primediffs */
     350         [ +  + ]:     259432 :   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                 :     257192 :     long off = cnt - ((start+(p>>1)) % p);
     357         [ +  + ]:  386191296 :     while (off >= 0) { data[off] = 1; off -= p; }
     358                 :            :   }
     359                 :       2240 : }
     360                 :            : 
     361                 :            : /* assume maxnum <= 436273289 < 2^29 */
     362                 :            : static void
     363                 :        948 : initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
     364                 :            : {
     365                 :        948 :   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                 :        948 :   maxnum |= 1; /* make it odd. */
     373                 :            :   /* base case */
     374         [ -  + ]:       1896 :   if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, lastp, p1); return; }
     375                 :            : 
     376                 :            :   /* Checked to be enough up to 40e6, attained at 155893 */
     377                 :        948 :   rootnum = usqrt(maxnum) | 1;
     378                 :        948 :   initprimes1(rootnum>>1, &psize, &last, p1);
     379                 :        948 :   end1 = p1 + psize - 1;
     380                 :        948 :   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                 :        948 :   asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
     385                 :            :                           &cache_model) - 1;
     386                 :            :   /* enough room on the stack ? */
     387                 :        948 :   alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
     388         [ -  + ]:        948 :   if (alloced)
     389                 :          0 :     p = (byteptr)pari_malloc(asize+1);
     390                 :            :   else
     391                 :        948 :     p = (byteptr)stack_malloc(asize+1);
     392                 :        948 :   end = p + asize; /* the 0 sentinel goes at end. */
     393                 :        948 :   curlow = last + 2; /* First candidate: know primes up to last (odd). */
     394                 :        948 :   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                 :        948 :   plast = p - 1;
     400                 :        948 :   p_prime_above = p1 + 2;
     401                 :        948 :   prime_above = 3;
     402         [ +  + ]:       3188 :   while (remains)
     403                 :            :   { /* cycle over arenas; performance not crucial */
     404                 :            :     unsigned char was_delta;
     405         [ +  + ]:       2240 :     if (asize > remains) { asize = remains; end = p + asize; }
     406                 :            :     /* Fake the upper limit appropriate for the given arena */
     407 [ +  + ][ +  - ]:     120740 :     while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
     408                 :     118500 :       prime_above += *p_prime_above++;
     409                 :       2240 :     was_delta = *p_prime_above;
     410                 :       2240 :     *p_prime_above = 0; /* sentinel for sieve_chunk */
     411                 :       2240 :     sieve_chunk(p1, curlow, p, asize);
     412                 :       2240 :     *p_prime_above = was_delta; /* restore */
     413                 :            : 
     414                 :       2240 :     p[asize] = 0; /* sentinel */
     415                 :       2240 :     for (q = p; ; plast = q++)
     416                 :            :     { /* q runs over addresses corresponding to primes */
     417         [ +  + ]:  236670440 :       while (*q) q++; /* use sentinel at end */
     418         [ +  + ]:   39260816 :       if (q >= end) break;
     419                 :   39258576 :       *curdiff++ = (unsigned char)(q-plast) << 1; /* < 255 for q < 436273291 */
     420                 :   39258576 :     }
     421                 :       2240 :     plast -= asize;
     422                 :       2240 :     remains -= asize;
     423                 :       2240 :     curlow += (asize<<1);
     424                 :            :   }
     425                 :        948 :   last = curlow - ((p - plast) << 1);
     426                 :        948 :   *curdiff++ = 0; /* sentinel */
     427                 :        948 :   *lenp = curdiff - p1;
     428                 :        948 :   *lastp = last;
     429         [ -  + ]:        948 :   if (alloced) pari_free(p); else avma = av;
     430                 :            : }
     431                 :            : 
     432                 :            : ulong
     433         [ +  - ]:    9093973 : 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                 :        948 : initprimes(ulong maxnum, long *lenp, ulong *lastp)
     444                 :            : {
     445                 :            :   byteptr t;
     446         [ -  + ]:        948 :   if (maxnum < 65537)
     447                 :          0 :     maxnum = 65537;
     448         [ -  + ]:        948 :   else if (maxnum > 436273289)
     449                 :          0 :     maxnum = 436273289;
     450                 :        948 :   t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
     451                 :        948 :   initprimes0(maxnum, lenp, lastp, t);
     452                 :        948 :   return (byteptr)pari_realloc(t, *lenp);
     453                 :            : }
     454                 :            : 
     455                 :            : void
     456                 :        948 : initprimetable(ulong maxnum)
     457                 :            : {
     458                 :            :   long len;
     459                 :            :   ulong last;
     460                 :        948 :   byteptr p = initprimes(maxnum, &len, &last), old = diffptr;
     461                 :        948 :   diffptrlen = minss(diffptrlen, len);
     462                 :        948 :   _maxprime  = minss(_maxprime,last); /*Protect against ^C*/
     463                 :        948 :   diffptr = p; diffptrlen = len; _maxprime = last;
     464         [ -  + ]:        948 :   if (old) free(old);
     465                 :        948 : }
     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                 :    4183558 : init_primepointer_lt(ulong a, byteptr *pd)
     480                 :            : {
     481                 :            :   ulong n, p;
     482                 :    4183558 :   prime_table_next_p(a, pd, &p, &n);
     483                 :    4183558 :   PREC_PRIME_VIADIFF(p, *pd);
     484                 :    4183558 :   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                 :         25 : boundfact(GEN n, ulong lim)
     507                 :            : {
     508      [ +  +  - ]:         25 :   switch(typ(n))
     509                 :            :   {
     510                 :         15 :     case t_INT: return Z_factor_limit(n,lim);
     511                 :            :     case t_FRAC: {
     512                 :         10 :       pari_sp av = avma;
     513                 :         10 :       GEN a = Z_factor_limit(gel(n,1),lim);
     514                 :         10 :       GEN b = Z_factor_limit(gel(n,2),lim);
     515                 :         10 :       gel(b,2) = ZC_neg(gel(b,2));
     516                 :         10 :       return gerepilecopy(av, merge_factor_i(a,b));
     517                 :            :     }
     518                 :            :   }
     519                 :          0 :   pari_err_TYPE("boundfact",n);
     520                 :         25 :   return NULL; /* not reached */
     521                 :            : }
     522                 :            : 
     523                 :            : /* NOT memory clean */
     524                 :            : GEN
     525                 :       9330 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
     526                 :            : {
     527                 :       9330 :   long i, j, l = lg(L);
     528                 :       9330 :   GEN e = new_chunk(l), P = new_chunk(l);
     529         [ +  + ]:      68947 :   for (i = j = 1; i < l; i++)
     530                 :            :   {
     531                 :      65888 :     ulong p = uel(L,i);
     532                 :      65888 :     long v = Z_lvalrem(N, p, &N);
     533 [ +  + ][ +  + ]:      65888 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
     534                 :            :   }
     535                 :       9330 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
     536                 :       9330 :   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                 :      15945 : factoru_pow(ulong n)
     547                 :            : {
     548                 :      15945 :   GEN f = cgetg(4,t_VEC);
     549                 :      15945 :   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                 :      15945 :   (void)new_chunk((15 + 1)*3);
     554                 :      15945 :   F = factoru(n);
     555                 :      15945 :   P = gel(F,1);
     556                 :      15945 :   E = gel(F,2); l = lg(P);
     557                 :      15945 :   avma = av;
     558                 :      15945 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
     559                 :      15945 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
     560                 :      15945 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
     561         [ +  + ]:      51285 :   for(i = 1; i < l; i++)
     562                 :            :   {
     563                 :      35340 :     p[i] = P[i];
     564                 :      35340 :     e[i] = E[i];
     565                 :      35340 :     c[i] = upowuu(p[i], e[i]);
     566                 :            :   }
     567                 :      15945 :   return f;
     568                 :            : }
     569                 :            : 
     570                 :            : static GEN
     571                 :      16950 : factorlim(GEN n, ulong lim)
     572         [ -  + ]:      16950 : { 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                 :       4730 : factor_pn_1_limit(GEN p, long n, ulong lim)
     577                 :            : {
     578                 :       4730 :   pari_sp av = avma;
     579                 :       4730 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
     580                 :       4730 :   long i, pp = itos_or_0(p);
     581         [ +  + ]:      13260 :   for(i=2; i<lg(d); i++)
     582                 :            :   {
     583                 :            :     GEN B;
     584 [ +  + ][ +  + ]:       8530 :     if (pp && d[i]%pp==0 && (
                 [ -  + ]
     585 [ #  # ][ +  + ]:       7655 :        ((pp&3)==1 && (d[i]&1)) ||
     586 [ +  + ][ +  + ]:       7655 :        ((pp&3)==3 && (d[i]&3)==2) ||
     587         [ +  + ]:       7600 :        (pp==2 && (d[i]&7)==4)))
     588                 :       3690 :     {
     589                 :       3690 :       GEN f=factor_Aurifeuille_prime(p,d[i]);
     590                 :       3690 :       B = factorlim(f, lim);
     591                 :       3690 :       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     592                 :       3690 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     593                 :            :     }
     594                 :            :     else
     595                 :       4840 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     596                 :       8530 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     597                 :            :   }
     598                 :       4730 :   return gerepilecopy(av, A);
     599                 :            : }
     600                 :            : GEN
     601                 :       4730 : factor_pn_1(GEN p, ulong n)
     602                 :       4730 : { 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                 :     222550 : RgV_is_ZVpos(GEN v)
     668                 :            : {
     669                 :     222550 :   long i, l = lg(v);
     670         [ +  + ]:     578677 :   for (i = 1; i < l; i++)
     671                 :            :   {
     672                 :     356132 :     GEN c = gel(v,i);
     673 [ +  - ][ +  + ]:     356132 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     674                 :            :   }
     675                 :     222550 :   return 1;
     676                 :            : }
     677                 :            : /* check whether v is a ZV with non-0 entries */
     678                 :            : static int
     679                 :        185 : RgV_is_ZVnon0(GEN v)
     680                 :            : {
     681                 :        185 :   long i, l = lg(v);
     682         [ +  + ]:        620 :   for (i = 1; i < l; i++)
     683                 :            :   {
     684                 :        470 :     GEN c = gel(v,i);
     685 [ +  - ][ +  + ]:        470 :     if (typ(c) != t_INT || !signe(c)) return 0;
     686                 :            :   }
     687                 :        185 :   return 1;
     688                 :            : }
     689                 :            : /* check whether v is a ZV with non-zero entries OR exactly [0] */
     690                 :            : static int
     691                 :        150 : RgV_is_ZV0(GEN v)
     692                 :            : {
     693                 :        150 :   long i, l = lg(v);
     694         [ +  + ]:        415 :   for (i = 1; i < l; i++)
     695                 :            :   {
     696                 :        305 :     GEN c = gel(v,i);
     697                 :            :     long s;
     698         [ -  + ]:        305 :     if (typ(c) != t_INT) return 0;
     699                 :        305 :     s = signe(c);
     700         [ +  + ]:        305 :     if (!s) return (l == 2);
     701                 :            :   }
     702                 :        150 :   return 1;
     703                 :            : }
     704                 :            : 
     705                 :            : static int
     706                 :     111450 : is_Z_factor_i(GEN f)
     707 [ +  - ][ +  + ]:     111450 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
                 [ +  + ]
     708                 :            : int
     709                 :     111115 : is_Z_factorpos(GEN f)
     710 [ +  + ][ +  - ]:     111115 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     711                 :            : int
     712                 :        150 : is_Z_factor(GEN f)
     713 [ +  - ][ +  - ]:        150 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     714                 :            : /* as is_Z_factorpos, also allow factor(0) */
     715                 :            : int
     716                 :        185 : is_Z_factornon0(GEN f)
     717 [ +  - ][ +  + ]:        185 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     718                 :            : GEN
     719                 :        120 : clean_Z_factor(GEN f)
     720                 :            : {
     721                 :        120 :   GEN P = gel(f,1);
     722                 :        120 :   long n = lg(P)-1;
     723 [ +  + ][ +  + ]:        120 :   if (n && equalim1(gel(P,1)))
     724                 :         40 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     725                 :        120 :   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                 :       1996 : check_arith_non0(GEN n, const char *f) {
     766   [ +  +  +  - ]:       1996 :   switch(typ(n))
     767                 :            :   {
     768                 :            :     case t_INT:
     769         [ +  + ]:       1846 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     770                 :       1806 :       return NULL;
     771                 :            :     case t_VEC:
     772 [ +  - ][ +  - ]:          5 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     773                 :          5 :       n = gel(n,2); /* fall through */
     774                 :            :     case t_MAT:
     775         [ +  + ]:        150 :       if (!is_Z_factornon0(n)) break;
     776                 :        115 :       return n;
     777                 :            :   }
     778                 :         35 :   pari_err_TYPE(f,n);
     779                 :       1921 :   return NULL;
     780                 :            : }
     781                 :            : /* n associated to a factorization of an integer */
     782                 :            : GEN
     783                 :     121035 : check_arith_all(GEN n, const char *f) {
     784   [ +  +  +  - ]:     121035 :   switch(typ(n))
     785                 :            :   {
     786                 :            :     case t_INT:
     787                 :     120885 :       return NULL;
     788                 :            :     case t_VEC:
     789 [ +  - ][ +  - ]:          5 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     790                 :          5 :       n = gel(n,2); /* fall through */
     791                 :            :     case t_MAT:
     792         [ -  + ]:        150 :       if (!is_Z_factor(n)) break;
     793                 :        150 :       return n;
     794                 :            :   }
     795                 :          0 :   pari_err_TYPE(f,n);
     796                 :     121035 :   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                 :     250535 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     808                 :            : {
     809                 :            :   GEN E, P;
     810         [ -  + ]:     250535 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     811                 :     250535 :   P = gel(F,1);
     812                 :     250535 :   E = gel(F,2);
     813                 :     250535 :   RgV_check_ZV(E, "divisors");
     814                 :     250535 :   *isint = RgV_is_ZV(P);
     815         [ +  + ]:     250535 :   if (*isint)
     816                 :            :   {
     817                 :     250525 :     long i, l = lg(P);
     818                 :            :     /* skip -1 */
     819 [ +  + ][ -  + ]:     250525 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     820                 :            :     /* test for 0 */
     821         [ +  + ]:     816525 :     for (i = 1; i < l; i++)
     822 [ +  + ][ +  - ]:     566005 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     823                 :          5 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     824                 :            :   }
     825                 :     250530 :   *pP = P;
     826                 :     250530 :   *pE = E;
     827                 :     250530 : }
     828                 :            : static void
     829                 :       7580 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     830                 :            : 
     831                 :            : int
     832                 :     258115 : divisors_init(GEN n, GEN *pP, GEN *pE)
     833                 :            : {
     834                 :            :   long i,l;
     835                 :            :   GEN E, P, e;
     836                 :            :   int isint;
     837                 :            : 
     838   [ +  +  +  + ]:     258115 :   switch(typ(n))
     839                 :            :   {
     840                 :            :     case t_INT:
     841         [ -  + ]:       7565 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     842                 :       7565 :       set_fact(absi_factor(n), &P,&E);
     843                 :       7565 :       isint = 1; break;
     844                 :            :     case t_VEC:
     845         [ -  + ]:          5 :       if (lg(n) != 3) pari_err_TYPE("divisors",n);
     846                 :          5 :       set_fact_check(gel(n,2), &P,&E, &isint);
     847                 :          5 :       break;
     848                 :            :     case t_MAT:
     849                 :     250530 :       set_fact_check(n, &P,&E, &isint);
     850                 :     250525 :       break;
     851                 :            :     default:
     852                 :         15 :       set_fact(factor(n), &P,&E);
     853                 :         15 :       isint = 0; break;
     854                 :            :   }
     855                 :     258110 :   l = lg(P);
     856                 :     258110 :   e = cgetg(l, t_VECSMALL);
     857         [ +  + ]:     845895 :   for (i=1; i<l; i++)
     858                 :            :   {
     859                 :     587790 :     e[i] = itos(gel(E,i));
     860         [ +  + ]:     587790 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     861                 :            :   }
     862                 :     258105 :   *pP = P; *pE = e; return isint;
     863                 :            : }
     864                 :            : 
     865                 :            : GEN
     866                 :     258085 : divisors(GEN n)
     867                 :            : {
     868                 :     258085 :   pari_sp av = avma;
     869                 :            :   long i, j, l;
     870                 :            :   ulong ndiv;
     871                 :            :   GEN *d, *t, *t1, *t2, *t3, P, E, e;
     872                 :     258085 :   int isint = divisors_init(n, &P, &E);
     873                 :            : 
     874                 :     258075 :   l = lg(E); e = cgetg(l, t_VECSMALL);
     875         [ +  + ]:     845750 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     876                 :     258075 :   ndiv = itou_or_0( zv_prod_Z(e) );
     877 [ +  - ][ -  + ]:     258075 :   if (!ndiv || ndiv & ~LGBITS) pari_err_OVERFLOW("divisors");
     878                 :     258075 :   d = t = (GEN*) cgetg(ndiv+1,t_VEC);
     879                 :     258075 :   *++d = gen_1;
     880         [ +  + ]:     258075 :   if (isint)
     881                 :            :   {
     882         [ +  + ]:     845710 :     for (i=1; i<l; i++)
     883         [ +  + ]:    1578520 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     884         [ +  + ]:    3825425 :         for (t2=d, t3=t1; t3<t2; ) *++d = mulii(*++t3, gel(P,i));
     885                 :     258065 :     e = ZV_sort((GEN)t);
     886                 :            :   } else {
     887         [ +  + ]:         40 :     for (i=1; i<l; i++)
     888         [ +  + ]:        130 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     889         [ +  + ]:        810 :         for (t2=d, t3=t1; t3<t2; ) *++d = gmul(*++t3, gel(P,i));
     890                 :         10 :     e = (GEN)t;
     891                 :            :   }
     892                 :     258075 :   return gerepileupto(av, e);
     893                 :            : }
     894                 :            : 
     895                 :            : GEN
     896                 :      21511 : divisorsu(ulong n)
     897                 :            : {
     898                 :      21511 :   pari_sp av = avma;
     899                 :            :   long i, j, l;
     900                 :            :   ulong nbdiv;
     901                 :      21511 :   GEN d, t, t1, t2, t3, P, e, fa = factoru(n);
     902                 :            : 
     903                 :      21511 :   P = gel(fa,1); l = lg(P);
     904                 :      21511 :   e = gel(fa,2);
     905                 :      21511 :   nbdiv = 1;
     906         [ +  + ]:      47317 :   for (i=1; i<l; i++) nbdiv *= 1+e[i];
     907                 :      21511 :   d = t = cgetg(nbdiv+1,t_VECSMALL);
     908                 :      21511 :   *++d = 1;
     909         [ +  + ]:      47317 :   for (i=1; i<l; i++)
     910         [ +  + ]:      58421 :     for (t1=t,j=e[i]; j; j--,t1=t2)
     911         [ +  + ]:      83557 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     912                 :      21511 :   vecsmall_sort(t);
     913                 :      21511 :   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                 :        110 : core2fa(GEN fa)
     927                 :            : {
     928                 :        110 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
     929                 :            :   long i;
     930         [ +  + ]:        295 :   for (i=1; i<lg(P); i++)
     931                 :            :   {
     932                 :        185 :     long e = itos(gel(E,i));
     933         [ +  + ]:        185 :     if (e & 1)  c = mulii(c, gel(P,i));
     934         [ +  + ]:        185 :     if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
     935                 :            :   }
     936                 :        110 :   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                 :        110 : core2partial(GEN n, long all)
     947                 :            : {
     948                 :        110 :   pari_sp av = avma;
     949         [ -  + ]:        110 :   if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
     950                 :        110 :   return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
     951                 :            : }
     952                 :            : static GEN
     953                 :         50 : core2_i(GEN n)
     954                 :            : {
     955                 :         50 :   GEN f = core(n);
     956         [ +  + ]:         50 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     957      [ +  +  + ]:         35 :   switch(typ(n))
     958                 :            :   {
     959                 :          5 :     case t_VEC: n = gel(n,1); break;
     960                 :         15 :     case t_MAT: n = factorback(n); break;
     961                 :            :   }
     962                 :         50 :   return mkvec2(f, sqrtint(diviiexact(n, f)));
     963                 :            : }
     964                 :            : GEN
     965                 :         45 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
     966                 :            : 
     967                 :            : GEN
     968         [ +  + ]:        190 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
     969                 :            : 
     970                 :            : static long
     971                 :         10 : _mod4(GEN c) {
     972                 :         10 :   long r, s = signe(c);
     973         [ -  + ]:         10 :   if (!s) return 0;
     974         [ -  + ]:         10 :   r = mod4(c); if (s < 0) r = 4-r;
     975                 :         10 :   return r;
     976                 :            : }
     977                 :            : 
     978                 :            : long
     979                 :        415 : corediscs(long D, ulong *f)
     980                 :            : {
     981                 :            :   /* D = f^2 dK */
     982         [ -  + ]:        415 :   long dK = D>=0 ? coreu(D) : -(long) coreu(-(ulong) D);
     983                 :        415 :   ulong dKmod4 = ((ulong)dK)&3UL;
     984 [ +  + ][ +  + ]:        415 :   if (dKmod4 == 2 || dKmod4 == 3)
     985                 :         70 :     dK *= 4;
     986         [ +  - ]:        415 :   if (f) *f = usqrt((ulong)(D/dK));
     987                 :        415 :   return D;
     988                 :            : }
     989                 :            : 
     990                 :            : GEN
     991                 :          5 : coredisc(GEN n)
     992                 :            : {
     993                 :          5 :   pari_sp av = avma;
     994                 :          5 :   GEN c = core(n);
     995         [ -  + ]:          5 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     996                 :          5 :   return gerepileuptoint(av, shifti(c,2));
     997                 :            : }
     998                 :            : 
     999                 :            : GEN
    1000                 :          5 : coredisc2(GEN n)
    1001                 :            : {
    1002                 :          5 :   pari_sp av = avma;
    1003                 :          5 :   GEN y = core2_i(n);
    1004                 :          5 :   GEN c = gel(y,1), f = gel(y,2);
    1005         [ -  + ]:          5 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
    1006                 :          5 :   y = cgetg(3,t_VEC);
    1007                 :          5 :   gel(y,1) = shifti(c,2);
    1008                 :          5 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
    1009                 :            : }
    1010                 :            : 
    1011                 :            : GEN
    1012         [ +  + ]:         10 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
    1013                 :            : 
    1014                 :            : long
    1015                 :        135 : omega(GEN n)
    1016                 :            : {
    1017                 :        135 :   pari_sp av = avma;
    1018                 :            :   GEN F, P;
    1019         [ +  + ]:        135 :   if ((F = check_arith_non0(n,"omega"))) {
    1020                 :            :     long n;
    1021                 :         15 :     P = gel(F,1); n = lg(P)-1;
    1022 [ +  + ][ +  + ]:         15 :     if (n && equalim1(gel(P,1))) n--;
    1023                 :         15 :     return n;
    1024                 :            :   }
    1025         [ +  + ]:        110 :   else if (lgefint(n) == 3)
    1026                 :            :   {
    1027         [ -  + ]:         81 :     if (n[2] == 1) return 0;
    1028                 :         81 :     F = factoru(n[2]);
    1029                 :            :   }
    1030                 :            :   else
    1031                 :         29 :     F = absi_factor(n);
    1032                 :        125 :   P = gel(F,1); avma = av; return lg(P)-1;
    1033                 :            : }
    1034                 :            : 
    1035                 :            : long
    1036                 :        135 : bigomega(GEN n)
    1037                 :            : {
    1038                 :        135 :   pari_sp av = avma;
    1039                 :            :   GEN F, E;
    1040         [ +  + ]:        135 :   if ((F = check_arith_non0(n,"bigomega")))
    1041                 :            :   {
    1042                 :         15 :     GEN P = gel(F,1);
    1043                 :         15 :     long n = lg(P)-1;
    1044                 :         15 :     E = gel(F,2);
    1045 [ +  + ][ +  + ]:         15 :     if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
    1046                 :         15 :     E = ZV_to_zv(E);
    1047                 :            :   }
    1048         [ +  + ]:        110 :   else if (lgefint(n) == 3)
    1049                 :            :   {
    1050         [ -  + ]:         85 :     if (n[2] == 1) return 0;
    1051                 :         85 :     F = factoru(n[2]);
    1052                 :         85 :     E = gel(F,2);
    1053                 :            :   }
    1054                 :            :   else
    1055                 :         25 :     E = ZV_to_zv(gel(absi_factor(n), 2));
    1056                 :        125 :   avma = av; return zv_sum(E);
    1057                 :            : }
    1058                 :            : 
    1059                 :            : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
    1060                 :            : ulong
    1061                 :     106591 : eulerphiu_fact(GEN f)
    1062                 :            : {
    1063                 :     106591 :   GEN P = gel(f,1), E = gel(f,2);
    1064                 :     106591 :   long i, m = 1, l = lg(P);
    1065         [ +  + ]:     422052 :   for (i = 1; i < l; i++)
    1066                 :            :   {
    1067                 :     315461 :     ulong p = P[i], e = E[i];
    1068         [ -  + ]:     315461 :     if (!e) continue;
    1069         [ +  + ]:     315461 :     if (p == 2)
    1070         [ +  + ]:      80606 :     { if (e > 1) m <<= e-1; }
    1071                 :            :     else
    1072                 :            :     {
    1073                 :     234855 :       m *= (p-1);
    1074         [ +  + ]:     234855 :       if (e > 1) m *= upowuu(p, e-1);
    1075                 :            :     }
    1076                 :            :   }
    1077                 :     106591 :   return m;
    1078                 :            : }
    1079                 :            : ulong
    1080                 :     106271 : eulerphiu(ulong n)
    1081                 :            : {
    1082                 :     106271 :   pari_sp av = avma;
    1083                 :            :   GEN F;
    1084         [ -  + ]:     106271 :   if (!n) return 2;
    1085                 :     106271 :   F = factoru(n);
    1086                 :     106271 :   avma = av; return eulerphiu_fact(F);
    1087                 :            : }
    1088                 :            : GEN
    1089                 :     102995 : eulerphi(GEN n)
    1090                 :            : {
    1091                 :     102995 :   pari_sp av = avma;
    1092                 :            :   GEN Q, F, P, E;
    1093                 :            :   long i, l;
    1094                 :            : 
    1095         [ +  + ]:     102995 :   if ((F = check_arith_all(n,"eulerphi")))
    1096                 :            :   {
    1097                 :         20 :     F = clean_Z_factor(F);
    1098         [ -  + ]:         20 :     n = (typ(n) == t_VEC)? gel(n,1): factorback(F);
    1099         [ +  + ]:         20 :     if (lgefint(n) == 3)
    1100                 :            :     {
    1101                 :            :       ulong e;
    1102                 :         15 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
    1103                 :         15 :       e = eulerphiu_fact(F);
    1104                 :         15 :       avma = av; return utoipos(e);
    1105                 :            :     }
    1106                 :            :   }
    1107         [ +  + ]:     102975 :   else if (lgefint(n) == 3) return utoipos(eulerphiu(uel(n,2)));
    1108                 :            :   else
    1109                 :         34 :     F = absi_factor(n);
    1110         [ +  + ]:         39 :   if (!signe(n)) return gen_2;
    1111                 :         29 :   P = gel(F,1);
    1112                 :         29 :   E = gel(F,2); l = lg(P);
    1113                 :         29 :   Q = cgetg(l, t_VEC);
    1114         [ +  + ]:        188 :   for (i = 1; i < l; i++)
    1115                 :            :   {
    1116                 :        159 :     GEN p = gel(P,i), q;
    1117                 :        159 :     ulong v = itou(gel(E,i));
    1118                 :        159 :     q = subiu(p,1);
    1119 [ +  + ][ +  + ]:        159 :     if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
    1120                 :        159 :     gel(Q,i) = q;
    1121                 :            :   }
    1122                 :     102995 :   return gerepileuptoint(av, ZV_prod(Q));
    1123                 :            : }
    1124                 :            : 
    1125                 :            : static GEN
    1126                 :         44 : numdiv_aux(GEN F)
    1127                 :            : {
    1128                 :         44 :   GEN x, E = gel(F,2);
    1129                 :         44 :   long i, l = lg(E);
    1130                 :         44 :   x = cgetg(l, t_VECSMALL);
    1131         [ +  + ]:        223 :   for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
    1132                 :         44 :   return x;
    1133                 :            : }
    1134                 :            : GEN
    1135                 :        135 : numdiv(GEN n)
    1136                 :            : {
    1137                 :        135 :   pari_sp av = avma;
    1138                 :            :   GEN F, E;
    1139                 :            :   long i, l;
    1140         [ +  + ]:        135 :   if ((F = check_arith_non0(n,"numdiv")))
    1141                 :            :   {
    1142                 :         15 :     F = clean_Z_factor(F);
    1143                 :         15 :     E = numdiv_aux(F);
    1144                 :            :   }
    1145         [ +  + ]:        110 :   else if (lgefint(n) == 3)
    1146                 :            :   {
    1147         [ -  + ]:         81 :     if (n[2] == 1) return gen_1;
    1148                 :         81 :     F = factoru(n[2]);
    1149                 :         81 :     E = gel(F,2); l = lg(E);
    1150         [ +  + ]:        462 :     for (i=1; i<l; i++) E[i]++;
    1151                 :            :   }
    1152                 :            :   else
    1153                 :         29 :     E = numdiv_aux(absi_factor(n));
    1154                 :        125 :   return gerepileuptoint(av, zv_prod_Z(E));
    1155                 :            : }
    1156                 :            : 
    1157                 :            : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
    1158                 :            : static GEN
    1159                 :        391 : u_euler_sumdiv(ulong p, long v)
    1160                 :            : {
    1161                 :        391 :   GEN u = utoipos(1 + p); /* can't overflow */
    1162         [ +  + ]:        401 :   for (; v > 1; v--) u = addsi(1, mului(p, u));
    1163                 :        391 :   return u;
    1164                 :            : }
    1165                 :            : /* 1 + q + ... + q^v */
    1166                 :            : static GEN
    1167                 :       1414 : euler_sumdiv(GEN q, long v)
    1168                 :            : {
    1169                 :       1414 :   GEN u = addui(1, q);
    1170         [ +  + ]:       1694 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
    1171                 :       1414 :   return u;
    1172                 :            : }
    1173                 :            : static GEN
    1174                 :        907 : u_euler_sumdivk(ulong p, long v, long k) { return euler_sumdiv(powuu(p,k), v); }
    1175                 :            : static GEN
    1176                 :        338 : euler_sumdivk(GEN p, long v, long k) { return euler_sumdiv(powiu(p,k), v); }
    1177                 :            : 
    1178                 :            : static GEN
    1179                 :         39 : sumdiv_aux(GEN F)
    1180                 :            : {
    1181                 :         39 :   GEN x, P = gel(F,1), E = gel(F,2);
    1182                 :         39 :   long i, l = lg(P);
    1183                 :         39 :   x = cgetg(l, t_VEC);
    1184         [ +  + ]:        208 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
    1185                 :         39 :   return x;
    1186                 :            : }
    1187                 :            : GEN
    1188                 :        135 : sumdiv(GEN n)
    1189                 :            : {
    1190                 :        135 :   pari_sp av = avma;
    1191                 :            :   GEN F, P, E;
    1192                 :            :   long i, l;
    1193                 :            : 
    1194         [ +  + ]:        135 :   if ((F = check_arith_non0(n,"sumdiv")))
    1195                 :            :   {
    1196                 :         15 :     F = clean_Z_factor(F);
    1197                 :         15 :     P = sumdiv_aux(F);
    1198                 :            :   }
    1199         [ +  + ]:        110 :   else if (lgefint(n) == 3)
    1200                 :            :   {
    1201         [ -  + ]:         86 :     if (n[2] == 1) return gen_1;
    1202                 :         86 :     F = factoru(n[2]);
    1203                 :         86 :     P = gel(F,1);
    1204                 :         86 :     E = gel(F,2); l = lg(P);
    1205         [ +  + ]:        477 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdiv(P[i], E[i]);
    1206                 :            :   }
    1207                 :            :   else
    1208                 :         24 :     P = sumdiv_aux(absi_factor(n));
    1209                 :        125 :   return gerepileuptoint(av, ZV_prod(P));
    1210                 :            : }
    1211                 :            : 
    1212                 :            : static GEN
    1213                 :         78 : sumdivk_aux(GEN F, long k)
    1214                 :            : {
    1215                 :         78 :   GEN x, P = gel(F,1), E = gel(F,2);
    1216                 :         78 :   long i, l = lg(P);
    1217                 :         78 :   x = cgetg(l, t_VEC);
    1218         [ +  + ]:        416 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdivk(gel(P,i), gel(E,i)[2], k);
    1219                 :         78 :   return x;
    1220                 :            : }
    1221                 :            : GEN
    1222                 :        480 : sumdivk(GEN n, long k)
    1223                 :            : {
    1224                 :        480 :   pari_sp av = avma;
    1225                 :            :   GEN E, F, P;
    1226                 :            :   long i, l, k1;
    1227                 :            : 
    1228         [ -  + ]:        480 :   if (!k) return numdiv(n);
    1229         [ +  + ]:        480 :   if (k == 1) return sumdiv(n);
    1230         [ -  + ]:        345 :   if (k ==-1) return gerepileupto(av, gdiv(sumdiv(n), n));
    1231                 :        345 :   k1 = k;
    1232         [ +  + ]:        345 :   if (k < 0)  k = -k;
    1233         [ +  + ]:        345 :   if ((F = check_arith_non0(n,"sumdivk")))
    1234                 :            :   {
    1235                 :         30 :     F = clean_Z_factor(F);
    1236                 :         30 :     P = sumdivk_aux(F,k);
    1237                 :            :   }
    1238         [ +  + ]:        295 :   else if (lgefint(n) == 3)
    1239                 :            :   {
    1240         [ +  + ]:        247 :     if (n[2] == 1) return gen_1;
    1241                 :        242 :     F = factoru(n[2]);
    1242                 :        242 :     P = gel(F,1);
    1243                 :        242 :     E = gel(F,2); l = lg(P);
    1244         [ +  + ]:       1119 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdivk(P[i], E[i], k);
    1245                 :            :   }
    1246                 :            :   else
    1247                 :         48 :     P = sumdivk_aux(absi_factor(n), k);
    1248                 :        320 :   P = ZV_prod(P);
    1249         [ +  + ]:        320 :   if (k1 > 0) return gerepileuptoint(av, P);
    1250                 :        450 :   return gerepileupto(av, gdiv(P, powiu(n,k)));
    1251                 :            : }
    1252                 :            : 
    1253                 :            : /* K t_VECSMALL of k >= 0 */
    1254                 :            : GEN
    1255                 :         20 : usumdivkvec(ulong n, GEN K)
    1256                 :            : {
    1257                 :         20 :   pari_sp av = avma;
    1258                 :         20 :   GEN F = factoru(n), P = gel(F,1), E = gel(F,2), Z, S;
    1259                 :         20 :   long i,j, l = lg(P), lK = lg(K);
    1260                 :         20 :   Z = cgetg(l, t_VEC);
    1261                 :         20 :   S = cgetg(lK, t_VEC);
    1262         [ +  + ]:         60 :   for (j=1; j<lK; j++)
    1263                 :            :   {
    1264                 :         40 :     long k = K[j];
    1265         [ +  + ]:         70 :     for (i=1; i<l; i++) gel(Z,i) = u_euler_sumdivk(P[i], E[i], k);
    1266                 :         40 :     gel(S,j) = ZV_prod(Z);
    1267                 :            :   }
    1268                 :         20 :   return gerepilecopy(av, S);
    1269                 :            : }
    1270                 :            : 
    1271                 :            : long
    1272                 :        135 : uissquarefree_fact(GEN f)
    1273                 :            : {
    1274                 :        135 :   GEN E = gel(f,2);
    1275                 :        135 :   long i, l = lg(E);
    1276         [ +  + ]:        265 :   for (i = 1; i < l; i++)
    1277         [ +  + ]:        200 :     if (E[i] > 1) return 0;
    1278                 :        135 :   return 1;
    1279                 :            : }
    1280                 :            : long
    1281                 :       5784 : uissquarefree(ulong n)
    1282                 :            : {
    1283         [ -  + ]:       5784 :   if (!n) return 0;
    1284                 :       5784 :   return moebiusu(n)? 1: 0;
    1285                 :            : }
    1286                 :            : long
    1287                 :        760 : Z_issquarefree(GEN n)
    1288                 :            : {
    1289      [ +  +  + ]:        760 :   switch(lgefint(n))
    1290                 :            :   {
    1291                 :          5 :     case 2: return 0;
    1292                 :          4 :     case 3: return uissquarefree(n[2]);
    1293                 :            :   }
    1294                 :        760 :   return moebius(n)? 1: 0;
    1295                 :            : }
    1296                 :            : long
    1297                 :         75 : issquarefree(GEN x)
    1298                 :            : {
    1299                 :            :   pari_sp av;
    1300                 :            :   GEN d;
    1301      [ +  +  - ]:         75 :   switch(typ(x))
    1302                 :            :   {
    1303                 :         10 :     case t_INT: return Z_issquarefree(x);
    1304                 :            :     case t_POL:
    1305         [ -  + ]:         65 :       if (!signe(x)) return 0;
    1306                 :         65 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
    1307                 :         65 :       avma = av; return (lg(d) == 3);
    1308                 :          0 :     default: pari_err_TYPE("issquarefree",x);
    1309                 :         75 :       return 0; /* not reached */
    1310                 :            :   }
    1311                 :            : }
    1312                 :            : 
    1313                 :            : /*********************************************************************/
    1314                 :            : /**                                                                 **/
    1315                 :            : /**                    DIGITS / SUM OF DIGITS                       **/
    1316                 :            : /**                                                                 **/
    1317                 :            : /*********************************************************************/
    1318                 :            : 
    1319                 :            : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
    1320                 :            : static void
    1321                 :      23825 : set_vexp(GEN v, long l)
    1322                 :            : {
    1323                 :            :   long m;
    1324         [ +  + ]:      33930 :   if (v[l]) return;
    1325                 :      10105 :   v[l] = 1; m = l>>1;
    1326                 :      10105 :   set_vexp(v, m);
    1327                 :      10105 :   set_vexp(v, l-m);
    1328                 :            : }
    1329                 :            : 
    1330                 :            : /* return all needed B^i for DAC algorithm, for lz digits */
    1331                 :            : static GEN
    1332                 :       3615 : get_vB(GEN T, long lz, void *E, struct bb_ring *r)
    1333                 :            : {
    1334                 :       3615 :   GEN vB, vexp = const_vecsmall(lz, 0);
    1335                 :       3615 :   long i, l = (lz+1) >> 1;
    1336                 :       3615 :   vexp[1] = 1;
    1337                 :       3615 :   vexp[2] = 1;
    1338                 :       3615 :   set_vexp(vexp, lz);
    1339                 :       3615 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
    1340                 :       3615 :   gel(vB, 1) = T;
    1341         [ +  + ]:      20655 :   for (i = 2; i <= l; i++)
    1342         [ +  + ]:      17040 :     if (vexp[i])
    1343                 :            :     {
    1344                 :      10105 :       long j = i >> 1;
    1345                 :      10105 :       GEN B2j = r->sqr(E, gel(vB,j));
    1346         [ +  + ]:      10105 :       gel(vB,i) = odd(i)? r->mul(E, B2j, T): B2j;
    1347                 :            :     }
    1348                 :       3615 :   return vB;
    1349                 :            : }
    1350                 :            : 
    1351                 :            : static void
    1352                 :      59875 : gen_digits_dac(GEN x, GEN vB, long l, GEN *z,
    1353                 :            :                void *E, GEN div(void *E, GEN a, GEN b, GEN *r))
    1354                 :            : {
    1355                 :            :   GEN q, r;
    1356                 :      59875 :   long m = l>>1;
    1357         [ +  + ]:      88552 :   if (l==1) { *z=x; return; }
    1358                 :      28677 :   q = div(E, x, gel(vB,m), &r);
    1359                 :      28677 :   gen_digits_dac(r, vB, m, z, E, div);
    1360                 :      28677 :   gen_digits_dac(q, vB, l-m, z+m, E, div);
    1361                 :            : }
    1362                 :            : 
    1363                 :            : static GEN
    1364                 :      13420 : gen_fromdigits_dac(GEN x, GEN vB, long i, long l, void *E,
    1365                 :            :                    GEN add(void *E, GEN a, GEN b),
    1366                 :            :                    GEN mul(void *E, GEN a, GEN b))
    1367                 :            : {
    1368                 :            :   GEN a, b;
    1369                 :      13420 :   long m = l>>1;
    1370         [ +  + ]:      13420 :   if (l==1) return gel(x,i);
    1371                 :       6170 :   a = gen_fromdigits_dac(x, vB, i, m, E, add, mul);
    1372                 :       6170 :   b = gen_fromdigits_dac(x, vB, i+m, l-m, E, add, mul);
    1373                 :      13420 :   return add(E, a, mul(E, b, gel(vB, m)));
    1374                 :            : }
    1375                 :            : 
    1376                 :            : static GEN
    1377                 :       2536 : gen_digits_i(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                 :            :   GEN z, vB;
    1381         [ +  + ]:       2536 :   if (n==1) retmkvec(gcopy(x));
    1382                 :       2521 :   vB = get_vB(B, n, E, r);
    1383                 :       2521 :   z = cgetg(n+1, t_VEC);
    1384                 :       2521 :   gen_digits_dac(x, vB, n, (GEN*)(z+1), E, div);
    1385                 :       2536 :   return z;
    1386                 :            : }
    1387                 :            : 
    1388                 :            : GEN
    1389                 :       2500 : gen_digits(GEN x, GEN B, long n, void *E, struct bb_ring *r,
    1390                 :            :                           GEN (*div)(void *E, GEN x, GEN y, GEN *r))
    1391                 :            : {
    1392                 :       2500 :   pari_sp av = avma;
    1393                 :       2500 :   return gerepilecopy(av, gen_digits_i(x, B, n, E, r, div));
    1394                 :            : }
    1395                 :            : 
    1396                 :            : GEN
    1397                 :       1080 : gen_fromdigits(GEN x, GEN B, void *E, struct bb_ring *r)
    1398                 :            : {
    1399                 :       1080 :   pari_sp av = avma;
    1400                 :       1080 :   long n = lg(x)-1;
    1401                 :       1080 :   GEN vB = get_vB(B, n, E, r);
    1402                 :       1080 :   GEN z = gen_fromdigits_dac(x, vB, 1, n, E, r->add, r->mul);
    1403                 :       1080 :   return gerepilecopy(av, z);
    1404                 :            : }
    1405                 :            : 
    1406                 :            : static GEN
    1407                 :       2400 : _addii(void *data /* ignored */, GEN x, GEN y)
    1408                 :       2400 : { (void)data; return addii(x,y); }
    1409                 :            : static GEN
    1410                 :        395 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
    1411                 :            : static GEN
    1412                 :       2585 : _mulii(void *data /* ignored */, GEN x, GEN y)
    1413                 :       2585 :  { (void)data; return mulii(x,y); }
    1414                 :            : static GEN
    1415                 :        312 : _dvmdii(void *data /* ignored */, GEN x, GEN y, GEN *r)
    1416                 :        312 :  { (void)data; return dvmdii(x,y,r); }
    1417                 :            : 
    1418                 :            : static struct bb_ring Z_ring = { _addii, _mulii, _sqri };
    1419                 :            : 
    1420                 :            : static GEN
    1421                 :        120 : check_basis(GEN B)
    1422                 :            : {
    1423         [ +  + ]:        120 :   if (!B) return utoipos(10);
    1424         [ -  + ]:        100 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
    1425         [ -  + ]:        100 :   if (cmpis(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
    1426                 :        120 :   return B;
    1427                 :            : }
    1428                 :            : 
    1429                 :            : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
    1430                 :            : static void
    1431                 :       2180 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
    1432                 :            : {
    1433                 :       2180 :   pari_sp av = avma;
    1434                 :            :   GEN q,r;
    1435                 :            :   long m;
    1436         [ +  + ]:       3263 :   if (l==1) { *z=itou(x); return; }
    1437                 :       1083 :   m=l>>1;
    1438                 :       1083 :   q = dvmdii(x, gel(vB,m), &r);
    1439                 :       1083 :   digits_dacsmall(q,vB,l-m,z);
    1440                 :       1083 :   digits_dacsmall(r,vB,m,z+l-m);
    1441                 :       1083 :   avma = av;
    1442                 :            : }
    1443                 :            : 
    1444                 :            : GEN
    1445                 :         35 : digits(GEN x, GEN B)
    1446                 :            : {
    1447                 :         35 :   pari_sp av=avma;
    1448                 :            :   long lz;
    1449                 :            :   GEN z, vB;
    1450         [ -  + ]:         35 :   if (typ(x)!=t_INT) pari_err_TYPE("digits",x);
    1451                 :         35 :   B = check_basis(B);
    1452         [ +  + ]:         35 :   if (!signe(x))       {avma = av; return cgetg(1,t_VEC); }
    1453         [ -  + ]:         30 :   if (absi_cmp(x,B)<0) {avma = av; retmkvec(absi(x)); }
    1454         [ +  + ]:         30 :   if (Z_ispow2(B))
    1455                 :            :   {
    1456                 :         10 :     long k = expi(B);
    1457         [ +  + ]:         10 :     if (k == 1) return binaire(x);
    1458         [ -  + ]:          5 :     if (k < BITS_IN_LONG)
    1459                 :            :     {
    1460                 :          0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
    1461                 :          0 :       z = binary_2k_zv(x, k);
    1462                 :          0 :       avma = av; return Flv_to_ZV(z);
    1463                 :            :     }
    1464                 :            :     else
    1465                 :            :     {
    1466                 :          5 :       avma = av; return binary_2k(x, k);
    1467                 :            :     }
    1468                 :            :   }
    1469         [ -  + ]:         20 :   if (signe(x) < 0) x = absi(x);
    1470                 :         20 :   lz = logint(x,B,NULL);
    1471         [ +  + ]:         20 :   if (lgefint(B)>3)
    1472                 :            :   {
    1473                 :          6 :     z = gerepileupto(av, gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii));
    1474                 :          6 :     vecreverse_inplace(z);
    1475                 :          6 :     return z;
    1476                 :            :   }
    1477                 :            :   else
    1478                 :            :   {
    1479                 :         14 :     vB = get_vB(B, lz, NULL, &Z_ring);
    1480                 :         14 :     (void)new_chunk(3*lz); /* HACK */
    1481                 :         14 :     z = zero_zv(lz);
    1482                 :         14 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1483                 :         35 :     avma = av; return vecsmall_to_vec(z);
    1484                 :            :   }
    1485                 :            : }
    1486                 :            : 
    1487                 :            : GEN
    1488                 :         30 : fromdigits(GEN x, GEN B)
    1489                 :            : {
    1490                 :         30 :   pari_sp av = avma;
    1491 [ +  - ][ -  + ]:         30 :   if (typ(x)!=t_VEC || !RgV_is_ZV(x)) pari_err_TYPE("fromdigits",x);
    1492                 :         30 :   B = check_basis(B);
    1493         [ +  + ]:         30 :   if (lg(x)==1) { avma = av; return gen_0; }
    1494                 :         25 :   x = vecreverse(x);
    1495                 :         30 :   return gerepileupto(av, gen_fromdigits(x, B, NULL, &Z_ring));
    1496                 :            : }
    1497                 :            : 
    1498                 :            : static ulong DS[] ={
    1499                 :            :   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,
    1500                 :            :   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,
    1501                 :            :   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,
    1502                 :            :   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,
    1503                 :            :   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,
    1504                 :            :   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,
    1505                 :            :   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,
    1506                 :            :   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,
    1507                 :            :   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,
    1508                 :            :   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,
    1509                 :            :   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,
    1510                 :            :   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,
    1511                 :            :   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,
    1512                 :            :   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,
    1513                 :            :   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,
    1514                 :            :   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,
    1515                 :            :   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,
    1516                 :            :   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,
    1517                 :            :   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,
    1518                 :            :   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,
    1519                 :            :   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,
    1520                 :            :   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,
    1521                 :            :   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,
    1522                 :            :   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,
    1523                 :            :   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,
    1524                 :            :   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,
    1525                 :            :   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,
    1526                 :            :   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,
    1527                 :            :   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,
    1528                 :            :   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,
    1529                 :            :   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,
    1530                 :            :   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,
    1531                 :            :   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,
    1532                 :            :   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,
    1533                 :            :   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,
    1534                 :            :   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,
    1535                 :            :   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,
    1536                 :            :   25,26,27
    1537                 :            : };
    1538                 :            : 
    1539                 :            : ulong
    1540                 :     253680 : sumdigitsu(ulong n)
    1541                 :            : {
    1542                 :     253680 :   ulong s = 0;
    1543         [ +  + ]:     972745 :   while (n) { s += DS[n % 1000]; n /= 1000; }
    1544                 :     253680 :   return s;
    1545                 :            : }
    1546                 :            : 
    1547                 :            : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1548                 :            : static ulong
    1549                 :         10 : sumdigits_block(ulong *res, long l)
    1550                 :            : {
    1551                 :         10 :   long s = sumdigitsu(*--res);
    1552         [ +  + ]:     253670 :   while (--l > 0) s += sumdigitsu(*--res);
    1553                 :         10 :   return s;
    1554                 :            : }
    1555                 :            : 
    1556                 :            : GEN
    1557                 :         25 : sumdigits(GEN n)
    1558                 :            : {
    1559                 :         25 :   pari_sp av = avma;
    1560                 :            :   ulong s, *res;
    1561                 :            :   long l;
    1562                 :            : 
    1563         [ -  + ]:         25 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1564                 :         25 :   l = lgefint(n);
    1565      [ +  +  + ]:         25 :   switch(l)
    1566                 :            :   {
    1567                 :          5 :     case 2: return gen_0;
    1568                 :         10 :     case 3: return utoipos(sumdigitsu(n[2]));
    1569                 :            :   }
    1570                 :         10 :   res = convi(n, &l);
    1571         [ +  - ]:         10 :   if ((ulong)l < ULONG_MAX / 81)
    1572                 :            :   {
    1573                 :         10 :     s = sumdigits_block(res, l);
    1574                 :         10 :     avma = av; return utoipos(s);
    1575                 :            :   }
    1576                 :            :   else /* Huge. Overflows ulong */
    1577                 :            :   {
    1578                 :          0 :     const long L = (long)(ULONG_MAX / 81);
    1579                 :          0 :     GEN S = gen_0;
    1580         [ #  # ]:          0 :     while (l > L)
    1581                 :            :     {
    1582                 :          0 :       S = addiu(S, sumdigits_block(res, L));
    1583                 :          0 :       res += L; l -= L;
    1584                 :            :     }
    1585         [ #  # ]:          0 :     if (l)
    1586                 :          0 :       S = addiu(S, sumdigits_block(res, l));
    1587                 :         25 :     return gerepileuptoint(av, S);
    1588                 :            :   }
    1589                 :            : }
    1590                 :            : 
    1591                 :            : GEN
    1592                 :         75 : sumdigits0(GEN x, GEN B)
    1593                 :            : {
    1594                 :         75 :   pari_sp av = avma;
    1595                 :            :   GEN z;
    1596                 :            :   long lz;
    1597                 :            : 
    1598         [ +  + ]:         75 :   if (!B) return sumdigits(x);
    1599         [ -  + ]:         55 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1600                 :         55 :   B = check_basis(B);
    1601         [ +  + ]:         55 :   if (Z_ispow2(B))
    1602                 :            :   {
    1603                 :         20 :     long k = expi(B);
    1604         [ +  + ]:         20 :     if (k == 1) { avma = av; return utoi(hammingweight(x)); }
    1605         [ +  + ]:         15 :     if (k < BITS_IN_LONG)
    1606                 :            :     {
    1607                 :         10 :       GEN z = binary_2k_zv(x, k);
    1608         [ -  + ]:         10 :       if (lg(z)-1 > 1<<(BITS_IN_LONG-k)) /* may overflow */
    1609                 :          0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1610                 :         10 :       avma = av; return utoi(zv_sum(z));
    1611                 :            :     }
    1612                 :          5 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1613                 :            :   }
    1614         [ -  + ]:         35 :   if (!signe(x))       { avma = av; return gen_0; }
    1615         [ -  + ]:         35 :   if (absi_cmp(x,B)<0) { avma = av; return absi(x); }
    1616         [ +  + ]:         35 :   if (equaliu(B,10))   { avma = av; return sumdigits(x); }
    1617                 :         30 :   lz = logint(x,B,NULL);
    1618                 :         30 :   z = gen_digits_i(x, B, lz, NULL, &Z_ring, _dvmdii);
    1619                 :         75 :   return gerepileuptoint(av, ZV_sum(z));
    1620                 :            : }

Generated by: LCOV version 1.9