题解 Luogu P2155 [SDOI2008]沙拉公主的困惑

Posted 2021-05-22 JustinRochester

传送门

发现好多人的做法并不对......

---

【分析】

先化简一波式子：

\\(\\quad Ans\\)
\\(\\displaystyle =\\sum_{i=1}^{N!}[\\gcd(i, M!)=1]\\)
\\(\\displaystyle =\\sum_{i=1}^{N!}\\sum_{d\\mid i\\wedge d\\mid (M!)}\\boldsymbol \\mu(d)\\)
\\(\\displaystyle =\\sum_{d\\mid (M!)}\\boldsymbol \\mu(d)\\sum_{i=1}^{N!}[d\\mid i]\\)

由于 \\(M\\leq N\\) 故 \\((M!)\\mid (N!)\\) 有 \\(d\\mid (N!)\\)

\\(\\therefore Ans\\)
\\(\\displaystyle =\\sum_{d\\mid (M!)}\\boldsymbol \\mu(d){N!\\over d}\\)
\\(\\displaystyle ={N!\\over M!}\\sum_{d\\mid (M!)}\\boldsymbol \\mu(d){M!\\over d}\\)
\\(\\displaystyle ={N!\\over M!}(\\boldsymbol \\mu*\\boldsymbol {id})(M!)\\)
\\(\\displaystyle ={N!\\over M!}\\boldsymbol \\varphi(M!)\\)
\\(\\displaystyle ={N!\\over M!}\\cdot M!\\cdot \\prod_{p\\mid (M!)}{p-1\\over p}\\)
\\(\\displaystyle =N!\\prod_{p\\leq M}{p-1\\over p}\\)

答案要求对 \\(R\\) 取模，但并没有表明 \\(R\\) 是否严格大于 \\(M\\)

故类似于 exLucas 的想法：设 \\(f_R(n)\\) 表示 \\(n\\) 去除 \\(R\\) 质因子后，剩余质因子的乘积；\\(g_R(n)\\) 表示 \\(n\\) 中 \\(R\\) 质因子的次数

故 \\(\\displaystyle n=f_R(n)\\cdot R^{\\displaystyle g_R(n)}\\)

代入所求式子得：\\(\\displaystyle Ans={f_R(N!)\\cdot f_R(\\prod_{p\\leq M}(p-1))\\over f_R(\\prod_{p\\leq M}p)}\\cdot R^{\\displaystyle g_R(N!)+g_R(\\prod_{p\\leq M}(p-1))-g_R(\\prod_{p\\leq M}p)}\\)

当 \\(\\displaystyle g_R(N!)+g_R(\\prod_{p\\leq M}(p-1))-g_R(\\prod_{p\\leq M}p)>0\\) 时，\\(Ans\\equiv 0\\pmod R\\)

当 \\(\\displaystyle g_R(N!)+g_R(\\prod_{p\\leq M}(p-1))-g_R(\\prod_{p\\leq M}p)=0\\) 时，\\(\\displaystyle Ans\\equiv f_R(N!)\\cdot f_R(\\prod_{p\\leq M}(p-1)) \\cdot f_R^{-1}(\\prod_{p\\leq M}p)\\pmod R\\)

---

现考虑如何线性求解三个数字的 \\(f_R\\) 和 \\(g_R\\)

对于 \\(N!\\) 有 \\(\\begin{cases} f_R(N!)\\equiv f_R((N-1)!)\\cdot f_R(N)\\pmod R \\\\ g_R(N!)\\equiv g_R((N-1)!)+g_R(N)\\pmod R \\end{cases}\\)

其中，\\(f_R(N)\\) 和 \\(g_R(N)\\) 可以通过对 \\(N\\) 暴力拆解得到，复杂度为 \\(T(N)=O(N)+O(N)+O({N\\over R})+O({N\\over R^2})+O({N\\over R^3})+\\cdots =O(N)+O({N\\over 1-{1\\over R}})=O(N)\\)

对于 \\(\\displaystyle \\prod_{p\\leq M} p\\) 有 \\(\\begin{cases} \\begin{cases} \\displaystyle f_R(\\prod_{p\\leq M}p)\\equiv f_R(\\prod_{p\\leq M-1}p)\\pmod R \\\\ \\displaystyle g_R(\\prod_{p\\leq M}p)\\equiv g_R(\\prod_{p\\leq M-1}p)\\pmod R \\end{cases}, M\\not\\in Prime \\\\\\ \\\\ \\begin{cases} \\displaystyle f_R(\\prod_{p\\leq M}p)\\equiv f_R(\\prod_{p\\leq M-1}p)\\cdot M^{[M\\neq R]}\\pmod R \\\\ \\displaystyle g_R(\\prod_{p\\leq M}p)\\equiv g_R(\\prod_{p\\leq M-1}p)+[M=R]\\pmod R \\end{cases}, M\\in Prime \\end{cases}\\)

复杂度为 \\(O(N)\\)

对于 \\(\\displaystyle \\prod_{p\\leq M}(p-1)\\) 有 \\(\\begin{cases} \\begin{cases} \\displaystyle f_R(\\prod_{p\\leq M}(p-1))\\equiv f_R(\\prod_{p\\leq M-1}(p-1))\\pmod R \\\\ \\displaystyle g_R(\\prod_{p\\leq M}(p-1))\\equiv g_R(\\prod_{p\\leq M-1}(p-1))\\pmod R \\end{cases}, M\\not\\in Prime \\\\\\ \\\\ \\begin{cases} \\displaystyle f_R(\\prod_{p\\leq M}(p-1))\\equiv f_R(\\prod_{p\\leq M-1}(p-1))\\cdot f_R(M-1)\\pmod R \\\\ \\displaystyle g_R(\\prod_{p\\leq M}(p-1))\\equiv g_R(\\prod_{p\\leq M-1}(p-1))+g_R(M-1)\\pmod R \\end{cases}, M\\in Prime \\end{cases}\\)

对于 \\(f_R(M-1)\\) 与 \\(g_R(M-1)\\) 同样是暴力拆解，复杂度不大于 \\(f_R(N!)\\) 和 \\(g_R(N!)\\) 的求解，故也是 \\(O(N)\\) 的

因此预处理总复杂度为 \\(O(N)+O(N)+O(N)=O(N)\\)

---

【代码】

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=1e7+10;
int R, prime[MAXN], cntprime, frac[MAXN], cntR[MAXN], prodP[MAXN], cntRP[MAXN], prodP1[MAXN], cntRP1[MAXN];
bool nprime[MAXN];
inline ll fpow(ll a,ll x) { ll ans=1; for(;x;x>>=1,a=a*a%R) if(x&1) ans=ans*a%R; return ans; }
inline int ans(int N, int M){
    if( cntR[N]+cntRP1[M]-cntRP[M] ) return 0;
    return (ll)frac[N]*prodP1[M]%R*fpow(prodP[M], R-2)%R;
}
inline void pre(){
    frac[0]=frac[1]=1;
    prodP[1]=prodP1[1]=1;
    for(int i=2;i<=1e7;++i){
        if(!nprime[i]) prime[++cntprime]=i;
        for(int j=1;j<=cntprime;++j)
            if((ll)prime[j]*i>1e7) break;
            else{
                nprime[prime[j]*i]=1;
                if(i%prime[j]==0) break;
            }
        cntR[i]=cntR[i-1];
        int Val=i;
        while(Val%R==0) ++cntR[i], Val/=R;
        frac[i]=(ll)frac[i-1]*Val%R;

        prodP[i]=prodP[i-1];
        cntRP[i]=cntRP[i-1];
        prodP1[i]=prodP1[i-1];
        cntRP1[i]=cntRP1[i-1];
        if(prime[cntprime]!=i) continue;
        if(i==R) ++cntRP[i];
        else prodP[i]=(ll)prodP[i]*i%R;
        Val=i-1;
        while(Val%R==0) ++cntRP1[i], Val/=R;
        prodP1[i]=(ll)prodP1[i]*Val%R;
    }
}
int main(){
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    int T, N, M; cin>>T>>R;
    pre();
    while(T--&&cin>>N>>M) cout<<ans(N, M)<<"\\n";
    cout.flush();
    return 0;
}