【BZOJ3817/UOJ42】Sum（类欧）

【BZOJ3817/UOJ42】Sum（类欧）

题面

BZOJ
UOJ

题解

令(x=sqrt r)，那么要求的式子是$$sum_{d=1}^n(-1)^{[dx]}$$
不难发现，对于每个(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(1-2*([dx]-2*[frac{dx}{2}])) ]

化简之后得到了

[n+sum_{d=1}^n (4*[frac{dx}{2}]-2*[dx]) ]

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

[egin{aligned}sum_{d=1}^n[kd]&=[d*(k-[frac{ax+b}{c}])]+d*[frac{ax+b}{c}]end{aligned} ]

显然后面那一半是可以直接计算的，而前面这一半令(k'=k-[frac{ax+b}{c}])
我们可以递归处理。
2.(k<1)
转化到一个坐标系上面去，我们要求的东西本质上就是有一条直线(y=kx)，要求解(xin[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-2*((n+1)/2));
			else printf("%lld
",n);
		}
		else printf("%lld
",n+Calc(1,0,2,n)*4-Calc(1,0,1,n)*2);
	}
	return 0;
}