[libmath-prime-util-perl] 06/35: Update partitions function

[Partha P. Mukherjee]

This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.33
in repository libmath-prime-util-perl.

commit 82894f097d2d4c4beb46633e4f9c32cd938bc5b5
Author: Dana Jacobsen <dana at acm.org>
Date:   Mon Oct 21 08:22:26 2013 -0700

    Update partitions function
---
 Changes                |  8 +++----
 lib/Math/Prime/Util.pm | 62 ++++++++++++++++++++++++++++----------------------
 2 files changed, 39 insertions(+), 31 deletions(-)

diff --git a/Changes b/Changes
index a9d976a..0b6784a 100644
--- a/Changes
+++ b/Changes
@@ -4,10 +4,10 @@ Revision history for Perl module Math::Prime::Util
 
     [API Changes]
 
-      - all_factors now includes 1 and n, making it identical to Pari's
-        divisors(n) function, but no longer identical to Math::Factor::XS's
-        factors(n) function.  This change allows consistency between
-        divisor_sum(n,0) and scalar all_factors(n).
+    - all_factors now includes 1 and n, making it identical to Pari's
+      divisors(n) function, but no longer identical to Math::Factor::XS's
+      factors(n) function.  This change allows consistency between
+      divisor_sum(n,0) and scalar all_factors(n).
 
     [ADDED]
 
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 8d128a3..d8301a3 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -1660,35 +1660,43 @@ sub liouville {
   return $l;
 }
 
-{
-  my @pent = (0, 1);
-  my @part = (1);
-  # TODO: make this faster.  Pari is almost instant for n < 1000000.
-  # We do seem to be faster than the stackoverflow or praxis solutions.
-  sub partitions {
-    my($n) = @_;
-    return 1 if defined $n && $n <= 0;
-    _validate_num($n) || _validate_positive_integer($n);
-    return abs($part[$n]) if defined $part[$n];
-    if ($n >= 128 && !defined $Math::BigInt::VERSION) {
-      eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
-      or do { croak "Cannot load Math::BigInt"; };
-    }
-    my $d = int(sqrt($n+1));
-    foreach my $i ( int((scalar @pent)/2) .. $d ) {
-      push @pent, int(($i*(3*$i+1))/2), int((($i+1)*(3*$i+2))/2);
-    }
-    foreach my $j (scalar @part .. $n) {
-      my $psum = ($n < 128) ? 0 : Math::BigInt->bzero;
-      foreach my $k (1 .. $n) {
-        last if $pent[$k] > $j;
-        my $gk = $part[ $j - $pent[$k] ];
-        $psum += (($k+1) & 2)  ?  $gk  :  -$gk;
-      }
-      $part[$j] = $psum;
+# This is faster than most regular implementations I've seen, e.g. the ones
+# on stackoverflow or praxis.  The same algorithm in C+GMP is 300x faster.
+# It's much slower than Pari's floating point Rademacher algorithm.
+# See 2011+ FLINT and Fredrik Johansson's work for state of the art.
+#   Perl-comb   partitions(10^5)  ~ 300 seconds
+#   GMP-comb    partitions(10^6)  ~ 120 seconds  ( ~300x faster than Perl-comb)
+#   Pari        partitions(10^8)  ~ 100 seconds  (~1000x faster than GMP-comb)
+#   Bober       partitions(10^9)  ~  20 seconds  (~  50x faster than Pari)
+#   Johansson   partitions(10^12) ~  10 seconds  (>1000x faster than Pari)
+sub partitions {
+  my($n) = @_;
+  return 1 if defined $n && $n <= 0;
+  _validate_num($n) || _validate_positive_integer($n);
+  if (!defined $Math::BigInt::VERSION) {
+    eval { require Math::BigInt; Math::BigInt->import(try=>'GMP,Pari'); 1; }
+    or do { croak "Cannot load Math::BigInt"; };
+  }
+  if ($_HAVE_GMP && defined &Math::Prime::Util::GMP::partitions) {
+    return Math::BigInt->new( '' . Math::Prime::Util::GMP::partitions($n) );
+  }
+  my $d = int(sqrt($n+1));
+  my @pent = (1);
+  foreach my $i ( 1 .. $d ) {
+    push @pent, int(($i*(3*$i+1))/2), int((($i+1)*(3*$i+2))/2);
+  }
+  my @part = (Math::BigInt->bone);
+  foreach my $j (scalar @part .. $n) {
+    my $psum = Math::BigInt->bzero;
+    my $k = 1;
+    foreach my $p (@pent) {
+      last if $p > $j;
+      if ((++$k) & 2) { $psum->badd( $part[ $j - $p ] ); }
+      else            { $psum->bsub( $part[ $j - $p ] ); }
     }
-    return $part[$n];
+    $part[$j] = $psum;
   }
+  return $part[$n];
 }
 
 sub chebyshev_theta {

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git