51Nod1220約數之和-數論

VictoryCzt發表於2018-12-19

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

【題目地址】

題意簡述

單組詢問，給定一個nn，求下面的式子的值在mod 109+7{\rm mod}\ 10^9+7意義下的值。

∑i=1n∑j=1nσ1(ij) \sum_{i=1}^n\sum_{j=1}^n\sigma_1(ij)

其中σ1(i)\sigma_1(i)表示ii的約數和，例如σ1(12)=1+2+3+4+6+12=28\sigma_1(12)=1+2+3+4+6+12=28

我們按照[SDOI2015]約數個數和的套路先考慮如何轉換σ1(ij)\sigma_1(ij)，其實也可以考慮對映，所以可以寫成如下形式：

σ1(ij)=∑x∣i∑y∣j[gcd(x,y)=1]xjy \sigma_1(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\frac{xj}{y}

簡略證明：

同樣，每個質因數直間互不影響，所以我們考慮，對於一個質因數pp，其中pa∣i,pb∣jp^a|i,p^b|j（這裡和上面連結的那道題一樣，ipa⊥pa,jpb⊥pb\frac{i}{p^a}\perp p^a,\frac{j}{p^b}\perp p^b），所以對於當x=p0x=p^0時，xjy\frac{xj}{y}的yy中肯定含有p0∼pbp^0\sim p^b，而xjxj中也只含有p0∼pbp^0\sim p^b（此時僅對於這一個質因數），所以此時可以將因子為p0∼pbp^0\sim p^b的取完，而另一部分則是在y=p0y=p^0時，xiy\frac{xi}{y}的y=1y=1了，而xixi肯定含有pb∼pa+bp^{b}\sim p^{a+b}，而算pbp^b時，是在x,yx,y同時為11的情況下，只會算一次，所以不用減去它。那麼此時對於因子pp的所有情況就算完了，而每次計算的時候都是一個沒有因子pp，一個有，所以要滿足gcd(x,y)=1gcd(x,y)=1，將其應用到所有的質因子上，就可以得到上面的式子。

那麼原式就可以轉化成：

∑i=1n∑j=1n∑x∣i∑y∣j[gcd(i,j)=1]xjy \sum_{i=1}^n\sum_{j=1}^n\sum_{x|i}\sum_{y|j}[gcd(i,j)=1]\frac{xj}{y}

我們用莫比烏斯反演的套路將其進行變形：

∑i=1n∑j=1n∑d=1nμ(d)∑x∣i∑y∣jxjy[d∣gcd(x,y)]∑d=1nμ(d)∑d∣x∑d∣yxy∑x∣in∑y∣jnj∑d=1nμ(d)∑d∣x∑d∣yxy∑i=1⌊nx⌋∑j=1⌊ny⌋jy∑d=1nμ(d)∑d∣x∑d∣yxy⌊nx⌋y∑j=1⌊ny⌋j∑d=1nμ(d)∑d∣x∑d∣yxy⌊nx⌋y⌊ny⌋×(⌊ny⌋+1)2∑d=1nμ(d)∑d∣x∑d∣yx⌊nx⌋⌊ny⌋×(⌊ny⌋+1)2∑d=1nμ(d)∑d∣xx⌊nx⌋∑d∣y⌊ny⌋×(⌊ny⌋+1)2 \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}^nj \\\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}\left\lfloor\frac{n}{x}\right\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}\left\lfloor\frac{n}{x}\right\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\left\lfloor\frac{n}{x}\right\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\left\lfloor\frac{n}{x}\right\rfloor\sum_{d|y} \frac{\lfloor\frac{n}{y}\rfloor\times(\lfloor\frac{n}{y}\rfloor+1)}{2}

而此時我們看後面那一大坨，分開考慮：

∑d∣xx⌊nx⌋=∑x=1⌊nd⌋xd⌊nxd⌋ \sum_{d|x}x\left\lfloor\frac{n}{x}\right\rfloor=\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}xd\left\lfloor\frac{n}{xd}\right\rfloor

我們把那一個dd拿出去，然後令f(n)=∑i=1ni⌊ni⌋f(n)=\sum\limits_{i=1}^ni\left\lfloor\frac{n}{i}\right\rfloor，其實根據前面連結的經驗，我們很容易知道：

f(n)=∑i=1ni⌊ni⌋=∑i=1nσ1(i) f(n)=\sum\limits_{i=1}^ni\left\lfloor\frac{n}{i}\right\rfloor=\sum_{i=1}^n\sigma_1(i)

證明：

首先可以看出每一個因子乘以出現次數，然後就是前面連結裡的證明了。

對於另一個，我們同樣轉換：

∑d∣y⌊ny⌋×(⌊ny⌋+1)2=∑y=1⌊nd⌋⌊nyd⌋×(⌊nyd⌋+1)2 \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}

其實我們令g(n)=∑i=1ni×(i+1)2g(n)=\sum_{i=1}^n\frac{i\times(i+1)}{2}，那麼可以得到：

g(n)=∑i=1n⌊ni⌋×(⌊ni⌋+1)2=∑i=1n∑j=1⌊ni⌋j 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

此時jj要有貢獻，條件為j≤⌊ni⌋j\leq\lfloor\frac{n}{i}\rfloor，我們將其變形可以得到：

j≤⌊ni⌋ j≤ni i≤nj i≤⌊nj⌋ j\leq\lfloor\frac{n}{i}\rfloor \\\ \\j\leq\frac{n}{i} \\\ \\i\leq\frac{n}{j} \\\ \\i\leq\lfloor\frac{n}{j}\rfloor

所以我們改變列舉順序，得到：
g(n)=∑j=1nj∑i=1⌊ni⌋1g(n)=∑i=1ni⌊ni⌋=f(n) g(n)=\sum_{j=1}^nj\sum_{i=1}^{\lfloor\frac{n}{i}\rfloor}1 \\ g(n)=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor=f(n)

所以原式變成了：

∑d=1nμ(d)d(f(nd))2 \sum_{d=1}^n\mu(d)d\left(f(\frac{n}{d})\right)^2

此時我們需要用杜教篩篩F(n)=μ(n)nF(n)=\mu(n)n，我們令G(n)=nG(n)=n，那麼F⨂G=1F\bigotimes G=1，所以直接套入杜教篩公式，得到：

S(n)=1−∑i=2niS(⌊ni⌋) S(n)=1-\sum_{i=2}^niS(\lfloor\frac{n}{i}\rfloor)

前面用杜教篩，後面預處理加O(n)O(\sqrt n)的算即可，複雜度O(n23)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%%%】

相關文章