大致题意： 求(1sim n)内所有满足(a>0)的(n)的复数约数(a+bi)的(a)之和。
解题思路
首先，我们设(x=(a+bi)(c+di)(1le xle n))，则：
由于(x)是一个实数，因此：
而由(②)式可得：
设(frac ab=frac pq(gcd(p,q)=1))，则(frac cd=-frac pq)，而原式就变成了：
然后可以发现(p+qi)与(p-qi)似乎可以利用平方差公式化简，即：
也就是说，只要(p^2+q^2)是(x)的约数，它就可以对答案造成贡献。
而造成的的贡献值是多少呢？
考虑每一个既为(p^2+q^2)的倍数，又为(x)的约数的数，都可以对答案造成贡献。
则设一个造成贡献的数为(w(p^2+q^2))。
由于(w(p^2+q^2))是(x)的约数，所以(w)是(frac x{p^2+q^2})的约数。
因此，(p^2+q^2)对答案造成的总贡献为：
如果我们将(p)提前，则可以发现剩下的部分恰好是一个(sigma)（约数和）函数，即：
那么枚举所有(p,q)可得总答案为：
考虑枚举(p^2+q^2=y)，则原式就相当于：
由于(sum_{y|x}sigma(frac xy))其实就是(sum_{i=1}^{lfloorfrac ny floor}sigma(i))，因此若我们设(D(i)=sum_{k=1}^isigma(k))，并设(f(i)=sum_{gcd(p,q)=1&&p^2+q^2=i}p)，则原式就可以转化为：
看到(D(lfloorfrac ny floor))，便可以想到用除法分块去搞。
设(F(i)=sum_{k=1}^if(i))，那么，我们现在的问题就是，如何快速求(F)和(D)的值。
如何求(F)
对于(F)，考虑杜教筛。
设(G(n)=sum_{p^2+q^2le n}p=sum_{p=1}^{lfloorsqrt n floor}pcdotlfloorsqrt{n-p^2} floor)，枚举(gcd(p,q))，可得：
然后就用杜教筛的经典转化套路，单独提出(d=1)，得到：
移项，得到式子为：
由于(G(n))可以通过(O(sqrt n))的时间计算出，然后我们除法分块递归求解(F)即可。
如何求(D)
(hl666)神仙说，他以前切过一道题就是求这道题里的(D)。
根据他的结论，我们得知：
这应该还是比较好懂的，就是枚举约数暴力算个数。
依然除法分块，(O(sqrt n))的时间复杂度内可以计算出。
代码
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define RL Reg LL
#define Con const
#define CI Con int&
#define CL Con LL&
#define I inline
#define W while
#define LL long long
#define X 1004535809
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
#define Dec(x,y) ((x-=(y))<0&&(x+=X))
using namespace std;
LL n;
I int gcd(CI x,CI y) {return y?gcd(y,x%y):x;}
I int XSum(CI x,CI y) {return x+y>=X?x+y-X:x+y;}
I int XSub(CI x,CI y) {return x-y<0?x-y+X:x-y;}
class DuSieve//杜教筛，虽然不该把求D也放在这里面。。。
{
private:
#define S 4641588
static Con int I2=X+1>>1;int Pc,P[S+5],d[S+5],f[S+5],D[S+5],F[S+5];
public:
I DuSieve()//初始化，线性筛出一部分D值，并求出一部分F值
{
RI i,j,k;for(d[1]=1,i=2;i<=S;++i)
for(!P[i]&&(d[P[++Pc]=i]=i+1),j=1;j<=Pc&&i*P[j]<=S;++j)
if(P[i*P[j]]=1,i%P[j]) d[i*P[j]]=1LL*d[i]*(P[j]+1)%X;
else {d[i*P[j]]=XSub(1LL*d[i]*(P[j]+1)%X,1LL*d[i/P[j]]*P[j]%X);break;}
for(i=1;i*i<=S;++i) for(k=i*i,j=1;k+j*j<=S;++j) gcd(i,j)==1&&Inc(f[k+j*j],i);
for(i=1;i<=S;++i) Inc(d[i],d[i-1]),Inc(f[i],f[i-1]);//统计前缀和
}
I int GetF(CL x)//求F
{
if(x<=S) return f[x];if(F[n/x]) return F[n/x];RI i,res=0;//对于小数据或已经求结果的部分直接返回答案
for(i=1;1LL*i*i<=x;++i) Inc(res,(LL)floor(sqrt(x-1LL*i*i))*i%X);//求G值和
for(i=2;1LL*i*i<=x;++i) Dec(res,1LL*GetF(x/(1LL*i*i))*i%X);return F[n/x]=res;//递归计算F
}
#define sum(x,y) (1LL*((x+y)%X)*((y-x+1)%X)%X*I2%X)
I int GetD(CL x)//求D
{
if(x<=S) return d[x];if(D[n/x]) return D[n/x];RL l,r;RI res=0;//对于小数据或已经求结果的部分直接返回答案
for(l=1;l<=x;l=r+1) r=x/(x/l),Inc(res,1LL*(x/r)%X*sum(l,r)%X);return D[n/x]=res;//除法分块
}
}D;
int main()
{
RL l,r,ans=0;for(scanf("%lld",&n),l=1;l<=n;l=r+1)//除法分块
r=n/(n/l),Inc(ans,1LL*XSub(D.GetF(r),D.GetF(l-1))*D.GetD(n/l)%X);//统计答案
return printf("%lld",XSum(XSum(ans,ans),D.GetD(n))),0;//输出答案
}