• + 1 comment

    Python3 solution

    import functools
    import sys
    
    class MathEngine():
        """A collection of algorithms for Modular Arithmetics.
    
           References:
               [1] Granville, A. "Binomial coefficients modulo prime powers".
               [2] Knuth, D. "The Art of Computer Programming", vol 1.
               [3] https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Polynomial_extended_Euclidean_algorithm
        """
        def __init__(self):
            self.pe = None
            self._cache_factorial = {}
    
        def binomial(self, n, m):
            """Find the binomial (n m) = n!/m!/(n-m)!."""
            if m == 0:
                if n == 0:
                    return 0
                else:
                    return 1
            b = self.factor_binomial(n, m)
            result = self.factors_to_integer(b)
            return result
    
        def binomial_prime_power(self, n, m, p):
            """Find the max power of `p` which divides a binomial (n m)."""
            result = self._number_of_carries(n-m,m,p)
            return result
    
        def factor_binomial(self, n, m):
            """Factor the binomial (n m) into primes.
    
            The algorithm is based on Kummer's theorem [1, page 7].
        
            Result:
                {d1:p1, ..., dn:pn}, so that (n m) = d1^p1 * ... * dn^pn.
            """
            result = {}
            for p in self.primes(n+1):
                power = self.binomial_prime_power(n,m,p)
                if power > 0:
                    result[p] = power
            return result
    
        def factor_integer(self, n):
            """Factor n into primes.
        
            Result:
                {d1:p1, ..., dn:pn}, so that n = d1^p1 * ... * dn^pn.
            """
            result = {}
            for p in self.primes(n+1):
                power = 0
                while n % p == 0 and n != 1:
                    n = n // p
                    power += 1
                if power > 0:
                    result[p] = power
                if n == 1:
                    break
            return result
    
        def factor_factorial(self, n):
            """Factor the factorial of n into primes.
    
            The algorithm is based on the Kummer's theorem formulated for binomial coefficients;
            for the case of factorials see [1, eq.(17)] or [2, eq.(1.2.5.8)].
        
            Result:
                {d1:p1, ..., dn:pn} such that n! = d1^p1 * ... * dn^pn.
            """
            result = {}
            for prime in self.primes(n+1):
                k = n
                result[prime] = 0
                while k > 0:
                    k = k // prime
                    result[prime] += k
            return result
    
        def factorial(self, n, p=0, mod=0):
            """Find the product of all integers <= n, i.e. the factorial of n;
            if p > 0 exclude integers divisible by p; if mod > 0 find result
            modulo mod."""
            _cache_key = (n,p,mod)
            if _cache_key in self._cache_factorial:
                return self._cache_factorial[_cache_key]
            result = 1
            if p == 0:
                if mod == 0:
                    for k in range(1,n+1):
                        result *= k
                else:
                    for k in range(1,n+1):
                        result = (result*k) % mod
            else:
                if mod == 0:
                    for k in range(1,n+1):
                        if k % p != 0:
                            result *= k
                else:
                    for k in range(1,n+1):
                        if k % p != 0:
                            result = (result*k) % mod
            self._cache_factorial[_cache_key] = result
            return result 
    
        def factors_to_integer(self, n):
            result = 1
            for prime,power in n.items():
                result *= prime**power
            return result
    
        def gcd_extended(self, a, b):
            """Find gcd(a,b) and coefficients s,t of the Bezout's identity, such that
            `a*s + b*t = gcd(a,b)`. See [3] for more details.
    
            Result:
                (gcd(a,b), s, t)
            """
            s, s_old = 0, 1
            t, t_old = 1, 0
            r, r_old = b, a
            while r != 0:
                k = r_old // r
                r_old, r = r, r_old - k * r
                s_old, s = s, s_old - k * s
                t_old, t = t, t_old - k * t
            return r_old, s_old, t_old
    
        def integer_digits(self, n, p):
            """Find digits of n in the p representation.
    
            Result:
                [a0,...,ak], such that n = a0 + a1*p + ... + ak*p^k, where ai < p.
            """
            result = []
            while n > 0:
                result.append(n % p)
                n //= p
            return result
    
        def integer_digits_len(self, n, p):
            result = 1
            while n > 0:
                result += 1
                n //= p
            return result
    
        def mod_binomial(self, n, m, mod):
            """Find the reminder from the division of the binomial (n m) by mod.
    
            The algorithm is simple and slow, but less prone to errors comparing
            to the faster and more complex generalized Lucas algorithm, for more
            details see `mod_binomial_lucas` method.
            """
            binomial_factors = self.factor_binomial(n, m)
            result = 1
            for prime, exp in binomial_factors.items():
                result = (result * pow(prime, exp, mod)) % mod
            return result
    
        def mod_binomial_lucas(self, n, m, p, q=1):
            """Find the reminder form the division of the binomial (n m) by p^q,
            where mod is a prime number. This routine employs Lucas' Theorem (for q=1)
            and its generalization for (q>1), see [1, Theorem 1].
            """
            e0 = self.binomial_prime_power(n,m,p)
            qq = q-e0
            if qq < 0:
                return 0
            p_e0, p_q, p_qq = p**e0, p**q, p**qq
            r = n - m
            d = self.integer_digits_len(n, p)
            _ni, _mi, _ri = n, m, r
            _least_positive_residue_n = []
            _least_positive_residue_m = []
            _least_positive_residue_r = []
            for i in range(0, d):
                _least_positive_residue_n.append(_ni % p_qq)
                _least_positive_residue_m.append(_mi % p_qq)
                _least_positive_residue_r.append(_ri % p_qq)
                _ni //= p
                _mi //= p
                _ri //= p
            _factorial = [1]
            curr = 1
            for x in range(1,p_qq):
                if x % p != 0:
                    curr = (curr*x) % p_qq
                _factorial.append(curr)
            result, i = 1, 0
            aaa, bbb, ccc = 1, 1, 1
            for i in range(0, d):
                aa = _least_positive_residue_n[i]
                bb = _least_positive_residue_m[i]
                cc = _least_positive_residue_r[i]
                a = _factorial[aa]
                b = _factorial[bb]
                c = _factorial[cc]
                aaa = (aaa*a) % p_qq
                bbb = (bbb*b) % p_qq
                ccc = (ccc*c) % p_qq
                i += 1
            result = aaa * self.mod_inverse(bbb*ccc % p_qq, p_qq) % p_qq
            if p != 2 or qq < 3:
                result *= (-1)**self._number_of_carries(m,r,p,qq-1)
            result = (result * p_e0) % p_q
            return result
    
        def _least_positive_residue(self, n, j, p, q):
            """Helper for mod_binomial_lucas method"""
            result = (n//p**j) % p**q
            return result
    
        def _number_of_carries(self, m, r, p, j=0):
            """Helper for mod_binomial_lucas method"""
            result, carry, jj = 0, 0, 0
            while (m > 0 or r > 0):
                if (m % p) + (r % p) + carry >= p:
                    carry = 1
                else:
                    carry = 0
                if jj >= j:
                    result += carry
                m, r = m // p, r // p
                jj += 1
            return result
    
        def mod_chinese(self, rs):
            n = functools.reduce(lambda a,b: a*b, rs.keys(), 1)
            result = 0
            for ni,ai in rs.items():
                nn = n // ni
                _, ri, si = self.gcd_extended(ni, nn)
                ei = si*nn
                result += ai*ei
            result = result % n
            return result
    
        def mod_inverse(self, n, mod):
            r, s, t = self.gcd_extended(n, mod)
            if r != 1:
                return None
            return s
    
        def primes(self, n):
            return self.pe.primes(n)
    
    class Solver(object):
        def __init__(self):
            self.me = MathEngine()
            self.mod = 142857
    
        def solution(self, n, r):
            """Find solution to the "nCr" problem."""
            mod_factors = {3: 3, 11: 1, 13: 1, 37: 1}
            mods = {p**q:self.me.mod_binomial_lucas(n, r, p, q) for p,q in mod_factors.items()}
            result = self.me.mod_chinese(mods)
            return result
    
    if __name__ == '__main__':
        T = int(sys.stdin.readline())
        solver = Solver()
        for i in range(T):
            n, r = [int(num) for num in sys.stdin.readline().split()]
            print(solver.solution(n, r))