51Nod1220约数之和-数论

推了一下午，写又花了一晚上，XD，自己还是太菜了

【题目地址】

题意简述

单组询问，给定一个$n$，求下面的式子的值在$\mathrm{m}\mathrm{o}\mathrm{d}10^9+7$意义下的值。

$$\sum _{i=1}^n\sum _{j=1}^n{\sigma }_1(ij)$$

其中${\sigma }_1(i)$表示$i$的约数和，例如${\sigma }_1(12)=1+2+3+4+6+12=28$

我们按照[SDOI2015]约数个数和的套路先考虑如何转换${\sigma }_1(ij)$，其实也可以考虑映射，所以可以写成如下形式：

$${\sigma }_1(ij)=\sum _{x|i}\sum _{y|j}[gcd(x,y)=1]\frac{xj}{y}$$

简略证明：

同样，每个质因数直间互不影响，所以我们考虑，对于一个质因数$p$，其中$p^a|i,p^b|j$（这里和上面链接的那道题一样，$\frac{i}{p^a}\perp p^a,\frac{j}{p^b}\perp p^b$），所以对于当$x=p^0$时，$\frac{xj}{y}$的$y$中肯定含有$p^0\sim p^b$，而$xj$中也只含有$p^0\sim p^b$（此时仅对于这一个质因数），所以此时可以将因子为$p^0\sim p^b$的取完，而另一部分则是在$y=p^0$时，$\frac{xi}{y}$的$y=1$了，而$xi$肯定含有$p^b\sim p^{a+b}$，而算$p^b$时，是在$x,y$同时为$1$的情况下，只会算一次，所以不用减去它。那么此时对于因子$p$的所有情况就算完了，而每次计算的时候都是一个没有因子$p$，一个有，所以要满足$gcd(x,y)=1$，将其应用到所有的质因子上，就可以得到上面的式子。

那么原式就可以转化成：

$$\sum _{i=1}^n\sum _{j=1}^n\sum _{x|i}\sum _{y|j}[gcd(i,j)=1]\frac{xj}{y}$$

我们用莫比乌斯反演的套路将其进行变形：

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^n\sum _{d=1}^n\mu (d)\sum _{x|i}\sum _{y|j}\frac{xj}{y}[d|gcd(x,y)] \\ \sum _{d=1}^n\mu (d)\sum _{d|x}\sum _{d|y}\frac{x}{y}\sum _{x|i}^n\sum _{y|j}^{n}j \\ \sum _{d=1}^n\mu (d)\sum _{d|x}\sum _{d|y}\frac{x}{y}\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{y}\rfloor }jy \\ \sum _{d=1}^n\mu (d)\sum _{d|x}\sum _{d|y}\frac{x}{y}\lfloor \frac{n}{x}\rfloor y\sum _{j=1}^{\lfloor \frac{n}{y}\rfloor }j \\ \sum _{d=1}^n\mu (d)\sum _{d|x}\sum _{d|y}\frac{x}{y}\lfloor \frac{n}{x}\rfloor y\frac{\lfloor \frac{n}{y}\rfloor \times (\lfloor \frac{n}{y}\rfloor +1)}{2} \\ \sum _{d=1}^n\mu (d)\sum _{d|x}\sum _{d|y}x\lfloor \frac{n}{x}\rfloor \frac{\lfloor \frac{n}{y}\rfloor \times (\lfloor \frac{n}{y}\rfloor +1)}{2} \\ \sum _{d=1}^n\mu (d)\sum _{d|x}x\lfloor \frac{n}{x}\rfloor \sum _{d|y}\frac{\lfloor \frac{n}{y}\rfloor \times (\lfloor \frac{n}{y}\rfloor +1)}{2} \end{gathered}$$

而此时我们看后面那一大坨，分开考虑：

$$\sum _{d|x}x\lfloor \frac{n}{x}\rfloor =\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }xd\lfloor \frac{n}{xd}\rfloor$$

我们把那一个$d$拿出去，然后令$f(n)=\sum _{i=1}^{n}i\lfloor \frac{n}{i}\rfloor$，其实根据前面链接的经验，我们很容易知道：

$$f(n)=\sum _{i=1}^{n}i\lfloor \frac{n}{i}\rfloor =\sum _{i=1}^n{\sigma }_1(i)$$

证明：

首先可以看出每一个因子乘以出现次数，然后就是前面链接里的证明了。

对于另一个，我们同样转换：

$$\begin{gathered} \sum _{d|y}\frac{\lfloor \frac{n}{y}\rfloor \times (\lfloor \frac{n}{y}\rfloor +1)}{2} \\ =\sum _{y=1}^{\lfloor \frac{n}{d}\rfloor }\frac{\lfloor \frac{n}{yd}\rfloor \times (\lfloor \frac{n}{yd}\rfloor +1)}{2} \end{gathered}$$

其实我们令$g(n)={\sum }_{i=1}^n\frac{i\times (i+1)}{2}$，那么可以得到：

$$g(n)=\sum _{i=1}^n\frac{\lfloor \frac{n}{i}\rfloor \times (\lfloor \frac{n}{i}\rfloor +1)}{2}=\sum _{i=1}^n\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }j$$

此时$j$要有贡献，条件为$j\le \lfloor \frac{n}{i}\rfloor$，我们将其变形可以得到：

$$\begin{gathered} j\le \lfloor \frac{n}{i}\rfloor \\ j\le \frac{n}{i} \\ i\le \frac{n}{j} \\ i\le \lfloor \frac{n}{j}\rfloor \end{gathered}$$

所以我们改变枚举顺序，得到：
$$\begin{gathered} g(n)=\sum _{j=1}^{n}j\sum _{i=1}^{\lfloor \frac{n}{i}\rfloor }1 \\ g(n)=\sum _{i=1}^{n}i\lfloor \frac{n}{i}\rfloor =f(n) \end{gathered}$$

所以原式变成了：

$$\sum _{d=1}^n\mu (d)d{\left(f(\frac{n}{d})\right)}^2$$

此时我们需要用杜教筛筛$F(n)=\mu (n)n$，我们令$G(n)=n$，那么$F⨂G=1$，所以直接套入杜教筛公式，得到：

$$S(n)=1-\sum _{i=2}^{n}iS(\lfloor \frac{n}{i}\rfloor )$$

前面用杜教筛，后面预处理加$O(\sqrt{n})$的算即可，复杂度$O(n^{\frac{2}{3}})$

#include<map>
#include<cstdio>
#include<cstring>
#include<algorithm>
#pragma GCC optimize(2)
#define ll long long
using namespace std;
const ll Mod=1e9+7;
const int M=2e6+10,MAX=2e6;
ll n,m,inv_2,cnt;
ll prime[M],mu[M],F[M],B[M];
ll &Rec(ll x){
	if(x<=MAX) return F[x];
	else return B[n/x];
}
bool vis[M];
ll s[M],p[M],sum[M],bum[M];
ll &Rec2(ll x){
	if(x<=MAX) return sum[x];
	else return bum[n/x];
}//空间储存技巧
ll fpow(ll a,ll b){
	ll res=1;
	for(;b;b>>=1,a=(a*a)%Mod){
		if(b&1)res=(res*a)%Mod;
	}
	return res;
}
void init(){
	memset(bum,-1,sizeof(bum));
	memset(B,-1,sizeof(B));
	s[1]=mu[1]=1;inv_2=fpow(2,Mod-2);
	for(int i=2;i<=MAX;i++){
		if(!vis[i]){
			prime[++cnt]=i;
			s[i]=p[i]=1+i;
			mu[i]=-1;
		}
		for(int j=1,v;j<=cnt&&i*prime[j]<=MAX;j++){
			v=i*prime[j];
			vis[v]=1;
			if(!(i%prime[j])){
				s[v]=s[i]/p[i]*(p[i]*prime[j]+1);
				p[v]=p[i]*prime[j]+1;
				break;
			}
			mu[v]=-mu[i];
			s[v]=s[i]*s[prime[j]];
			p[v]=1+prime[j];
		}
	}
	for(int i=1;i<=MAX;i++){
		sum[i]=(sum[i-1]+s[i]%Mod)%Mod;
		F[i]=(F[i-1]+(i*mu[i])%Mod)%Mod;
		if(F[i]<0)F[i]+=Mod;
	}
}
ll S(ll x){if(x>=Mod)x%=Mod;return (x*(x+1)%Mod)*inv_2%Mod;}
ll Sarea(ll L,ll R){return ((S(R)-S(L-1))%Mod+Mod)%Mod;}
ll GetS(ll x){
	ll &ans=Rec2(x);
	if(ans!=-1) return ans;
	ans=0;
	for(ll i=1,j;i<=x;i=j+1){
		j=(x/(x/i));
		ans=(ans+Sarea(i,j)*(x/i)%Mod)%Mod;
	}
	return ans;
}
ll calc(ll x){
	ll &ans=Rec(x);
	if(ans!=-1) return ans;
	ans=1;
	for(ll i=2,j;i<=x;i=j+1){
		j=(x/(x/i));
		ans=(ans-Sarea(i,j)*calc(x/i)%Mod)%Mod;
	}
	if(ans<0)ans+=Mod;
	return ans;
}
ll Carea(ll L,ll R){return ((calc(R)-calc(L-1))%Mod+Mod)%Mod;}
ll Sqr(ll a){return a*a%Mod;}
ll solve(ll x){
	ll ans=0;
	for(ll i=1,j;i<=x;i=j+1){
		j=(x/(x/i));
		ans=(ans+Carea(i,j)*Sqr(GetS(x/i))%Mod)%Mod;
	}
	return ans;
}
int main(){
	init();
	scanf("%lld",&n);
	printf("%lld\n",solve(n));
	return 0;
}

另外一种类似的推导方式【IN-orz%%%】