【模板】拓展卢卡斯

最新推荐文章于 2019-11-01 22:18:32 发布

转载最新推荐文章于 2019-11-01 22:18:32 发布 · 120 阅读

0 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/LiGuanlin1124/p/11128840.html

本文深入探讨了ExLucas算法的实现细节，包括快速幂、扩展欧几里得算法、组合数计算等核心模块，并通过具体实例展示了算法的优化过程。文章还提供了完整的代码示例，便于读者理解和应用。

粘代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
template<typename T>
inline void read(T&x)
{
    T f = 1,c = 0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){c=c*10+ch-'0';ch=getchar();}
    x = f*c;
}
ll n,m,P;
ll fastpow(ll x,ll y,ll MOD)
{
    ll ret = 1;
    while(y)
    {
        if(y&1)ret=ret*x%MOD;
        x=x*x%MOD;y>>=1;
    }
    return ret;
}
void exgcd(ll a,ll b,ll&x,ll&y)
{
    if(!b){x=1,y=0;return ;}
    exgcd(b,a%b,y,x);y-=a/b*x;
}
ll inv(ll a,ll MOD){ll x,y;exgcd(a,MOD,x,y);return (x%MOD+MOD)%MOD;}
ll p[30],pk[30],pnt,f[30][1000050];
ll jc[1000050],jny[1000050];
void init()
{
    ll now = P;
    for(ll i=2;i*i<=now;i++)if(now%i==0)
    {
        p[++pnt]=i,pk[pnt]=1;
        while(now%i==0)now/=i,pk[pnt]*=i;
    }
    if(now!=1)
    {
        p[++pnt]=now,pk[pnt]=now;
    }
    if(p[1]==MOD)
    {
        jc[0] = 1;
        for(ll i=1;i<P;i++)jc[i]=jc[i-1]*i%P;
        jny[P-1] = inv(jc[P-1],P);
        for(ll i=P-1;i;i--)jny[i-1]=jny[i]*i%P;
        return ;
    }
    for(ll i=1;i<=pnt;i++)
    {
        f[i][0]=1;
        for(ll j=1;j<pk[i];j++)
            if(j%p[i])f[i][j]=f[i][j-1]*j%pk[i];
            else f[i][j]=f[i][j-1];
    }
}
ll calp(ll x,ll P){ll ret=0;for(x/=P;x;x/=P)ret+=x;return ret;}
ll calc(ll x,ll i)
{
    if(!x)return 1;
    return fastpow(f[i][pk[i]-1],x/pk[i],pk[i])*f[i][x%pk[i]]%pk[i]*calc(x/p[i],i);
}
ll Lucas(ll x,ll y)
{
    if(x<y)return 0;
    if(x<P)return jc[x]*jny[y]%P*jny[x-y]%P;
    return Lucas(x/P,y/P)*Lucas(x%P,y%P)%P;
}
ll eLucas(ll x,ll y,int i)
{
    ll k = calp(x,p[i])-calp(y,p[i])-calp(x-y,p[i]);
    return fastpow(p[i],k,pk[i])
            *calc(x,i)%pk[i]
            *inv(calc(y,i),pk[i])%pk[i]
            *inv(calc(x-y,i),pk[i])%pk[i];
}
ll ExLucas(ll x,ll y)
{
    if(p[1]==MOD)return Lucas(x,y);
    ll ans = 0;
    for(ll i=1;i<=pnt;i++)
        ans= (ans + eLucas(x,y,i) * (P/pk[i])% P * inv(P/pk[i],pk[i]) % P)%P;
    return ans;
}
int main()
{
    read(n),read(m),read(P);
    init();
    printf("%lld\n",ExLucas(n,m));
    return 0;
}

View Code

upd：有一个地方：

ll calc(ll x,ll i)
{
    if(!x)return 1;
    return fastpow(f[i][pk[i]-1],x/pk[i],pk[i])*f[i][x%pk[i]]%pk[i]*calc(x/p[i],i);
}

曾经的代码

他可以这样变强：

ll mul(ll x,int i)
{
    if(!x)return 1;ll bas = pk[i]-1;
    if(p[i]==2)bas=(pk[i]==4)?3:1;
    return 1ll*fastpow(bas,x/pk[i],pk[i])*f[i][x%pk[i]]%pk[i]*mul(x/p[i],i)%pk[i];
}

变强的代码

所以说对于质数$p$（$p$不为$2$）与他的整数次幂$pk$存在$\prod _{i=1,i⊥p}^{pk}i ≡ pk-1 \; (mod\;pk)$？

先证一下这个：$$(p-2)!≡1\;(mod\;p)$$

这个可好证了……

先特判$3$时显然成立。

因为$1$的逆元是自己，$-1*-1=1$即$(p-1)*(p-1)≡1\;(mod\;p)$，所以这两个的逆元都是自己。

剩下的$2,3……p-2$共$p-3$个数，在$p$不为$2$时能两两一对组成逆元。

所以这一堆乘一起为$1$。

所以$(p-2)!≡1\;(mod\;p)$。

现在我们来证$\prod _{i=1,i⊥p}^{pk}i ≡ pk-1 \; (mod\;pk)$。

证法和上面那个差不多。

所有$p$的倍数没有逆元，$1$和$pk-1$的逆元是自己，剩下两两配对都是$1$。

upd:刚才证完了非$2$部分，现在来证一下$2$的部分。

刚才的伪证是在只有$1$和$pk-1$的逆元是自己的前提下，而当$p=2$且$k>2$时，$2^{n-1}-1$和$2^{n-1}+1$的逆元也是自己。

很巧的是这四个数乘一起得$1$。

（应该没有其他的锅了吧……）

转载于:https://www.cnblogs.com/LiGuanlin1124/p/11128840.html