题目链接

https://loj.ac/p/6052

题解

挺有趣的数论题，我的做法复杂度是 \(O(n^{\frac{3}{4}})\) 的，似乎比题解劣一些……
另外本题并不需要用到有关高斯整数的一些知识。（我也不会）

首先试图把有关复数的奇怪题意转化成我们比较熟悉的东西。
为了方便把题目中的 \(d\) 改成 \(-d\)，考虑 \((a+b{\rm i})(c-d{\rm i})=(ac+bd)+(ad-bc){\rm i}\)，那么由结果是一个整数可得 \(ad-bc=0\)，也就是 \((c,d)\) 与 \((a,b)\) 要成比例。
先不考虑 \(b\le 0\) 的情况。假设枚举 \(a,b\)，考虑哪些 \((c,d)\) 能够与之对答案产生贡献。
其需要满足两个条件：
(1) \((c,d)\) 与 \((a,b)\) 成比例。
设 \(g=\gcd(a,b),a=a'g,b=b'g\)，则该条件等价于存在 \(k\) 使 \(c=a'k,d=b'k\).
(2) \(ac+bd\le n\).
\(ac+bd=kg(a'^2+b'^2)\le n\)，因此满足上述两个条件的 \((c,d)\) 的对数为 \(\lfloor\frac{n}{g(a'^2+b'^2)}\rfloor\).
枚举 \(g,a',b'\)（为了方便使用 \(a,b\) 代替），可得答案等于

\[\sum^n_{g=1}g\sum^{\sqrt{n/g}}_{a=1,b=1}[\gcd(a,b)=1]a\cdot \lfloor\frac{n}{g(a^2+b^2)}\rfloor \]

于是我们把问题化成了我们熟悉的数论函数求和形式。

下一步是套路的莫比乌斯反演：

\[Ans=\sum^{\sqrt n}_{d=1}d\mu(d)\sum^{n/d^2}_{g=1}g\sum^{\sqrt{n/gd^2}}_{a=1,b=1}a\cdot \lfloor\frac{n/gd^2}{(a^2+b^2)}\rfloor \]

现在需要在低于线性的时间复杂度内求出它。
尽管有这么多求和符号，但是我们可以从中观察出一些特点：

  1. 看起来最棘手的是最内层，但是内层的求和结果似乎只和 \(n/gd^2\) 有关。那么假设我们把 \(\sum^{\sqrt m}_{a=1,b=1}\lfloor\frac{m}{a^2+b^2}\rfloor\) 看成关于 \(m\) 的函数 \(f(m)\)，那么外层枚举 \(d,g\) 后所有的求和项都形如 \(f(\lfloor\frac{n}{k}\rfloor)\)，其中 \(k\) 是整数。而这是一个我们相当喜欢的形式——它总共只有 \(O(\sqrt n)\) 种取值，而且假设我们能在 \(O(\sqrt m)\) 的时间内求出 \(f(m)\)，那么对每个 \(k\) 求出 \(f(\lfloor\frac{n}{k}\rfloor)\) 的总时间复杂度为 \(\sum^{\sqrt n}_{k=1}\sqrt\frac{n}{k}+\sum^{\sqrt n}_{k=1}\sqrt k=n^{\frac{1}{2}}\int^{n^{\frac{1}{2}}}_{0}t^{-\frac{1}{2}}{\rm d}t+\int^{n^{\frac{1}{2}}}_{0}t^{\frac{1}{2}}{\rm d}t=O(n^{\frac{3}{4}})\).
    进一步观察 \(\lfloor\frac{n}{k}\rfloor\) 这些数（称作“关键点”）形成的结构，其实相当于它们把 \([1,n]\) 分成了 \(O(\sqrt n)\) 段，每一段内 \(\lfloor\frac{n}{i}\rfloor\) 的值相等，关键点位于每一段的右端点。我们把这种划分方式称为“\(n\) 的划分”。设 \(x,y\) 位于同一段，\(m\) 是一个关键点，则有 \(\lfloor\frac{m}{x}\rfloor=\lfloor\frac{m}{y}\rfloor\). 换句话说，对于关键点 \(m\)，\(m\) 的划分中的关键点一定是 \(n\) 的划分中的关键点的子集，因此如果对某种信息我们对 \(n\) 的划分的每一段维护出了前缀和，那么就可以快速地查询 \(m\) 的划分中的一段的和。
    由此我们得到 \(O(\sqrt m)\) 时间求出 \(f(m)\) 的方法：设 \(c(m)\) 表示满足 \(a^2+b^2\le m\) 的 \(a\) 之和。显然我们可以做到 \(O(\sqrt m)\) 时间求出 \(c(m)\)——直接枚举 \(a\le \sqrt m\) 然后用 sqrt 函数 \(O(1)\) 累加进答案即可。现在我们对 \(m=\lfloor\frac{n}{k}\rfloor\) 的每种可能取值都求出了 \(c\) 函数的值，那么对相应的 \(f(m)\)，枚举 \(a^2+b^2\) 的取值在 \(m\) 的划分（注意不是 \(n\)）中所在的段 \(i\)，对结果的贡献是 \(\lfloor\frac{m}{i}\rfloor\cdot S(i)\)，\(S(i)\) 就是 \(a^2+b^2\) 落在这一段里的 \(a\) 值之和，可以由 \(n\) 的划分中对应的两段 \(c\) 值相减得到。
  2. 外层的枚举具有极其相似的结构：假设固定了 \(d\)，那么剩下的就是 \(\sum^{n/d^2}_{g=1}g\cdot f(\lfloor\frac{n/d^2}{g}\rfloor)\). 设函数 \(h(m)=\sum_{g=1}g\cdot f(\lfloor\frac{m}{g}\rfloor)\)，那我们需要的依然只是 \(h\) 函数在所有 \(\lfloor\frac{n}{k}\rfloor\) (\(k\) 为整数)处的点值。在求 \(h(m)\) 时使用整除分块，和前面不同的是此处求 \(h\) 的过程是递归的，像杜教筛一样写一个记忆化搜索，递归到所有 \(\lfloor\frac{m}{k}\rfloor\) 并通过它们的 \(f\) 和 \(h\) 函数值求出 \(h(m)\).
    而最后枚举 \(d\)，将 \(d\mu(d)h(\frac{n}{d^2})\) 累加进答案就可以了。

最后的答案乘以 \(2\)，再加上 \(b=0\) 的答案（就是 \(1\) 到 \(n\) 所有数的约数和）。

总时间复杂度 \(O(n^{\frac{3}{4}})\)，空间复杂度 \(O(\sqrt n)\)，或许可以用预处理优化到更快（？）

实际上这道题给我们的启发就是，对于任何奇形怪状的函数，如果我们只需要并保存它 \(\lfloor\frac{n}{k}\rfloor\) 处的取值，就可以在 \(O(n^{\frac{3}{4}})\) 的时间内完成某种类似于卷积的操作。设 \(f\)、\(g\) 是两个这样的函数，\(h(n)=\sum^n_{i=1}f(i)g(\lfloor\frac{n}{i}\rfloor)\)，那么就可以求出 \(h(n)\) 在所有关键点处的取值。或许这在某种意义上更接近于杜教筛的本质。

代码

#include<bits/stdc++.h>
#define llong long long
#define mkpr make_pair
#define x first
#define y second
#define iter iterator
#define riter reverse_iterator
#define y1 Lorem_ipsum_
#define tm dolor_sit_amet_
using namespace std;

inline int read()
{
	int x = 0,f = 1; char ch = getchar();
	for(;!isdigit(ch);ch=getchar()) {if(ch=='-') f = -1;}
	for(; isdigit(ch);ch=getchar()) {x = x*10+ch-48;}
	return x*f;
}

const int mxN = 1e5;
const int P = 1004535809;
llong n,TH;
int mu[mxN+3],pri[mxN+3]; bool isp[mxN+3]; int prin;
llong c[mxN*2+3],f[mxN*2+3],g[mxN*2+3];

void EulerSieve()
{
	isp[1] = true; mu[1] = 1;
	for(int i=2; i<=mxN; i++)
	{
		if(!isp[i]) {pri[++prin] = i; mu[i] = -1;}
		for(int j=1; j<=prin&&i*pri[j]<=mxN; j++)
		{
			isp[i*pri[j]] = true;
			if(i%pri[j]==0) {mu[i*pri[j]] = 0; break;}
			mu[i*pri[j]] = -mu[i];
		}
	}
}

llong getid(llong x) {return x<=TH?x:TH+n/x;}
llong getx(llong x) {return x<=TH?x:n/(x-TH);}

llong sum2(llong x) {x%=P; return (x*(x+1ll)/2ll)%P;}

llong getf(llong n)
{
	int id = getid(n); if(n==1ll) {return 0ll;} if(f[id]) return f[id];
	llong tmp = 0ll,tmp2 = 0ll;
	for(llong i=1ll; i<=n; )
	{
		llong j = n/(n/i),k = n/i; int idj = getid(j);
		f[id] = (f[id]+((n/i)%P)*(c[idj]-tmp+P))%P;
		if(i>1ll) {g[id] = (g[id]+getf(k)*(sum2(j)-tmp2))%P;}
		tmp = c[idj],tmp2 = sum2(j);
		i = j+1;
	}
	g[id] = (g[id]+f[id])%P;
	return f[id];
}

int main()
{
	EulerSieve();
	scanf("%lld",&n); TH = sqrt(n);
	for(llong i=1ll; i<=n; )
	{
		llong j = n/(n/i),k = n/i; int id = getid(i);
		for(llong a=1ll; a*a<=j; a++)
		{
			c[id] += a*(llong)sqrt(j-a*a);
			c[id] %= P;
		}
		i = j+1;
	}
	getf(n);
	llong ans = 0ll;
	for(llong i=1; i*i<=n; i++)
	{
		int id = getid(n/(i*i));
		ans = (ans+mu[i]*i*g[id]%P+P)%P;
	}
	ans = ans*2%P;
	for(llong i=1; i<=n; )
	{
		llong j = n/(n/i);
		ans = (ans+(sum2(j)-sum2(i-1)+P)*((n/i)%P))%P;
		i = j+1;
	}
	printf("%lld\n",ans);
	return 0;
}