类欧几里得算法学习笔记

来更了

题意

求类似

\[ \sum_{i=0}^n i^{k_1}\left\lfloor \frac{ai+b}{c} \right\rfloor ^{k_2} \]

的东西

简单情况

\[ f(a,b,c,n)=\sum_{i=0}^n \left\lfloor \frac{ai+b}{c} \right\rfloor \]

概述

这其实相当于求一条直线下的整点数

考虑到直线斜率非常大的时候,可以快速转化为斜率 \(<1\) 的情况,并且规模可以缩小到原先的至多 \(\frac{1}{2}\)

而直线斜率小的时候可以用对称不改变规模转化为大的情况

这样总次数就是 \(\mathcal O(\log c)\)

做法

特判

\(a=0\) 时,直接算

\(a\ge c\)\(b\ge c\)

首先有

\[ \begin{align} f(a,b,c,n) = & \sum_{i=0}^n \left\lfloor \frac{ai+b}{c} \right\rfloor \\ = & \sum_{i=0}^n \left( \left\lfloor \frac{a}{c} \right\rfloor i + \left\lfloor \frac{b}{c} \right\rfloor +\left\lfloor \frac{(a\bmod c)i+(b\bmod c)}{c} \right\rfloor\right) \\ = & \frac{n(n+1)}{2} \left\lfloor \frac{a}{c} \right\rfloor + (n+1) \left\lfloor \frac{b}{c} \right\rfloor + f(a\bmod c, b\bmod c, c, n) \end{align} \]

前面的两项可以直接计算,接下来只考虑 \(a<c,b<c\) 的情况

\(a,b<c\)

枚举一个 \(j\)

\[ \sum_{i=0}^n \left\lfloor \frac{ai+b}{c} \right\rfloor = \sum_{i=0}^n \sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1} 1 \]

\(m=\left\lfloor \frac{a\times n+b}{c} \right\rfloor\)

于是

\[ \begin{align} & \sum_{i=0}^n \left\lfloor \frac{ai+b}{c} \right\rfloor \\ = & \sum_{i=0}^n \sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1} 1 \\ = & \sum_{i=0}^n \sum_{j=0}^{m-1} [j+1\le \left\lfloor \frac{ai+b}{c} \right\rfloor] \\ = & \sum_{i=0}^n \sum_{j=0}^{m-1} [i\ge\left\lceil \frac{cj+c-b}{a} \right\rceil] \\ = & \sum_{i=0}^n \sum_{j=0}^{m-1} [i\ge\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor+1] \\ = & \sum_{j=0}^{m-1} \sum_{i=0}^n [i\ge\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor+1] \\ = & \sum_{j=0}^{m-1} n-\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor \\ = & mn-f(c,c-b-1,a,m-1) \end{align} \]

里面转化为关于 \(i\) 的条件自己推一推就好了,注意直接做会得到上取整

主要套路就是交换两个 \(\sum\)

复杂度

观察第 \(a,c\) 两个参数,\((a,c)\to(a\bmod c,c)\to(c, a\bmod c)\)

于是复杂度和欧几里得算法相同

代码

可以自然溢出,最终答案在范围内即可

听说自然溢出要 unsigned,那代码就仅供参考了

1
2
3
4
5
6
7
ll f(int a, int b, int c, int n){
ll ans=(ll)n*(n+1)/2*(a/c)+(ll)(n+1)*(b/c);
a%=c, b%=c;
if(!a) return ans;
int m=((ll)n*a+b)/c;
return ans+(ll)m*n-f(c, c-b-1, a, m-1);
}

另一种情况

UOJ #42. 【清华集训2014】Sum

\[ \sum_{d = 1}^{n}{(-1)^{\left\lfloor \sqrt{d\times r\times d} \right\rfloor}} \]

做法

显然只要求 \(\sum_{d=1}^n (\left\lfloor d \sqrt r \right\rfloor \bmod 2)\) 就好了

上式即 \(\sum_{d=1}^n \left\lfloor d \sqrt r \right\rfloor -2 \sum_{d=1}^n \left\lfloor \frac{d \sqrt r}{2} \right\rfloor\)

\[ f(a,b,c,n) = \sum_{i=1}^n \left\lfloor i \times \frac{a\sqrt r + b}{c}\right\rfloor \]

\(\frac{a\sqrt r +b}{c}\ge 1\)

\[ \begin{align} f(a,b,c,n) & = \sum_{i=1}^n \left\lfloor i \times \frac{a\sqrt r + b}{c}\right\rfloor \\ & = \sum_{i=1}^n \left\lfloor i \left\lfloor \frac{a\sqrt r + b}{c}\right\rfloor - i\left(\frac{a\sqrt r + b}{c} - \left\lfloor \frac{a\sqrt r + b}{c} \right\rfloor\right) \right\rfloor \end{align} \]

第一部分显然可以提出直接计算,于是我们只考虑 \(\frac{a\sqrt r + b}{c}<1\) 的情况

\(\frac{a\sqrt r + b}{c}<1\)

\(m=\left\lfloor n\times \frac{a\sqrt r + b}{c} \right\rfloor\)

\[ \begin{align} f(a,b,c,n) & = \sum_{i=1}^n \sum_{j=1}^m [j\le \left\lfloor i\times \frac{a\sqrt r + b}{c} \right\rfloor] \\ & = \sum_{j=1}^m \sum_{i=1}^n [i\ge \left\lceil j\times \frac{ac\sqrt r-bc}{a^2r-b^2} \right\rceil] \\ & = \sum_{j=1}^m n-\left\lceil j\times \frac{ac\sqrt r-bc}{a^2r-b^2} \right\rceil+1 \end{align} \]

其中第二行可以通过移项和分母有理化得到,不展开说明了

注意这里上取整解决不了,可以发现上取整和下取整 \(+1\) 只有在整数时不同,这里要是整数必然有 \(r\) 是完全平方数

  • 可以发现当 \(\frac{a^2r-b^2}{\gcd(ac\sqrt r-bc, a^2r-b^2)} \mid j\) 时会出现整数,特判即可

  • 另一种解决方法是当 \(\sqrt r\) 是整数的时候可以很轻易地计算原式,在一开始判掉

复杂度

这里 \(f(a,b,c,n)\to f(ac,-bc,a^2r-b^2,m)\),不妨令 \(x=\frac{a \sqrt r+b}{c}<1\),这可以看做 \((n,x)\to (n', \frac{1}{x})\)

其中 \(n'=\lfloor n\times x \rfloor\)显然每次严格减小,但是没什么用

\(\alpha=\frac{\sqrt{5}-1}{2}\)

  • \(0<x\le \alpha\),有 \(n'=\lfloor n\times x \rfloor \le \alpha n\)

  • \(\alpha<x<1\),有\(1<\frac{1}{x}<\frac{1}{\alpha}=\alpha+1\),于是 \(n''=\lfloor n'\times (\frac{1}{x}-1)\rfloor<\alpha n'<\alpha n\)

于是递归至多两层后 \(n\) 会变成至多 \(\alpha\) 倍,层数至多是 \(2\log_{\frac{1}{\alpha}} n+\mathcal O(1)\)

总复杂度 \(\mathcal O(\log^2 n)\)

但是这里数字会变得很大,每次约分一下就好了,我也不知道为什么能过,实测int可过

代码

不太懂自然溢出,反正能过

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>

using namespace std;

int T, r, n;
double sr;
inline int gcd(int a, int b){ return b?gcd(b, a%b):a;}
int f(int a, int b, int c, int n){
if(n==0) return 0;
int g=gcd(gcd(a, b), c);
a/=g, b/=g, c/=g;
int t=(a*sr+b)/c;
b-=t*c;
int m=n*(a*sr+b)/c, ssr=sr, ans=0;
if(ssr*ssr==r && a*a*r-b*b) ans=m/((a*a*r-b*b)/gcd(a*c*ssr-b*c, a*a*r-b*b));
return n*(n+1)/2*t+m*n-f(a*c, -b*c, a*a*r-b*b, m)+ans;
}
int main() {
scanf("%d", &T);
while(T--){
scanf("%d%d", &n, &r), sr=sqrt(r);
printf("%d\n", n-2*(f(1, 0, 1, n)-2*f(1, 0, 2, n)));
}
return 0;
}

一般的普通情况

LOJ #138. 类欧几里得算法

\[ \sum_{i=0}^n i^{k_1}\left\lfloor \frac{ai+b}{c} \right\rfloor ^{k_2} \]

做法

咕咕咕

这里钦定 \(0^0=1\)

特判

首先,若 \(a=0\)\(k_2=0\)

问题转化为 \(\lambda\sum\limits_{i=0}^n i^{k_1}\) 的形式,其中 \(\lambda\) 是一个易求的常数

这可以直接插值

\(n<0\) 时返回 \(0\)

但是鲁棒性较好的代码不需要某些特判

\(a\ge c\)\(b\ge c\)

\[ \begin{align} f(a,b,c,n,k_1,k_2) &= \sum_{i=0}^n i^{k_1}\left\lfloor \frac{ai+b}{c} \right\rfloor ^{k_2} \\ &= \sum_{i=0}^n i^{k_1}\left( \left\lfloor \frac{a}{c} \right\rfloor i + \left\lfloor \frac{b}{c} \right\rfloor +\left\lfloor \frac{(a\bmod c)i+(b\bmod c)}{c} \right\rfloor \right) ^{k_2} \end{align} \]

前面的两项可以插值

三项式暴力展开,形式是若干个 \(\lambda \times f(a\bmod c, b\bmod c, c, n, k_1', k_2')\)其中 \(\lambda\) 是一个易求的常数

可以发现其中 \(k_1'+k_2'\le k_1+k_2\)

\(a,b<c\)

我们可以把 \(x^{k_2}\) 转化为 \(\sum\limits_{j=0}^{x-1}\left((j+1)^{k_2}-x^{k_2}\right)\)

同样令 \(m=\left\lfloor \frac{a\times n+b}{c} \right\rfloor\)

\[ \begin{align} f(a,b,c,n,k_1,k_2) &= \sum_{i=0}^n i^{k_1}\left\lfloor \frac{ai+b}{c} \right\rfloor ^{k_2} \\ &= \sum_{i=0}^n i^{k_1} \sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1} \left((j+1)^{k_2}-j^{k_2}\right) \\ &= \sum_{i=0}^n i^{k_1}\sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right) [j\le \left\lfloor \frac{ai+b}{c} \right\rfloor-1] \\ &= \sum_{i=0}^n i^{k_1}\sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right) [i>\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor] \\ &= \sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right) \sum_{i=0}^n i^{k_1}[i>\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor] \\ &= \sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right) \left(\sum_{i=0}^n i^{k_1} -\sum_{i=0}^{\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor}i^{k_1}\right)\\ &= \left(\sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right)\right) \times \left( \sum_{i=0}^n i^{k_1} \right) - \sum_{j=0}^{m-1} \left((j+1)^{k_2}-j^{k_2}\right) \sum_{i=0}^{\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor}i^{k_1} \end{align} \]

前面部分可以插值,后面的 \(\left((j+1)^{k_2}-j^{k_2}\right)\) 是一个关于 \(j\)\(k_2-1\) 次多项式 \(A(x)\)\(\sum\limits_{i=0}^{\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor}i^{k_1}\) 是一个关于 \(\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor\)\(\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor\)\(k_1+1\) 次多项式 \(B(x)\)

转化为 \(\sum\limits_{j=0}^{m-1} A(j) \times B\left(\left\lfloor \frac{cj+c-b-1}{a} \right\rfloor\right)\)

枚举两个多项式的项 \(x^{k_1'}\)\(x^{k_2'}\) 有贡献是 \([x^{k_1'}]A(x) \times [x^{k_2'}] B(x) \times f(c, c-b-1, a, m-1, k_1', k_2')\);

并且也有 \(k_1'+k_2'\le k_1+k_2\)

复杂度

可以发现每层的 \(n,a,b,c\) 都是相同的,记忆化 \(k1,k2\) 两维

注意到 \((a,c)\to (a\bmod c,c) \to(c,a\bmod c)\),层数是 \(\mathcal O(\log c)\)

每次计算是 \(\mathcal O\left((k_1+k_2)^2\right)\)

所以总复杂度就是 \(\mathcal O\left((k_1+k_2)^4\log c\right)\)

跑得很快就是了

代码

代码中的多项式 \(A_{k_2}(x) = \sum\limits_{i=0}^x \left((i+1)^{k_2}-i^{k_2}\right)\),和上面的 \(A(x)\) 不同

\(B_{k_1}(x) = \sum\limits_{i=0}^x i^{k_1}\)

  • update on 2019-1-7 19:24:20:

之前的代码记忆化层数开了60,在Luogu 上的模板题里被 Sooke 卡到了66层qaq..

还被卡常了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include<string.h>

using namespace std;
#define ll long long

const int N = 12, P = 1000000007;
int n, a, b, c, k1, k2, T, Y[N], C[N][N], F[75][N][N];
struct Poly{
int a[N];
inline Poly(){ memset(a, 0, sizeof a);}
inline Poly operator +(const Poly &rhs)const{
Poly ans;
for(int i=0; i<N; ++i) ans.a[i]=(a[i]+rhs.a[i])%P;
return ans;
}
inline Poly operator *(const Poly &rhs)const{
Poly ans;
for(int i=0; i<N; ++i) for(int j=0; i+j<N; ++j)
ans.a[i+j]=(ans.a[i+j]+(ll)a[i]*rhs.a[j])%P;
return ans;
}
inline int calc(int n, int k){
int ans=a[k+1];
for(int i=k; ~i; --i) ans=((ll)ans*n+a[i])%P;
return ans;
}
} A[N], B[N];
inline int Pow(ll x, int y=P-2){
int ans=1;
for(x+=P; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
return ans;
}
void solve(Poly &A, int n){
for(int i=0; i<=n; ++i){
Poly t;
t.a[0]=Y[i];
for(int j=0; j<=n; ++j) if(j!=i) t.a[0]=(ll)t.a[0]*Pow(i-j)%P;
for(int j=0; j<=n; ++j) if(j!=i){
Poly d;
d.a[0]=P-j, d.a[1]=1;
t=t*d;
}
A=A+t;
}
}
int f(int n, int a, int b, int c, int k1, int k2, int dep=0){
if(~F[dep][k1][k2]) return F[dep][k1][k2];
int ans=0;
if(!a) ans=(ll)B[k1].calc(n, k1)*Pow(b/c, k2)%P;
else if(!k2) ans=B[k1].calc(n, k1);
else if(a>=c || b>=c){
int x=a/c, y=b/c;
a%=c, b%=c;
for(int i=0, I=1; i<=k2; ++i, I=(ll)I*x%P)
for(int j=0, J=I; i+j<=k2; ++j, J=(ll)J*y%P)
ans=(ans+(ll)f(n, a, b, c, k1+i, k2-i-j, dep+1)*C[k2][i]%P*C[k2-i][j]%P*J)%P;
}
else{
int m=((ll)a*n+b)/c;
ans=(ll)A[k2].calc(m-1, k2)*B[k1].calc(n, k1)%P;
for(int j=0; j<k2; ++j) for(int i=0; i<=k1+1; ++i)
ans=(ans-(ll)f(m-1, c, c-b-1, a, j, i, dep+1)*C[k2][j]%P*B[k1].a[i])%P;
}
return F[dep][k1][k2]=(ans+P)%P;
}
int main() {
C[0][0]=1;
for(int i=1; i<=10; ++i) for(int j=0; j<=i; ++j) C[i][j]=(C[i-1][j]+(j?C[i-1][j-1]:0))%P;
for(int i=0; i<=10; ++i){
Y[0]=1-!i;
for(int j=1; j<=i+1; ++j) Y[j]=((ll)P+Y[j-1]+Pow(j+1, i)-Pow(j, i))%P;
solve(A[i], i+1);
Y[0]=!i;
for(int j=1; j<=i+1; ++j) Y[j]=(Y[j-1]+Pow(j, i))%P;
solve(B[i], i+1);
}

scanf("%d", &T);
while(T--){
scanf("%d%d%d%d%d%d", &n, &a, &b, &c, &k1, &k2);
memset(F, -1, sizeof F);
printf("%d\n", f(n, a, b, c, k1, k2));
}
return 0;
}
0%