类欧几里得算法

问题

求以下形式表达式的值。

$$ \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;
}