【BZOJ3817/UOJ42】Sum（类欧）

题面

BZOJ
UOJ

题解

令\(x=\sqrt r\)，那么要求的式子是$$\sum_{d=1}^n(-1)$$
不难发现，对于每个\(d\)而言的取值只和\([dx]\)的奇偶性相关。
如果\(x\)是个整数，也就是\(r\)是完全平方数的时候，显然是可以直接算答案的。
计算答案的时候显然之和有几个奇数或者几个偶数相关（只要求一个另外一个就是补集）
比如说我们来求有几个是偶数，那么要满足的条件就是\([dx]=2*[\frac{dx}{2}]\)。
先重新写下式子，我们写成这个样子$$\sum_{d=1}^n(1-2*([dx]\ mod\ 2)$$
这个显然成立，当\([dx]\)是偶数的时候贡献是\(1\)，奇数的时候贡献是\(0\)。
而\([dx]\ mod\ 2\)可以写成减法的形式。所以原式可以写成

化简之后得到了

现在把模型转化一下，变得更加一般一点，那么要求解的东西是$$\sum _{d=1}^n[\frac{ax+b}{c}d]$$
听说这个玩意叫做类欧？类欧几里得算法。
令\(k=\frac{ax+b}{c}\)，那么式子可以化简为\(\sum[kd]\)
这里进行分类讨论
1.当\(k>=1\)的时候

显然后面那一半是可以直接计算的，而前面这一半令\(k'=k-[\frac{ax+b}{c}]\)
我们可以递归处理。
2.\(k<1\)
转化到一个坐标系上面去，我们要求的东西本质上就是有一条直线\(y=kx\)，要求解\(x\in[1,n]\)时，与\(x=n\)、\(x\)正半轴围成的三角形内部整点的个数。
我们把这个三角形逆时针旋转\(90\)度，再沿着\(y\)轴翻转过来，让长度为\(n\)的那条边靠着\(y\)轴，这样子翻转过来之后，\(k\)变成了倒数，\(n\)变成了\([nk]\)，然后拿矩形减去多出来的部分，不就是矩形减去\(k>=1\)的情况了吗？
再讨论一下取倒数的结果，斜率\(k=\frac{ax+b}{c}\)，取倒数之后\(k=\frac{c}{ax+b}\)，然后分母有理化一下\(k=\frac{c(ax-b)}{a^2r-b^2}\)。
总的时间复杂度是\(log\)级别的。

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
inline int read()
{
	int x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
double x;ll n,r;
ll Calc(ll a,ll b,ll c,ll n)
{
	if(!n)return 0;
	ll d=__gcd(a,__gcd(b,c));a/=d;b/=d;c/=d;
	ll k=(a*x+b)/c;
	if(!k)
	{
		ll N=(a*x+b)/c*n;
		return n*N-Calc(a*c,-b*c,a*a*r-b*b,N);
	}
	else
		return k*n*(n+1)/2+Calc(a,b-c*k,c,n);
}
int main()
{
	int T=read();
	while(T--)
	{
		n=read();r=read();x=sqrt(r);
		ll k=x;
		if(k*k==r)
		{
			if(k&1)printf("%lld\n",n-2*((n+1)/2));
			else printf("%lld\n",n);
		}
		else printf("%lld\n",n+Calc(1,0,2,n)*4-Calc(1,0,1,n)*2);
	}
	return 0;
}