【模板】类欧几里得算法

强烈不建议看洛谷的模板,不是很好。

【模板】类欧几里得算法(洛谷)
【模板】类欧几里得算法(LOJ)

类欧几里得算法可以干什么

其实类欧几里得算法跟欧几里得算法一点关系都没有,只是复杂度相似。

他可以解决如下问题:

已知 n , a , b , c , k 1 , k 2 n,a,b,c,k_1,k_2 n,a,b,c,k1,k2,求 ∑ i = 0 n i k 1 ⌊ a n + b c ⌋ k 2 \sum\limits_{i=0}^ni^{k_1}\left\lfloor\dfrac{an+b}{c}\right\rfloor^{k_2} i=0nik1can+bk2,我们规定 0 0 = 1 0^0=1 00=1

你会发现洛谷的模板仅仅只有 k 1 + k 2 ⩽ 2 k_1+k_2\leqslant 2 k1+k22 的情况,所以建议配合 LOJ 模板学习。

类欧几里得算法的推导

我们首先观察 a = 0 a=0 a=0 的情况,式子变为 ⌊ b c ⌋ k 2 ∑ i = 0 n i k 1 \left\lfloor\dfrac{b}{c}\right\rfloor^{k_2}\sum\limits_{i=0}^n i^{k_1} cbk2i=0nik1,用拉格朗日插值计算 ∑ i = 0 n i k 1 \sum\limits_{i=0}^ni^{k_1} i=0nik1 即可,对于 ⌊ a n + b c ⌋ = 0 \left\lfloor\dfrac{an+b}{c}\right\rfloor=0 can+b=0 的情况同理。

接下来考虑 a ⩾ c a\geqslant c ac b ⩾ c b\geqslant c bc 的情况,我们要尽可能缩小问题的规模。

a ′ = a   m o d   c a^\prime=a\bmod c a=amodc,使用二项式定理有:
∑ i = 0 n i k 1 ( ⌊ a c ⌋ i + ⌊ a ′ i + b c ⌋ ) k 2 = ∑ i = 0 n i k 1 ∑ j = 0 k 2 ( k 2 j ) ( ⌊ a c ⌋ i ) j ⌊ a ′ i + b c ⌋ k 2 − j = ∑ j = 0 k 2 ( k 2 j ) ⌊ a c ⌋ j ∑ i = 0 n i k 1 + j ⌊ a ′ i + b c ⌋ k 2 − j \begin{aligned} &\sum\limits_{i=0}^ni^{k_1}\left(\left\lfloor\dfrac{a}{c}\right\rfloor i+\left\lfloor\dfrac{a^\prime i+b}{c}\right\rfloor\right)^{k_2}\\ =&\sum\limits_{i=0}^ni^{k_1}\sum\limits_{j=0}^{k_2} \binom{k_2}{j}\left(\left\lfloor\dfrac{a}{c}\right\rfloor i\right)^j\left\lfloor\dfrac{a^\prime i+b}{c}\right\rfloor^{k_2-j}\\ =&\sum\limits_{j=0}^{k_2}\binom{k_2}{j}\left\lfloor\dfrac{a}{c}\right\rfloor^j\sum\limits_{i=0}^ni^{k_1+j}\left\lfloor\dfrac{a^\prime i+b}{c}\right\rfloor^{k_2-j} \end{aligned} ==i=0nik1(cai+cai+b)k2i=0nik1j=0k2(jk2)(cai)jcai+bk2jj=0k2(jk2)caji=0nik1+jcai+bk2j

如果我们定义 φ ( n , a , b , c , k 1 , k 2 ) = ∑ i = 0 n i k 1 ⌊ a i + b c ⌋ k 2 \varphi(n,a,b,c,k_1,k_2)=\sum\limits_{i=0}^ni^{k_1}\left\lfloor\dfrac{ai+b}{c}\right\rfloor^{k_2} φ(n,a,b,c,k1,k2)=i=0nik1cai+bk2,则上面推导出来的即为:
φ ( n , a , b , c , k 1 , k 2 ) = ∑ j = 0 k 2 ( k 2 j ) ⌊ a c ⌋ j φ ( n , a   m o d   c , b , c , k 1 + j , k 2 − j ) \varphi(n,a,b,c,k_1,k_2)=\sum\limits_{j=0}^{k_2}\binom{k_2}{j}\left\lfloor\dfrac{a}{c}\right\rfloor^j\varphi(n,a\bmod c,b,c,k_1+j,k_2-j) φ(n,a,b,c,k1,k2)=j=0k2(jk2)cajφ(n,amodc,b,c,k1+j,k2j)

接下来令 b ′ = b   m o d   c b^\prime =b\bmod c b=bmodc,继续考虑:
∑ i = 0 n i k 1 ( ⌊ b c ⌋ + ⌊ a i + b ′ c ⌋ ) k 2 = ∑ i = 0 n i k 1 ∑ j = 0 k 2 ( k 2 j ) ⌊ b c ⌋ j ⌊ a i + b ′ c ⌋ k 2 − j = ∑ j = 0 k 2 ( k 2 j ) ⌊ b c ⌋ j ∑ i = 0 n i k 1 ⌊ a i + b ′ c ⌋ k 2 − j \begin{aligned} &\sum\limits_{i=0}^ni^{k_1}\left(\left\lfloor\dfrac{b}{c}\right\rfloor+\left\lfloor\dfrac{ai+b^\prime}{c}\right\rfloor\right)^{k_2}\\ =&\sum\limits_{i=0}^ni^{k_1}\sum\limits_{j=0}^{k_2}\binom{k_2}{j}\left\lfloor\dfrac{b}{c}\right\rfloor^j\left\lfloor\dfrac{ai+b^\prime}{c}\right\rfloor^{k_2-j}\\ =&\sum\limits_{j=0}^{k_2}\binom{k_2}{j}\left\lfloor\dfrac{b}{c}\right\rfloor^j\sum\limits_{i=0}^ni^{k_1}\left\lfloor\dfrac{ai+b^\prime}{c}\right\rfloor^{k_2-j} \end{aligned} ==i=0nik1(cb+cai+b)k2i=0nik1j=0k2(jk2)cbjcai+bk2jj=0k2(jk2)cbji=0nik1cai+bk2j

所以有:
φ ( n , a , b , c , k 1 , k 2 ) = ∑ j = 0 k 2 ( k 2 j ) ⌊ b c ⌋ j φ ( n , a , b   m o d   c , c , k 1 , k 2 − j ) \varphi(n,a,b,c,k_1,k_2)=\sum\limits_{j=0}^{k_2}\binom{k_2}{j}\left\lfloor\dfrac{b}{c}\right\rfloor^j\varphi(n,a,b\bmod c,c,k_1,k_2-j) φ(n,a,b,c,k1,k2)=j=0k2(jk2)cbjφ(n,a,bmodc,c,k1,k2j)

最后我们考虑 a < c    and    b < c a<c\;\text{and}\;b<c a<candb<c 的情况,在此之前我们要学习一个人类智慧式子(这式子看上去就很废话文学):
x k = ∑ i = 1 x i k − ( i − 1 ) k x^k=\sum\limits_{i=1}^{x}i^k-(i-1)^k xk=i=1xik(i1)k

定义 m = ⌊ a n + b c ⌋ m=\left\lfloor\dfrac{an+b}{c}\right\rfloor m=can+b,此时我们直接将 ⌊ a i + b c ⌋ \left\lfloor\dfrac{ai+b}{c}\right\rfloor cai+b 尝试变形。
∑ i = 0 n i k 1 ∑ j = 0 ⌊ a i + b c ⌋ − 1 [ ( j + 1 ) k 2 − j k 2 ] = ∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 n i k 1 [ ⌊ a i + b c ⌋ ⩾ j + 1 ] = ∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 n i k 1 [ i > ⌊ c j + c − b − 1 a ⌋ ] = ∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 n i k 1 − ∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 ⌊ c j + c − b − 1 a ⌋ i k 1 \begin{aligned} &\sum\limits_{i=0}^ni^{k_1}\sum\limits_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}[(j+1)^{k_2}-j^{k_2}]\\ =&\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^ni^{k_1}\left[\left\lfloor\dfrac{ai+b}{c}\right\rfloor\geqslant j+1\right]\\ =&\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^ni^{k_1}\left[i>\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor\right]\\ =&\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^ni^{k_1}-\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^{\lfloor\frac{cj+c-b-1}{a}\rfloor}i^{k_1} \end{aligned} ===i=0nik1j=0cai+b1[(j+1)k2jk2]j=0m1[(j+1)k2jk2]i=0nik1[cai+bj+1]j=0m1[(j+1)k2jk2]i=0nik1[i>acj+cb1]j=0m1[(j+1)k2jk2]i=0nik1j=0m1[(j+1)k2jk2]i=0acj+cb1ik1

前面一半很好解决,我们直接看后面一半。令 f ( x ) = ( x + 1 ) k 2 − x k 2 , g ( x ) = ∑ i = 0 x i k 1 f(x)=(x+1)^{k_2}-x^{k_2},g(x)=\sum\limits_{i=0}^x i^{k_1} f(x)=(x+1)k2xk2,g(x)=i=0xik1。所以后半部分推导为:
∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 ⌊ c j + c − b − 1 a ⌋ i k 1 = ∑ i = 0 m − 1 ∑ j = 0 k 2 − 1 f ( i ) i j ∑ k = 0 k 1 + 1 g ( k ) ⌊ c j + c − b − 1 a ⌋ k = ∑ j = 0 k 2 − 1 f ( i ) ∑ k = 0 k 1 + 1 g ( k ) ∑ i = 0 m − 1 i j ⌊ c j + c − b − 1 a ⌋ k \begin{aligned} &\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^{\lfloor\frac{cj+c-b-1}{a}\rfloor}i^{k_1}\\ =&\sum\limits_{i=0}^{m-1}\sum\limits_{j=0}^{k_2-1}f(i)i^j\sum\limits_{k=0}^{k_1+1}g(k)\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor^k\\ =&\sum\limits_{j=0}^{k_2-1}f(i)\sum\limits_{k=0}^{k_1+1}g(k)\sum\limits_{i=0}^{m-1}i^j\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor^k \end{aligned} ==j=0m1[(j+1)k2jk2]i=0acj+cb1ik1i=0m1j=0k21f(i)ijk=0k1+1g(k)acj+cb1kj=0k21f(i)k=0k1+1g(k)i=0m1ijacj+cb1k

所以有:
φ ( n , a , b , c , k 1 , k 2 ) = ∑ j = 0 m − 1 [ ( j + 1 ) k 2 − j k 2 ] ∑ i = 0 n i k 1 − ∑ j = 0 k 2 − 1 f ( i ) ∑ k = 0 k 1 + 1 g ( k ) φ ( m − 1 , c , c − b − 1 , a , j , k ) \varphi(n,a,b,c,k_1,k_2)=\sum\limits_{j=0}^{m-1}[(j+1)^{k_2}-j^{k_2}]\sum\limits_{i=0}^ni^{k_1}-\sum\limits_{j=0}^{k_2-1}f(i)\sum\limits_{k=0}^{k_1+1}g(k)\varphi(m-1,c,c-b-1,a,j,k) φ(n,a,b,c,k1,k2)=j=0m1[(j+1)k2jk2]i=0nik1j=0k21f(i)k=0k1+1g(k)φ(m1,c,cb1,a,j,k)

正确使用插值可以在 O ( T k 4 log ⁡ n ) \mathcal{O}(Tk^4\log n) O(Tk4logn) 的时间内解决这个问题。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=15;
const int mod=1e9+7;
struct node{
    int a[N][N];
    void init(){memset(a,0,sizeof(a));}
};
int read(){
    char ch=getchar();
    int w=0,s=1;
    while(ch>'9'||ch<'0'){
        if(ch=='-') s=-1;
        ch=getchar();
    }
    while(ch<='9'&&ch>='0'){
        w=(w<<1)+(w<<3)+ch-'0';
        ch=getchar();
    }
    return w*s;
}
int Mul(int x,int y){return (x%mod*y%mod)%mod;}
int Add(int x,int y){return (x+y)%mod;}
int Dec(int x,int y){return (x-y+mod)%mod;}
int Pow(int a,int k){
    if(k==0) return 1;
    int ans=1;
    while(k){
        if(k&1) ans=Mul(ans,a);
        a=Mul(a,a);
        k>>=1;
    }
    return ans;
}
int inv(int x){return Pow(x,mod-2);}
int x[N],y[N],sum[N][N],C[N][N];
node solve(int n,int a,int b,int c){
    node ans;
    ans.init();
    if(a==0||(a*n+b)/c==0){
        for(register int k1=0;k1<=10;k1++){
            int tmp=0,now=1;
            for(register int i=0;i<=k1+1;i++){
                tmp=Add(tmp,Mul(now,sum[k1][i]));
                now=Mul(now,n);
            }
            int lst=(a*n+b)/c;
            now=1;
            for(register int k2=0;k1+k2<=10;k2++){
                ans.a[k1][k2]=Mul(now,tmp);
                now=Mul(now,lst);
            }
        }
        return ans;
    }
    if(a>=c){
        node tmp=solve(n,a%c,b,c);
        for(register int k1=0;k1<=10;k1++){
            for(register int k2=0;k1+k2<=10;k2++){
                int now=1;
                for(register int i=0;i<=k2;i++){
                    ans.a[k1][k2]=Add(ans.a[k1][k2],Mul(C[k2][i],Mul(now,tmp.a[k1+i][k2-i])));
                    now=Mul(now,a/c);
                }
            }
        }
        return ans;
    }
    if(b>=c){
        node tmp=solve(n,a,b%c,c);
        for(register int k1=0;k1<=10;k1++){
            for(register int k2=0;k1+k2<=10;k2++){
                int now=1;
                for(register int i=0;i<=k2;i++){
                    ans.a[k1][k2]=Add(ans.a[k1][k2],Mul(C[k2][i],Mul(now,tmp.a[k1][k2-i])));
                    now=Mul(now,b/c);
                }
            }
        }
        return ans;
    }
    int m=(a*n+b)/c;
    node tmp=solve(m-1,c,c-b-1,a);
    for(register int k1=0;k1<=10;k1++){
        int lst=0,now=1;
        for(register int i=0;i<=k1+1;i++){
            lst=Add(lst,Mul(now,sum[k1][i]));
            now=Mul(now,n);
        }
        for(register int k2=0;k1+k2<=10;k2++){
            ans.a[k1][k2]=Mul(Pow(m,k2),lst);
            for(register int i=0;i<=k2-1;i++){
                for(register int j=0;j<=k1+1;j++){
                    ans.a[k1][k2]=Add(ans.a[k1][k2],Dec(mod,Mul(C[k2][i],Mul(sum[k1][j],tmp.a[i][j]))));
                }
            }
        }
    }
    return ans;
}
void lagrange(int n,int k){
    int p[N],q[N];
    memset(p,0,sizeof(p));
    p[0]=1;
    for(register int i=1;i<=n;i++){
        for(register int j=i-1;j>=0;j--){
            p[j+1]=Add(p[j+1],p[j]);
            p[j]=Dec(mod,Mul(p[j],x[i]));
        }
    }
    for(register int i=1;i<=n;i++){
        memset(q,0,sizeof(q));
        for(register int j=n-1;j>=0;j--) q[j]=Add(p[j+1],Mul(q[j+1],x[i]));
        int now=1;
        for(register int j=1;j<=n;j++) if(j!=i) now=Mul(now,Dec(x[i],x[j]));
        now=inv(Add(now,mod));
        for(register int j=0;j<=n;j++) q[j]=Mul(q[j],now);
        for(register int j=0;j<=n;j++) sum[k][j]=Add(sum[k][j],Mul(q[j],y[i]));
    }
}
signed main(){
    for(register int i=0;i<=10;i++){
        memset(y,0,sizeof(y));
        memset(x,0,sizeof(x));
        y[0]=Pow(0,i);
        for(register int j=1;j<=i+2;j++){
            y[j]=Add(y[j-1],Pow(j,i));
            x[j]=j;
        }
        lagrange(i+2,i);
    }
    for(register int i=0;i<=10;i++){
        C[i][0]=1;
        for(register int j=1;j<=i;j++){
            C[i][j]=C[i-1][j-1]+C[i-1][j];
        }
    }
    int T=read();
    while(T--){
        int n=read(),a=read(),b=read(),c=read(),k1=read(),k2=read();
        printf("%lld\n",solve(n,a,b,c).a[k1][k2]);
    }
}

本文参考:

【LOJ138】类欧几里得算法
【LOJ】#138. 类欧几里得算法
WeLikeStudying 的博客 - P5710

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值