完了我是数学姿势越来越弱了，感觉这种CXRdalao秒掉的题我都要做好久
一些前置推导
首先我们很容易得出((a+bi)(c+di)=k Leftrightarrow ac-bd=k,ad+bc=0)
我们可以直接(ad+bc=0Rightarrow ad=-bcRightarrow 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^2le n} p=sum_{p=1}^{lfloorsqrt n floor} pcdot lfloor sqrt{n-p^2} floor)
考虑直接枚举(gcd(p,q))，那么即有：
根据杜教筛的套路，我们直接把(d=1)的情况提出来，那么就有：
即得到(F(n)=G(n)-sum_{dge 2} dcdot F(lfloorfrac{n}{d^2} floor))
这里由于求解单个(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;
}