from math import comb, perm, factorial as fact TableSize = 121 Modulus = 2 ** 64 # 2 ** 64 2 ** 128 Cmax = 67 # 67 132 Pmax = 25 # 25 41 Fmax = Pmax F = [] # odd factorial S = [] # shift back to factorial Finv = [] # multiplicative inverse of odd fact for n in range(TableSize): f = fact(n) s = (f & -f).bit_length() - 1 odd_f = (f >> s) % Modulus inv_f = pow(odd_f, -1, Modulus) assert odd_f * inv_f % Modulus == 1 assert (odd_f << s) % Modulus == f % Modulus F.append(odd_f) S.append(s) Finv.append(inv_f) def fact_small(n): return F[n] << S[n] def perm_small(n, k): return (F[n] * Finv[n-k] % Modulus) << (S[n] - S[n-k]) def comb_small(n, k): return (F[n] * Finv[k] * Finv[n-k] % Modulus) << (S[n] - S[k] - S[n-k]) assert all(fact_small(n) == fact(n) for n in range(Fmax+1)) assert all(perm_small(n, k) == perm(n, k) for n in range(Pmax+1) for k in range(0, n+1)) assert all(comb_small(n, k) == comb(n, k) for n in range(Cmax+1) for k in range(0, n+1)) limits = bytearray(TableSize) # Assumes k <= n-k for n in range(0, TableSize): for k in range(0, n+1): if comb(n, k) != comb_small(n, k) % Modulus: break else: k += 1 limits[n] = k print(list(limits)) for n in range(TableSize): for k in range(0, n+1): if k < limits[n]: assert comb_small(n, k) == comb(n, k) # Example: print(comb(100, 15)) print(comb_small(100, 15))