Optimizations for is_powerfree, squarefree_count. Optimation and bigi…

…nt fix for PP powerfree_count. Add nth_powerfree(n,k) with refining estimate instead of interpolation.
danaj · May 20, 2024 · 95f1922 · 95f1922
1 parent 37d97d9
commit 95f1922
Show file tree

Hide file tree

Showing 11 changed files with 269 additions and 59 deletions.
diff --git a/Changes b/Changes
@@ -79,6 +79,7 @@ Revision history for Perl module Math::Prime::Util
     - powerful_count(n[,k])           count of k-powerful numbers <= n
     - nth_powerful(n[,k])             the nth k-powerful number (default k=2)
     - powerfree_count(n[,k])          count of k-powerfree numbers <= n
+    - nth_powerfree(n[,k])            The nth k-powerfree number (default k=2)
     - powerfree_sum(n[,k])            sum of k-powerfree numbers <= n
     - powerfree_part(n[,k])           remove excess powers from n
     - powerfree_part_sum(n[,k])       sum of k-powerfree parts for 1 to n

diff --git a/XS.xs b/XS.xs
@@ -1801,22 +1801,25 @@ void is_square_free(IN SV* svn)
 void is_powerfree(IN SV* svn, IN int k = 2)
   ALIAS:
     powerfree_count = 1
-    powerfree_sum = 2
-    powerfree_part = 3
-    powerfree_part_sum = 4
+    nth_powerfree = 2
+    powerfree_sum = 3
+    powerfree_part = 4
+    powerfree_part_sum = 5
   PREINIT:
     int status;
     UV n, res;
   PPCODE:
-    status = _validate_and_set(&n, aTHX_ svn, IFLAG_ABS);
+    status = _validate_and_set(&n, aTHX_ svn, (ix==2) ? IFLAG_POS : IFLAG_ABS);
     res = 0;
     if (status == 1) {
       switch (ix) {
         case 0:  res = is_powerfree(n,k);    break;
         case 1:  res = powerfree_count(n,k); break;
-        case 2:  res = powerfree_sum(n,k);   break;
-        case 3:  res = powerfree_part(n,k);  break;
-        case 4:
+        case 2:  if (n == 0 || k < 2)  XSRETURN_UNDEF;
+                 res = nth_powerfree(n,k);   break;
+        case 3:  res = powerfree_sum(n,k);   break;
+        case 4:  res = powerfree_part(n,k);  break;
+        case 5:
         default: res = powerfree_part_sum(n,k);  break;
       }
       if (ix == 0)
@@ -1830,9 +1833,10 @@ void is_powerfree(IN SV* svn, IN int k = 2)
     switch (ix) {
       case  0: _vcallsub_with_gmp(0.00,"is_powerfree"); break;
       case  1: _vcallsub_with_pp("powerfree_count"); break;
-      case  2: _vcallsub_with_pp("powerfree_sum"); break;
-      case  3: _vcallsub_with_pp("powerfree_part"); break;
-      case  4:
+      case  2: _vcallsub_with_pp("nth_powerfree"); break;
+      case  3: _vcallsub_with_pp("powerfree_sum"); break;
+      case  4: _vcallsub_with_pp("powerfree_part"); break;
+      case  5:
       default: _vcallsub_with_pp("powerfree_part_sum"); break;
     }
     return;

diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
@@ -76,7 +76,7 @@ our @EXPORT_OK =
       perfect_power_count_lower perfect_power_count_upper
       nth_perfect_power nth_perfect_power_approx
       nth_perfect_power_lower nth_perfect_power_upper
-      powerfree_count powerfree_sum
+      nth_powerfree powerfree_count powerfree_sum
       powerfree_part powerfree_part_sum
       smooth_count rough_count powersum
       lucky_count lucky_count_approx lucky_count_lower lucky_count_upper
@@ -3122,6 +3122,18 @@ L<OEIS A013928|http://oeis.org/A013928>.
 With C<k == 3> this produces the sequence
 L<OEIS A060431|http://oeis.org/A060431>.
 
+=head2 nth_powerfree
+
+Given a non-negative integer C<n> and an optional non-negative integer C<k>,
+returns the C<n>-th k-powerfree number.
+If C<k> is omitted, C<k=2> is used.
+Returns undef if C<k> is less than 2 or C<n=0>.  Returns 1 for C<n=1>.
+
+With C<k == 2> this produces the sequence
+L<OEIS A005117|http://oeis.org/A005117>.
+With C<k == 3> this produces the sequence
+L<OEIS A004709|http://oeis.org/A004709>.
+
 =head2 powerfree_sum
 
 Given an integer C<n> and an optional non-negative integer C<k>, returns
@@ -4974,7 +4986,7 @@ integer C<n>, efficiently compute C<lucasu(P,Q,k) mod |n|>.
 
 This corresponds to gmpy2's C<lucasu_mod> function.
 
-When C<(P,Q) = (1,-1)> this returns the modular fibonacci sequence.  This
+When C<(P,Q) = (1,-1)> this returns the modular Fibonacci sequence.  This
 corresponds to Sidef's C<fibmod> function.
 
 =head2 lucasvmod

diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
@@ -130,6 +130,20 @@ BEGIN {
 *Mforprimes = \&Math::Prime::Util::forprimes;
 *MLi = \&Math::Prime::Util::LogarithmicIntegral;
 
+if (defined $Math::Prime::Util::GMP::VERSION && $Math::Prime::Util::GMP::VERSION >= 0.53) {
+  *Saddint = \&Math::Prime::Util::GMP::addint;
+  *Ssubint = \&Math::Prime::Util::GMP::subint;
+  *Smulint = \&Math::Prime::Util::GMP::mulint;
+  *Sdivint = \&Math::Prime::Util::GMP::divint;
+  *Spowint = \&Math::Prime::Util::GMP::powint;
+} else {
+  *Saddint = \&Math::Prime::Util::addint;
+  *Ssubint = \&Math::Prime::Util::subint;
+  *Smulint = \&Math::Prime::Util::mulint;
+  *Sdivint = \&Math::Prime::Util::divint;
+  *Spowint = \&Math::Prime::Util::powint;
+}
+
 my $_precalc_size = 0;
 sub prime_precalc {
   my($n) = @_;
@@ -1739,29 +1753,128 @@ sub powerfree_count {
 
   return (($n >= 1) ? 1 : 0)  if $k < 2 || $n <= 1;
 
-  my $count = $n;
+  my $count = 0;
   my $nk = Mrootint($n, $k);
 
-  if ($nk < 100 || $nk > 1e8) {
-    Math::Prime::Util::forsquarefree(
-      sub {
-        $count += ((scalar(@_) & 1) ? -1 : 1) * Mdivint($n, Mpowint($_, $k));
-      },
-      2, $nk
-    );
-  } else {
-    my @mu = (0, Mmoebius(1, $nk));
+  # If we can do everything native, do that.
+  if ($n < SINTMAX && $nk < 20000) {
+    use integer;
+    my @mu = Mmoebius(0, $nk);
     foreach my $i (2 .. $nk) {
-      next if $mu[$i] == 0;
-      $count += $mu[$i] * Mdivint($n, Mpowint($i, $k));
+      $count += $mu[$i] * $n/($i**$k) if $mu[$i];
+    }
+    return Maddint($count,$n);
+  } elsif ($n < SINTMAX && $nk < 1e8) {
+    # Split out the trailing n/i^k = 1, saves memory and time if large enough.
+    use integer;
+    my $L1 = Mrootint($n/2,$k);
+    my @mu = Mmoebius(0, $L1);
+    foreach my $i (2 .. $L1) {
+      $count += $mu[$i] * $n/($i**$k) if $mu[$i];
+    }
+    #@mu = Mmoebius($L1+1, $nk);   my $c1 = 0;   $c1 += $_ for @mu;
+    my $c1 = Math::Prime::Util::mertens($nk) - Math::Prime::Util::mertens($L1);
+    return Mvecsum($count,$c1,$n);
+  }
+
+  # Simple way.  All the bigint math kills performance.
+  # Math::Prime::Util::forsquarefree(
+  #   sub {
+  #     my $t = Mdivint($n, Mpowint($_, $k));
+  #     $count = (scalar(@_) & 1) ? Msubint($count,$t) : Maddint($count,$t);
+  #   },
+  #   2, $nk
+  # );
+
+  # Optimization 1:  pull out all the ranges at the end with small constant
+  #                  multiplications.
+  # Optimization 2:  Use GMP basic arithmetic functions if possible, saving
+  #                  all the bigint object overhead.  Can be 10x faster.
+
+  my $A = Msqrtint($nk);
+  my @L = (0, $nk, map { Mrootint(Mdivint($n,$_),$k) } 2..$A);
+  my @C;
+
+  Math::Prime::Util::forsquarefree(
+    sub {
+      $count = (scalar(@_) & 1)
+             ? Ssubint($count, Sdivint($n, Spowint($_, $k)))
+             : Saddint($count, Sdivint($n, Spowint($_, $k)));
+    },
+    2, $L[$A]
+  );
+  for my $i (2 .. $A) {
+    my($c, $lo, $hi) = (0, $L[$i], $L[$i-1]);
+    if ($i < 15) {
+      $c = Math::Prime::Util::mertens($hi) - Math::Prime::Util::mertens($lo);
+    } else {
+      $c += $_ for Mmoebius( Maddint($lo,1), $hi );
     }
+    push @C, $c * ($i-1);
+    @C = (Mvecsum(@C)) if scalar(@C) > 100000;  # Save/restrict memory.
   }
-  $count;
+  my $ctot = Mvecsum(@C); # Can typically be done in native math.
+  Mvecsum($count, $n, $ctot);
+}
+
+sub nth_powerfree {
+  my($n, $k) = @_;
+  _validate_positive_integer($n);
+  if (defined $k) { _validate_positive_integer($k); }
+  else            { $k = 2; }
+
+  return undef if $n == 0 || $k < 2;
+  return $n if $n < 4;
+
+  # 1. zm is the zeta multiplier, qk is the expected value.
+  my($zm, $qk);
+  if ($n <= 2**52) {
+    $zm = ($k == 2) ? 1.644934066848226 : 1.0 + RiemannZeta($k);
+    $qk = $zm * "$n";
+    $qk = int("$qk");
+  } else {
+    do { require Math::BigFloat; Math::BigFloat->import(); }
+      if !defined $Math::BigFloat::VERSION;
+    require Math::Prime::Util::ZetaBigFloat;
+    my $acc = length("$n")+10;
+    my $bk = Math::BigFloat->new($k);  $bk->accuracy($acc);
+    $zm = Math::Prime::Util::ZetaBigFloat::RiemannZeta($bk)->badd(1);
+    $qk = $zm->copy->bmul("$n");
+    $qk = Math::BigInt->new($qk->bfloor->bstr);
+  }
+  $zm = $zm->numify() if ref($zm);
+
+  my($count, $diff);
+  # In practice this converges very rapidly, usually needing only one iteration.
+  for (1 .. 10) {
+    # 2. Get the actual count at qk and the difference from our goal.
+    $count = Math::Prime::Util::powerfree_count($qk,$k);
+    $diff = ($count >= $n)  ?  $count-$n  :  $n-$count;
+    # print "qk $qk  count $count  diff $diff\n";
+
+    # 3. If not close, update the estimate using the expected density zm.
+    last if $diff <= 300;   # Threshold could be improved.
+    if ($count > $n) {
+      $qk -= int("$diff" * $zm);
+    } else {
+      $qk += int("$diff" * $zm);
+    }
+  }
+
+  # 4. Make sure we're on a powerfree number.
+  $qk-- while !Math::Prime::Util::is_powerfree($qk,$k);
+
+  # 5. Walk forward or backward to next/prev powerfree number.
+  my $adder = ($count < $n) ? 1 : -1;
+  while ($count != $n) {
+    do { $qk += $adder; } while !Math::Prime::Util::is_powerfree($qk,$k);
+    $count += $adder;
+  }
+  $qk;
 }
 
 sub powerfree_sum {
   my($n, $k) = @_;
-  $n = -$n if defined $n && $n < 0;
   _validate_positive_integer($n);
   if (defined $k) { _validate_positive_integer($k); }
   else            { $k = 2; }
@@ -7815,7 +7928,7 @@ sub is_frobenius_khashin_pseudoprime {
       $k = kronecker($c, $n);
     } while $k == 1;
   }
-  return 0 if $k == 0 || ($k == 2 && !($n % 3));;
+  return 0 if $k == 0 || ($k == 2 && !($n % 3));
 
   my $ea = ($k == 2) ? 2 : 1;
   my($ra,$rb,$a,$b,$d) = ($ea,1,$ea,1,$n-1);

diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
@@ -165,6 +165,7 @@ sub entropy_bytes {
 *perfect_power_count = \&Math::Prime::Util::PP::perfect_power_count;
 *is_powerfree = \&Math::Prime::Util::PP::is_powerfree;
 *powerfree_count = \&Math::Prime::Util::PP::powerfree_count;
+*nth_powerfree = \&Math::Prime::Util::PP::nth_powerfree;
 *powerfree_sum = \&Math::Prime::Util::PP::powerfree_sum;
 *powerfree_part = \&Math::Prime::Util::PP::powerfree_part;
 *powerfree_part_sum = \&Math::Prime::Util::PP::powerfree_part_sum;

diff --git a/lib/ntheory.pm b/lib/ntheory.pm
@@ -380,6 +380,7 @@ Tags:
   smooth_count(n,k)                   count of k-smooth numbers <= n
   rough_count(n,k)                    count of k-rough numbers <= n
   powerfree_count(n[,k])              count of k-powerfree numbers <= n
+  nth_powerfree(n[,k])                the nth k-powerfree number
   powerfree_sum(n[,k])                sum of k-powerfree numbers <= n
   powerfree_part(n[,k])               remove excess powers so n is k-free
   powerfree_part_sum(n[,k])           sum of k-powerfree parts for 1 to n