divisors(0) returns 0/[]. divisors(n,k) limits to divisors <= k. Fix …

…kronecker overflow with 64-bit signed.

Changes

@@ -11,6 +11,10 @@ Revision history for Perl module Math::Prime::Util
     - is_strong_pseudoprime passes for bases equal to 0 mod n.  This matches
       what the GMP code always did, and means primes return 1 for any base.
 
+    - divisor_sum(0) returns 0.  divisors(0) returns an empty list.
+
+    - divisors takes an optional second argument to prune the returned list.
+
     [ADDED]
 
     - powersum(n,k)                   sum of k-th powers of positive ints <= n
@@ -187,6 +191,8 @@ Revision history for Perl module Math::Prime::Util
     - is_pseudoprime and is_euler_pseudoprime have consistent XS, PP, and GMP
       behaviour.  Specifically, single digit inputs and bases a multiple of n.
 
+    - kronecker in XS with 63-bit signed a and 64-bit unsigned n fixed.
+
     [FUNCTIONALITY AND PERFORMANCE]
 
     - LMO prime count is 1.2-1.5x faster.

MANIFEST

@@ -265,6 +265,7 @@ xt/check-nth-bounds.pl
 xt/chinese.pl
 xt/create-pc-tables.pl
 xt/division.pl
+xt/divisors.pl
 xt/divtest.pl
 xt/foralmostprimes.pl
 xt/moebius-mertens.pl
@@ -274,6 +275,7 @@ xt/primality-aks.pl
 xt/primality-proofs.pl
 xt/small-is-next-prev.pl
 xt/factor-holf.pl
+xt/kronecker.pl
 xt/legendre_phi.t
 xt/lucasseq.pl
 xt/lucky.c

XS.xs

@@ -2340,79 +2340,107 @@ _pidigits(IN int digits)
     XPUSHs(sv_2mortal(newSVpvn(out, digits+1)));
     Safefree(out);
 
+void inverse_totient(IN SV* svn)
+  PREINIT:
+    U32 gimme_v;
+    int status, it_overflow;
+    UV i, n, ntotients;
+  PPCODE:
+    gimme_v = GIMME_V;
+    status = _validate_and_set(&n, aTHX_ svn, IFLAG_POS);
+    it_overflow = (status == 1 && gimme_v == G_ARRAY && n > UV_MAX/7.5);
+    if (status == 1 && !it_overflow) {
+      if (gimme_v == G_SCALAR) {
+        XSRETURN_UV( inverse_totient_count(n) );
+      } else if (gimme_v == G_ARRAY) {
+        UV* tots = inverse_totient_list(&ntotients, n);
+        EXTEND(SP, (IV)ntotients);
+        for (i = 0; i < ntotients; i++)
+          PUSHs(sv_2mortal(newSVuv( tots[i] )));
+        Safefree(tots);
+      }
+    } else {
+      _vcallsubn(aTHX_ gimme_v, VCALL_GMP|VCALL_PP, "inverse_totient", 1, 0);
+      return; /* skip implicit PUTBACK */
+    }
+
 void
 factor(IN SV* svn)
   ALIAS:
     factor_exp = 1
-    divisors = 2
-    inverse_totient = 3
   PREINIT:
     U32 gimme_v;
-    int status, i, nfactors, it_overflow;
+    int status, i, nfactors;
     UV n;
   PPCODE:
     gimme_v = GIMME_V;
     status = _validate_and_set(&n, aTHX_ svn, IFLAG_POS);
-    it_overflow = (status == 1 && ix==3 && gimme_v == G_ARRAY && n > UV_MAX/7.5 );
-    if (status == 1 && !it_overflow) {
+    if (status == 1) {
       UV factors[MPU_MAX_FACTORS+1];
       UV exponents[MPU_MAX_FACTORS+1];
       if (gimme_v == G_SCALAR) {
-        UV res;
-        switch (ix) {
-          case 0:  res = factor(n, factors);        break;
-          case 1:  res = factor_exp(n, factors, 0); break;
-          case 2:  res = divisor_sum(n, 0);         break;
-          default: res = inverse_totient_count(n);  break;
-        }
-        PUSHs(sv_2mortal(newSVuv( res )));
+        UV res = (ix == 0)  ?  factor(n, factors)  :  factor_exp(n, factors, 0);
+        XSRETURN_UV(res);
       } else if (gimme_v == G_ARRAY) {
-        switch (ix) {
-          case 0:  nfactors = factor(n, factors);
-                   EXTEND(SP, nfactors);
-                   for (i = 0; i < nfactors; i++)
-                     PUSHs(sv_2mortal(newSVuv( factors[i] )));
-                   break;
-          case 1:  nfactors = factor_exp(n, factors, exponents);
-                   /* if (n == 1)  XSRETURN_EMPTY; */
-                   EXTEND(SP, nfactors);
-                   for (i = 0; i < nfactors; i++) {
-                     AV* av = newAV();
-                     av_push(av, newSVuv(factors[i]));
-                     av_push(av, newSVuv(exponents[i]));
-                     PUSHs( sv_2mortal(newRV_noinc( (SV*) av )) );
-                   }
-                   break;
-          case 2: {
-                     UV ndivisors;
-                     UV* divs = _divisor_list(n, &ndivisors);
-                     EXTEND(SP, (IV)ndivisors);
-                     for (i = 0; (UV)i < ndivisors; i++)
-                       PUSHs(sv_2mortal(newSVuv( divs[i] )));
-                     Safefree(divs);
-                   }
-                   break;
-          default: {
-                     UV ntotients;
-                     UV* tots = inverse_totient_list(&ntotients, n);
-                     EXTEND(SP, (IV)ntotients);
-                     for (i = 0; (UV)i < ntotients; i++)
-                       PUSHs(sv_2mortal(newSVuv( tots[i] )));
-                     Safefree(tots);
-                   }
-                   break;
+        if (ix == 0) {
+          nfactors = factor(n, factors);
+          EXTEND(SP, nfactors);
+          for (i = 0; i < nfactors; i++)
+            PUSHs(sv_2mortal(newSVuv( factors[i] )));
+        } else {
+          nfactors = factor_exp(n, factors, exponents);
+          /* if (n == 1)  XSRETURN_EMPTY; */
+          EXTEND(SP, nfactors);
+          for (i = 0; i < nfactors; i++) {
+            AV* av = newAV();
+            av_push(av, newSVuv(factors[i]));
+            av_push(av, newSVuv(exponents[i]));
+            PUSHs( sv_2mortal(newRV_noinc( (SV*) av )) );
+          }
         }
       }
     } else {
-      switch (ix) {
-        case 0:  _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor", 1, 0);     break;
-        case 1:  _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor_exp", 1, 0); break;
-        case 2:  _vcallsubn(aTHX_ gimme_v, VCALL_GMP|VCALL_PP, "divisors", 1, 0);   break;
-        default: _vcallsubn(aTHX_ gimme_v, VCALL_GMP|VCALL_PP, "inverse_totient", 1, 0);   break;
-      }
+      if (ix == 0)
+        _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor", 1, 0);
+      else
+        _vcallsubn(aTHX_ gimme_v, VCALL_ROOT, "_generic_factor_exp", 1, 0);
       return; /* skip implicit PUTBACK */
     }
 
+void divisors(IN SV* svn, IN SV* svk = 0)
+  PREINIT:
+    U32 gimme_v;
+    int status;
+    UV n, k, i, ndivisors, *divs;
+  PPCODE:
+    gimme_v = GIMME_V;
+    status = _validate_and_set(&n, aTHX_ svn, IFLAG_POS);
+    k = n;
+    if (status == 1 && svk != 0) {
+      status = _validate_and_set(&k, aTHX_ svk, IFLAG_POS);
+      if (k > n)  k = n;
+    }
+    if (status != 1) {
+      _vcallsubn(aTHX_ gimme_v, VCALL_GMP|VCALL_PP, "divisors", items, 53);
+      return; /* skip implicit PUTBACK */
+    }
+    if (GIMME_V == G_VOID) {
+      /* Nothing */
+    } else if (GIMME_V == G_SCALAR && k >= n) {
+      ndivisors = divisor_sum(n, 0);
+      PUSHs(sv_2mortal(newSVuv( ndivisors )));
+    } else {
+      divs = divisor_list(n, &ndivisors, k);
+      if (GIMME_V == G_SCALAR) {
+        PUSHs(sv_2mortal(newSVuv( ndivisors )));
+      } else {
+        EXTEND(SP, (IV)ndivisors);
+        for (i = 0; i < ndivisors; i++)
+          PUSHs(sv_2mortal(newSVuv( divs[i] )));
+      }
+      Safefree(divs);
+    }
+
 void
 trial_factor(IN UV n, ...)
   ALIAS:
@@ -2475,19 +2503,20 @@ divisor_sum(IN SV* svn, ...)
   PREINIT:
     UV n, k, sigma;
   PPCODE:
-    sigma = 0;
-    if (items == 1) {
-      if ( _validate_and_set(&n, aTHX_ svn, IFLAG_POS))
-        sigma = divisor_sum(n, 1);
+    if (items == 1 && _validate_and_set(&n, aTHX_ svn, IFLAG_POS)) {
+      sigma = divisor_sum(n, 1);
+      if (n <= 1 || sigma != 0)
+        XSRETURN_UV(sigma);
     } else {
       SV* svk = ST(1);
       if ( (!SvROK(svk) || (SvROK(svk) && SvTYPE(SvRV(svk)) != SVt_PVCV)) &&
            _validate_and_set(&n, aTHX_ svn, IFLAG_POS) &&
-           _validate_and_set(&k, aTHX_ svk, IFLAG_POS) )
+           _validate_and_set(&k, aTHX_ svk, IFLAG_POS) ) {
         sigma = divisor_sum(n, k);
+        if (n <= 1 || sigma != 0)
+          XSRETURN_UV(sigma);
+      }
     }
-    if (sigma != 0)   /* sigma 0 means overflow */
-      XSRETURN_UV(sigma);
     _vcallsub_with_pp("divisor_sum");
     return; /* skip implicit PUTBACK */
 
@@ -4337,7 +4366,7 @@ fordivisors (SV* block, IN SV* svn)
       return;
     }
 
-    divs = _divisor_list(n, &ndivisors);
+    divs = divisor_list(n, &ndivisors, UV_MAX);
 
     START_FORCOUNT;
     SAVESPTR(GvSV(PL_defgv));

factor.c

@@ -367,41 +367,51 @@ int trial_factor(UV n, UV *factors, UV f, UV last)
 }
 
 
-static void _divisors_from_factors(UV nfactors, UV* fp, UV* fe, UV* res) {
-  UV s, count = 1;
+static UV _divisors_from_factors(UV nfactors, UV* fp, UV* fe, UV* res, UV k) {
+  UV s, t, count = 1;
 
   res[0] = 1;
   for (s = 0; s < nfactors; s++) {
     UV i, j, scount = count, p = fp[s], e = fe[s], mult = 1;
     for (j = 0; j < e; j++) {
       mult *= p;
-      for (i = 0; i < scount; i++)
-        res[count++] = res[i] * mult;
+      for (i = 0; i < scount; i++) {
+        t = res[i] * mult;
+        if (t <= k)
+          res[count++] = t;
+      }
     }
   }
+  return count;
 }
 
-UV* _divisor_list(UV n, UV *num_divisors)
+UV* divisor_list(UV n, UV *num_divisors, UV maxd)
 {
   UV factors[MPU_MAX_FACTORS+1];
   UV exponents[MPU_MAX_FACTORS+1];
   UV* divs;
   int i, nfactors, ndivisors;
 
-  if (n <= 1) {
-    New(0, divs, 2, UV);
-    if (n == 0) {  divs[0] = 0;  divs[1] = 1;  *num_divisors = 2;  }
-    if (n == 1) {  divs[0] = 1;                *num_divisors = 1;  }
+  if (n == 0 || maxd == 0) {
+    *num_divisors = 0;
+    return 0;
+  } else if (n == 1 || maxd == 1) {
+    New(0, divs, 1, UV);
+    divs[0] = 1;
+    *num_divisors = 1;
     return divs;
   }
+
+  if (maxd > n) maxd = n;
+
   /* Factor and convert to factor/exponent pair */
   nfactors = factor_exp(n, factors, exponents);
   /* Calculate number of divisors, allocate space, fill with divisors */
   ndivisors = exponents[0] + 1;
   for (i = 1; i < nfactors; i++)
     ndivisors *= (exponents[i] + 1);
   New(0, divs, ndivisors, UV);
-  _divisors_from_factors(nfactors, factors, exponents, divs);
+  ndivisors = _divisors_from_factors(nfactors, factors, exponents, divs, maxd);
   /* Sort divisors (numeric ascending) */
   qsort(divs, ndivisors, sizeof(UV), _numcmp);
   /* Return number of divisors and list */
@@ -430,8 +440,8 @@ UV divisor_sum(UV n, UV k)
   UV product = 1;
 
   if (k > 11 || (k > 0 && n >= sigma_overflow[k-1])) return 0;
-  if (n <= 1)                               /* n=0  divisors are [0,1] */
-    return (n == 1) ? 1 : (k == 0) ? 2 : 1; /* n=1  divisors are [1]   */
+  /*   divisors(0) = []   divisors(1) = [1]  */
+  if (n <= 1)  return n;
   nfac = factor(n,factors);
   if (k == 0) {
     for (i = 0; i < nfac; i++) {

factor.h

@@ -21,7 +21,7 @@ extern int squfof_factor(UV n, UV *factors, UV rounds);
 extern int lehman_factor(UV n, UV *factors, int dotrial);
 extern int cheb_factor(UV n, UV *factors, UV B, UV initx);
 
-extern UV* _divisor_list(UV n, UV *num_divisors);
+extern UV* divisor_list(UV n, UV *num_divisors, UV maxd);
 
 extern int prime_omega(UV n);     /* number of distinct prime factors */
 extern int prime_bigomega(UV n);  /* number of prime factors w/ multiplicity */

lib/Math/Prime/Util.pm

@@ -4155,7 +4155,7 @@ C<chebyshev_theta(n^(1/k))>.
   say "Number of divisors: sigma_0($n) = ", divisor_sum($n, 0);
 
 Given a single non-negative integer C<n>, returns the sum of the
-divisors of C<n>, including 1 and itself.
+divisors of C<n>, including 1 and itself.  We return 0 for C<n=0>.
 
 An optional second non-negative integer C<k> may be given, indicating
 the sum should use the C<k-th> powers of the divisors.
@@ -5695,6 +5695,13 @@ C<DivisorSigma[0,n]> functions.
 
 Also see the L</for_divisors> functions for looping over the divisors.
 
+When C<n=0> we return the empty set (zero in scalar context).
+
+An optional second positive integer argument C<k> indicates that the results
+should not include any value larger than C<k>.  This is especially useful
+when the number has thousands of divisors and we may only be interested in
+the small ones.
+
 
 =head2 trial_factor