modular_binomial_fast.sf

#!/usr/bin/ruby

# Translated by: Trizen
# Date: 27 April 2022
# https://github.com/trizen

# Fast algorithm for computing the binomial coefficient modulo some integer m.

# The implementation is based on Lucas' Theorem and its generalization given in the paper
# Andrew Granville "The Arithmetic Properties of Binomial Coefficients", In Proceedings of
# the Organic Mathematics Workshop, Simon Fraser University, December 12-14, 1995.

# Translation of binomod.gp v1.5 by Max Alekseyev, with some minor optimizations.

# See also:
#   https://home.gwu.edu/~maxal/gpscripts/

func factorial_without_prime(n,p,pk,from,res) {
    return 1 if (n <= 1)

    if (p > n) {
        return factorialmod(n, pk)
    }

    # Simple implementation:
    # return (1..n -> reduce {|a,b| (p `divides` b) ? a : mulmod(a, b, pk) })

    if (*from == n) {
        return *res
    }

    if (*from > n) {
        *from = 0
        *res  = 1
    }

    var r = *res;

    for v in ((*from + 1) .. n) {
        (p `divides` v) || do {
            r = mulmod(r, v, pk)
        }
    }

    *res   = r
    *from  = n

    return r
}

func lucas_theorem(n,k,p) {     # p is prime

    var r = 1

    while (k) {

        var (nqp, np) = divmod(n,p)
        var (kqp, kp) = divmod(k,p)

        if (kp > np) { return 0 }

        r = mulmod(r, divmod(
            factorialmod(np, p),
            factorialmod(kp, p) *
            factorialmod(np-kp, p), p),
        p)

        (n,k) = (nqp,kqp)
    }

    return r
}

func modular_binomial(n,k,m) {

    if (m.is_zero)               { return NaN }
    if (m.is_one)                { return 0 }
    if (k.is_neg)                { k = (n - k) }
    if (k.is_neg)                { return 0 }
    if (n.is_neg)                { return (((-1)**k * __FUNC__(-n + k - 1, k, m)) % m) }
    if (k > n)                   { return 0 }
    if (k.is_zero || (k == n))   { return 1%m }
    if (k.is_one  || (k == n-1)) { return n%m }
    if (n-k < k)                 { k = (n - k) }

    var F = []

    for p,q in (m.abs.factor_exp) {

        if (q == 1) {
            F << Mod(lucas_theorem(n,k,p), p)
            next
        }

        var d = n.ilog(p)+1

        var np = {|i| idiv(n, p**i) % p }.map(0 .. d)
        var kp = {|i| idiv(k, p**i) % p }.map(0 .. d)

        var e = []

        for i in (0 .. d) {
            e[i] = ((np[i] < (kp[i]+(i>0 ? e[i-1] : 0))) ? 1 : 0)
        }

        for i in (d-1 `downto` 0) {
            e[i] += e[i+1]
        }

        if (e[0] >= q) {
            F << Mod(0, p**q)
            next
        }

        var rq = q-e[0]

        var pq  = p**q
        var prq = p**rq

        var N = {|i| idiv(n,   p**i) % prq }.map(0 .. d)
        var K = {|i| idiv(k,   p**i) % prq }.map(0 .. d)
        var R = {|i| idiv(n-k, p**i) % prq }.map(0 .. d)

        var NKR = N.range.map {|i|
            [N[i], K[i], R[i], N[i]+K[i]+R[i]]
        }.sort_by { .tail }

        N = NKR.map{ _[0] }
        K = NKR.map{ _[1] }
        R = NKR.map{ _[2] }

        var acc = [1]

        for k in (1 .. min(max(N+K+R), 1e3)) {
            if (p `divides` k) {
                acc.push(acc[-1])
            } else {
                acc.push(mulmod(acc[-1], k, prq))
            }
        }

        var from = 0
        var res  = 1

        var v = (p**e[0] * prod(0 .. d, {|j|

            var pairs = []

            var a = (acc[N[j]] \\ pairs.push([\a, N[j]]))
            var b = (acc[K[j]] \\ pairs.push([\b, K[j]]))
            var c = (acc[R[j]] \\ pairs.push([\c, R[j]]))

            pairs.sort_by { .tail }.each_2d {|x,y|
                *x = factorial_without_prime(y, p, prq, \from, \res)
            }

            Mod(a, pq) / (b * c)
        }))

        if ((p>2 || rq<3) && (rq <= e.len)) {
            v *= (-1)**e[rq-1]
        }

        F << v
    }

    F || return 0
    chinese(F...).lift % m
}

say ("binomial(10^10, 10^5) mod 13! = ", modular_binomial(1e10, 1e5, 13!))

#
## Run some tests
#

assert_eq(modular_binomial(10, 2, 43), 2)
assert_eq(modular_binomial(10, 8, 43), 2)

assert_eq(modular_binomial(10, 2, 24), 21)
assert_eq(modular_binomial(10, 8, 24), 21)

assert_eq(modular_binomial(100, 42, -127), binomial(100, 42) % -127)

assert_eq(modular_binomial(12,    5,   100000),      792)
assert_eq(modular_binomial(16,    4,   100000),      1820)
assert_eq(modular_binomial(100,   50,  139),         71)
assert_eq(modular_binomial(1000,  10,  1243),        848)
assert_eq(modular_binomial(124,   42,  1234567),     395154)
assert_eq(modular_binomial(1e9,   1e4, 1234567),     833120)
assert_eq(modular_binomial(1e10,  1e5, 1234567),     589372)
assert_eq(modular_binomial(-1e10, 1e5, 4233330243),  2865877173)

assert_eq(modular_binomial(1e10, 1e4, 13!), 1845043200)
assert_eq(modular_binomial(1e10, 1e5, 13!), 1556755200)
assert_eq(modular_binomial(1e10, 1e6, 13!), 5748019200)

assert_eq(modular_binomial(-1e10, 1e4, 13!), 4151347200)
assert_eq(modular_binomial(-1e10, 1e5, 13!), 1037836800)
assert_eq(modular_binomial(-1e10, 1e6, 13!), 2075673600)

assert_eq(modular_binomial(3, 1,  9), binomial(3, 1) % 9)
assert_eq(modular_binomial(4, 1, 16), binomial(4, 1) % 16)

assert_eq(modular_binomial(1e9, 1e5, 43*97*503), 585492)
assert_eq(modular_binomial(1e9,1e6,5041689707), 15262431)
assert_eq(modular_binomial(1e7,1e5, 43**2 * 97**3 * 13**4), 1778017500428)
assert_eq(modular_binomial(1e7,1e5, 42**2 * 97**3 * 13**4), 10015143223176)
assert_eq(modular_binomial(1e9,1e5,12345678910), 4517333900)
assert_eq(modular_binomial(1e9,1e6,13**2 * 5**6), 2598375)
assert_eq(modular_binomial(1e10, 1e5, 1234567), 589372)

assert_eq(modular_binomial(1e5, 1e3, 43), binomial(1e5,1e3)%43)
assert_eq(modular_binomial(1e5, 1e3, 43*97), binomial(1e5,1e3) % (43*97))
assert_eq(modular_binomial(1e5, 1e3, 43*97*43), binomial(1e5,1e3) % (43*97*43))
assert_eq(modular_binomial(1e5, 1e3, 43*97*(5**5)), binomial(1e5,1e3) % (43*97*(5**5)))
assert_eq(modular_binomial(1e5,1e3, next_prime(1e4)**2), binomial(1e5, 1e3) % next_prime(1e4)**2)
assert_eq(modular_binomial(1e5,1e3, next_prime(1e4)), binomial(1e5, 1e3) % next_prime(1e4))
assert_eq(modular_binomial(1e6,1e3, next_prime(1e5)), binomial(1e6, 1e3) % next_prime(1e5))
assert_eq(modular_binomial(1e6,1e3, next_prime(1e7)), binomial(1e6, 1e3) % next_prime(1e7))
assert_eq(modular_binomial(1234567, 1e3, 20!), binomial(1234567, 1e3)%20!)
assert_eq(modular_binomial(1234567, 1e4, 20!), binomial(1234567, 1e4)%20!)

# --- 1. Edge Cases & Base Constraints ---
assert_eq(modular_binomial(10, 5, 1), 0)
assert_eq(modular_binomial(10, -5, 7), 0)
assert_eq(modular_binomial(-5, -2, 7), 0)
assert_eq(modular_binomial(-10, 5, 17), 4)
assert_eq(modular_binomial(-10, 4, 17), 1)
assert_eq(modular_binomial(5, 10, 7), 0)
assert_eq(modular_binomial(10, 0, 13), 1)
assert_eq(modular_binomial(10, 10, 13), 1)
assert_eq(modular_binomial(10, 1, 13), 10)
assert_eq(modular_binomial(10, 9, 13), 10)
assert_eq(modular_binomial(20, 17, 11), binomial(20, 3) % 11)

# --- 2. Algorithmic Fast Paths ---
# k <= 10^4 (bypasses prime factorization)
assert_eq(modular_binomial(10**6, 5000, 10**9 + 7), binomial(10**6, 5000) % (10**9 + 7))

# Square-free m (triggers Lucas Theorem directly)
assert_eq(modular_binomial(10**7, 10**5 + 1, 2310), binomial(10**7, 10**5 + 1) % 2310)

# n < p^q (triggers small_k routine directly on prime power)
assert_eq(modular_binomial(5000, 2500, 101**3), binomial(5000, 2500) % 101**3)

# e_0 >= q (correctly evaluates to 0 mod p^q)
assert_eq(modular_binomial(343, 7, 49), 0)

# --- 3. Prime Power Generalizations (Granville) ---
# p^2 optimization path
assert_eq(modular_binomial(10**7, 10**4 + 5, 1009**2), binomial(10**7, 10**4 + 5) % 1009**2)

# p^e (e >= 3) polynomial fast path
assert_eq(modular_binomial(10**6, 10**4 + 3, 11**3), binomial(10**6, 10**4 + 3) % 11**3)
assert_eq(modular_binomial(10**7, 10**4 + 1, 5**5), binomial(10**7, 10**4 + 1) % 5**5)

# Deep p^e (e = 10) polynomial expansion holding precision
assert_eq(modular_binomial(10**8, 10**4 + 7, 3**10), binomial(10**8, 10**4 + 7) % 3**10)

# --- 4. Composite Modulus / Chinese Remainder Theorem ---
# CRT successfully recombines Lucas, p^2, and p^e results
assert_eq(modular_binomial(10**7, 10**4 + 11, 11**3 * 13**2 * 17 * 19), binomial(10**7, 10**4 + 11) % (11**3 * 13**2 * 17 * 19))

# --- 5. The "Buggy Block" Edge Case Check ---
# Safely handles sign carry adjustments when q is vastly larger than the base-p representation of n
assert_eq(modular_binomial(100, 50, 2**40), binomial(100, 50) % 2**40)

# Tests from Math::Prime::Util
assert_eq(binomialmod(0,0,7), 1)
assert_eq(modular_binomial(0,1,7), 0)
assert_eq(modular_binomial(0,2,7), 0)
assert_eq(modular_binomial(3,0,7), 1)
assert_eq(modular_binomial(7,5,11), 10)
assert_eq(modular_binomial(950,100,123456), 24942)
assert_eq(modular_binomial(950,100,7), 2)
assert_eq(modular_binomial(8100,4000,1155), 924)
assert_eq(modular_binomial(950,100,1000000007), 640644226)
assert_eq(modular_binomial(189,34,877), 81)
assert_eq(modular_binomial(189,34,253009), 47560)
assert_eq(modular_binomial(189,34,36481), 14169)
assert_eq(modular_binomial(1900,17,41), 0)
assert_eq(modular_binomial(5000,654,101223721), 59171352)
assert_eq(modular_binomial(-112,5,351), 313)
assert_eq(modular_binomial(-189,34,877), 141)
assert_eq(modular_binomial(-23,-29,377), 117)
assert_eq(modular_binomial(189,-34,877), 0)
assert_eq(modular_binomial(100000000,87654321,1005973), 937361)
assert_eq(modular_binomial(100000000,7654321,1299709), 582708)
assert_eq(modular_binomial(100000000,7654321,12345678), 4152168)
assert_eq(modular_binomial(100000,7654,32768), 12288)
assert_eq(modular_binomial(100000,7654,196608), 110592)
assert_eq(modular_binomial(100000,7654,101223721), 5918452)
assert_eq(modular_binomial(100000000,7654321,32768), 24576)
assert_eq(modular_binomial(100000000,7654321,196608), 122880)
assert_eq(modular_binomial(100000000,7654321,101223721), 5463123)

__END__
var upto = 10
for n in (-upto .. upto), k in (-upto .. upto), m in (-upto .. upto) {
    if (m != 0) {
        say "Testing: binomial(#{n}, #{k}, #{m})"
        assert_eq(modular_binomial(n, k, m), binomial(n, k) % m)
    }
}