[libmath-prime-util-perl] 40/55: Extend binomial to negative arguments

[Partha P. Mukherjee]

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

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

commit 8f60e4120025c38686c43c760baed82a845be542
Author: Dana Jacobsen <dana at acm.org>
Date:   Tue May 13 19:45:50 2014 -0700

    Extend binomial to negative arguments
---
 XS.xs                       | 22 ++++++++++++++++------
 lib/Math/Prime/Util.pm      | 10 ++++++++--
 lib/Math/Prime/Util/PP.pm   | 19 ++++++++++++++++---
 lib/Math/Prime/Util/PPFE.pm |  4 ++--
 t/19-moebius.t              |  8 ++++++++
 5 files changed, 50 insertions(+), 13 deletions(-)

diff --git a/XS.xs b/XS.xs
index 3cdacc0..4c589dc 100644
--- a/XS.xs
+++ b/XS.xs
@@ -766,18 +766,27 @@ znorder(IN SV* sva, IN SV* svn)
   PREINIT:
     int astatus, nstatus;
   PPCODE:
-    astatus = _validate_int(aTHX_ sva, 0);
-    nstatus = _validate_int(aTHX_ svn, 0);
-    if (astatus == 1 && nstatus == 1) {
+    astatus = _validate_int(aTHX_ sva, (ix==1) ? 2 : 0);
+    nstatus = _validate_int(aTHX_ svn, (ix==1) ? 2 : 0);
+    if (astatus != 0 && nstatus != 0) {
       UV a = my_svuv(sva);
       UV n = my_svuv(svn);
       UV ret;
       switch (ix) {
         case 0:  ret = znorder(a, n);
                  break;
-        case 1:  ret = binomial(a, n);
-                 if (ret == 0 && n <= a)
+        case 1:  if (nstatus == -1) {
+                   ret = 0;
+                 } else if (astatus == -1) {
+                   ret = binomial( -my_sviv(sva)+n-1, n );
+                   if (ret > 0 && ret <= (UV)IV_MAX)
+                     XSRETURN_IV( (IV)ret * ((n&1) ? -1 : 1) );
                    goto overflow;
+                 } else {
+                   ret = binomial(a, n);
+                   if (ret == 0 && n <= a)
+                     goto overflow;
+                 }
                  break;
         case 2:  ret = jordan_totient(a, n);
                  if (ret == 0 && n > 1)
@@ -844,7 +853,8 @@ kronecker(IN SV* sva, IN SV* svb)
       } else if (ix == 1) {
         UV n = (astatus == -1) ? (UV)(-(my_sviv(sva))) : my_svuv(sva);
         UV k = (bstatus == -1) ? (UV)(-(my_sviv(svb))) : my_svuv(svb);
-        XSRETURN_UV( valuation(n, k) );
+        /* valuation of 0-2 is very common, so return a constant if possible */
+        RETURN_NPARITY( valuation(n, k) );
       } else {
         UV a, n, ret = 0;
         n = (bstatus != -1) ? my_svuv(svb) : (UV)(-(my_sviv(svb)));
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 77945ff..660b81c 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -2201,11 +2201,17 @@ C<mpz_legendre(a,n)> functions.
 
 =head2 binomial
 
-Given non-negative arguments C<n> and C<k>, returns the binomial
-coefficient C<n*(n-1)*...*(n-k+1)/k!>, also known as the choose function.
+Given integer arguments C<n> and C<k>, returns the binomial coefficient
+C<n*(n-1)*...*(n-k+1)/k!>, also known as the choose function.  Negative
+arguments use the L<Kronenburg extensions|http://arxiv.org/abs/1105.3689/>.
 This corresponds to Pari's C<binomial(n,k)> function, Mathematica's
 C<Binomial[n,k]> function, and GMP's C<mpz_bin_ui> function.
 
+The rules for negative C<n> follow the same rules as the functions mentioned,
+which does not match Math::BigInt.  Additionally, Pari does not use the
+extension for C<n E<lt>0 ; k E<lt>= n>, while GMP's API doesn't allow this
+case at all.
+
 
 =head2 znorder
 
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index d795a7e..4e5d244 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1787,9 +1787,22 @@ sub kronecker {
 }
 sub binomial {
   my($n, $k) = @_;
-  my $r = Math::BigInt->new("$n")->bnok("$k");
-  $r = _bigint_to_int($r) if $r->bacmp(''.~0) <= 0;
-  return $r;
+  if ($k <= 0) { return ($k == 0) ? 1 : 0; }
+  my $r;
+  if ($n >= 0) {
+    $r = Math::BigInt->new("$n")->bnok("$k");
+    $r = _bigint_to_int($r) if $r->bacmp(''.~0) <= 0;
+  } else {
+    # Math::BigInt is incorrect for negative n
+    $r = Math::BigInt->new( ''.(-$n+$k-1) )->bnok("$k");
+    if ($k & 1) {
+      $r->bneg;
+      $r = _bigint_to_int($r) if $r->bacmp(''.(~0>>1)) <= 0;
+    } else {
+      $r = _bigint_to_int($r) if $r->bacmp(''.~0) <= 0;
+    }
+  }
+  $r;
 }
 
 sub _is_perfect_square {
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index 5b0f901..cf0ecb9 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -212,8 +212,8 @@ sub kronecker {
 
 sub binomial {
   my($n, $k) = @_;
-  _validate_positive_integer($n);
-  _validate_positive_integer($k);
+  _validate_integer($n);
+  _validate_integer($k);
   return Math::Prime::Util::PP::binomial($n, $k);
 }
 
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 5be67ca..b9409b0 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -410,6 +410,14 @@ my @binomials = (
  [ 67,31, "11923179284862717872" ], # ...and this
  [ 228,12, "30689926618143230620" ],# But the result of this is too big.
  [ 177,78, "3314450882216440395106465322941753788648564665022000" ],
+ [ -10,5, -2002 ],
+ [ -11,22, 64512240 ],
+ [ -12,23, -286097760 ],
+ [ -12,-23, 0 ],
+ [  12,-23, 0 ],
+ [  12,-12, 0 ],
+ [ -12,0, 1 ],
+ [  0,-1, 0 ],
 );
 
 # These are slow with XS, and *really* slow with PP.

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