强烈不建议看洛谷的模板,不是很好。
【模板】类欧几里得算法(洛谷)
【模板】类欧几里得算法(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=0∑nik1⌊can+b⌋k2,我们规定 0 0 = 1 0^0=1 00=1。
你会发现洛谷的模板仅仅只有 k 1 + k 2 ⩽ 2 k_1+k_2\leqslant 2 k1+k2⩽2 的情况,所以建议配合 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} ⌊cb⌋k2i=0∑nik1,用拉格朗日插值计算 ∑ i = 0 n i k 1 \sum\limits_{i=0}^ni^{k_1} i=0∑nik1 即可,对于 ⌊ a n + b c ⌋ = 0 \left\lfloor\dfrac{an+b}{c}\right\rfloor=0 ⌊can+b⌋=0 的情况同理。
接下来考虑 a ⩾ c a\geqslant c a⩾c 或 b ⩾ c b\geqslant c b⩾c 的情况,我们要尽可能缩小问题的规模。
令
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=0∑nik1(⌊ca⌋i+⌊ca′i+b⌋)k2i=0∑nik1j=0∑k2(jk2)(⌊ca⌋i)j⌊ca′i+b⌋k2−jj=0∑k2(jk2)⌊ca⌋ji=0∑nik1+j⌊ca′i+b⌋k2−j
如果我们定义
φ
(
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=0∑nik1⌊cai+b⌋k2,则上面推导出来的即为:
φ
(
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=0∑k2(jk2)⌊ca⌋jφ(n,amodc,b,c,k1+j,k2−j)
接下来令
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=0∑nik1(⌊cb⌋+⌊cai+b′⌋)k2i=0∑nik1j=0∑k2(jk2)⌊cb⌋j⌊cai+b′⌋k2−jj=0∑k2(jk2)⌊cb⌋ji=0∑nik1⌊cai+b′⌋k2−j
所以有:
φ
(
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=0∑k2(jk2)⌊cb⌋jφ(n,a,bmodc,c,k1,k2−j)
最后我们考虑
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=1∑xik−(i−1)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=0∑nik1j=0∑⌊cai+b⌋−1[(j+1)k2−jk2]j=0∑m−1[(j+1)k2−jk2]i=0∑nik1[⌊cai+b⌋⩾j+1]j=0∑m−1[(j+1)k2−jk2]i=0∑nik1[i>⌊acj+c−b−1⌋]j=0∑m−1[(j+1)k2−jk2]i=0∑nik1−j=0∑m−1[(j+1)k2−jk2]i=0∑⌊acj+c−b−1⌋ik1
前面一半很好解决,我们直接看后面一半。令
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)k2−xk2,g(x)=i=0∑xik1。所以后半部分推导为:
∑
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=0∑m−1[(j+1)k2−jk2]i=0∑⌊acj+c−b−1⌋ik1i=0∑m−1j=0∑k2−1f(i)ijk=0∑k1+1g(k)⌊acj+c−b−1⌋kj=0∑k2−1f(i)k=0∑k1+1g(k)i=0∑m−1ij⌊acj+c−b−1⌋k
所以有:
φ
(
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=0∑m−1[(j+1)k2−jk2]i=0∑nik1−j=0∑k2−1f(i)k=0∑k1+1g(k)φ(m−1,c,c−b−1,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