【BZOJ1951】古代猪文（SDOI2010）-数论大集合

本文链接：https://blog.youkuaiyun.com/Maxwei_wzj/article/details/80002341

本文解析了一道涉及费马小定理、中国剩余定理等数论知识的复杂题目，通过枚举约数和使用Lucas定理计算组合数对质数取模的结果，最终给出了O(nlogn)的时间复杂度解决方案。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

测试地址：古代猪文
做法：本题需要用到的数论知识有：费马小定理，枚举约数，中国剩余定理/合并模线性方程，扩展欧几里得，Lucas定理。可以说是一道很复杂的数论题了。
注意到题目要求的式子是：

G \sum k | n C n / k n mod p

$G^{\sum_{k|n}C_n^{n/k}}\mod p$
注意到

p=999911659 p = 999911659 $p=999911659$ 是一个质数，所以根据费马小定理，答案等同于：

G \sum k | n C n / k n mod (p - 1) mod p

$G^{\sum_{k|n}C_n^{n/k}\mod (p-1)}\mod p$
所以我们只需计算里面那个和式，然后快速幂算出最后的答案即可。枚举约数

k k $k$ 是

O (\sqrt{n})

$O(\sqrt n)$ 的，而计算组合数的话，首先要注意

p−1 p − 1 $p-1$ 不是质数，它质因数分解的结果为

2×3×4679×35617 2 × 3 × 4679 × 35617 $2\times 3\times 4679\times 35617$ ，我们没有办法，只能分别求出结果对这四个质数取模的取值，然后用中国剩余定理或者合并模线性方程来算出最后的答案，这里我用的是合并模线性方程的做法，即如有两个同余方程：

M≡r1(moda1) M ≡ r 1 ( mod a 1 ) $M\equiv r_1(\mod a_1)$ ，

M≡r2(moda2) M ≡ r 2 ( mod a 2 ) $M\equiv r_2(\mod a_2)$ ，那么它们合并之后的结果为：

M≡r1+a1x0(moda1a2gcd(a1,a2)) M ≡ r 1 + a 1 x 0 ( mod a 1 a 2 gcd ( a 1 , a 2 ) ) $M\equiv r_1+a_1x_0(\mod \frac{a_1a_2}{\gcd(a_1,a_2)})$ ，其中

x0 x 0 $x_0$ 为方程

a1x1−a2x2=r2−r1 a 1 x 1 − a 2 x 2 = r 2 − r 1 $a_1x_1-a_2x_2=r_2-r_1$ 的一个解，可以用扩展欧几里得算法求出。
而求组合数对一个质数取模的结果可以用Lucas定理来解决，即

Cmn%p=Cm/pn/pCm%pn%p%p C n m % p = C n / p m / p C n % p m % p % p $C_n^m\% p=C_{n/p}^{m/p}C_{n\%p}^{m\%p}\%p$ 。这样我们就解决了这道题，时间复杂度为

O(n−−√logn) O ( n log ⁡ n ) $O(\sqrt n\log n)$ 。
还有一点要注意，有一个点是

G=p G = p $G=p$ 的，而若是你算出来的幂数模

p−1 p − 1 $p-1$ 的结果是

0 0 $0$ ，那么你会输出

1

$1$ ，而正确答案是

0 0 $0$ ，要特判一下。
我傻逼的地方：一大堆公式记错，例如线性求逆元的状态转移方程是

i n v (i) = (p - p / i) \times i n v (p % i) % p

$inv(i)=(p-p/i)\times inv(p\%i)\% p$ ，以及合并模线性方程的公式也记错了……看来以后要经常背……
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int p;
ll n,g,mod[10]={999911659,2,3,4679,35617};
ll fac[40010],inv[40010],fi[40010],ans=0,tot,lasta=1,x;

ll C(ll n,ll m)
{
    if (n<m) return 0;
    if (n>=mod[p]||m>=mod[p])
        return C(n/mod[p],m/mod[p])*C(n%mod[p],m%mod[p])%mod[p];
    return fac[n]*fi[m]%mod[p]*fi[n-m]%mod[p];
}

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=(s*ss)%mod[0];
        ss=(ss*ss)%mod[0];b>>=1;
    }
    return s;
}

ll exgcd(ll a,ll b)
{
    ll x0=1,y0=0,x1=0,y1=1;
    while(b)
    {
        ll tmp,q;
        q=a/b;
        tmp=x0,x0=x1,x1=tmp-q*x1;
        tmp=y0,y0=y1,y1=tmp-q*y1;
        tmp=a,a=b,b=tmp%b;
    }
    return x0;
}

int main()
{
    scanf("%lld%lld",&n,&g);

    for(p=1;p<=4;p++)
    {
        fac[0]=fac[1]=inv[0]=inv[1]=fi[0]=fi[1]=1;
        for(ll i=2;i<=mod[p];i++)
        {
            fac[i]=fac[i-1]*i%mod[p];
            inv[i]=(mod[p]-mod[p]/i)*inv[mod[p]%i]%mod[p];
            fi[i]=fi[i-1]*inv[i]%mod[p];
        }
        tot=0;
        for(ll i=1;i*i<=n;i++)
            if (n%i==0)
            {
                tot+=C(n,i);
                if (i*i!=n) tot+=C(n,n/i);
            }
        x=exgcd(lasta,mod[p]);
        x=(x*(tot-ans)%mod[p]+mod[p])%mod[p];
        ans=(ans+x*lasta)%(lasta*mod[p]);
        lasta*=mod[p];
    }
    ans=power(g,ans);
    if (g==mod[0]) printf("0");
    else printf("%lld",ans);

    return 0;
}