fafa法数 - 莫比乌斯反演 - 数论

要说神仙数论题
老爷爷最强无敌
（咦怎么这么押韵的说）
题目大意：
给你$L,R$，求有多少$n\in [L,R]$，满足存在$1\le a<b,c\ge 1$，使得$n=a{b}^c$。$1\le L\le R\le 8\times 10^{16}$
题解：
首先当c确定时，为了避免重复计数，我们要保证a中的质因子的指数小于c。
其次，从这个角度出发，$c\ge 4$都是没有意义的。
因为若$2|c$，则$n=a(b^2{)}^{\frac{c}{2}}$，否则$n=ab(b^2{)}^{\frac{c-1}{2}}$，可以验证这两种情况都是合法的。
若$c=2$，则有$n=a{b}^2>a^3$有$a<\sqrt[3]{n}$，因此（令$p_2(n),p_3(n)$分别表示n中是否含平方、立方因子）对答案的贡献是：
$$\sum _{1\le a<\sqrt[3]{n}}{p}_2(a)\left(\lfloor \sqrt{\frac{n}{a}}\rfloor -a\right)$$
若$c=3$，则这么计算是有问题的，因为例如$1\times 4^3=1\times 8^2$，可能某些$c=3$会被$c=2$计算。
记$n=b{y}^3$，我们知道其中一种可行但可能不合法的$c=2$的情况是$n=by\times y^2$。若想使得其可能合法，就要将$by$中的平方质因子取尽，即要找到一个最大的$k$，满足$k^2|by$，并变换$n=\frac{by}{k^2}\times (ky{)}^2$。
这样仍然可能是没有被计算的，因为$c=2$的情况只会计算$\frac{by}{k^2}<ky$的部分，也就是$\frac{by}{k^2}\ge ky$即$k\le \sqrt[3]{b}$的部分是属于$c=3$的答案但是不会被$c=2$统计到。
因此我们可以类似的枚举$b,k,y$来计算答案：
$$\sum _{1\le b<\sqrt[4]{n}}{p}_3(b)\sum _{b<y\le \sqrt[3]{\frac{n}{b}}}\sum _{1\le k\le \sqrt[3]{b}}\left[k^2|by\right]p_2\left(\frac{by}{k^2}\right)$$
进一步简单的数论转化，令：
$$g=\left(k^2,b\right),k^{\prime }=\frac{k^2}{g},b^{\prime }=\frac{b}{g},y^{\prime }=\frac{y}{k^{\prime }}$$
则：
$$\left[k^2|by\right]=[k^{\prime }g|b^{\prime }gy]=[k^{\prime }|b^{\prime }y]=[k^{\prime }|y]$$
同时：
$$p_2\left(\frac{by}{k^2}\right)=p_2\left(b^{\prime }{y}^{\prime }\right)=p_2\left(b^{\prime }\right)p_2\left(y^{\prime }\right)\left[b^{\prime }\perp y^{\prime }\right]$$
对后面的$\left[b^{\prime }\perp y^{\prime }\right]$做莫比乌斯反演并化简整个式子（省略一些简单推导）：
$$\sum _{1\le b<\sqrt[4]{n}}{p}_3(b)\sum _{1\le k\le \sqrt[3]{b}}{p}_2\left(b^{\prime }\right)\sum _{d|b^{\prime }}\mu (d)\sum _{\frac{b}{k^{\prime }}<y^{\prime }\le \frac{1}{k^{\prime }}\sqrt[3]{\frac{n}{b}}}{p}_2\left(y^{\prime }\right)\left[d|y^{\prime }\right]$$
记$cn{t}_d(m)={\sum }_{i=1}^m{p}_2(id)$，则上式等于：
$$\sum _{1\le b<\sqrt[4]{n}}{p}_3(b)\sum _{1\le k\le \sqrt[3]{b}}{p}_2\left(b^{\prime }\right)\sum _{d|b^{\prime }}\mu (d)\left(cn{t}_d\left(\lfloor \frac{1}{k^{\prime }d}\sqrt[3]{\frac{n}{b}}\rfloor \right)-cn{t}_d\left(\lfloor \frac{b}{k^{\prime }d}\rfloor \right)\right)$$
这样暴力枚举$b,k,d$，可以大概算出复杂度不会并且远远不超过$O\left(n^{\frac{5}{12}}\right)$，预处理前$O\left(n^{\frac{1}{3}}\right)$的$p_2(n)$以及前$O\left(n^{\frac{1}{4}}\right)$的$p_3$和$\mu$，还有$cn{t}_d(m)$即可，其中$O(d)=O\left(n^{\frac{1}{4}}\right),O(m)=O\left(\frac{\sqrt[3]{n}}{d}\right)$。
哇终于写(chao)完了$\mathrm{Q}\mathrm{w}\mathrm{Q}$

#include<bits/stdc++.h>
#define gc getchar()
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define N3 440000
#define N4 100000
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
vector<int> dvs[N3],cnt[N4];
int p2[N3],p3[N3],mu[N3],pri[N3],mnp[N3],np[N3];
inline int gcd(int a,int b) { return a?gcd(b%a,a):b; }
inline int Mysqrt(lint x)
{
	int L=1,R=1000000000,mid=(L+R)>>1;//1e9
	while(L<=R)
	{
		if((lint)mid*mid<=x) L=mid+1;
		else R=mid-1;mid=(L+R)>>1;
	}
	return R;
}
inline int Mycbrt(lint x)
{
	int L=1,R=1000000,mid=(L+R)>>1;//1e6
	while(L<=R)
	{
		if((lint)mid*mid*mid<=x) L=mid+1;
		else R=mid-1;mid=(L+R)>>1;
	}
	return R;
}
inline lint calc(lint n)
{
	if(!n) return 0;lint ans=0;
	for(int a=1;(lint)a*a*a<n;a++)
		if(p2[a]) ans+=Mysqrt(n/a)-a;
	for(int b=1;(lint)b*b*b*b<n;b++) if(p3[b])
		for(int k=1;k*k*k<=b;k++)
		{
			int k2dg=gcd(k*k,b),kp=k*k/k2dg,bp=b/k2dg;
			if(!p2[bp]) continue;
			Rep(i,dvs[bp])
			{
				int d=dvs[bp][i],s=b/kp,t=Mycbrt(n/b)/kp;
				ans+=mu[d]*(cnt[d][t/d]-cnt[d][s/d]);
			}
		}
	return ans;
}
inline int get_mnp_mu(int n)
{
	mnp[1]=1,mu[1]=1;
	for(int i=2,cnt=0;i<=n;i++)
	{
		if(!np[i]) pri[++cnt]=i,mnp[i]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&pri[j]<=n/i;j++)
		{
			int x=pri[j]*i;np[x]=1,mnp[x]=pri[j],mu[x]=-mu[i];
			if(i%pri[j]==0) { mu[x]=0;break; }
		}
	}
	return 0;
}
inline int prelude(lint n)
{
	int n3=Mycbrt(n)+1,n4=Mysqrt(Mysqrt(n))+1;
	get_mnp_mu(n3);
	rep(i,1,n3)
	{
		p2[i]=p3[i]=1;int x=i;
		while(x>1)
		{
			int t=mnp[x],cnt=0;while(x%t==0) x/=t,cnt++;
			if(cnt>1) p2[i]=0;if(cnt>2) p3[i]=0;
			if(!p2[i]&&!p3[i]) break;
		}
	}
	rep(i,1,n4) rep(j,1,n4/i) if(mu[j]) dvs[i*j].pb(j);
	rep(d,1,n4)
	{
		cnt[d].resize(n3/d+1),cnt[d][0]=0;
		rep(i,1,(int)cnt[d].size()-1)
			cnt[d][i]=cnt[d][i-1]+p2[i*d];
	}
	return 0;
}
int main()
{
	lint L,R;cin>>L>>R;prelude(R);
	return cout<<calc(R)-calc(L-1)<<endl,0;
}