Gym - 102460M DivModulo

题意：

定义 $ndmodd=m\times d^{x}dmodd=mmodd$，其中 $n=m\times d^x$ 且 $d$ 不是 $m$ 的因子。先给定 $M,N,D$，求 $C(M,N)dmodD$。$(N\le M\le 4\times 10^{18},2\le D\le 1.6\times 10^7)$

链接：

https://vjudge.net/problem/Gym-102460M

解题思路：

由于习惯，将题目里的 $NM$ 对调一下，即求 $C(N,M)dmodD$。记 $N/D$ 为 $N$ 去掉所有 $D$ 因子后所得数，当 $D$ 为质数时，
$$C(N,M)=\frac{N!}{M!\ast (N-M)!}=\frac{N!/D}{M!/D\ast (N-M)!/D}\ast D^K$$
很显然，答案是前面那部分，由于分母不含 $D$ 因子，分母与 $D$ 互质，逆元存在，
$$ans=N!/D\ast inv(M!/D)\ast inv(N-M)!/D)modD$$
而
n ! ≡ ∏ d ∣ i i ∗ ∏ d ∤ i i ≡ {   d ⌊ n d ⌋ ∗ ⌊ n d ⌋ !   } ∗ {   [   ∏ d ∤ i i ∈ [ 1 , p ] i   ] ⌊ n p ⌋ ∗ ∏ d ∤ i i ∈ [ 1 , n % p ] i   } ( m o d   p ) n! \equiv \prod\limits_{d\mid i}i * \prod\limits_{d \nmid i}i \equiv \{~d^{\lfloor\frac{n}{d}\rfloor} * \lfloor\cfrac{n}{d}\rfloor!~\} * \{~[~{\prod\limits_{d\nmid i}^{i\in [1,p]} i~]}^{\lfloor\frac{n}{p}\rfloor} * \prod\limits_{d \nmid i}^{i\in[1,n\%p]}i~\} \quad (mod~p) n!≡d∣i∏​i∗d∤i∏​i≡{ d⌊dn​⌋∗⌊dn​⌋! }∗{ [ d∤i∏i∈[1,p]​i ]⌊pn​⌋∗d∤i∏i∈[1,n%p]​i }(mod p)
记
$$f(n)=n!/dmodp,g(n)=\prod _{d\nmid i}^{i\in [1,n]}imodp$$
有
$$f(n)=f(n/d)\ast g(p{)}^{\lfloor \frac{n}{p}\rfloor }\ast g(n\%p)modp$$
$O(p)$ 预处理出 $g(n)$，则可以快速求得 $f(n)$。
 
当 $D$ 不是质数时，分母与 $D$ 不一定互质，将 $D$ 唯一分解：
$$D=D_1\ast D_2\ast ...\ast D_s=p_1^{{\alpha }_1}\ast p_2^{{\alpha }_2}\ast ...\ast p_s^{{\alpha }_s}$$
依然可以求出
$$C(N,M)=\frac{N!}{M!\ast (N-M)!}=\frac{N!/p_i}{M!/p_i\ast (N-M)!/p_i}\ast p_i^{k_i}$$
进一步可以得到 $D^K$ 中 $K=\underset{i=1}{\overset{s}{\min }}\frac{k_i}{{\alpha }_i}$，
记 $C(N,M)=A$，有
$$ans\equiv A/D\equiv \frac{A}{D^K}\equiv \frac{A/p_i\ast p_i^{k_i}}{D^K}\equiv \frac{A/p_i\ast p_i^{k_i-K\ast {\alpha }_i}}{(\frac{D}{D_i}{)}^K}(mod{D}_i)$$
此时分母与模数互质，可以求出 $ans\%D_i$ 的结果，再用 $crt$ 求解即可。

参考代码：

#include<bits/stdc++.h>
 
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;
#define sz(a) ((int)a.size())
#define pb push_back
#define lson (rt << 1)
#define rson (rt << 1 | 1)
#define gmid (l + r >> 1)
const int maxn = 1e5 + 5;
const int inf = 0x3f3f3f3f;
const int mod = 1e9 + 7;

ll tp[16000005];
ll m[maxn], a[maxn], pi[maxn], num[maxn], ki[maxn];
ll n, N, M, D, K;

ll qpow(ll a, ll b, ll p){

    ll ret = 1;
    while(b){

        if(b & 1) ret = ret * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}

void exgcd(ll a, ll b, ll &x, ll &y, ll &d){

    if(!b) d = a, x = 1, y = 0;
    else exgcd(b, a % b, y, x, d), y -= a / b * x;
}

ll inv(ll a, ll p){

    ll x, y, d;
    exgcd(a, p, x, y, d);
    x = (x % p + p) % p;
    return x;
}

ll crt(){

    ll M = 1, ret = 0, x, y, d;
    for(int i = 1; i <= n; ++i) M *= m[i];
    for(int i = 1; i <= n; ++i){

        exgcd(M / m[i], m[i], x, y, d);
        x = (x % m[i] + m[i]) % m[i];
        ret = (ret + M / m[i] * x % M * a[i] % M) % M;
    }
    return ret;
}

ll cal(ll x, ll p){

    ll ret = 0;
    while(x){

        ret += x / p;
        x /= p;
    }
    return ret;
}

ll cal_fac(ll p){

    return cal(N, p) - cal(M, p) - cal(N - M, p);
}

ll dfs(ll n, ll pi, ll p){

    if(!n) return 1;
    ll ret = qpow(tp[p], n / p, p) * tp[n % p] % p;
    return ret * dfs(n / pi, pi, p) % p;
}

ll cal_rem(ll pi, ll p){

    tp[0] = 1;
    for(ll i = 1; i <= p; ++i){

        tp[i] = tp[i - 1];
        if(i % pi) tp[i] = tp[i] * i % p;
    }
    return dfs(N, pi, p) * inv(dfs(M, pi, p), p) % p * inv(dfs(N - M, pi, p), p) % p;
}

void init(){

    ll x = D;
    for(ll i = 2; i * i <= x; ++i){

        if(x % i) continue;
        num[++n] = 0, pi[n] = i, m[n] = 1;
        while(x % i == 0) ++num[n], x /= i, m[n] *= i;
    }
    if(x > 1) num[++n] = 1, pi[n] = m[n] = x;
}

int main(){
 
    ios::sync_with_stdio(0); cin.tie(0);
    cin >> N >> M >> D;
    init();
    ki[1] = cal_fac(pi[1]), K = ki[1] / num[1];
    for(int i = 2; i <= n; ++i){
        
        ki[i] = cal_fac(pi[i]);
        K = min(K, ki[i] / num[i]);
    }
    for(int i = 1; i <= n; ++i){

        a[i] = cal_rem(pi[i], m[i]) * qpow(pi[i], ki[i] - K * num[i], m[i]) % m[i] * qpow(inv(D / m[i], m[i]), K, m[i]) % m[i];
    }
    ll ret = crt();
    cout << ret << endl;
    return 0;
}