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 16393-29b9383) Lines: 575 675 85.2 %
Date: 2014-04-24 Functions: 68 76 89.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 324 461 70.3 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /*********************************************************************/
      15                 :            : /**                                                                 **/
      16                 :            : /**                     ARITHMETIC FUNCTIONS                        **/
      17                 :            : /**                        (second part)                            **/
      18                 :            : /**                                                                 **/
      19                 :            : /*********************************************************************/
      20                 :            : #include "pari.h"
      21                 :            : #include "paripriv.h"
      22                 :            : 
      23                 :            : static ulong _maxprime = 0;
      24                 :            : static ulong diffptrlen;
      25                 :            : 
      26                 :            : /* Building/Rebuilding the diffptr table. The actual work is done by the
      27                 :            :  * following two subroutines;  the user entry point is the function
      28                 :            :  * initprimes() below.  initprimes1() is the old algorithm, called when
      29                 :            :  * maxnum (size) is moderate. Must be called after pari_init_stack() )*/
      30                 :            : static void
      31                 :        832 : initprimes1(ulong size, long *lenp, ulong *lastp, byteptr p1)
      32                 :            : {
      33                 :        832 :   pari_sp av = avma;
      34                 :            :   long k;
      35                 :        832 :   byteptr q, r, s, p = (byteptr)stack_calloc(size+2), fin = p + size;
      36                 :            : 
      37         [ +  + ]:       8320 :   for (r=q=p,k=1; r<=fin; )
      38                 :            :   {
      39         [ +  + ]:      11648 :     do { r+=k; k+=2; r+=k; } while (*++q);
      40         [ +  + ]:     264576 :     for (s=r; s<=fin; s+=k) *s = 1;
      41                 :            :   }
      42                 :        832 :   r = p1; *r++ = 2; *r++ = 1; /* 2 and 3 */
      43                 :        832 :   for (s=q=p+1; ; s=q)
      44                 :            :   {
      45         [ +  + ]:     293696 :     do q++; while (*q);
      46         [ +  + ]:     104000 :     if (q > fin) break;
      47                 :     103168 :     *r++ = (unsigned char) ((q-s) << 1);
      48                 :     103168 :   }
      49                 :        832 :   *r++ = 0;
      50                 :        832 :   *lenp = r - p1;
      51                 :        832 :   *lastp = ((s - p) << 1) + 1;
      52                 :        832 :   avma = av;
      53                 :        832 : }
      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                 :        832 : good_arena_size(ulong slow2_size, ulong total, ulong fixed_to_cache,
     198                 :            :                 cache_model_t *cache_model)
     199                 :            : {
     200                 :        832 :   ulong asize, cache_arena = cache_model->arena;
     201                 :            :   double Xmin, Xmax, A, B, C1, C2, D, V;
     202                 :        832 :   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         [ -  + ]:        832 :   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         [ +  - ]:        832 :   if (cache_arena - fixed_to_cache > 10 * slow2_size) {
     237                 :        832 :       asize = cache_arena - fixed_to_cache;
     238         [ -  + ]:        832 :       if (asize > total) asize = total; /* Automatically false... */
     239                 :        832 :       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                 :        832 :   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                 :       1984 : sieve_chunk(byteptr known_primes, ulong s, byteptr data, ulong n)
     342                 :            : {
     343                 :       1984 :   ulong p, cnt = n-1, start = s, delta = 1;
     344                 :            :   byteptr q;
     345                 :            : 
     346                 :       1984 :   memset(data, 0, n);
     347                 :       1984 :   start >>= 1;  /* (start - 1)/2 */
     348                 :       1984 :   start += n; /* Corresponds to the end */
     349                 :            :   /* data corresponds to start, q runs over primediffs */
     350         [ +  + ]:     229216 :   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                 :     227232 :     long off = cnt - ((start+(p>>1)) % p);
     357         [ +  + ]:  338868000 :     while (off >= 0) { data[off] = 1; off -= p; }
     358                 :            :   }
     359                 :       1984 : }
     360                 :            : 
     361                 :            : /* assume maxnum <= 436273290 < 2^29 */
     362                 :            : static void
     363                 :        832 : initprimes0(ulong maxnum, long *lenp, ulong *lastp, byteptr p1)
     364                 :            : {
     365                 :        832 :   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                 :        832 :   maxnum |= 1; /* make it odd. */
     373                 :            :   /* base case */
     374         [ -  + ]:       1664 :   if (maxnum < 1ul<<17) { initprimes1(maxnum>>1, lenp, lastp, p1); return; }
     375                 :            : 
     376                 :            :   /* Checked to be enough up to 40e6, attained at 155893 */
     377                 :        832 :   rootnum = usqrt(maxnum) | 1;
     378                 :        832 :   initprimes1(rootnum>>1, &psize, &last, p1);
     379                 :        832 :   end1 = p1 + psize - 1;
     380                 :        832 :   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                 :        832 :   asize = good_arena_size((ulong)(rootnum * slow2_in_roots), remains+1, 0,
     385                 :            :                           &cache_model) - 1;
     386                 :            :   /* enough room on the stack ? */
     387                 :        832 :   alloced = (((byteptr)avma) <= ((byteptr)bot) + asize);
     388         [ -  + ]:        832 :   if (alloced)
     389                 :          0 :     p = (byteptr)pari_malloc(asize+1);
     390                 :            :   else
     391                 :        832 :     p = (byteptr)stack_malloc(asize+1);
     392                 :        832 :   end = p + asize; /* the 0 sentinel goes at end. */
     393                 :        832 :   curlow = last + 2; /* First candidate: know primes up to last (odd). */
     394                 :        832 :   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                 :        832 :   plast = p - 1;
     400                 :        832 :   p_prime_above = p1 + 2;
     401                 :        832 :   prime_above = 3;
     402         [ +  + ]:       2816 :   while (remains)
     403                 :            :   { /* cycle over arenas; performance not crucial */
     404                 :            :     unsigned char was_delta;
     405         [ +  + ]:       1984 :     if (asize > remains) { asize = remains; end = p + asize; }
     406                 :            :     /* Fake the upper limit appropriate for the given arena */
     407 [ +  + ][ +  - ]:     105984 :     while (prime_above*prime_above <= curlow + (asize << 1) && *p_prime_above)
     408                 :     104000 :       prime_above += *p_prime_above++;
     409                 :       1984 :     was_delta = *p_prime_above;
     410                 :       1984 :     *p_prime_above = 0; /* sentinel for sieve_chunk */
     411                 :       1984 :     sieve_chunk(p1, curlow, p, asize);
     412                 :       1984 :     *p_prime_above = was_delta; /* restore */
     413                 :            : 
     414                 :       1984 :     p[asize] = 0; /* sentinel */
     415                 :       1984 :     for (q = p; ; plast = q++)
     416                 :            :     { /* q runs over addresses corresponding to primes */
     417         [ +  + ]:  207710784 :       while (*q) q++; /* use sentinel at end */
     418         [ +  + ]:   34456768 :       if (q >= end) break;
     419                 :   34454784 :       *curdiff++ = (unsigned char)(q-plast) << 1; /* < 255 for q < 436273291 */
     420                 :   34454784 :     }
     421                 :       1984 :     plast -= asize;
     422                 :       1984 :     remains -= asize;
     423                 :       1984 :     curlow += (asize<<1);
     424                 :            :   }
     425                 :        832 :   last = curlow - ((p - plast) << 1);
     426                 :        832 :   *curdiff++ = 0; /* sentinel */
     427                 :        832 :   *lenp = curdiff - p1;
     428                 :        832 :   *lastp = last;
     429         [ -  + ]:        832 :   if (alloced) pari_free(p); else avma = av;
     430                 :            : }
     431                 :            : 
     432                 :            : ulong
     433         [ +  - ]:    9614526 : 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 <= 436273290: the LHS ensures modular function
     439                 :            :  * have enough fast primes to work, the RHS ensures that p_{n+1} - p_n < 255 */
     440                 :            : byteptr
     441                 :        832 : initprimes(ulong maxnum, long *lenp, ulong *lastp)
     442                 :            : {
     443                 :            :   byteptr t;
     444         [ -  + ]:        832 :   if (maxnum < 65537)
     445                 :          0 :     maxnum = 65537;
     446         [ -  + ]:        832 :   else if (maxnum > 436273290)
     447                 :          0 :     maxnum = 436273290;
     448                 :        832 :   t = (byteptr)pari_malloc((size_t) (1.09 * maxnum/log((double)maxnum)) + 146);
     449                 :        832 :   initprimes0(maxnum, lenp, lastp, t);
     450                 :        832 :   return (byteptr)pari_realloc(t, *lenp);
     451                 :            : }
     452                 :            : 
     453                 :            : void
     454                 :        832 : initprimetable(ulong maxnum)
     455                 :            : {
     456                 :            :   long len;
     457                 :            :   ulong last;
     458                 :        832 :   byteptr p = initprimes(maxnum, &len, &last), old = diffptr;
     459                 :        832 :   diffptrlen = minss(diffptrlen, len);
     460                 :        832 :   _maxprime  = minss(_maxprime,last); /*Protect against ^C*/
     461                 :        832 :   diffptr = p; diffptrlen = len; _maxprime = last;
     462         [ -  + ]:        832 :   if (old) free(old);
     463                 :        832 : }
     464                 :            : 
     465                 :            : /* all init_primepointer_xx routines set *ptr to the corresponding place
     466                 :            :  * in prime table */
     467                 :            : /* smallest p >= a */
     468                 :            : ulong
     469                 :          0 : init_primepointer_geq(ulong a, byteptr *pd)
     470                 :            : {
     471                 :            :   ulong n, p;
     472                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     473                 :          0 :   return p;
     474                 :            : }
     475                 :            : /* largest p < a */
     476                 :            : ulong
     477                 :    4494746 : init_primepointer_lt(ulong a, byteptr *pd)
     478                 :            : {
     479                 :            :   ulong n, p;
     480                 :    4494746 :   prime_table_next_p(a, pd, &p, &n);
     481                 :    4494746 :   PREC_PRIME_VIADIFF(p, *pd);
     482                 :    4494746 :   return p;
     483                 :            : }
     484                 :            : /* largest p <= a */
     485                 :            : ulong
     486                 :          0 : init_primepointer_leq(ulong a, byteptr *pd)
     487                 :            : {
     488                 :            :   ulong n, p;
     489                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     490         [ #  # ]:          0 :   if (p != a) PREC_PRIME_VIADIFF(p, *pd);
     491                 :          0 :   return p;
     492                 :            : }
     493                 :            : /* smallest p > a */
     494                 :            : ulong
     495                 :          0 : init_primepointer_gt(ulong a, byteptr *pd)
     496                 :            : {
     497                 :            :   ulong n, p;
     498                 :          0 :   prime_table_next_p(a, pd, &p, &n);
     499         [ #  # ]:          0 :   if (p == a) NEXT_PRIME_VIADIFF(p, *pd);
     500                 :          0 :   return p;
     501                 :            : }
     502                 :            : 
     503                 :            : GEN
     504                 :         35 : boundfact(GEN n, ulong lim)
     505                 :            : {
     506      [ +  +  - ]:         35 :   switch(typ(n))
     507                 :            :   {
     508                 :         25 :     case t_INT: return Z_factor_limit(n,lim);
     509                 :            :     case t_FRAC: {
     510                 :         10 :       pari_sp av = avma;
     511                 :         10 :       GEN a = Z_factor_limit(gel(n,1),lim);
     512                 :         10 :       GEN b = Z_factor_limit(gel(n,2),lim);
     513                 :         10 :       gel(b,2) = ZC_neg(gel(b,2));
     514                 :         10 :       return gerepilecopy(av, merge_factor_i(a,b));
     515                 :            :     }
     516                 :            :   }
     517                 :          0 :   pari_err_TYPE("boundfact",n);
     518                 :         35 :   return NULL; /* not reached */
     519                 :            : }
     520                 :            : 
     521                 :            : /* NOT memory clean */
     522                 :            : GEN
     523                 :       7881 : Z_smoothen(GEN N, GEN L, GEN *pP, GEN *pe)
     524                 :            : {
     525                 :       7881 :   long i, j, l = lg(L);
     526                 :       7881 :   GEN e = new_chunk(l), P = new_chunk(l);
     527         [ +  + ]:      59284 :   for (i = j = 1; i < l; i++)
     528                 :            :   {
     529                 :      56881 :     ulong p = (ulong)L[i];
     530                 :      56881 :     long v = Z_lvalrem(N, p, &N);
     531 [ +  + ][ +  + ]:      56881 :     if (v) { P[j] = p; e[j] = v; j++; if (is_pm1(N)) { N = NULL; break; } }
     532                 :            :   }
     533                 :       7881 :   P[0] = evaltyp(t_VECSMALL) | evallg(j); *pP = P;
     534                 :       7881 :   e[0] = evaltyp(t_VECSMALL) | evallg(j); *pe = e; return N;
     535                 :            : }
     536                 :            : /***********************************************************************/
     537                 :            : /**                                                                   **/
     538                 :            : /**                    SIMPLE FACTORISATIONS                          **/
     539                 :            : /**                                                                   **/
     540                 :            : /***********************************************************************/
     541                 :            : /* Factor n and output [p,e,c] where
     542                 :            :  * p, e and c are vecsmall with n = prod{p[i]^e[i]} and c[i] = p[i]^e[i] */
     543                 :            : GEN
     544                 :      18031 : factoru_pow(ulong n)
     545                 :            : {
     546                 :      18031 :   GEN f = cgetg(4,t_VEC);
     547                 :      18031 :   pari_sp av = avma;
     548                 :            :   GEN F, P, E, p, e, c;
     549                 :            :   long i, l;
     550                 :            :   /* enough room to store <= 15 * [p,e,p^e] (OK if n < 2^64) */
     551                 :      18031 :   (void)new_chunk((15 + 1)*3);
     552                 :      18031 :   F = factoru(n);
     553                 :      18031 :   P = gel(F,1);
     554                 :      18031 :   E = gel(F,2); l = lg(P);
     555                 :      18031 :   avma = av;
     556                 :      18031 :   gel(f,1) = p = cgetg(l,t_VECSMALL);
     557                 :      18031 :   gel(f,2) = e = cgetg(l,t_VECSMALL);
     558                 :      18031 :   gel(f,3) = c = cgetg(l,t_VECSMALL);
     559         [ +  + ]:      58397 :   for(i = 1; i < l; i++)
     560                 :            :   {
     561                 :      40366 :     p[i] = P[i];
     562                 :      40366 :     e[i] = E[i];
     563                 :      40366 :     c[i] = upowuu(p[i], e[i]);
     564                 :            :   }
     565                 :      18031 :   return f;
     566                 :            : }
     567                 :            : 
     568                 :            : static GEN
     569                 :      16875 : factorlim(GEN n, ulong lim)
     570         [ -  + ]:      16875 : { return lim? Z_factor_limit(n, lim): Z_factor(n); }
     571                 :            : /* factor p^n - 1, assuming p prime. If lim != 0, limit factorization to
     572                 :            :  * primes <= lim */
     573                 :            : GEN
     574                 :       4718 : factor_pn_1_limit(GEN p, long n, ulong lim)
     575                 :            : {
     576                 :       4718 :   pari_sp av = avma;
     577                 :       4718 :   GEN A = factorlim(subiu(p,1), lim), d = divisorsu(n);
     578                 :       4718 :   long i, pp = itos_or_0(p);
     579         [ +  + ]:      13208 :   for(i=2; i<lg(d); i++)
     580                 :            :   {
     581                 :            :     GEN B;
     582 [ +  + ][ +  + ]:       8490 :     if (pp && d[i]%pp==0 && (
                 [ -  + ]
     583 [ #  # ][ +  + ]:       7623 :        ((pp&3)==1 && (d[i]&1)) ||
     584 [ +  + ][ +  + ]:       7623 :        ((pp&3)==3 && (d[i]&3)==2) ||
     585         [ +  + ]:       7558 :        (pp==2 && (d[i]&7)==4)))
     586                 :       3667 :     {
     587                 :       3667 :       GEN f=factor_Aurifeuille_prime(p,d[i]);
     588                 :       3667 :       B = factorlim(f, lim);
     589                 :       3667 :       A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     590                 :       3667 :       B = factorlim(diviiexact(polcyclo_eval(d[i],p), f), lim);
     591                 :            :     }
     592                 :            :     else
     593                 :       4823 :       B = factorlim(polcyclo_eval(d[i],p), lim);
     594                 :       8490 :     A = merge_factor(A, B, (void*)&cmpii, cmp_nodata);
     595                 :            :   }
     596                 :       4718 :   return gerepilecopy(av, A);
     597                 :            : }
     598                 :            : GEN
     599                 :       4718 : factor_pn_1(GEN p, ulong n)
     600                 :       4718 : { return factor_pn_1_limit(p, n, 0); }
     601                 :            : 
     602                 :            : #if 0
     603                 :            : static GEN
     604                 :            : to_mat(GEN p, long e) {
     605                 :            :   GEN B = cgetg(3, t_MAT);
     606                 :            :   gel(B,1) = mkcol(p);
     607                 :            :   gel(B,2) = mkcol(utoipos(e)); return B;
     608                 :            : }
     609                 :            : /* factor phi(n) */
     610                 :            : GEN
     611                 :            : factor_eulerphi(GEN n)
     612                 :            : {
     613                 :            :   pari_sp av = avma;
     614                 :            :   GEN B = NULL, A, P, E, AP, AE;
     615                 :            :   long i, l, v = vali(n);
     616                 :            : 
     617                 :            :   l = lg(n);
     618                 :            :   /* result requires less than this: at most expi(n) primes */
     619                 :            :   (void)new_chunk(bit_accuracy(l) * (l /*p*/ + 3 /*e*/ + 2 /*vectors*/) + 3+2);
     620                 :            :   if (v) { n = shifti(n, -v); v--; }
     621                 :            :   A = Z_factor(n); P = gel(A,1); E = gel(A,2); l = lg(P);
     622                 :            :   for(i = 1; i < l; i++)
     623                 :            :   {
     624                 :            :     GEN p = gel(P,i), q = subis(p,1), fa;
     625                 :            :     long e = itos(gel(E,i)), w;
     626                 :            : 
     627                 :            :     w = vali(q); v += w; q = shifti(q,-w);
     628                 :            :     if (! is_pm1(q))
     629                 :            :     {
     630                 :            :       fa = Z_factor(q);
     631                 :            :       B = B? merge_factor(B, fa, (void*)&cmpii, cmp_nodata): fa;
     632                 :            :     }
     633                 :            :     if (e > 1) {
     634                 :            :       if (B) {
     635                 :            :         gel(B,1) = shallowconcat(gel(B,1), p);
     636                 :            :         gel(B,2) = shallowconcat(gel(B,2), utoipos(e-1));
     637                 :            :       } else
     638                 :            :         B = to_mat(p, e-1);
     639                 :            :     }
     640                 :            :   }
     641                 :            :   avma = av;
     642                 :            :   if (!B) return v? to_mat(gen_2, v): trivial_fact();
     643                 :            :   A = cgetg(3, t_MAT);
     644                 :            :   P = gel(B,1); E = gel(B,2); l = lg(P);
     645                 :            :   AP = cgetg(l+1, t_COL); gel(A,1) = AP; AP++;
     646                 :            :   AE = cgetg(l+1, t_COL); gel(A,2) = AE; AE++;
     647                 :            :   /* prepend "2^v" */
     648                 :            :   gel(AP,0) = gen_2;
     649                 :            :   gel(AE,0) = utoipos(v);
     650                 :            :   for (i = 1; i < l; i++)
     651                 :            :   {
     652                 :            :     gel(AP,i) = icopy(gel(P,i));
     653                 :            :     gel(AE,i) = icopy(gel(E,i));
     654                 :            :   }
     655                 :            :   return A;
     656                 :            : }
     657                 :            : #endif
     658                 :            : 
     659                 :            : /***********************************************************************/
     660                 :            : /**                                                                   **/
     661                 :            : /**         CHECK FACTORIZATION FOR ARITHMETIC FUNCTIONS              **/
     662                 :            : /**                                                                   **/
     663                 :            : /***********************************************************************/
     664                 :            : static int
     665                 :     218788 : RgV_is_ZVpos(GEN v)
     666                 :            : {
     667                 :     218788 :   long i, l = lg(v);
     668         [ +  + ]:     579966 :   for (i = 1; i < l; i++)
     669                 :            :   {
     670                 :     361183 :     GEN c = gel(v,i);
     671 [ +  - ][ +  + ]:     361183 :     if (typ(c) != t_INT || signe(c) <= 0) return 0;
     672                 :            :   }
     673                 :     218788 :   return 1;
     674                 :            : }
     675                 :            : /* check whether v is a ZV with non-0 entries */
     676                 :            : static int
     677                 :        215 : RgV_is_ZVnon0(GEN v)
     678                 :            : {
     679                 :        215 :   long i, l = lg(v);
     680         [ +  + ]:        740 :   for (i = 1; i < l; i++)
     681                 :            :   {
     682                 :        560 :     GEN c = gel(v,i);
     683 [ +  - ][ +  + ]:        560 :     if (typ(c) != t_INT || !signe(c)) return 0;
     684                 :            :   }
     685                 :        215 :   return 1;
     686                 :            : }
     687                 :            : /* check whether v is a ZV with non-zero entries OR exactly [0] */
     688                 :            : static int
     689                 :        110 : RgV_is_ZV0(GEN v)
     690                 :            : {
     691                 :        110 :   long i, l = lg(v);
     692         [ +  + ]:        260 :   for (i = 1; i < l; i++)
     693                 :            :   {
     694                 :        190 :     GEN c = gel(v,i);
     695                 :            :     long s;
     696         [ -  + ]:        190 :     if (typ(c) != t_INT) return 0;
     697                 :        190 :     s = signe(c);
     698         [ +  + ]:        190 :     if (!s) return (l == 2);
     699                 :            :   }
     700                 :        110 :   return 1;
     701                 :            : }
     702                 :            : 
     703                 :            : static int
     704                 :     109564 : is_Z_factor_i(GEN f)
     705 [ +  - ][ +  + ]:     109564 : { return typ(f) == t_MAT && lg(f) == 3 && RgV_is_ZVpos(gel(f,2)); }
                 [ +  + ]
     706                 :            : int
     707                 :     109239 : is_Z_factorpos(GEN f)
     708 [ +  + ][ +  - ]:     109239 : { return is_Z_factor_i(f) && RgV_is_ZVpos(gel(f,1)); }
     709                 :            : int
     710                 :        110 : is_Z_factor(GEN f)
     711 [ +  - ][ +  - ]:        110 : { return is_Z_factor_i(f) && RgV_is_ZV0(gel(f,1)); }
     712                 :            : /* as is_Z_factorpos, also allow factor(0) */
     713                 :            : int
     714                 :        215 : is_Z_factornon0(GEN f)
     715 [ +  - ][ +  + ]:        215 : { return is_Z_factor_i(f) && RgV_is_ZVnon0(gel(f,1)); }
     716                 :            : GEN
     717                 :        115 : clean_Z_factor(GEN f)
     718                 :            : {
     719                 :        115 :   GEN P = gel(f,1);
     720                 :        115 :   long n = lg(P)-1;
     721 [ +  + ][ +  + ]:        115 :   if (n && equalim1(gel(P,1)))
     722                 :         35 :     return mkmat2(vecslice(P,2,n), vecslice(gel(f,2),2,n));
     723                 :        115 :   return f;
     724                 :            : }
     725                 :            : 
     726                 :            : /* n associated to a factorization of a positive integer: either N (t_INT)
     727                 :            :  * a factorization matrix faN, or a t_VEC: [N, faN] */
     728                 :            : GEN
     729                 :          0 : check_arith_pos(GEN n, const char *f) {
     730   [ #  #  #  # ]:          0 :   switch(typ(n))
     731                 :            :   {
     732                 :            :     case t_INT:
     733         [ #  # ]:          0 :       if (signe(n) <= 0 ) pari_err_DOMAIN(f, "argument", "<=", gen_0, gen_0);
     734                 :          0 :       return NULL;
     735                 :            :     case t_VEC:
     736 [ #  # ][ #  # ]:          0 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     737                 :          0 :       n = gel(n,2); /* fall through */
     738                 :            :     case t_MAT:
     739         [ #  # ]:          0 :       if (!is_Z_factorpos(n)) break;
     740                 :          0 :       return n;
     741                 :            :   }
     742                 :          0 :   pari_err_TYPE(f,n);
     743                 :          0 :   return NULL;
     744                 :            : }
     745                 :            : /* n associated to a factorization of a non-0 integer */
     746                 :            : GEN
     747                 :       2157 : check_arith_non0(GEN n, const char *f) {
     748   [ +  +  +  - ]:       2157 :   switch(typ(n))
     749                 :            :   {
     750                 :            :     case t_INT:
     751         [ +  + ]:       2007 :       if (!signe(n)) pari_err_DOMAIN(f, "argument", "=", gen_0, gen_0);
     752                 :       1967 :       return NULL;
     753                 :            :     case t_VEC:
     754 [ +  - ][ +  - ]:          5 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     755                 :          5 :       n = gel(n,2); /* fall through */
     756                 :            :     case t_MAT:
     757         [ +  + ]:        150 :       if (!is_Z_factornon0(n)) break;
     758                 :        115 :       return n;
     759                 :            :   }
     760                 :         35 :   pari_err_TYPE(f,n);
     761                 :       2082 :   return NULL;
     762                 :            : }
     763                 :            : /* n associated to a factorization of an integer */
     764                 :            : GEN
     765                 :     105490 : check_arith_all(GEN n, const char *f) {
     766   [ +  +  +  - ]:     105490 :   switch(typ(n))
     767                 :            :   {
     768                 :            :     case t_INT:
     769                 :     105380 :       return NULL;
     770                 :            :     case t_VEC:
     771 [ +  - ][ +  - ]:          5 :       if (lg(n) != 3 || typ(gel(n,1)) != t_INT) break;
     772                 :          5 :       n = gel(n,2); /* fall through */
     773                 :            :     case t_MAT:
     774         [ -  + ]:        110 :       if (!is_Z_factor(n)) break;
     775                 :        110 :       return n;
     776                 :            :   }
     777                 :          0 :   pari_err_TYPE(f,n);
     778                 :     105490 :   return NULL;
     779                 :            : }
     780                 :            : 
     781                 :            : /***********************************************************************/
     782                 :            : /**                                                                   **/
     783                 :            : /**                MISCELLANEOUS ARITHMETIC FUNCTIONS                 **/
     784                 :            : /**                (ultimately depend on Z_factor())                  **/
     785                 :            : /**                                                                   **/
     786                 :            : /***********************************************************************/
     787                 :            : /* set P,E from F. Check whether F is an integer and kill "factor" -1 */
     788                 :            : static void
     789                 :     250535 : set_fact_check(GEN F, GEN *pP, GEN *pE, int *isint)
     790                 :            : {
     791                 :            :   GEN E, P;
     792         [ -  + ]:     250535 :   if (lg(F) != 3) pari_err_TYPE("divisors",F);
     793                 :     250535 :   P = gel(F,1);
     794                 :     250535 :   E = gel(F,2);
     795                 :     250535 :   RgV_check_ZV(E, "divisors");
     796                 :     250535 :   *isint = RgV_is_ZV(P);
     797         [ +  + ]:     250535 :   if (*isint)
     798                 :            :   {
     799                 :     250525 :     long i, l = lg(P);
     800                 :            :     /* skip -1 */
     801 [ +  + ][ -  + ]:     250525 :     if (l>1 && signe(gel(P,1)) < 0) { E++; P = vecslice(P,2,--l); }
     802                 :            :     /* test for 0 */
     803         [ +  + ]:     816525 :     for (i = 1; i < l; i++)
     804 [ +  + ][ +  - ]:     566005 :       if (!signe(gel(P,i)) && signe(gel(E,i)))
     805                 :          5 :         pari_err_DOMAIN("divisors", "argument", "=", gen_0, F);
     806                 :            :   }
     807                 :     250530 :   *pP = P;
     808                 :     250530 :   *pE = E;
     809                 :     250530 : }
     810                 :            : static void
     811                 :       7590 : set_fact(GEN F, GEN *pP, GEN *pE) { *pP = gel(F,1); *pE = gel(F,2); }
     812                 :            : 
     813                 :            : int
     814                 :     258125 : divisors_init(GEN n, GEN *pP, GEN *pE)
     815                 :            : {
     816                 :            :   long i,l;
     817                 :            :   GEN E, P, e;
     818                 :            :   int isint;
     819                 :            : 
     820   [ +  +  +  + ]:     258125 :   switch(typ(n))
     821                 :            :   {
     822                 :            :     case t_INT:
     823         [ -  + ]:       7575 :       if (!signe(n)) pari_err_DOMAIN("divisors", "argument", "=", gen_0, gen_0);
     824                 :       7575 :       set_fact(absi_factor(n), &P,&E);
     825                 :       7575 :       isint = 1; break;
     826                 :            :     case t_VEC:
     827         [ -  + ]:          5 :       if (lg(n) != 3) pari_err_TYPE("divisors",n);
     828                 :          5 :       set_fact_check(gel(n,2), &P,&E, &isint);
     829                 :          5 :       break;
     830                 :            :     case t_MAT:
     831                 :     250530 :       set_fact_check(n, &P,&E, &isint);
     832                 :     250525 :       break;
     833                 :            :     default:
     834                 :         15 :       set_fact(factor(n), &P,&E);
     835                 :         15 :       isint = 0; break;
     836                 :            :   }
     837                 :     258120 :   l = lg(P);
     838                 :     258120 :   e = cgetg(l, t_VECSMALL);
     839         [ +  + ]:     845965 :   for (i=1; i<l; i++)
     840                 :            :   {
     841                 :     587850 :     e[i] = itos(gel(E,i));
     842         [ +  + ]:     587850 :     if (e[i] < 0) pari_err_TYPE("divisors [denominator]",n);
     843                 :            :   }
     844                 :     258115 :   *pP = P; *pE = e; return isint;
     845                 :            : }
     846                 :            : 
     847                 :            : GEN
     848                 :     258095 : divisors(GEN n)
     849                 :            : {
     850                 :     258095 :   pari_sp av = avma;
     851                 :            :   long i, j, l;
     852                 :            :   ulong ndiv;
     853                 :            :   GEN *d, *t, *t1, *t2, *t3, P, E, e;
     854                 :     258095 :   int isint = divisors_init(n, &P, &E);
     855                 :            : 
     856                 :     258085 :   l = lg(E); e = cgetg(l, t_VECSMALL);
     857         [ +  + ]:     845820 :   for (i=1; i<l; i++) e[i] = E[i]+1;
     858                 :     258085 :   ndiv = itou_or_0( zv_prod_Z(e) );
     859 [ +  - ][ -  + ]:     258085 :   if (!ndiv || ndiv & ~LGBITS) pari_err_OVERFLOW("divisors");
     860                 :     258085 :   d = t = (GEN*) cgetg(ndiv+1,t_VEC);
     861                 :     258085 :   *++d = gen_1;
     862         [ +  + ]:     258085 :   if (isint)
     863                 :            :   {
     864         [ +  + ]:     845780 :     for (i=1; i<l; i++)
     865         [ +  + ]:    1578705 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     866         [ +  + ]:    3826665 :         for (t2=d, t3=t1; t3<t2; ) *++d = mulii(*++t3, gel(P,i));
     867                 :     258075 :     e = ZV_sort((GEN)t);
     868                 :            :   } else {
     869         [ +  + ]:         40 :     for (i=1; i<l; i++)
     870         [ +  + ]:        130 :       for (t1=t,j=E[i]; j; j--,t1=t2)
     871         [ +  + ]:        810 :         for (t2=d, t3=t1; t3<t2; ) *++d = gmul(*++t3, gel(P,i));
     872                 :         10 :     e = (GEN)t;
     873                 :            :   }
     874                 :     258085 :   return gerepileupto(av, e);
     875                 :            : }
     876                 :            : 
     877                 :            : GEN
     878                 :      21706 : divisorsu(ulong n)
     879                 :            : {
     880                 :      21706 :   pari_sp av = avma;
     881                 :            :   long i, j, l;
     882                 :            :   ulong nbdiv;
     883                 :      21706 :   GEN d, t, t1, t2, t3, P, e, fa = factoru(n);
     884                 :            : 
     885                 :      21706 :   P = gel(fa,1); l = lg(P);
     886                 :      21706 :   e = gel(fa,2);
     887                 :      21706 :   nbdiv = 1;
     888         [ +  + ]:      48156 :   for (i=1; i<l; i++) nbdiv *= 1+e[i];
     889                 :      21706 :   d = t = cgetg(nbdiv+1,t_VECSMALL);
     890                 :      21706 :   *++d = 1;
     891         [ +  + ]:      48156 :   for (i=1; i<l; i++)
     892         [ +  + ]:      60251 :     for (t1=t,j=e[i]; j; j--,t1=t2)
     893         [ +  + ]:      89094 :       for (t2=d, t3=t1; t3<t2; ) *(++d) = *(++t3) * P[i];
     894                 :      21706 :   vecsmall_sort(t);
     895                 :      21706 :   return gerepileupto(av, t);
     896                 :            : }
     897                 :            : 
     898                 :            : static GEN
     899                 :          0 : corefa(GEN fa)
     900                 :            : {
     901                 :          0 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1;
     902                 :            :   long i;
     903         [ #  # ]:          0 :   for (i=1; i<lg(P); i++)
     904         [ #  # ]:          0 :     if (mod2(gel(E,i))) c = mulii(c,gel(P,i));
     905                 :          0 :   return c;
     906                 :            : }
     907                 :            : static GEN
     908                 :        120 : core2fa(GEN fa)
     909                 :            : {
     910                 :        120 :   GEN P = gel(fa,1), E = gel(fa,2), c = gen_1, f = gen_1;
     911                 :            :   long i;
     912         [ +  + ]:        315 :   for (i=1; i<lg(P); i++)
     913                 :            :   {
     914                 :        195 :     long e = itos(gel(E,i));
     915         [ +  + ]:        195 :     if (e & 1)  c = mulii(c, gel(P,i));
     916         [ +  + ]:        195 :     if (e != 1) f = mulii(f, powiu(gel(P,i), e >> 1));
     917                 :            :   }
     918                 :        120 :   return mkvec2(c,f);
     919                 :            : }
     920                 :            : GEN
     921                 :          0 : corepartial(GEN n, long all)
     922                 :            : {
     923                 :          0 :   pari_sp av = avma;
     924         [ #  # ]:          0 :   if (typ(n) != t_INT) pari_err_TYPE("corepartial",n);
     925                 :          0 :   return gerepileuptoint(av, corefa(Z_factor_limit(n,all)));
     926                 :            : }
     927                 :            : GEN
     928                 :        120 : core2partial(GEN n, long all)
     929                 :            : {
     930                 :        120 :   pari_sp av = avma;
     931         [ -  + ]:        120 :   if (typ(n) != t_INT) pari_err_TYPE("core2partial",n);
     932                 :        120 :   return gerepilecopy(av, core2fa(Z_factor_limit(n,all)));
     933                 :            : }
     934                 :            : static GEN
     935                 :         60 : core2_i(GEN n)
     936                 :            : {
     937                 :         60 :   GEN f = core(n);
     938         [ +  + ]:         60 :   if (!signe(f)) return mkvec2(gen_0, gen_1);
     939      [ +  +  + ]:         45 :   switch(typ(n))
     940                 :            :   {
     941                 :          5 :     case t_VEC: n = gel(n,1); break;
     942                 :         15 :     case t_MAT: n = factorback(n); break;
     943                 :            :   }
     944                 :         60 :   return mkvec2(f, sqrtint(diviiexact(n, f)));
     945                 :            : }
     946                 :            : GEN
     947                 :         50 : core2(GEN n) { pari_sp av = avma; return gerepilecopy(av, core2_i(n)); }
     948                 :            : 
     949                 :            : GEN
     950         [ +  + ]:        190 : core0(GEN n,long flag) { return flag? core2(n): core(n); }
     951                 :            : 
     952                 :            : static long
     953                 :         20 : _mod4(GEN c) {
     954                 :         20 :   long r, s = signe(c);
     955         [ -  + ]:         20 :   if (!s) return 0;
     956         [ -  + ]:         20 :   r = mod4(c); if (s < 0) r = 4-r;
     957                 :         20 :   return r;
     958                 :            : }
     959                 :            : 
     960                 :            : GEN
     961                 :         10 : coredisc(GEN n)
     962                 :            : {
     963                 :         10 :   pari_sp av = avma;
     964                 :         10 :   GEN c = core(n);
     965         [ -  + ]:         10 :   if (_mod4(c)<=1) return c; /* c = 0 or 1 mod 4 */
     966                 :         10 :   return gerepileuptoint(av, shifti(c,2));
     967                 :            : }
     968                 :            : 
     969                 :            : GEN
     970                 :         10 : coredisc2(GEN n)
     971                 :            : {
     972                 :         10 :   pari_sp av = avma;
     973                 :         10 :   GEN y = core2_i(n);
     974                 :         10 :   GEN c = gel(y,1), f = gel(y,2);
     975         [ -  + ]:         10 :   if (_mod4(c)<=1) return gerepilecopy(av, y);
     976                 :         10 :   y = cgetg(3,t_VEC);
     977                 :         10 :   gel(y,1) = shifti(c,2);
     978                 :         10 :   gel(y,2) = gmul2n(f,-1); return gerepileupto(av, y);
     979                 :            : }
     980                 :            : 
     981                 :            : GEN
     982         [ +  + ]:         10 : coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
     983                 :            : 
     984                 :            : long
     985                 :        140 : omega(GEN n)
     986                 :            : {
     987                 :        140 :   pari_sp av = avma;
     988                 :            :   GEN F, P;
     989         [ +  + ]:        140 :   if ((F = check_arith_non0(n,"omega"))) {
     990                 :            :     long n;
     991                 :         15 :     P = gel(F,1); n = lg(P)-1;
     992 [ +  + ][ +  + ]:         15 :     if (n && equalim1(gel(P,1))) n--;
     993                 :         15 :     return n;
     994                 :            :   }
     995         [ +  + ]:        115 :   else if (lgefint(n) == 3)
     996                 :            :   {
     997         [ -  + ]:         81 :     if (n[2] == 1) return 0;
     998                 :         81 :     F = factoru(n[2]);
     999                 :            :   }
    1000                 :            :   else
    1001                 :         34 :     F = absi_factor(n);
    1002                 :        130 :   P = gel(F,1); avma = av; return lg(P)-1;
    1003                 :            : }
    1004                 :            : 
    1005                 :            : long
    1006                 :        140 : bigomega(GEN n)
    1007                 :            : {
    1008                 :        140 :   pari_sp av = avma;
    1009                 :            :   GEN F, E;
    1010         [ +  + ]:        140 :   if ((F = check_arith_non0(n,"bigomega")))
    1011                 :            :   {
    1012                 :         15 :     GEN P = gel(F,1);
    1013                 :         15 :     long n = lg(P)-1;
    1014                 :         15 :     E = gel(F,2);
    1015 [ +  + ][ +  + ]:         15 :     if (n && equalim1(gel(P,1))) E = vecslice(E,2,n);
    1016                 :         15 :     E = ZV_to_zv(E);
    1017                 :            :   }
    1018         [ +  + ]:        115 :   else if (lgefint(n) == 3)
    1019                 :            :   {
    1020         [ -  + ]:         89 :     if (n[2] == 1) return 0;
    1021                 :         89 :     F = factoru(n[2]);
    1022                 :         89 :     E = gel(F,2);
    1023                 :            :   }
    1024                 :            :   else
    1025                 :         26 :     E = ZV_to_zv(gel(absi_factor(n), 2));
    1026                 :        130 :   avma = av; return zv_sum(E);
    1027                 :            : }
    1028                 :            : 
    1029                 :            : /* assume f = factoru(n), possibly with 0 exponents. Return phi(n) */
    1030                 :            : ulong
    1031                 :     106216 : eulerphiu_fact(GEN f)
    1032                 :            : {
    1033                 :     106216 :   GEN P = gel(f,1), E = gel(f,2);
    1034                 :     106216 :   long i, m = 1, l = lg(P);
    1035         [ +  + ]:     421207 :   for (i = 1; i < l; i++)
    1036                 :            :   {
    1037                 :     314991 :     ulong p = P[i], e = E[i];
    1038         [ -  + ]:     314991 :     if (!e) continue;
    1039         [ +  + ]:     314991 :     if (p == 2)
    1040         [ +  + ]:      80261 :     { if (e > 1) m <<= e-1; }
    1041                 :            :     else
    1042                 :            :     {
    1043                 :     234730 :       m *= (p-1);
    1044         [ +  + ]:     234730 :       if (e > 1) m *= upowuu(p, e-1);
    1045                 :            :     }
    1046                 :            :   }
    1047                 :     106216 :   return m;
    1048                 :            : }
    1049                 :            : ulong
    1050                 :     105986 : eulerphiu(ulong n)
    1051                 :            : {
    1052                 :     105986 :   pari_sp av = avma;
    1053                 :            :   GEN F;
    1054         [ -  + ]:     105986 :   if (!n) return 2;
    1055                 :     105986 :   F = factoru(n);
    1056                 :     105986 :   avma = av; return eulerphiu_fact(F);
    1057                 :            : }
    1058                 :            : GEN
    1059                 :     103000 : eulerphi(GEN n)
    1060                 :            : {
    1061                 :     103000 :   pari_sp av = avma;
    1062                 :            :   GEN Q, F, P, E;
    1063                 :            :   long i, l;
    1064                 :            : 
    1065         [ +  + ]:     103000 :   if ((F = check_arith_all(n,"eulerphi")))
    1066                 :            :   {
    1067                 :         20 :     F = clean_Z_factor(F);
    1068         [ -  + ]:         20 :     n = (typ(n) == t_VEC)? gel(n,1): factorback(F);
    1069         [ +  + ]:         20 :     if (lgefint(n) == 3)
    1070                 :            :     {
    1071                 :            :       ulong e;
    1072                 :         15 :       F = mkmat2(ZV_to_nv(gel(F,1)), ZV_to_nv(gel(F,2)));
    1073                 :         15 :       e = eulerphiu_fact(F);
    1074                 :         15 :       avma = av; return utoipos(e);
    1075                 :            :     }
    1076                 :            :   }
    1077         [ +  + ]:     102980 :   else if (lgefint(n) == 3) return utoipos(eulerphiu((ulong)n[2]));
    1078                 :            :   else
    1079                 :         34 :     F = absi_factor(n);
    1080         [ +  + ]:         39 :   if (!signe(n)) return gen_2;
    1081                 :         29 :   P = gel(F,1);
    1082                 :         29 :   E = gel(F,2); l = lg(P);
    1083                 :         29 :   Q = cgetg(l, t_VEC);
    1084         [ +  + ]:        188 :   for (i = 1; i < l; i++)
    1085                 :            :   {
    1086                 :        159 :     GEN p = gel(P,i), q;
    1087                 :        159 :     ulong v = itou(gel(E,i));
    1088                 :        159 :     q = subiu(p,1);
    1089 [ +  + ][ +  + ]:        159 :     if (v != 1) q = mulii(q, v == 2? p: powiu(p, v-1));
    1090                 :        159 :     gel(Q,i) = q;
    1091                 :            :   }
    1092                 :     103000 :   return gerepileuptoint(av, ZV_prod(Q));
    1093                 :            : }
    1094                 :            : 
    1095                 :            : static GEN
    1096                 :         49 : numdiv_aux(GEN F)
    1097                 :            : {
    1098                 :         49 :   GEN x, E = gel(F,2);
    1099                 :         49 :   long i, l = lg(E);
    1100                 :         49 :   x = cgetg(l, t_VECSMALL);
    1101         [ +  + ]:        238 :   for (i=1; i<l; i++) x[i] = itou(gel(E,i))+1;
    1102                 :         49 :   return x;
    1103                 :            : }
    1104                 :            : GEN
    1105                 :        140 : numdiv(GEN n)
    1106                 :            : {
    1107                 :        140 :   pari_sp av = avma;
    1108                 :            :   GEN F, E;
    1109                 :            :   long i, l;
    1110         [ +  + ]:        140 :   if ((F = check_arith_non0(n,"numdiv")))
    1111                 :            :   {
    1112                 :         15 :     F = clean_Z_factor(F);
    1113                 :         15 :     E = numdiv_aux(F);
    1114                 :            :   }
    1115         [ +  + ]:        115 :   else if (lgefint(n) == 3)
    1116                 :            :   {
    1117         [ -  + ]:         81 :     if (n[2] == 1) return gen_1;
    1118                 :         81 :     F = factoru(n[2]);
    1119                 :         81 :     E = gel(F,2); l = lg(E);
    1120         [ +  + ]:        462 :     for (i=1; i<l; i++) E[i]++;
    1121                 :            :   }
    1122                 :            :   else
    1123                 :         34 :     E = numdiv_aux(absi_factor(n));
    1124                 :        130 :   return gerepileuptoint(av, zv_prod_Z(E));
    1125                 :            : }
    1126                 :            : 
    1127                 :            : /* 1 + p + ... + p^v, p != 2^BIL - 1 */
    1128                 :            : static GEN
    1129                 :        401 : u_euler_sumdiv(ulong p, long v)
    1130                 :            : {
    1131                 :        401 :   GEN u = utoipos(1 + p); /* can't overflow */
    1132         [ +  + ]:        421 :   for (; v > 1; v--) u = addsi(1, mului(p, u));
    1133                 :        401 :   return u;
    1134                 :            : }
    1135                 :            : /* 1 + q + ... + q^v */
    1136                 :            : static GEN
    1137                 :       1434 : euler_sumdiv(GEN q, long v)
    1138                 :            : {
    1139                 :       1434 :   GEN u = addui(1, q);
    1140         [ +  + ]:       1734 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
    1141                 :       1434 :   return u;
    1142                 :            : }
    1143                 :            : static GEN
    1144                 :        927 : u_euler_sumdivk(ulong p, long v, long k) { return euler_sumdiv(powuu(p,k), v); }
    1145                 :            : static GEN
    1146                 :        338 : euler_sumdivk(GEN p, long v, long k) { return euler_sumdiv(powiu(p,k), v); }
    1147                 :            : 
    1148                 :            : static GEN
    1149                 :         39 : sumdiv_aux(GEN F)
    1150                 :            : {
    1151                 :         39 :   GEN x, P = gel(F,1), E = gel(F,2);
    1152                 :         39 :   long i, l = lg(P);
    1153                 :         39 :   x = cgetg(l, t_VEC);
    1154         [ +  + ]:        208 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdiv(gel(P,i), itou(gel(E,i)));
    1155                 :         39 :   return x;
    1156                 :            : }
    1157                 :            : GEN
    1158                 :        140 : sumdiv(GEN n)
    1159                 :            : {
    1160                 :        140 :   pari_sp av = avma;
    1161                 :            :   GEN F, P, E;
    1162                 :            :   long i, l;
    1163                 :            : 
    1164         [ +  + ]:        140 :   if ((F = check_arith_non0(n,"sumdiv")))
    1165                 :            :   {
    1166                 :         15 :     F = clean_Z_factor(F);
    1167                 :         15 :     P = sumdiv_aux(F);
    1168                 :            :   }
    1169         [ +  + ]:        115 :   else if (lgefint(n) == 3)
    1170                 :            :   {
    1171         [ -  + ]:         91 :     if (n[2] == 1) return gen_1;
    1172                 :         91 :     F = factoru(n[2]);
    1173                 :         91 :     P = gel(F,1);
    1174                 :         91 :     E = gel(F,2); l = lg(P);
    1175         [ +  + ]:        492 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdiv(P[i], E[i]);
    1176                 :            :   }
    1177                 :            :   else
    1178                 :         24 :     P = sumdiv_aux(absi_factor(n));
    1179                 :        130 :   return gerepileuptoint(av, ZV_prod(P));
    1180                 :            : }
    1181                 :            : 
    1182                 :            : static GEN
    1183                 :         78 : sumdivk_aux(GEN F, long k)
    1184                 :            : {
    1185                 :         78 :   GEN x, P = gel(F,1), E = gel(F,2);
    1186                 :         78 :   long i, l = lg(P);
    1187                 :         78 :   x = cgetg(l, t_VEC);
    1188         [ +  + ]:        416 :   for (i=1; i<l; i++) gel(x,i) = euler_sumdivk(gel(P,i), gel(E,i)[2], k);
    1189                 :         78 :   return x;
    1190                 :            : }
    1191                 :            : GEN
    1192                 :        490 : sumdivk(GEN n, long k)
    1193                 :            : {
    1194                 :        490 :   pari_sp av = avma;
    1195                 :            :   GEN E, F, P;
    1196                 :            :   long i, l, k1;
    1197                 :            : 
    1198         [ -  + ]:        490 :   if (!k) return numdiv(n);
    1199         [ +  + ]:        490 :   if (k == 1) return sumdiv(n);
    1200         [ -  + ]:        355 :   if (k ==-1) return gerepileupto(av, gdiv(sumdiv(n), n));
    1201                 :        355 :   k1 = k;
    1202         [ +  + ]:        355 :   if (k < 0)  k = -k;
    1203         [ +  + ]:        355 :   if ((F = check_arith_non0(n,"sumdivk")))
    1204                 :            :   {
    1205                 :         30 :     F = clean_Z_factor(F);
    1206                 :         30 :     P = sumdivk_aux(F,k);
    1207                 :            :   }
    1208         [ +  + ]:        305 :   else if (lgefint(n) == 3)
    1209                 :            :   {
    1210         [ +  + ]:        257 :     if (n[2] == 1) return gen_1;
    1211                 :        252 :     F = factoru(n[2]);
    1212                 :        252 :     P = gel(F,1);
    1213                 :        252 :     E = gel(F,2); l = lg(P);
    1214         [ +  + ]:       1149 :     for (i=1; i<l; i++) gel(P,i) = u_euler_sumdivk(P[i], E[i], k);
    1215                 :            :   }
    1216                 :            :   else
    1217                 :         48 :     P = sumdivk_aux(absi_factor(n), k);
    1218                 :        330 :   P = ZV_prod(P);
    1219         [ +  + ]:        330 :   if (k1 > 0) return gerepileuptoint(av, P);
    1220                 :        460 :   return gerepileupto(av, gdiv(P, powiu(n,k)));
    1221                 :            : }
    1222                 :            : 
    1223                 :            : /* K t_VECSMALL of k >= 0 */
    1224                 :            : GEN
    1225                 :         20 : usumdivkvec(ulong n, GEN K)
    1226                 :            : {
    1227                 :         20 :   pari_sp av = avma;
    1228                 :         20 :   GEN F = factoru(n), P = gel(F,1), E = gel(F,2), Z, S;
    1229                 :         20 :   long i,j, l = lg(P), lK = lg(K);
    1230                 :         20 :   Z = cgetg(l, t_VEC);
    1231                 :         20 :   S = cgetg(lK, t_VEC);
    1232         [ +  + ]:         60 :   for (j=1; j<lK; j++)
    1233                 :            :   {
    1234                 :         40 :     long k = K[j];
    1235         [ +  + ]:         70 :     for (i=1; i<l; i++) gel(Z,i) = u_euler_sumdivk(P[i], E[i], k);
    1236                 :         40 :     gel(S,j) = ZV_prod(Z);
    1237                 :            :   }
    1238                 :         20 :   return gerepilecopy(av, S);
    1239                 :            : }
    1240                 :            : 
    1241                 :            : long
    1242                 :         75 : uissquarefree_fact(GEN f)
    1243                 :            : {
    1244                 :         75 :   GEN E = gel(f,2);
    1245                 :         75 :   long i, l = lg(E);
    1246         [ +  + ]:        135 :   for (i = 1; i < l; i++)
    1247         [ +  + ]:        105 :     if (E[i] > 1) return 0;
    1248                 :         75 :   return 1;
    1249                 :            : }
    1250                 :            : long
    1251                 :       5378 : uissquarefree(ulong n)
    1252                 :            : {
    1253         [ -  + ]:       5378 :   if (!n) return 0;
    1254                 :       5378 :   return moebiusu(n)? 1: 0;
    1255                 :            : }
    1256                 :            : long
    1257                 :        765 : Z_issquarefree(GEN n)
    1258                 :            : {
    1259      [ +  +  + ]:        765 :   switch(lgefint(n))
    1260                 :            :   {
    1261                 :          5 :     case 2: return 0;
    1262                 :          8 :     case 3: return uissquarefree(n[2]);
    1263                 :            :   }
    1264                 :        765 :   return moebius(n)? 1: 0;
    1265                 :            : }
    1266                 :            : long
    1267                 :         25 : issquarefree(GEN x)
    1268                 :            : {
    1269                 :            :   pari_sp av;
    1270                 :            :   GEN d;
    1271      [ +  +  - ]:         25 :   switch(typ(x))
    1272                 :            :   {
    1273                 :         15 :     case t_INT: return Z_issquarefree(x);
    1274                 :            :     case t_POL:
    1275         [ -  + ]:         10 :       if (!signe(x)) return 0;
    1276                 :         10 :       av = avma; d = RgX_gcd(x, RgX_deriv(x));
    1277                 :         10 :       avma = av; return (lg(d) == 3);
    1278                 :          0 :     default: pari_err_TYPE("issquarefree",x);
    1279                 :         25 :       return 0; /* not reached */
    1280                 :            :   }
    1281                 :            : }
    1282                 :            : 
    1283                 :            : /*********************************************************************/
    1284                 :            : /**                                                                 **/
    1285                 :            : /**                    DIGITS / SUM OF DIGITS                       **/
    1286                 :            : /**                                                                 **/
    1287                 :            : /*********************************************************************/
    1288                 :            : 
    1289                 :            : static ulong DS[] ={
    1290                 :            :   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,
    1291                 :            :   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,
    1292                 :            :   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,
    1293                 :            :   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,
    1294                 :            :   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,
    1295                 :            :   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,
    1296                 :            :   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,
    1297                 :            :   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,
    1298                 :            :   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,
    1299                 :            :   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,
    1300                 :            :   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,
    1301                 :            :   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,
    1302                 :            :   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,
    1303                 :            :   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,
    1304                 :            :   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,
    1305                 :            :   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,
    1306                 :            :   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,
    1307                 :            :   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,
    1308                 :            :   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,
    1309                 :            :   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,
    1310                 :            :   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,
    1311                 :            :   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,
    1312                 :            :   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,
    1313                 :            :   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,
    1314                 :            :   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,
    1315                 :            :   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,
    1316                 :            :   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,
    1317                 :            :   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,
    1318                 :            :   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,
    1319                 :            :   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,
    1320                 :            :   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,
    1321                 :            :   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,
    1322                 :            :   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,
    1323                 :            :   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,
    1324                 :            :   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,
    1325                 :            :   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,
    1326                 :            :   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,
    1327                 :            :   25,26,27
    1328                 :            : };
    1329                 :            : 
    1330                 :            : ulong
    1331                 :     253680 : sumdigitsu(ulong n)
    1332                 :            : {
    1333                 :     253680 :   ulong s = 0;
    1334         [ +  + ]:     972745 :   while (n) { s += DS[n % 1000]; n /= 1000; }
    1335                 :     253680 :   return s;
    1336                 :            : }
    1337                 :            : 
    1338                 :            : /* res=array of 9-digits integers, return \sum_{0 <= i < l} sumdigits(res[i]) */
    1339                 :            : static ulong
    1340                 :         10 : sumdigits_block(ulong *res, long l)
    1341                 :            : {
    1342                 :         10 :   long s = sumdigitsu(*--res);
    1343         [ +  + ]:     253670 :   while (--l > 0) s += sumdigitsu(*--res);
    1344                 :         10 :   return s;
    1345                 :            : }
    1346                 :            : 
    1347                 :            : GEN
    1348                 :         25 : sumdigits(GEN n)
    1349                 :            : {
    1350                 :         25 :   pari_sp av = avma;
    1351                 :            :   ulong s, *res;
    1352                 :            :   long l;
    1353                 :            : 
    1354         [ -  + ]:         25 :   if (typ(n) != t_INT) pari_err_TYPE("sumdigits", n);
    1355                 :         25 :   l = lgefint(n);
    1356      [ +  +  + ]:         25 :   switch(l)
    1357                 :            :   {
    1358                 :          5 :     case 2: return gen_0;
    1359                 :         10 :     case 3: return utoipos(sumdigitsu(n[2]));
    1360                 :            :   }
    1361                 :         10 :   res = convi(n, &l);
    1362         [ +  - ]:         10 :   if ((ulong)l < ULONG_MAX / 81)
    1363                 :            :   {
    1364                 :         10 :     s = sumdigits_block(res, l);
    1365                 :         10 :     avma = av; return utoipos(s);
    1366                 :            :   }
    1367                 :            :   else /* Huge. Overflows ulong */
    1368                 :            :   {
    1369                 :          0 :     const long L = (long)(ULONG_MAX / 81);
    1370                 :          0 :     GEN S = gen_0;
    1371         [ #  # ]:          0 :     while (l > L)
    1372                 :            :     {
    1373                 :          0 :       S = addiu(S, sumdigits_block(res, L));
    1374                 :          0 :       res += L; l -= L;
    1375                 :            :     }
    1376         [ #  # ]:          0 :     if (l)
    1377                 :          0 :       S = addiu(S, sumdigits_block(res, l));
    1378                 :         25 :     return gerepileuptoint(av, S);
    1379                 :            :   }
    1380                 :            : }
    1381                 :            : 
    1382                 :            : static GEN
    1383                 :         90 : check_basis(GEN B)
    1384                 :            : {
    1385         [ +  + ]:         90 :   if (!B) return utoipos(10);
    1386         [ -  + ]:         80 :   if (typ(B)!=t_INT) pari_err_TYPE("digits",B);
    1387         [ -  + ]:         80 :   if (cmpis(B,2)<0) pari_err_DOMAIN("digits","B","<",gen_2,B);
    1388                 :         90 :   return B;
    1389                 :            : }
    1390                 :            : 
    1391                 :            : static void
    1392                 :        660 : digits_dac(GEN x, GEN vB, long l, GEN* z)
    1393                 :            : {
    1394                 :            :   GEN q,r;
    1395                 :            :   long m;
    1396         [ +  + ]:        972 :   if (l==1) { *z=x; return; }
    1397                 :        312 :   m=l>>1;
    1398                 :        312 :   q = dvmdii(x, gel(vB,m), &r);
    1399                 :        312 :   digits_dac(q,vB,l-m,z);
    1400                 :        312 :   digits_dac(r,vB,m,z+l-m);
    1401                 :            : }
    1402                 :            : 
    1403                 :            : /* x has l digits in base B, write them to z[0..l-1], vB[i] = B^i */
    1404                 :            : static void
    1405                 :       2180 : digits_dacsmall(GEN x, GEN vB, long l, ulong* z)
    1406                 :            : {
    1407                 :       2180 :   pari_sp av = avma;
    1408                 :            :   GEN q,r;
    1409                 :            :   long m;
    1410         [ +  + ]:       3263 :   if (l==1) { *z=itou(x); return; }
    1411                 :       1083 :   m=l>>1;
    1412                 :       1083 :   q = dvmdii(x, gel(vB,m), &r);
    1413                 :       1083 :   digits_dacsmall(q,vB,l-m,z);
    1414                 :       1083 :   digits_dacsmall(r,vB,m,z+l-m);
    1415                 :       1083 :   avma = av;
    1416                 :            : }
    1417                 :            : 
    1418                 :            : /* set v[i] = 1 iff B^i is needed in the digits_dac algorithm */
    1419                 :            : static void
    1420                 :        490 : set_vexp(GEN v, long l)
    1421                 :            : {
    1422                 :            :   long m;
    1423         [ +  + ]:        710 :   if (v[l]) return;
    1424                 :        220 :   v[l] = 1; m = l>>1;
    1425                 :        220 :   set_vexp(v, m);
    1426                 :        220 :   set_vexp(v, l-m);
    1427                 :            : }
    1428                 :            : /* return all needed B^i for DAC algorithm, for lz digits */
    1429                 :            : static GEN
    1430                 :         50 : get_vB(GEN B, long lz)
    1431                 :            : {
    1432                 :         50 :   GEN vB, vexp = const_vecsmall(lz, 0);
    1433                 :         50 :   long i, l = (lz+1) >> 1;
    1434                 :         50 :   vexp[1] = 1;
    1435                 :         50 :   vexp[2] = 1;
    1436                 :         50 :   set_vexp(vexp, lz);
    1437                 :         50 :   vB = zerovec(lz); /* unneeded entries remain = 0 */
    1438                 :         50 :   gel(vB, 1) = B;
    1439         [ +  + ]:        735 :   for (i = 2; i <= l; i++)
    1440         [ +  + ]:        685 :     if (vexp[i])
    1441                 :            :     {
    1442                 :        220 :       long j = i >> 1;
    1443                 :        220 :       GEN B2j = sqri(gel(vB,j));
    1444         [ +  + ]:        220 :       gel(vB,i) = odd(i)? mulii(B2j, B): B2j;
    1445                 :            :     }
    1446                 :         50 :   return vB;
    1447                 :            : }
    1448                 :            : 
    1449                 :            : GEN
    1450                 :         35 : digits(GEN x, GEN B)
    1451                 :            : {
    1452                 :         35 :   pari_sp av=avma;
    1453                 :            :   long lz;
    1454                 :            :   GEN z, vB;
    1455         [ -  + ]:         35 :   if (typ(x)!=t_INT) pari_err_TYPE("digits",x);
    1456                 :         35 :   B = check_basis(B);
    1457         [ +  + ]:         35 :   if (!signe(x))       {avma = av; return cgetg(1,t_VEC); }
    1458         [ -  + ]:         30 :   if (absi_cmp(x,B)<0) {avma = av; retmkvec(absi(x)); }
    1459         [ +  + ]:         30 :   if (Z_ispow2(B))
    1460                 :            :   {
    1461                 :         10 :     long k = expi(B);
    1462         [ +  + ]:         10 :     if (k == 1) return binaire(x);
    1463         [ -  + ]:          5 :     if (k < BITS_IN_LONG)
    1464                 :            :     {
    1465                 :          0 :       (void)new_chunk(4*(expi(x) + 2)); /* HACK */
    1466                 :          0 :       z = binary_2k_zv(x, k);
    1467                 :          0 :       avma = av; return Flv_to_ZV(z);
    1468                 :            :     }
    1469                 :            :     else
    1470                 :            :     {
    1471                 :          5 :       avma = av; return binary_2k(x, k);
    1472                 :            :     }
    1473                 :            :   }
    1474         [ -  + ]:         20 :   if (signe(x) < 0) x = absi(x);
    1475                 :         20 :   lz = logint(x,B,NULL);
    1476                 :         20 :   vB = get_vB(B, lz);
    1477         [ +  + ]:         20 :   if (lgefint(B)>3)
    1478                 :            :   {
    1479                 :          6 :     z = zerovec(lz);
    1480                 :          6 :     digits_dac(x,vB,lz,(GEN*)(z+1));
    1481                 :          6 :     return gerepilecopy(av,z);
    1482                 :            :   }
    1483                 :            :   else
    1484                 :            :   {
    1485                 :         14 :     (void)new_chunk(3*lz); /* HACK */
    1486                 :         14 :     z = zero_zv(lz);
    1487                 :         14 :     digits_dacsmall(x,vB,lz,(ulong*)(z+1));
    1488                 :         35 :     avma = av; return vecsmall_to_vec(z);
    1489                 :            :   }
    1490                 :            : }
    1491                 :            : 
    1492                 :            : GEN
    1493                 :         75 : sumdigits0(GEN x, GEN B)
    1494                 :            : {
    1495                 :         75 :   pari_sp av = avma;
    1496                 :            :   GEN vB, z;
    1497                 :            :   long lz;
    1498                 :            : 
    1499         [ +  + ]:         75 :   if (!B) return sumdigits(x);
    1500         [ -  + ]:         55 :   if (typ(x) != t_INT) pari_err_TYPE("sumdigits", x);
    1501                 :         55 :   B = check_basis(B);
    1502         [ +  + ]:         55 :   if (Z_ispow2(B))
    1503                 :            :   {
    1504                 :         20 :     long k = expi(B);
    1505         [ +  + ]:         20 :     if (k == 1) { avma = av; return utoi(hammingweight(x)); }
    1506         [ +  + ]:         15 :     if (k < BITS_IN_LONG)
    1507                 :            :     {
    1508                 :         10 :       GEN z = binary_2k_zv(x, k);
    1509         [ -  + ]:         10 :       if (lg(z)-1 > 1<<(BITS_IN_LONG-k)) /* may overflow */
    1510                 :          0 :         return gerepileuptoint(av, ZV_sum(Flv_to_ZV(z)));
    1511                 :         10 :       avma = av; return utoi(zv_sum(z));
    1512                 :            :     }
    1513                 :          5 :     return gerepileuptoint(av, ZV_sum(binary_2k(x, k)));
    1514                 :            :   }
    1515         [ -  + ]:         35 :   if (!signe(x))       { avma = av; return gen_0; }
    1516         [ -  + ]:         35 :   if (absi_cmp(x,B)<0) { avma = av; return absi(x); }
    1517         [ +  + ]:         35 :   if (equaliu(B,10))   { avma = av; return sumdigits(x); }
    1518                 :         30 :   lz = logint(x,B,NULL);
    1519                 :         30 :   vB = get_vB(B, lz);
    1520                 :         30 :   z = zerovec(lz);
    1521                 :         30 :   digits_dac(x,vB,lz,(GEN*)(z+1));
    1522                 :         75 :   return gerepileuptoint(av, ZV_sum(z));
    1523                 :            : }

Generated by: LCOV version 1.9