【LOJ6052】「雅礼集训 2017 Day11」DIV（杜教筛）
大致题意： 求\(1\sim n\)内所有满足\(a>0\)的\(n\)的复数约数\(a+bi\)的\(a\)之和。
解题思路
首先，我们设\(x=(a+bi)(c+di)(1\le x\le 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}^{\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\)，则原式就可以转化为：
看到\(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)\)，可得：
然后就用杜教筛的经典转化套路，单独提出\(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;//输出答案
}
