123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301 |
- #!/usr/bin/perl
- # Efficient algorithm for computing `binomial(n, k) mod m`, based on the factorization of `m`.
- # Algorithm by Andrew Granville:
- # https://www.scribd.com/document/344759427/BinCoeff-pdf
- # Algorithm translated from (+some optimizations):
- # https://github.com/hellman/libnum/blob/master/libnum/modular.py
- # Translated by: Trizen
- # Date: 29 September 2017
- # Edit: 28 April 2022
- # https://github.com/trizen
- use 5.020;
- use strict;
- use warnings;
- use ntheory qw(:all);
- use experimental qw(signatures);
- sub modular_binomial ($n, $k, $m) {
- if ($m == 0) {
- return undef;
- }
- if ($m == 1) {
- return 0;
- }
- if ($k < 0) {
- $k = subint($n, $k);
- }
- if ($k < 0) {
- return 0;
- }
- if ($n < 0) {
- return modint(mulint(powint(-1, $k), __SUB__->(subint($k, $n) - 1, $k, $m)), $m);
- }
- if ($k > $n) {
- return 0;
- }
- if ($k == 0 or $k == $n) {
- return modint(1, $m);
- }
- if ($k == 1 or $k == subint($n, 1)) {
- return modint($n, $m);
- }
- my @congruences;
- foreach my $pair (factor_exp(absint($m))) {
- my ($p, $e) = @$pair;
- if ($e == 1) {
- push @congruences, [lucas_theorem($n, $k, $p), $p];
- }
- else {
- push @congruences, [modular_binomial_prime_power($n, $k, $p, $e), powint($p, $e)];
- }
- }
- modint(chinese(@congruences), $m);
- }
- #<<<
- #~ sub factorial_prime_pow ($n, $p) {
- #~ divint(subint($n, sumdigits($n, $p)), subint($p, 1));
- #~ }
- #>>>
- sub factorial_prime_pow ($n, $p) {
- my $count = 0;
- my $ppow = $p;
- while ($ppow <= $n) {
- $count = addint($count, divint($n, $ppow));
- $ppow = mulint($ppow, $p);
- }
- return $count;
- }
- sub binomial_prime_pow ($n, $k, $p) {
- #<<<
- factorial_prime_pow($n, $p)
- - factorial_prime_pow($k, $p)
- - factorial_prime_pow(subint($n, $k), $p);
- #>>>
- }
- sub factorial_without_prime ($n, $p, $pk) {
- return 1 if ($n <= 1);
- if ($p > $n) {
- return factorialmod($n, $pk);
- }
- my $r = 1;
- my $t = 0;
- foreach my $v (1 .. $n) {
- if (++$t == $p) {
- $t = 0;
- }
- else {
- $r = mulmod($r, $v, $pk);
- }
- }
- return $r;
- }
- sub lucas_theorem ($n, $k, $p) { # p is prime
- my $r = 1;
- while ($k) {
- my $np = modint($n, $p);
- my $kp = modint($k, $p);
- if ($kp > $np) { return 0 }
- my $rp = subint($np, $kp);
- my $x = factorialmod($np, $p);
- my $y = factorialmod($kp, $p);
- my $z = factorialmod($rp, $p);
- $y = mulmod($y, $z, $p);
- $x = divmod($x, $y, $p);
- $r = mulmod($r, $x, $p);
- $n = divint($n, $p);
- $k = divint($k, $p);
- }
- return $r;
- }
- sub binomial_non_prime_part ($n, $k, $p, $e) {
- my $pe = powint($p, $e);
- my $r = subint($n, $k);
- my $acc = 1;
- my @fact_pe = (1);
- if ($pe < ~0 and $p < $n) {
- my $count = 0;
- foreach my $x (1 .. vecmin(1e4, $pe - 1)) {
- if (++$count == $p) {
- $count = 0;
- }
- else {
- $acc = mulmod($acc, $x, $pe);
- }
- push @fact_pe, $acc;
- }
- }
- my $top = 1;
- my $bottom = 1;
- my $is_negative = 0;
- my $digits = 0;
- while ($n) {
- if ($digits >= $e) {
- $is_negative ^= modint($n, 2);
- $is_negative ^= modint($r, 2);
- $is_negative ^= modint($k, 2);
- }
- my $np = modint($n, $pe);
- my $rp = modint($r, $pe);
- my $kp = modint($k, $pe);
- #<<<
- $top = mulmod($top, ($fact_pe[$np] // factorial_without_prime($np, $p, $pe)), $pe);
- $bottom = mulmod($bottom, ($fact_pe[$rp] // factorial_without_prime($rp, $p, $pe)), $pe);
- $bottom = mulmod($bottom, ($fact_pe[$kp] // factorial_without_prime($kp, $p, $pe)), $pe);
- #>>>
- $n = divint($n, $p);
- $r = divint($r, $p);
- $k = divint($k, $p);
- ++$digits;
- }
- my $res = divmod($top, $bottom, $pe);
- if ($is_negative and ($p != 2 or $e < 3)) {
- $res = subint($pe, $res);
- }
- return $res;
- }
- sub modular_binomial_prime_power ($n, $k, $p, $e) {
- my $pow = binomial_prime_pow($n, $k, $p);
- if ($pow >= $e) {
- return 0;
- }
- my $er = $e - $pow;
- my $r = modint(binomial_non_prime_part($n, $k, $p, $er), powint($p, $er));
- my $pe = powint($p, $e);
- return mulmod(powmod($p, $pow, $pe), $r, $pe);
- }
- use Test::More tests => 44;
- is(modular_binomial(10, 2, 43), 2);
- is(modular_binomial(10, 8, 43), 2);
- is(modular_binomial(10, 2, 24), 21);
- is(modular_binomial(10, 8, 24), 21);
- is(modular_binomial(100, 42, -127), binomial(100, 42) % -127);
- is(modular_binomial(12, 5, 100000), 792);
- is(modular_binomial(16, 4, 100000), 1820);
- is(modular_binomial(100, 50, 139), 71);
- is(modular_binomial(1000, 10, 1243), 848);
- is(modular_binomial(124, 42, 1234567), 395154);
- is(modular_binomial(1e9, 1e4, 1234567), 833120);
- is(modular_binomial(1e10, 1e5, 1234567), 589372);
- is(modular_binomial(1e10, 1e5, 4233330243), 3403056024);
- is(modular_binomial(-1e10, 1e5, 4233330243), 2865877173);
- is(modular_binomial(1e10, 1e4, factorial(13)), 1845043200);
- is(modular_binomial(1e10, 1e5, factorial(13)), 1556755200);
- is(modular_binomial(1e10, 1e6, factorial(13)), 5748019200);
- is(modular_binomial(-1e10, 1e4, factorial(13)), 4151347200);
- is(modular_binomial(-1e10, 1e5, factorial(13)), 1037836800);
- is(modular_binomial(-1e10, 1e6, factorial(13)), 2075673600);
- is(modular_binomial(3, 1, 9), binomial(3, 1) % 9);
- is(modular_binomial(4, 1, 16), binomial(4, 1) % 16);
- is(modular_binomial(1e9, 1e5, 43 * 97 * 503), 585492);
- is(modular_binomial(1e9, 1e6, 5041689707), 15262431);
- is(modular_binomial(1e7, 1e5, 43**2 * 97**3 * 13**4), 1778017500428);
- is(modular_binomial(1e7, 1e5, 42**2 * 97**3 * 13**4), 10015143223176);
- is(modular_binomial(1e9, 1e5, 12345678910), 4517333900);
- is(modular_binomial(1e9, 1e6, 13**2 * 5**6), 2598375);
- is(modular_binomial(1e10, 1e5, 1234567), 589372);
- is(modular_binomial(1e5, 1e3, 43), binomial(1e5, 1e3) % 43);
- is(modular_binomial(1e5, 1e3, 43 * 97), binomial(1e5, 1e3) % (43 * 97));
- is(modular_binomial(1e5, 1e3, 43 * 97 * 43), binomial(1e5, 1e3) % (43 * 97 * 43));
- is(modular_binomial(1e5, 1e3, 43 * 97 * (5**5)), binomial(1e5, 1e3) % (43 * 97 * (5**5)));
- is(modular_binomial(1e5, 1e3, next_prime(1e4)**2), binomial(1e5, 1e3) % next_prime(1e4)**2);
- is(modular_binomial(1e5, 1e3, next_prime(1e4)), binomial(1e5, 1e3) % next_prime(1e4));
- is(modular_binomial(1e6, 1e3, next_prime(1e5)), binomial(1e6, 1e3) % next_prime(1e5));
- is(modular_binomial(1e6, 1e3, next_prime(1e7)), binomial(1e6, 1e3) % next_prime(1e7));
- is(modular_binomial(1234567, 1e3, factorial(20)), binomial(1234567, 1e3) % factorial(20));
- is(modular_binomial(1234567, 1e4, factorial(20)), binomial(1234567, 1e4) % factorial(20));
- is(modular_binomial(1e6, 1e3, powint(2, 128) + 1), binomial(1e6, 1e3) % (powint(2, 128) + 1));
- is(modular_binomial(1e6, 1e3, powint(2, 128) - 1), binomial(1e6, 1e3) % (powint(2, 128) - 1));
- is(modular_binomial(1e6, 1e4, (powint(2, 128) + 1)**2), binomial(1e6, 1e4) % ((powint(2, 128) + 1)**2));
- is(modular_binomial(1e6, 1e4, (powint(2, 128) - 1)**2), binomial(1e6, 1e4) % ((powint(2, 128) - 1)**2));
- is(modular_binomial(-10, -9, -10), binomial(-10, -9) % -10);
- say("binomial(10^10, 10^5) mod 13! = ", modular_binomial(1e10, 1e5, factorial(13)));
- say modular_binomial(12, 5, 100000); #=> 792
- say modular_binomial(16, 4, 100000); #=> 1820
- say modular_binomial(100, 50, 139); #=> 71
- say modular_binomial(1000, 10, 1243); #=> 848
- say modular_binomial(124, 42, 1234567); #=> 395154
- say modular_binomial(1e9, 1e4, 1234567); #=> 833120
- say modular_binomial(1e10, 1e5, 1234567); #=> 589372
- __END__
- my $upto = 10;
- foreach my $n (-$upto .. $upto) {
- foreach my $k (-$upto .. $upto) {
- foreach my $m (-$upto .. $upto) {
- next if ($m == 0);
- say "Testing: binomial($n, $k, $m)";
- is(modular_binomial($n, $k, $m), binomial($n, $k) % $m);
- }
- }
- }
|