[libmath-prime-util-perl] 35/55: Add binomial(n,k)
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:41 UTC 2015
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 935293c49e37b6c75e23f94257ca4ed6bc5f3b05
Author: Dana Jacobsen <dana at acm.org>
Date: Mon May 12 16:49:21 2014 -0700
Add binomial(n,k)
---
Changes | 1 +
TODO | 2 +-
XS.xs | 18 ++++++++++++------
lib/Math/Prime/Util.pm | 12 ++++++++++--
lib/Math/Prime/Util/PP.pm | 6 ++++++
lib/Math/Prime/Util/PPFE.pm | 7 +++++++
t/19-moebius.t | 22 +++++++++++++++++++++-
util.c | 20 ++++++++++++++++++++
util.h | 1 +
9 files changed, 79 insertions(+), 10 deletions(-)
diff --git a/Changes b/Changes
index d1e2367..a3d99d5 100644
--- a/Changes
+++ b/Changes
@@ -8,6 +8,7 @@ Revision history for Perl module Math::Prime::Util
- invmod(a,n) inverse of a modulo n
- forpart { ... } n[,{...}] loop over partitions (Pari 2.6.x)
- vecsum(...) sum list of integers
+ - binomial(n,k) binomial coefficient
[FUNCTIONALITY AND PERFORMANCE]
diff --git a/TODO b/TODO
index 2f0bb1c..71bd5c1 100644
--- a/TODO
+++ b/TODO
@@ -75,6 +75,6 @@
- XS alias forprimes, forcomposites, forpart
-- More Pari: vecsum, parforprime, binomial
+- More Pari: parforprime, vecmin, vecmax
- is_power should take a third argument that gets set to the root
diff --git a/XS.xs b/XS.xs
index b4ec692..41bcd91 100644
--- a/XS.xs
+++ b/XS.xs
@@ -760,8 +760,9 @@ divisor_sum(IN SV* svn, ...)
void
znorder(IN SV* sva, IN SV* svn)
ALIAS:
- jordan_totient = 1
- legendre_phi = 2
+ binomial = 1
+ jordan_totient = 2
+ legendre_phi = 3
PREINIT:
int astatus, nstatus;
PPCODE:
@@ -774,11 +775,15 @@ znorder(IN SV* sva, IN SV* svn)
switch (ix) {
case 0: ret = znorder(a, n);
break;
- case 1: ret = jordan_totient(a, n);
+ case 1: ret = binomial(a, n);
+ if (ret == 0 && n <= a)
+ goto overflow;
+ break;
+ case 2: ret = jordan_totient(a, n);
if (ret == 0 && n > 1)
goto overflow;
break;
- case 2:
+ case 3:
default: ret = legendre_phi(a, n);
break;
}
@@ -788,8 +793,9 @@ znorder(IN SV* sva, IN SV* svn)
overflow:
switch (ix) {
case 0: _vcallsub_with_gmp("znorder"); break;
- case 1: _vcallsub_with_pp("jordan_totient"); break;
- case 2:
+ case 1: _vcallsub_with_gmp("binomial"); break;
+ case 2: _vcallsub_with_pp("jordan_totient"); break;
+ case 3:
default: _vcallsub_with_pp("legendre_phi"); break;
}
return; /* skip implicit PUTBACK */
diff --git a/lib/Math/Prime/Util.pm b/lib/Math/Prime/Util.pm
index 403c605..77945ff 100644
--- a/lib/Math/Prime/Util.pm
+++ b/lib/Math/Prime/Util.pm
@@ -44,8 +44,8 @@ our @EXPORT_OK =
moebius mertens euler_phi jordan_totient exp_mangoldt liouville
partitions
chebyshev_theta chebyshev_psi
- divisor_sum
- carmichael_lambda kronecker znorder znprimroot znlog legendre_phi
+ divisor_sum carmichael_lambda
+ kronecker binomial znorder znprimroot znlog legendre_phi
ExponentialIntegral LogarithmicIntegral RiemannZeta RiemannR
);
our %EXPORT_TAGS = (all => [ @EXPORT_OK ]);
@@ -2199,6 +2199,14 @@ function, and GMP's C<mpz_kronecker(a,n)>, C<mpz_jacobi(a,n)>, and
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.
+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.
+
+
=head2 znorder
$order = znorder(2, next_prime(10**16)-6);
diff --git a/lib/Math/Prime/Util/PP.pm b/lib/Math/Prime/Util/PP.pm
index 61d038f..48153bc 100644
--- a/lib/Math/Prime/Util/PP.pm
+++ b/lib/Math/Prime/Util/PP.pm
@@ -1785,6 +1785,12 @@ sub kronecker {
}
return ($b == 1) ? $k : 0;
}
+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;
+}
sub _is_perfect_square {
my($n) = @_;
diff --git a/lib/Math/Prime/Util/PPFE.pm b/lib/Math/Prime/Util/PPFE.pm
index ca96288..5b0f901 100644
--- a/lib/Math/Prime/Util/PPFE.pm
+++ b/lib/Math/Prime/Util/PPFE.pm
@@ -210,6 +210,13 @@ sub kronecker {
return Math::Prime::Util::PP::kronecker(@_);
}
+sub binomial {
+ my($n, $k) = @_;
+ _validate_positive_integer($n);
+ _validate_positive_integer($k);
+ return Math::Prime::Util::PP::binomial($n, $k);
+}
+
sub znorder {
my($a, $n) = @_;
_validate_positive_integer($a);
diff --git a/t/19-moebius.t b/t/19-moebius.t
index 9fd42f4..cba20f2 100644
--- a/t/19-moebius.t
+++ b/t/19-moebius.t
@@ -7,7 +7,7 @@ use Math::Prime::Util
qw/moebius mertens euler_phi jordan_totient divisor_sum exp_mangoldt
chebyshev_theta chebyshev_psi carmichael_lambda znorder liouville
znprimroot znlog kronecker legendre_phi gcd lcm is_power valuation
- invmod vecsum
+ invmod vecsum binomial
/;
my $extra = defined $ENV{EXTENDED_TESTING} && $ENV{EXTENDED_TESTING};
@@ -397,6 +397,20 @@ my @vecsums = (
[ "55340232221128654848", "18446744073709551616","18446744073709551616","18446744073709551616" ],
);
+my @binomials = (
+ [ 0,0, 1 ],
+ [ 0,1, 0 ],
+ [ 1,0, 1 ],
+ [ 1,1, 1 ],
+ [ 1,2, 0 ],
+ [ 13,13, 1 ],
+ [ 13,14, 0 ],
+ [ 40,19, "131282408400" ],
+ [ 67,31, "11923179284862717872" ],
+ [ 228,12, "30689926618143230620" ],
+ [ 177,78, "3314450882216440395106465322941753788648564665022000" ],
+);
+
# These are slow with XS, and *really* slow with PP.
if (!$usexs) {
%big_mertens = map { $_ => $big_mertens{$_} }
@@ -432,6 +446,7 @@ plan tests => 0 + 1
+ scalar(@valuations)
+ 3 + scalar(@invmods)
+ scalar(@vecsums)
+ + scalar(@binomials)
+ scalar(keys %powers)
+ scalar(keys %primroots) + 2
+ scalar(keys %jordan_totients)
@@ -656,6 +671,11 @@ foreach my $r (@vecsums) {
my($exp, @vals) = @$r;
is( vecsum(@vals), $exp, "vecsum(@vals) = $exp" );
}
+###### binomial
+foreach my $r (@binomials) {
+ my($n, $k, $exp) = @$r;
+ is( binomial($n,$k), $exp, "binomial($n,$k)) = $exp" );
+}
sub cmp_closeto {
my $got = shift;
diff --git a/util.c b/util.c
index 0dcdc03..1702d7a 100644
--- a/util.c
+++ b/util.c
@@ -1324,6 +1324,26 @@ int kronecker_ss(IV a, IV b) {
return kronecker_su(a, -b) * ((a < 0) ? -1 : 1);
}
+/* Thanks to MJD and RosettaCode */
+UV binomial(UV n, UV k) {
+ UV d, r = 1;
+ if (k >= n) return (k == n);
+ if (k > n/2) k = n-k;
+ for (d = 1; d <= k; d++) {
+ if (r >= UV_MAX/n) { /* Possible overflow */
+ UV g = gcd_ui(r, d);
+ r /= g;
+ if (r >= UV_MAX/n) return 0; /* Unavoidable overflow */
+ r *= n--;
+ r /= (d/g);
+ } else {
+ r *= n--;
+ r /= d;
+ }
+ }
+ return r;
+}
+
UV totient(UV n) {
UV i, nfacs, totient, lastf, facs[MPU_MAX_FACTORS+1];
if (n <= 1) return n;
diff --git a/util.h b/util.h
index 52529b7..15ea85e 100644
--- a/util.h
+++ b/util.h
@@ -44,6 +44,7 @@ extern int kronecker_uu(UV a, UV b);
extern int kronecker_su(IV a, UV b);
extern int kronecker_ss(IV a, IV b);
+extern UV binomial(UV n, UV k);
extern UV modinverse(UV a, UV p); /* Returns 1/a mod p */
extern UV divmod(UV a, UV b, UV n); /* Returns a/b mod n */
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git
More information about the Pkg-perl-cvs-commits
mailing list