大致题意： 求\(1\sim n\)内所有满足\(a>0\)的约数\(a+bi\)的\(a\)之和。
解题思路
首先，我们设\(x=(a+bi)(c+di)(1\le x\le n)\)，则：
\[x=(ac-bd)+(ad+bc)i\]
由于\(x\)是一个实数，因此：
\[①②\begin{cases}ac-bd=x,①\\ad+bc=0.②\end{cases}\]
而由\(②②\)式可得：
\[ad+bc=0⇒ad=-bc⇒\frac ab=-\frac cd\]
设\(\frac ab=\frac pq(gcd(p,q)=1)\)，则\(\frac cd=-\frac pq\)，而原式就变成了：
\[x=s(p+qi)\cdot t(p-qi)\]
然后可以发现\(p+qi\)与\(p-qi\)似乎可以利用平方差公式化简，即：
\[x=s\cdot t\cdot(p^2+q^2)\]
也就是说，只要\(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\)对答案造成的总贡献为：
\[\sum_{w|\frac x{p^2+q^2}}p\cdot w\]
如果我们将\(p\)提前，则可以发现剩下的部分恰好是一个\(\sigma\)（约数和）函数，即：
\[p\cdot\sigma(\frac x{p^2+q^2})\]
那么枚举所有\(p,q\)可得总答案为：
\[\sum_{gcd(p,q)=1}p\sum_{(p^2+q^2)|x}\sigma(\frac x{p^2+q^2})\]
考虑枚举\(p^2+q^2=y\)，则原式就相当于：
\[\sum_y\sum_{gcd(p,q)=1\&\&p^2+q^2=y}p\sum_{y|x}\sigma(\frac xy)\]
由于\(\sum_{y|x}\sigma(\frac xy)\)其实就是\(\sum_{i=1}^{\lfloor\frac ny\rfloor}\sigma(i)\)，因此若我们设\(D(i)=\sum_{k=1}^i\sigma(k)\)，并设\(f(i)=\sum_{gcd(p,q)=1\&\&p^2+q^2=i}p\)，则原式就可以转化为：
\[\sum_yf(y)\cdot D(\lfloor\frac ny\rfloor)\]
看到\(D(\lfloor\frac ny\rfloor)\)，便可以想到用去搞。
设\(F(i)=\sum_{k=1}^if(i)\)，那么，我们现在的问题就是，如何快速求\(F\)和\(D\)的值。
如何求\(F\)
对于\(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)\)，可得：
\[G(n)=\sum_{d\ge1}d\cdot F(\lfloor\frac n {d^2}\rfloor)\]
然后就用杜教筛的经典转化套路，单独提出\(d=1\)，得到：
\[G(n)=F(n)+\sum_{d\ge2}d\cdot F(\lfloor\frac n {d^2}\rfloor)\]
移项，得到式子为：
\[F(n)=G(n)-\sum_{d\ge2}d\cdot F(\lfloor\frac n {d^2}\rfloor)\]
由于\(G(n)\)可以通过\(O(\sqrt n)\)的时间计算出，然后我们除法分块递归求解\(F\)即可。
如何求\(D\)
\(hl666\)神仙说，他以前切过一道题就是求这道题里的\(D\)。
根据他的结论，我们得知：
\[D(n)=\sum_{i=1}^n\sigma(i)=\sum_{d=1}^nd\cdot\lfloor\frac nd\rfloor\]
这应该还是比较好懂的，就是枚举约数暴力算个数。
依然除法分块，\(O(\sqrt n)\)的时间复杂度内可以计算出。
代码
#include#define Tp template #define Ts template #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;//输出答案}