类欧几里得算法

问题

求以下形式表达式的值。

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

简单情况

求

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

概述

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

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

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

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

做法

一

\(a=0\) 时，直接算。

二

若 \(a\ge c\) 或 \(b\ge c\)：

首先有

\[ \begin{aligned} f(a,b,c,n) = & \sum_{i=0}^n \left\lfloor \frac{a\times i+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{aligned} \]

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

三

若 \(a,b<c\)：

枚举一个 \(j\)

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

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

于是

\[ \begin{aligned} \sum_{i=0}^n \left\lfloor \frac{a\times i+b}{c} \right\rfloor = & \sum_{i=0}^n \sum_{j=0}^{\left\lfloor \frac{a\times i+b}{c} \right\rfloor-1} 1 \\ = & \sum_{i=0}^n \sum_{j=0}^{m-1} [j+1\le \left\lfloor \frac{a\times i+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{aligned} \]

其中第二行到第三行的变换自己推一推就好了，注意直接做会得到上取整。

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

复杂度

观察第 \(a,c\) 两个变量，\((a,c)\to(a\bmod c,c)\to(c, a\bmod c)\)

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

代码

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

听说自然溢出要 unsigned，代码仅供参考。

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);
}

一般情况

LOJ #138. 类欧几里得算法

\[ \sum_{i=0}^n i^{k_1}\left\lfloor \frac{a\times i+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\) 的情况。

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

一

若 \(a\ge c\) 或 \(b\ge c\)：

令

\[ \begin{aligned} f(a,b,c,n,k_1,k_2) &= \sum_{i=0}^n i^{k_1}\left\lfloor \frac{a\times i+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{aligned} \]

前面的两项可以插值。

三项式暴力展开，形式是若干个 \(\lambda \times f(a\bmod c, b\bmod c, c, n, k_1', k_2')\)，且其中 \(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}-j^{k_2}\right)\)。

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

\[ \begin{aligned} f(a,b,c,n,k_1,k_2) &= \sum_{i=0}^n i^{k_1}\left\lfloor \frac{a\times i+b}{c} \right\rfloor ^{k_2} \\ &= \sum_{i=0}^n i^{k_1} \sum_{j=0}^{\left\lfloor \frac{a\times i+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{a\times i+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{aligned} \]

前面部分可以插值，后面的 \(\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\) 都是相同的，记忆化 \(k_1,k_2\) 两维。

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

每次计算是 \(O\left((k_1+k_2)^2\right)\)，所以总复杂度就是 \(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。

还被卡常了。

#include<cstdio>
#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;
}

一种特殊情况

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{aligned} 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{aligned} \]

第一部分可以提出直接计算，于是我们只考虑 \(\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{aligned} 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{aligned} \]

第二行可以通过移项和分母有理化得到。

注意这里上取整解决不了，可以发现上取整和下取整 \(+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+O(1)\)。

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

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

代码

不太懂自然溢出，反正能过。

#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;
}