完了我是数学姿势越来越弱了，感觉这种CXRdalao秒掉的题我都要做好久
一些前置推导
首先我们很容易得出\((a+bi)(c+di)=k \Leftrightarrow ac-bd=k,ad+bc=0\)
我们可以直接\(ad+bc=0\Rightarrow ad=-bc\Rightarrow \frac{a}b=-\frac{c}{d}\)
考虑把这个分数化为最简的形式，那么就意味着我们要把\(\gcd\)拿出来
我们令\(\frac{a}b=\frac{p}{q}(\gcd(p,q)=1)\)，那么\(\frac{c}d=-\frac{p}q\)
把这个代回去就有\(x(p+qi)\cdot q(p-qi)=k\)
然后直接把式子乘起来就有\(xy(p^2+q^2)=k\)
那么我们可以发现，如果\(p^2+q^2|k\)，那么它就可以对答案产生贡献
然后考虑求贡献和，即所有是\(p^2+q^2\)倍数以及\(k\)的倍数的数的个数
假设这个数\(M=w(p^2+q^2)\)，那么我们知道\(w|\frac{k}{p^2+q^2}\)，可以列出贡献的式子：
\]
然后我们枚举\(p,q\)后统计答案，发现不好维护，因此可以直接枚举\(t=p^2+q^2\)，则原式等于：
\]
然后我们考虑简化这个式子，首先\(\sum_{y|k}\sigma(\frac{k}t)\)其实就是\(\sum_{i=1}^{\frac{n}t}\sigma i\)
所以我们记\(\sigma\)的前缀和为\(D\)，然后为了方便把\(\sum_{\gcd(p,q)=1\&\&p^2+q^2=t} p\)设为\(F\)，这样原式即为：
\]
是我们熟悉的除法分块形式，所以考虑分别求出\(D,F\)的值，由于这里的数据范围比较大所以我们考虑用杜教筛
求解\(D\)
先讲比较简单的\(D\)的求解，首先如果是小范围我们可以直接用线性筛筛出单个的\(\sigma\)然后做前缀和
然后有一个很简单的结论，我们可以直接暴力枚举约数算个数，即：
\]
这个直接除法分块一下，然后总体就是\(O(n^{\frac{2}3})\)的
求解\(F\)
首先还是小范围答案，我们可以线性筛出素数的时候直接枚举\(p,q\)然后算贡献即可
我们考虑对\(F\)做前缀和，即令\(G(n)=\sum_{p^2+q^2\le n} p=\sum_{p=1}^{\lfloor\sqrt n\rfloor} p\cdot \lfloor \sqrt{n-p^2}\rfloor\)
考虑直接枚举\(\gcd(p,q)\)，那么即有：
\]
根据杜教筛的套路，我们直接把\(d=1\)的情况提出来，那么就有：
\]
即得到\(F(n)=G(n)-\sum_{d\ge 2} d\cdot F(\lfloor\frac{n}{d^2}\rfloor)\)
这里由于求解单个\(G\)是\(\sqrt n\)的，因此总体复杂度还是\(n^{\frac{2}3}\)的
综上，我们总算是把这道杜教筛的练手题做完了，然后我把一个\(x\)打成\(n\)调了一晚上233
CODE
#include<cstdio>
#include<map>
#include<cmath>
#define RI register int
#define CI const int&
using namespace std;
typedef long long LL;
const int N=5000000,mod=1004535809,inv2=502267905;
int prime[N+5],cnt,ans; bool vis[N+5];
LL n,ds[N+5],fs[N+5]; map <LL,int> _ds,_fs;
inline void inc(LL& x,const LL y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void inc(int& x,CI y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void dec(int& x,CI y)
{
if ((x-=y)<0) x+=mod;
}
inline int sum(CI x,CI y)
{
int t=x+y; return t>=mod?t-mod:t;
}
inline int sub(CI x,CI y)
{
int t=x-y; return t<0?t+mod:t;
}
inline int gcd(CI x,CI y)
{
return y?gcd(y,x%y):x;
}
inline int Sum(const LL& l,const LL& r)
{
return ((l+r)%mod)*((r-l+1)%mod)%mod*inv2%mod;
}
#define Pi prime[j]
inline void init(CI n)
{
RI i,j; ds[1]=vis[1]=1; for (i=2;i<=n;++i)
{
if (!vis[i]) prime[++cnt]=i,ds[i]=i+1;
for (j=1;j<=cnt&&i*Pi<=n;++j)
{
vis[i*Pi]=1; if (i%Pi) ds[i*Pi]=ds[i]*(Pi+1);
else { ds[i*Pi]=ds[i]*(Pi+1)-Pi*ds[i/Pi]; break; }
}
}
for (i=1;i*i<=n;++i)
{
int t=i*i; for (j=1;j*j+t<=n;++j) if (gcd(i,j)==1) fs[j*j+t]+=i;
}
for (i=1;i<=n;++i) ds[i]%=mod,fs[i]%=mod,inc(ds[i],ds[i-1]),inc(fs[i],fs[i-1]);
}
#undef Pi
inline int Ds(const LL& x)
{
if (x<=N) return ds[x]; if (_ds.count(x)) return _ds[x]; int ret=0;
for (LL l=1,r;l<=x;l=r+1) r=x/(x/l),inc(ret,1LL*Sum(l,r)*((x/r)%mod)%mod); return _ds[x]=ret;
}
inline int Fs(const LL& x)
{
if (x<=N) return fs[x]; if (_fs.count(x)) return _fs[x]; int ret=0; register LL i;
for (i=1;i*i<=x;++i) inc(ret,i*((LL)floor(sqrt(x-i*i))%mod)%mod);
for (i=2;i*i<=x;++i) dec(ret,i*Fs(x/(i*i))%mod); return _fs[x]=ret;
}
int main()
{
scanf("%lld",&n); init(N); for (LL l=1,r;l<=n;l=r+1)
r=n/(n/l),inc(ans,1LL*sub(Fs(r),Fs(l-1))*Ds(n/l)%mod);
return printf("%d",sum(sum(ans,ans),Ds(n))),0;
}