BZOJ-1951 [Sdoi2010]古代猪文(CRT+Lucas+欧拉降幂)

题目描述

  给定整数 \(n,q(1\leq n,q\leq 10^9)\)，计算 \(q^{\sum_{d|n}C_{n}^{d}}\mod {999911659}\)。

分析

  若 \(q=999911659\)，则上式为 \(0\)，否则因为 \(q,n\) 互质，根据扩展欧拉定理：

\[q^{\sum_{d|n}C_{n}^{d}}\equiv q^{\sum_{d|n}C_{n}^{d}\mod {999911658}}\pmod{999911659} \]

  所以本题的关键是如何计算 \(\sum_{d|n}C_{n}^{d}\mod {999911658}\)。

  分解质因数，发现 \(999911658=2\times 3\times 4679\times 35617\)。

  枚举 \(n\) 的因子 \(d\)，用 \(\text{Lucas}\) 定理求组合数 \(C_{n}^{d}\)，分别计算出 \(\sum_{d|n}C_{n}^{d}\) 对 \(2,3,4679,35617\) 四个质数取模的结果，记为 \(a_1,a_2,a_3,a_4\)。在计算过程中，对于一个质数 \(p\)，预处理 \(p\) 以内的所有阶乘以及阶乘的模 \(p\) 乘法逆元，就能快速计算按照 \(p\) 进制表示之后，每一位对应的组合数。

  最后用中国剩余定理求解线性同余方程组：

\[ \begin{cases}x\equiv a_1\pmod {2}\\ x\equiv a_2\pmod {3}\\ x\equiv a_3\pmod {4679}\\ x\equiv a_4\pmod {35617}\end{cases} \]

  即可得到 \(\sum_{d|n}C_{n}^{d}\mod {999911658}\) 的最小非负整数解 \(x\)。再用快速幂求 \(q^x\) 即可。

代码

#include<bits/stdc++.h>
using namespace std;
long long n,q;
const int mod=999911658;
long long b[]={0,2,3,4679,35617};
long long quick_pow(long long a,long long b,long long mod)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
long long fac[50010],a[10];
void init(long long p)
{
    fac[0]=1;
    for(int i=1;i<=p;i++)
        fac[i]=fac[i-1]*i%p;
}
long long C(long long n,long long m,long long p)
{
    if(m>n)
        return 0;
    return fac[n]*quick_pow(fac[m],p-2,p)%p*quick_pow(fac[n-m],p-2,p)%p;
}
long long Lucas(long long n,long long m,long long p)
{
    if(m>n)
        return 0;
    if(n==0)
        return 1;
    return Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return a;
    }
    long long gcd=exgcd(b,a%b,x,y);
    int temp=x;x=y;y=temp-a/b*y;
    return gcd;
}
long long CRT()
{
    long long ans=0,lcm=mod,x,y;
    for(int i=1;i<=4;i++)
    {
        long long temp=lcm/b[i];
        exgcd(temp,b[i],x,y);
        x=(x%b[i]+b[i])%b[i];
        ans=(ans+temp*x%lcm*a[i]%lcm)%lcm;
    }
    return (ans+lcm)%lcm;
}
int main()
{
    cin>>n>>q;
    if(q%(mod+1)==0)
    {
        puts("0");
        return 0;
    }
    for(int j=1;j<=4;j++)
    {
        init(b[j]);
        for(long long i=1;i*i<=n;i++)
        {
            if(n%i==0)
            {
                a[j]=(a[j]+Lucas(n,i,b[j]))%b[j];
                if(i*i!=n)
                    a[j]=(a[j]+Lucas(n,n/i,b[j]))%b[j];
            }
        }
    }
    long long x=CRT();
    cout<<quick_pow(q,x,mod+1)<<endl;
    return 0;
}