【LOJ6052】「雅礼集训 2017 Day11」DIV（杜教筛）

点此看题面

大致题意： 求(1sim n)内所有满足(a>0)的(n)的复数约数(a+bi)的(a)之和。

解题思路

首先，我们设(x=(a+bi)(c+di)(1le xle n))，则：

[x=(ac-bd)+(ad+bc)i ]

由于(x)是一个实数，因此：

[egin{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=scdot tcdot(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}}pcdot w ]

如果我们将(p)提前，则可以发现剩下的部分恰好是一个(sigma)（约数和）函数，即：

[pcdotsigma(frac x{p^2+q^2}) ]

那么枚举所有(p,q)可得总答案为：

[sum_{gcd(p,q)=1}psum_{(p^2+q^2)|x}sigma(frac x{p^2+q^2}) ]

考虑枚举(p^2+q^2=y)，则原式就相当于：

[sum_ysum_{gcd(p,q)=1&&p^2+q^2=y}psum_{y|x}sigma(frac xy) ]

由于(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)，则原式就可以转化为：

[sum_yf(y)cdot D(lfloorfrac ny floor) ]

看到(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))，可得：

[G(n)=sum_{dge1}dcdot F(lfloorfrac n {d^2} floor) ]

然后就用杜教筛的经典转化套路，单独提出(d=1)，得到：

[G(n)=F(n)+sum_{dge2}dcdot F(lfloorfrac n {d^2} floor) ]

移项，得到式子为：

[F(n)=G(n)-sum_{dge2}dcdot F(lfloorfrac n {d^2} floor) ]

由于(G(n))可以通过(O(sqrt n))的时间计算出，然后我们除法分块递归求解(F)即可。

如何求(D)

(hl666)神仙说，他以前切过一道题就是求这道题里的(D)。

根据他的结论，我们得知：

[D(n)=sum_{i=1}^nsigma(i)=sum_{d=1}^ndcdotlfloorfrac nd floor ]

这应该还是比较好懂的，就是枚举约数暴力算个数。

依然除法分块，(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;//输出答案
}