[51Nod 1227] 平均最小公倍数 (杜教筛)

题目描述

i=abj=1ilcm(i,j)ilargesum_{i=a}^bsum_{j=1}^ifrac{lcm(i,j)}i
1<=a<=b<=1091 <= a <= b <= 10^9

题目分析

这道题其实是[51Nod 1238] 最小公倍数之和题解的简化版,或者说是本质…
就直接上公式了

f(n)=i=1nlcm(n,i)n=i=1ni(n,i)large f(n)=sum_{i=1}^nfrac{lcm(n,i)}n=sum_{i=1}^nfrac i{(n,i)},则Ans=i=abf(i)Ans=sum_{i=a}^bf(i)
f(n)=dndiid[(i,n)==d]=dndiid[(id,nd)==1]=dni=1ndi[(i,nd)==1]=dni=1di[(i,d)==1]large f(n)=sum_{d|n}sum_{d|i}frac id[(i,n)==d]\=sum_{d|n}sum_{d|i}frac id[(frac id,frac nd)==1]\=sum_{d|n}sum_{i=1}^{lfloorfrac nd floor}i[(i,frac nd)==1]\=sum_{d|n}sum_{i=1}^di[(i,d)==1]
此处有一个常识
i=1ni[(i,n)==1]=φ(n)n+[n==1]2sum_{i=1}^ni[(i,n)==1]=frac {varphi(n)n+[n==1]}2

  • 证明如下
    • n>1n>1时,若(i,n)=1    (ni,n)=1(i,n)=1iff(n-i,n)=1,所以与nn互质的数是成对出现,且他们的和为nn
    • 再加之n=1n=1的特殊情况,可得
      i=1ni[(i,n)==1]=φ(n)n+[n==1]2sum_{i=1}^ni[(i,n)==1]=frac {varphi(n)n+[n==1]}2

继续
f(n)=dnφ(d)d+[d==1]2=1+dnφ(d)d2large herefore f(n)=sum_{d|n}frac {varphi(d)d+[d==1]}2\=frac{1+sum_{d|n}varphi(d)d}2
i=1nf(i)=n+i=1ndnφ(d)d2=n+d=1nφ(d)ddi12=n+d=1nφ(d)dnd2large herefore sum_{i=1}^nf(i)=frac{n+sum_{i=1}^nsum_{d|n}varphi(d)d}2\=\frac{n+sum_{d=1}^nvarphi(d)dsum_{d|i}1}2\=frac{n+sum_{d=1}^nvarphi(d)dlfloorfrac nd floor}2
此时就可以用整除分块优化+杜教筛计算d=1nφ(d)dlarge sum_{d=1}^nvarphi(d)d
h(n)=φ(d)d,g(n)=i=1nh(d)large h(n)=varphi(d)d,g(n)=sum_{i=1}^nh(d)
n=dnφ(d)n2=dnφ(d)n=dnφ(d)dnd=dnh(d)ndi=1ni2=i=1ndih(d)id=d=1nh(d)diid=d=1nh(d)i=1ndilargeecause n=sum_{d|n}varphi(d)\ herefore n^2=sum_{d|n}varphi(d)n=sum_{d|n}varphi(d)dcdotfrac nd=sum_{d|n}h(d)frac nd\ herefore sum_{i=1}^ni^2=sum_{i=1}^nsum_{d|i}h(d)frac id\=sum_{d=1}^nh(d)sum_{d|i}frac id\=sum_{d=1}^nh(d)sum_{i=1}^{lfloorfrac nd floor}i
注意我们是在杜教筛,不能到这里就把i=1ndisum_{i=1}^{lfloorfrac nd floor}i看做Θ(1)Theta(1)可求的式子而之后再也不做变换,那样往往会陷入更麻烦的方法或者死胡同里去,接着往下
i=1ni2=i=1nid=1nih(d)=i=1nig(ni)g(n)=i=1ni2i=2nig(ni)large herefore sum_{i=1}^ni^2=sum_{i=1}^nisum_{d=1}^{lfloorfrac ni floor}h(d)=sum_{i=1}^nicdot g({lfloorfrac ni floor})\ herefore g(n)= sum_{i=1}^ni^2-sum_{i=2}^nicdot g({lfloorfrac ni floor})
然后套杜教筛即可

虽然在外面套了一层整除分块优化,但由于记忆化的原因,不影响时间复杂度,预处理出一部分gg后总复杂度为Θ(n23)largeTheta(n^{frac 23})

…有兴趣的可以去了解一下[51Nod 1238] 最小公倍数之和题解,比这道题恶心点

AC code
#include <cstdio>
#include <cstring>
#include <cmath>
#include <tr1/unordered_map>
#include <algorithm>
using namespace std;
using namespace tr1;
typedef long long LL;
const int inv6 = 166666668;
const int inv2 = 500000004;
const int MAXN = 5e6 + 1;
const int mod = 1e9 + 7;
int Prime[MAXN], phi[MAXN], Cnt;
bool IsnotPrime[MAXN];
LL g[MAXN];
void init()
{
	phi[1] = 1;
	for(int i = 2; i < MAXN; ++i)
	{
		if(!IsnotPrime[i])
			Prime[++Cnt] = i, phi[i] = i-1;
		for(int j = 1; j <= Cnt && i * Prime[j] < MAXN; ++j)
		{
			IsnotPrime[i * Prime[j]] = 1;
			if(i % Prime[j] == 0) { phi[i * Prime[j]] = phi[i] * Prime[j]; break; }
			phi[i * Prime[j]] = phi[i] * phi[Prime[j]];
		}
	}
	for(int i = 1; i < MAXN; ++i)
		g[i] = (g[i-1] + 1ll * phi[i] * i % mod) % mod;
}
unordered_map<LL, LL> G;
inline LL sum1(LL i, LL j) { return ((i+j)%mod) * ((j-i+1)%mod) % mod * inv2 % mod; }
inline LL sum(LL n)
{
	if(n < MAXN) return g[n];
	if(G.count(n)) return G[n];
	LL ret = (n%mod) * ((n+1)%mod) % mod * ((2*n+1)%mod) % mod * inv6 % mod;
	for(LL i = 2, j; i <= n; i=j+1)
	{
		j = n/(n/i);
		ret = (ret - sum(n/i) * sum1(i, j) % mod) % mod;
	}
	return G[n]=ret;
}

LL solve(LL n)
{
	LL ret = n%mod;
	for(LL i = 1, j; i <= n; i=j+1)
	{
		j = n/(n/i);
		ret = (ret + ((sum(j)-sum(i-1))%mod) * ((n/i)%mod) % mod) % mod;
	}
	return ret * inv2 % mod;
}

int main ()
{
	LL a, b; init();
	scanf("%lld%lld", &a, &b);
	printf("%lld
", ((solve(b)-solve(a-1)) % mod + mod) % mod);
}