It's a little late, but I had a thought that code could be made more general, slightly faster, and much easier to understand if the popcount logic were to be replaced with a table that records how many bits of the factorial were shifted out to make it odd.

from math import comb, perm, factorial as fact

Modulus = 2 ** 64
Cmax = 67
Pmax = 25
Fmax = Pmax

F = []          # odd factorial
S = []          # shift back to factorial
Finv = []       # multiplicative inverse of odd fact
for n in range(Cmax+1):
    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))