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 <d...@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 _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits