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))