原文链接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(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) = Oleft(sum_{p=1}^{sqrt n} sum_{i=1}^{sqrt[3]{frac{n}{p^2}}} sqrt{leftlfloor frac{n}{p^2i} ight floor} ight) $$
$$ecause Oleft(sum_{i=1}^{sqrt[3]{n}}sqrt{ frac n i} ight) = Oleft(int _{1}^{sqrt[3]{n}} sqrt{frac n x} { m d} x ight)\=Oleft(sqrt{n}int _{1}^{sqrt[3]{n}} x^{-0.5} { m d} x ight)=Oleft(sqrt{n}cdot 0.5 left(sqrt[3]{n} ight) ^{0.5} ight)=Oleft(n^{frac 2 3 } ight)$$
$$ herefore T(n) = Oleft(sum_{p=1}^{sqrt n} sum_{i=1}^{sqrt[3]{frac{n}{p^2}}} sqrt{lfloor frac{n}{p^2i} floor} ight)=Oleft(sum_{p=1}^{sqrt{n}} left(frac n {p^2} ight) ^{frac 2 3 } ight)=Oleft(n^{frac 2 3}int _1^sqrt n x^{-frac 4 3} { m d} x ight) \=Oleft( n^frac 23 cdot -1cdot (left (sqrt n ight)^{-frac 13} - 1) ight) = Oleft(n^frac 23 (1-n^{-frac 16}) ight) = Oleft(n^frac 23 ight)$$
没想到这个“暴力”的复杂度居然是对的！代码短小精悍，吊打前几种做法了！
代码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;
}
代码2 - 整除分块
#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 N=11000005;
const LL nb=1e11,un=N-5;
int p[N/5],pcnt;
bool vis[N],_u2[N];
char u[N];
int u2[N],w[N];
unordered_map <LL,LL> U2,W;
LL A;
void get_prime(int n){
pcnt=0,u[1]=1;
for (register int i=2;i<=n;i++){
if (!vis[i])
p[++pcnt]=i,u[i]=-1,w[i]=1;
for (register int j=1;j<=pcnt&&p[j]*i<=n;j++){
vis[i*p[j]]=1;
if (i%p[j])
u[i*p[j]]=-u[i],w[i*p[j]]=w[i]+1;
else {
u[i*p[j]]=0,w[i*p[j]]=w[i];
break;
}
}
}
}
void prework(int n){
get_prime(n);
for (int i=1;i<=n;i++){
_u2[i]=u[i]*u[i];
// u[i]+=u[i-1];
u2[i]=u2[i-1]+_u2[i];
w[i]=1LL<<w[i];
w[i]+=w[i-1];
}
}
LL getU2(LL n){
if (n<=un)
return u2[n];
if (U2.count(n))
return U2[n];
LL ans=0;
ans=0;
for (LL i=1;i*i<=n;i++)
ans+=u[i]?u[i]*(n/(i*i)):0;
return U2[n]=ans;
}
LL getW(LL n){
if (n<=un)
return w[n];
if (W.count(n))
return W[n];
LL ans=0;
LL pr=0,nw,sq=(LL)sqrt(n);
for (LL i=1;i<=sq;i++)
ans+=_u2[i]?n/i:0;
pr=u2[sq];
for (LL i=sq+1,j;i<=n;i=j+1){
j=n/(n/i);
nw=j<=un?u2[j]:getU2(j);
ans+=(nw-pr)*(n/i);
pr=nw;
}
return W[n]=ans;
}
LL getS(LL n){
if (!n)
return 0;
LL ans=0,pr=0,nw;
for (LL i=1,j;i<=n;i=j+1){
j=n/(n/i);
nw=j<=un?w[j]:getW(j);
ans+=(nw-pr)*(n/i);
pr=nw;
}
return ans;
}
LL solve(LL n){
if (!n)
return 0;
A=n;
return (getS(n)+n)/2;
}
int main(){
prework(N-5);
LL a=read(),b=read();
cout<<solve(b)-solve(a-1)<<endl;
return 0;
}
代码3 - 优秀的暴力
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N=5e5+10;
int mu[N],l,r;
int calc(int x){
int re=0;
for (int i=1;i*i<=x;i++)
if (mu[i]){
int tmp=0,m=x/(i*i);
for (int j=1;j*j*j<=m;j++){
tmp+=(m/(j*j)-j)*3+1;
for (int k=j+1;j*k*k<=m;k++)
tmp+=(m/(j*k)-k)*6+3;
}
re+=mu[i]*tmp;
}
return (re+x)/2;
}
signed main(){
mu[1]=1;
for (int i=1;i<N;i++)
for (int j=i+i;j<N;j+=i)
mu[j]-=mu[i];
scanf("%lld%lld",&l,&r);
printf("%lld",calc(r)-calc(l-1));
return 0;
}