/*
2024.2.12
常见求组合数4种模型 :
1. 数据范围支持预处理出所有组合数 C[n][m] 表示 C(n,m) 0 <= n , m <= 4000 O(N^2)
2. 从公式出发 C(n , m) = n! / ( m! * (n - m) !) 若数据范围不支持第一种预处理可以考虑若允许条件下考虑预处理出阶乘和模p下的逆元(保证精度)
∴ C(n , m) = fact[n] * infact[m] * infact[n - m] % p O(NlogN)
3. 当求的组合数组数很小,n\m巨大的时候无法进行预处理出发,则可以考虑从 lucas定理出发 有 C(a , b) = C(a % p , b % p) * C(a / p , b / p) (mod p)同余
计算的时候直接从定义出发时间复杂度在O(logp * logN * 20 * 20) ( 1 <= b <= a <= 1e18 , 1 <= p <= 1e5 , 1 <= n <= 20 约等于 4e7)
4. 没有取模的条件时且 0 <= b <= a <= 5000 涉及高精度计算
*/
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std ;
const int p = 1e9 + 7;
int p ;
int qpow(int a,int b) {
int ret = 1 ;
while ( b ) {
if ( b & 1 ) ret = (long long)ret * a % p ;
a = (long long) a * a % p ;
b >>= 1 ;
}
return ret ;
}
int inv(int x) {
return qpow(x , p - 2) ;
}
struct get_comb_1
{
#define MAXN 2010
int C[MAXN][MAXN] ;
void comb()
{
for(int i = 0 ; i < MAXN ; i++)
for(int j = 0 ; j <= i ; j++)
if ( !j ) C[i][j] = 1;
else C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % p ;
}
};
struct get_comb_2
{
#define MAXN 100010
int fact[MAXN] , infact[MAXN] ;
void comb() {
fact[0] = infact[0] = 1 ;
for(int i = 1 ; i < MAXN ; i++) {
fact[i] = (long long)fact[i - 1] * i % p ;
infact[i] = (long long)infact[i - 1] * inv(i) % p ;
}
}
int gets(int a,int b) {
return (__int128)fact[a] * infact[b] * infact[a - b] % p ;
}
};
struct get_comb_3
{
int C(int a, int b) {
int ret = 1 ;
for(int i = 1 , j = a ; i <= b ; i++ , j-- )
{
ret = (long long) ret * j % p ;
ret = (long long) ret * inv(i) % p ;
}
return ret ;
}
int lucas(long long a,long long b) {
if ( a < p && b < p ) return C(a , b) ;
else return (long long) C(a % p , b % p) * lucas(a / p , b / p) % p ;
}
};
struct get_comb_4
{
#define MAX_N 5010
int primes[MAX_N] ;
int sum[MAX_N] ;
bool st[MAX_N] ;
void init() {
for(int i = 0 ; i < MAX_N ; i++) primes[i] = sum[i] = 0 , st[i] = false ;
}
void get_prime() {
for(int i = 2 ; i < MAX_N ; i++)
{
if ( !st[i] ) {
primes[++primes[0]] = i ;
}
for(int j = 1 ; primes[j] <= MAX_N / i ; j++) {
st[i * primes[j]] = true ;
if ( i % primes[j] == 0 ) break ;
}
}
}
vector<int> mul(vector<int> A , int B ) {
vector<int> C ;
int t = 0 ;
for(int i = 0 ; i < (int)A.size() ; i++) {
t += A[i] * B ;
C.push_back(t % 10) ;
t /= 10 ;
}
while ( t ) C.push_back( t % 10 ) , t /= 10 ;
while ( C.size() > 1 && C.back() == 0 ) C.pop_back() ;
return C ;
}
int gets(int x,int p) {
int ret = 0 ;
while ( x ) {
ret += x / p ;
x /= p ;
}
return ret ;
}
vector<int> Comb(int a, int b) {
init() ;
get_prime() ;
for(int i = 1 ; i <= primes[0] ; i++) {
int p = primes[i] ;
sum[i] = gets(a , p) - gets(b , p) - gets(a - b , p) ;
}
vector<int> C(1 , 1) ;
for(int i = 1 ; i <= primes[0] ; i++)
for(int j = 1 ; j <= sum[i] ; j++)
C = mul(C , primes[i]) ;
return C ;
}
};