bzoj 1951(Lucas定理+中国剩余定理)

传送门
新技能get:为没有mod的组合数创造mod,然后使用Lucas定理
这里写图片描述

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD=999911659;
ll p[4]={2,3,4679,35617};
ll fac[4][36000],g;
ll d[36000],a[36000];int tot=0;
int n;
inline void init(int j) {
    fac[j][0]=fac[j][1]=1;
    for (int i=2;i<p[j];++i)
        fac[j][i]=fac[j][i-1]*i%p[j];
}
ll fpow(ll a,ll b,ll p) {
    ll ret=1;
    while (b) {
        if (b&1) ret=ret*a%p;
        b>>=1,a=a*a%p;
    }
    return ret;
}
ll lucas(int i,ll n,ll m) {
    ll ret=1;
    while (n&&m) {
        ll nn=n%p[i],mm=m%p[i];
        if (nn<mm) return 0;
        ret=ret*fac[i][nn]%p[i]*fpow(fac[i][nn-mm]*fac[i][mm]%p[i],p[i]-2,p[i])%p[i];
        n/=p[i],m/=p[i];
    }
    return ret;
}
ll exgcd(ll a,ll b,ll &x,ll &y) {
    if (!b) {
        x=1,y=0;
        return a;
    }
    ll d=exgcd(b,a%b,y,x);
    y-=(a/b)*x;
    return d;
}
ll CRT(int n,ll a[],ll p[]) {
    ll M=1;
    for (int i=0;i<n;++i) M*=p[i];
    ll ret=0;
    for (int i=0;i<n;++i) {
        ll w=M/p[i],x,y;
        ll d=exgcd(w,p[i],x,y);
        ret=(ret+x*w%M*a[i]%M)%M;
    }
    return (ret+M)%M;
}
int main() {
    scanf("%d%lld",&n,&g);
    if (g==MOD) {puts("0");return 0;}
    for (int i=1;i*i<=n;++i) 
        if (n%i==0) {
            if (i*i==n) d[++tot]=i;
            else d[++tot]=i,d[++tot]=n/i;
        }
    for (int j=0;j<4;++j) init(j);
    for (int j=1;j<=tot;++j) {
        for (int i=0;i<4;++i)
            a[i]+=lucas(i,n,d[j]),a[i]%=p[i];
    }
    ll ans=CRT(4,a,p);
    printf("%lld\n",fpow(g,ans,MOD));
    return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值