Processing math: 100%

51nod1238 最小公倍数之和 V3 莫比乌斯函数 杜教筛

题意:求ni=1nj=1lcm(i,j).
  题解:虽然网上很多题解说用mu卡不过去,,,不过试了一下貌似时间还挺充足的,。。也许有时间用phi试试?
  因为是用的莫比乌斯函数求的,所以推导比大部分题解多。。。而且我写式子一般都比较详细,所以可能看上去很多式子,实际上是因为每一步都写了,几乎没有跳过的。所以应该都可以看懂的。
  末尾的e函数是指的e[1]=1,e[x]=0(x!=1)这样一个函数
  ni=1nj=1lcm(i,j)
  ni=1ni=1ijgcd(i,j)
  枚举gcd
  nd=1ndi=1ndj=1[gcd(i,j)==d]ijd
  因为(ijd2d=ijd),所以:
  nd=1ndi=1ndj=1[gcd(i,j)==d]ijd
  nd=1dndi=1ndj=1ij[gcd(i,j)==1]
  nd=1dndi=1ndj=1ijk|gcd(i,j)μ(k)
  枚举k,再枚举k的倍数。
  nd=1dndk=1μ(k)ndki=1ikndkj=1jk
  设S(n)=ni=1i
  nd=1dndk=1μ(k)k2S(ndk)
  枚举T=dk
  nT=1S(nT)2k|Tμ(k)k2Tk
  nT=1S(nT)2k|Tμ(k)kT
  nT=1S(nk)2Tk|Tμ(k)k
  设f(T)=Tk|Tμ(k)k,卷上id2,因为S(nk)可以数论分块,所以我们只需要快速求出区间[l,r]内的f之和即可,显然求出f的前缀和即可解决问题

  (fid2)(n)=i|nf(i)n2i2=i|nik|iμ(k)kn2i2
  i|nk|iμ(k)kn2i=ni|nk|iμ(k)kni
  设h(i)=k|iμ(k)k,则原式:
  ni|nh(i)ni=n(hid)(n)
  (fid2)(n)=n(hid)(n)
  h(n)=k|nμ(k)k=(μid)1
  fid2=n[(μid)1id]=n[(μid)id1]
  其中(μid)id=i|nμ(i)ini=ni|nμ(i)=e
  所以
  n[(μid)id1]=n[e1]=n
  带入杜教筛的式子:
  g(1)S(n)=ni=1(fg)(i)ni=2g(i)S(ni)
  =ni=1ini=2i2S(ni)
  然后直接上杜教就可以了.
  其实还有一个问题。。。一开始预处理的前缀和怎么求?
  要知道前缀和,首先要求出f.
  因为f(T)=Tk|Tμ(k)k,所以如果我们可以快速求出k|Tμ(k)k,然后就只需要再O(n)的乘上T就可以了.
  我们先预处理出μ(k),然后对于每一个k,枚举它的倍数,统计贡献。那么复杂度为 n1+n2+...+nn=nlogn(此处的n为原题面的23次方,即要预处理的f个数)

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define LL long long
#define RL register LL
#define AC 3000
#define ac 5000000
#define p 1000000007LL
//#define h(x) ((x <= block) ? sum[x] : S[n / x])

LL n, ans, block; 
LL mu[ac], S[AC], sum[ac], inv[AC];
int pri[ac], tot;
bool z[ac], vis[AC];

inline LL h(LL x)
{
	return ((x <= block) ? sum[x] : S[n / x]);
}

inline void up(LL & a, LL b)
{
	a += b;
	if(a >= p) a -= p;
	if(a <= -p) a += p;
}

LL count(LL l, LL r){
	return (r - l + 1) % p * ((r + l) % p) % p * inv[2] % p;
}

void pre()
{
	scanf("%lld", &n), block = pow(n, 0.66666);
	mu[1] = 1;
	for(R i = 2; i <= block; i ++)
	{
		if(!z[i]) pri[++ tot] = i, mu[i] = -1;
		for(R j = 1; j <= tot; j ++)
		{
			int now = pri[j];
			if(i * now > block) break;
			z[i * now] = true;
			if(!(i % now)) break;
			mu[i * now] = - mu[i];
		}
	}
	inv[1] = 1;
	for(R i = 2; i <= 10; i ++) inv[i] = (p - p / i) * inv[p % i] % p;
	for(R i = 1; i <= block; i ++)//枚举mu(i)
		for(R j = 1; j; j ++)//枚举i的倍数
		{
			if(j * i > block) break;
			up(sum[i * j], mu[i] * i % p);
		}
	for(R i = 1; i <= block; i ++) sum[i] = sum[i] * i % p; 
	for(R i = 1; i <= block; i ++) up(sum[i], sum[i - 1]);//算出f数组后还要统计前缀和
}

LL get(LL x)
{
	x %= p;
	return x * (x + 1) % p * (2 * x + 1) % p * inv[6] % p;
}

void cal(LL x)
{
	if(x <= block || vis[n / x]) return ;
	LL rnt = count(1, x);
	for(RL i = 2, lim, now; i <= x; i = lim + 1)
	{
		lim = x / (x / i), now = x / i, cal(now);
		up(rnt, - ((get(lim) - get(i - 1)) % p * h(now) % p));
	}
	S[n / x] = rnt, vis[n / x] = true;
}

void work()
{
	for(RL i = 1, lim, now, x; i <= n; i = lim + 1)
	{
		lim = n / (n / i), now = n / i, x = count(1, now);
		up(ans, (x % p * x % p * ((h(lim) - h(i - 1)) % p) % p));
	}
	printf("%lld\n", (ans + p) % p);
}

int main()
{
	//freopen("in.in", "r", stdin);
	pre();
	cal(n);
	work();
//	fclose(stdin);
	return 0;
}
posted @   ww3113306  阅读(348)  评论(1编辑  收藏  举报
知识共享许可协议
本作品采用知识共享署名-非商业性使用-禁止演绎 3.0 未本地化版本许可协议进行许可。
点击右上角即可分享
微信分享提示