51Nod1222 最小公倍数计数 数论 Min_25 筛

Posted 2021-02-11 zhouzhendong

原文链接https://www.cnblogs.com/zhouzhendong/p/51Nod1222.html

题意

　　给定 $a,b$， 求

$$sum_{n=a}^b sum_{i=1}^n sum_{j=1}^i [{ m lcm } (i,j) = n]$$

$$a,bleq 10^{11}$$

$${ m Time Limit } = 6s$$

题解

　　本题做法很多。

Min_25 筛

　　先差分一下，转化成求前缀和。

　　先把原题的统计无序数对转化成统计有序数对，最终 $ans‘ = (ans+n)/2$ 即可。

　　设集合 $P$ 表示素数集合。

　　设 $c(n,p)$ 表示最大的使得 $p^{c(n,p)}|n$ 的数。

　　若 ${ m lcm } (i,j) = n$ ，则

$$forall p in P, c(n,p)=max(c(i,p),c(j,p))$$

　　所以，$forall pin P$ ，$c(i,p)$ 和 $c(j,p)$ 共有 $2c(n,p) +1 $ 种取值方法。

　　所以，设

$$n=prod_i p_i^{k_i} (p_iin P)$$

　　则

$$ sum_{i=1}^n sum_{j=1}^i [{ m lcm } (i,j) = n] = prod_t (2k_t+1) $$

　　显然这个式子满足 Min_25 筛的条件，直接筛就好了。

整除分块

$$S(k) = sum_{n=1}^k sigma _0 (n^2) =S(n) = sum_{i=1}^{n} sigma_0(i^2)\\=sum_{i=1}^n sum_{d|i} 2^{omega(d)}\\=sum_{d=1}^nlfloor frac nd floor 2^{omega(d)}$$

整除分块一下，考虑

$$sum_{i=1}^n 2^{omega (i)} = sum_{i=1}^n sum_{d|i} mu(d)^2\\=sum_{d=1}^nlfloor frac nd   floor mu(d) ^2$$

再整除分块，再推

$$sum_{i=1}^n mu (i)^2 = sum_{i=1}^n mu(i) lfloor frac n {i^2} floor$$

于是只要求个 $mu$ 的前缀和。

看起来要杜教筛，但是……

由于这里 $i^2leq n$，所以只需要暴力预处理前缀和就好了。

注意这个做法可能会被卡常数！

复杂度正确的暴力

$$S(k) = sum_{n=1}^ksum_{i=1}^{n}sum_{j=1}^n [{ m lcm}(i,j) = n]$$

$$=sum_{n=1}^ksum_{d=1}^nsum_{i=1}^nsum_{j=1}^n[ijd=n] [gcd(i,j) = 1]$$

$$=sum_{p=1}^kmu(p) sum_{i,j,d} [ijdleq frac{k} {p^2}]$$

先枚举个 $p$ ，再强制 $i,j,d$ 单调不降，枚举 $i$ 再枚举 $j$ ，暴力计算。

证明一下时间复杂度：

$$T(n) = O(sum_{p=1}^{sqrt n} sum_{i=1}^{sqrt[3]{frac{n}{p^2}}} lfloor frac{n}{p^2i} floor) $$

$$O(sum_{i=1}^{sqrt[3]{n}} frac n i) = O(int _{1}^{sqrt[3]{n}} frac n x { m d} x)=O(nint _{1}^{sqrt[3]{n}} x^{-1} { m d} x)$$

代码1 - Min_25 筛

#pragma GCC optimize("Ofast","inline")
#include <bits/stdc++.h>
#define clr(x) memset(x,0,sizeof (x))
using namespace std;
typedef long long LL;
LL read(){
	LL x=0,f=0;
	char ch=getchar();
	while (!isdigit(ch))
		f|=ch==‘-‘,ch=getchar();
	while (isdigit(ch))
		x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
	return f?-x:x;
}
const int Base=1000005,N=Base*2+5;
LL n,cn,a,b,base;
LL h[N],ps[N],cnt;
LL p[N],pcnt;
#define ID(i) ((i)<=base?i:cnt-cn/(i)+1)
LL f(int e){
	return e*2+1;
}
LL g(LL n,LL m){
	LL ans=max(0LL,h[ID(n)]-h[ID(p[m-1])]);
	for (int i=m;i<=pcnt&&p[i]*p[i]<=n;i++){
		LL nn=n/p[i];
		for (int e=1;nn>0;e++,nn/=p[i])
			ans+=f(e)*((e>1)+g(nn,i+1));
	}
	return ans;
}
LL _solve(LL _n){
	cn=n=_n,base=(LL)sqrt(n),cnt=pcnt=0;
	for (LL i=1;i<=n;i=ps[cnt]+1)
		ps[++cnt]=n/(n/i),h[cnt]=ps[cnt]-1;
	p[0]=1;
	for (LL i=2;i<=base;i++)
		if (h[i]!=h[i-1]){
			p[++pcnt]=i;
			LL i2=i*i;
			for (LL j=cnt;ps[j]>=i2;j--)
				h[j]-=h[ID(ps[j]/i)]-(pcnt-1);
		}
	for (LL i=1;i<=cnt;i++)
		h[i]*=3;
	return g(n,1)+1;
}
LL solve(LL n){
	return (_solve(n)+n)/2;
}
int main(){
	a=read(),b=read();
	cout<<solve(b)-solve(a-1)<<endl;
	return 0;
}