In [1]:
import numpy as np
In [2]:
N = 10**6
MOD = 10**9+7
In [3]:
Nsq = 10**3
fact = np.arange(N,dtype=np.int64).reshape(Nsq,Nsq); fact[0,0] = 1
for n in range(1,Nsq):
    fact[:,n] *= fact[:,n-1]; fact[:,n] %= MOD
for n in range(1,Nsq):
    fact[n] *= fact[n-1,-1]; fact[n] %= MOD
fact = fact.ravel()
In [4]:
fact_inv = np.arange(1,N+1,dtype=np.int64)[::-1].reshape(Nsq,Nsq); fact_inv[0,0] = pow(int(fact[N-1]),MOD-2,MOD)
for n in range(1,Nsq):
    fact_inv[:,n] *= fact_inv[:,n-1]; fact_inv[:,n] %= MOD
for n in range(1,Nsq):
    fact_inv[n] *= fact_inv[n-1,-1]; fact_inv[n] %= MOD
fact_inv = fact_inv.ravel()[::-1]
In [5]:
print(fact[:20])
print(fact_inv[:20])

[        1         1         2         6        24       120       720
      5040     40320    362880   3628800  39916800 479001600 227020758
 178290591 674358851 789741546 425606191 660911389 557316307]
[        1         1 500000004 166666668  41666667 808333339 301388891
 900198419 487524805 831947206 283194722 571199524 380933296 490841026
 320774361 821384963 738836565 514049213 639669405 402087866]
In [6]:
n = 10000
comb = fact[n] * fact_inv[:n+1] % MOD * fact_inv[n::-1] % MOD
print(comb[:20])

[        1     10000  49995000 616668838 709582588 797500005   2082363
 115876526 494264530 965683250 814128327  12902709 490429945 878019727
 270210871  21715928  52096223 772275699 812680514 219833368]