题目链接
题解
挺有趣的数论题，我的做法复杂度是 (O(n^{frac{3}{4}})) 的，似乎比题解劣一些……
另外本题并不需要用到有关高斯整数的一些知识。（我也不会）
首先试图把有关复数的奇怪题意转化成我们比较熟悉的东西。
为了方便把题目中的 (d) 改成 (-d)，考虑 ((a+b{
m i})(c-d{
m i})=(ac+bd)+(ad-bc){
m i})，那么由结果是一个整数可得 (ad-bc=0)，也就是 ((c,d)) 与 ((a,b)) 要成比例。
先不考虑 (ble 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+bdle n).
(ac+bd=kg(a'^2+b'^2)le n)，因此满足上述两个条件的 ((c,d)) 的对数为 (lfloorfrac{n}{g(a'^2+b'^2)}
floor).
枚举 (g,a',b')（为了方便使用 (a,b) 代替），可得答案等于
于是我们把问题化成了我们熟悉的数论函数求和形式。
下一步是套路的莫比乌斯反演：
现在需要在低于线性的时间复杂度内求出它。
尽管有这么多求和符号，但是我们可以从中观察出一些特点：
- 看起来最棘手的是最内层，但是内层的求和结果似乎只和 (n/gd^2) 有关。那么假设我们把 (sum^{sqrt m}_{a=1,b=1}lfloorfrac{m}{a^2+b^2}
floor) 看成关于 (m) 的函数 (f(m))，那么外层枚举 (d,g) 后所有的求和项都形如 (f(lfloorfrac{n}{k}
floor))，其中 (k) 是整数。而这是一个我们相当喜欢的形式——它总共只有 (O(sqrt n)) 种取值，而且假设我们能在 (O(sqrt m)) 的时间内求出 (f(m))，那么对每个 (k) 求出 (f(lfloorfrac{n}{k}
floor)) 的总时间复杂度为 (sum^{sqrt n}_{k=1}sqrtfrac{n}{k}+sum^{sqrt n}_{k=1}sqrt k=n^{frac{1}{2}}int^{n^{frac{1}{2}}}_{0}t^{-frac{1}{2}}{
m d}t+int^{n^{frac{1}{2}}}_{0}t^{frac{1}{2}}{
m d}t=O(n^{frac{3}{4}})).
进一步观察 (lfloorfrac{n}{k} floor) 这些数（称作“关键点”）形成的结构，其实相当于它们把 ([1,n]) 分成了 (O(sqrt n)) 段，每一段内 (lfloorfrac{n}{i} floor) 的值相等，关键点位于每一段的右端点。我们把这种划分方式称为“(n) 的划分”。设 (x,y) 位于同一段，(m) 是一个关键点，则有 (lfloorfrac{m}{x} floor=lfloorfrac{m}{y} floor). 换句话说，对于关键点 (m)，(m) 的划分中的关键点一定是 (n) 的划分中的关键点的子集，因此如果对某种信息我们对 (n) 的划分的每一段维护出了前缀和，那么就可以快速地查询 (m) 的划分中的一段的和。
由此我们得到 (O(sqrt m)) 时间求出 (f(m)) 的方法：设 (c(m)) 表示满足 (a^2+b^2le m) 的 (a) 之和。显然我们可以做到 (O(sqrt m)) 时间求出 (c(m))——直接枚举 (ale sqrt m) 然后用sqrt函数 (O(1)) 累加进答案即可。现在我们对 (m=lfloorfrac{n}{k} floor) 的每种可能取值都求出了 (c) 函数的值，那么对相应的 (f(m))，枚举 (a^2+b^2) 的取值在 (m) 的划分（注意不是 (n)）中所在的段 (i)，对结果的贡献是 (lfloorfrac{m}{i} floorcdot S(i))，(S(i)) 就是 (a^2+b^2) 落在这一段里的 (a) 值之和，可以由 (n) 的划分中对应的两段 (c) 值相减得到。 - 外层的枚举具有极其相似的结构：假设固定了 (d)，那么剩下的就是 (sum^{n/d^2}_{g=1}gcdot f(lfloorfrac{n/d^2}{g}
floor)). 设函数 (h(m)=sum_{g=1}gcdot f(lfloorfrac{m}{g}
floor))，那我们需要的依然只是 (h) 函数在所有 (lfloorfrac{n}{k}
floor) ((k) 为整数)处的点值。在求 (h(m)) 时使用整除分块，和前面不同的是此处求 (h) 的过程是递归的，像杜教筛一样写一个记忆化搜索，递归到所有 (lfloorfrac{m}{k}
floor) 并通过它们的 (f) 和 (h) 函数值求出 (h(m)).
而最后枚举 (d)，将 (dmu(d)h(frac{n}{d^2})) 累加进答案就可以了。
最后的答案乘以 (2)，再加上 (b=0) 的答案（就是 (1) 到 (n) 所有数的约数和）。
总时间复杂度 (O(n^{frac{3}{4}}))，空间复杂度 (O(sqrt n))，或许可以用预处理优化到更快（？）
实际上这道题给我们的启发就是，对于任何奇形怪状的函数，如果我们只需要并保存它 (lfloorfrac{n}{k} floor) 处的取值，就可以在 (O(n^{frac{3}{4}})) 的时间内完成某种类似于卷积的操作。设 (f)、(g) 是两个这样的函数，(h(n)=sum^n_{i=1}f(i)g(lfloorfrac{n}{i} floor))，那么就可以求出 (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
",ans);
return 0;
}