TopCoder - 12584 SRM 582 Div 1 SemiPerfectPower (莫比乌斯反演)

题目大意：

给定$L,R$，求$[L,R]$中能够表示为$a\cdot b^c(1\leq a1)$的数(SemiPerfect数)的个数

$R\leq 8\cdot 10^{16}$

解题思路

首先显然可以通过作差转化为求$[1,n]$内的个数

接下来考虑简化$c$的情形

推论：任何一个SemiPerfect数可以表示为$c\leq 3$的形式

证明：若$c>3$，对于$n=a\cdot b^c(c>3)$

当$2|c$时，显然存在形如$n=a\cdot (b^{\frac{c} {2} })^{2}$的表示

当$2\not |c$时，可以表示为$n=(a\cdot c)\cdot (b^{\frac{c-1} {2} })^2$同样合法

接下来考虑对于两种情况分类讨论

为了便于叙述，令$F(n)=\max \lbrace k\in \N^+|\ k^2|n\rbrace$

$G(n)=[\nexists k>1,k^3|n]$

Part1 $c=2$

为了避免重复，强制每一个数$n$的唯一表示为

$n=a\cdot b^2(F(a)=1,a<b)$

由于$aa^3$，即$a<n^{\frac{1} {3} }$

暴力枚举$a$，预处理$n^{\frac{1} {3} }$中所有的$G(i)$即可

Part $c=3$

同样的，限制条件为$n=a\cdot b^3,(G(a)=1,a<b)$，得到$a<n^{\frac{1} {4} }$

但是由于$c=2,3$两部分有重复，还需额外考虑强制$n$不存在形如$n=a’\cdot b’^2$的表示

假设已知$n=a\cdot b^3$，不存在$n=a’\cdot b’^2$的判定条件是

$a’=\frac{n} {F^2(n)}\ge F(n)$，即$F(n)\leq n^{\frac{1} {3} }$

同时由于$F(n)=F(a\cdot b)b$

得到$F(a,b)\leq a^{\frac{1} {3} }$

由于等号右边包含$a$，考虑枚举$a$，易求出$L=F(a),d=\frac{a} {L^2}$，得到$F(a\cdot b)$的另一种表达形式

$F(a\cdot b)=L \cdot \gcd (d,b)\cdot F(\frac{b} {\gcd(d,b)})\leq a^{\frac{1} {3} }$

上面的转化意为：$L$为$a$中已经成对的部分自然取出，然后优先考虑为$D$匹配$b$中的因数成对，对于剩下的部分再重新计算答案

化简该式，得到$L\cdot \gcd(d,b)F(\frac{b} {\gcd(d,b)})\leq a^{\frac{1} {3} }$

式子包含$\gcd $，似乎具有莫比乌斯反演的性质

考虑计算$b\in [a+1,(\frac{n} {a})^{\frac{1} {3} }]$的数量

观察到$a^{\frac{1} {3} }\leq n^{\frac{1} {12} }$，上限只有$25$左右，可以考虑直接枚举$F(\frac{b} {\gcd(b,d)})$

令枚举的$g=\gcd(b,d),F(\frac{b} {g})=f$，计算$\gcd(\frac{b} {g},\frac{d} {g})=1,g\cdot f\cdot L\leq a^{\frac{1} {3} }$的$b$的数量

考虑直接从$\frac{b} {g}$中把$f^2$的因数提取出来，令$\begin{aligned} L'=\lceil \frac{a+1} {gf^2}\rceil,R'=\lfloor \frac{(\frac{n} {a})^{\frac{1} {3} }} {gf^2}\rfloor \end{aligned}$，令$i=\frac{b} {gf^2},x=\frac{d} {g}$，得到新的限制条件式子为

$\gcd(x,f)=1,\gcd(i,x)=1,F(i)=1,i\in[L’,R’]$

在确定了$g,f$之后，需要考虑的限制就是$\gcd(i,x)=1,F(i)=1,i\in[L’,R’]$

由于包含$\gcd$，不妨用莫比乌斯反演计算该式，得到表达式为

$Sum=\sum_{k|x}\mu(k)\sum_{i\in [L’,R’]} [k|i\ \text{and}\ F(i)=1]$

对于$k$，计算$\sum_{i\in[L’,R’]}[k|i\ \text{and}\ F(i)=1]$可以归纳为计算

$\sum_{i=\frac{n} {k} } [F(ik)=1]$，一共有$n^{\frac{1} {3} }\ln n$种不同的权值，可以暴力预处理得到

枚举$d$的因数对于所有的上面的式子计算，可能的$g,f$并不多，可以直接枚举

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define pb push_back
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)

class SemiPerfectPower {
public:
	static const int N=450000,M=20000;
	int w[N],notpri[N],pri[N],pc,F[N],G[N];
	vector <int> S[N],Fac[N];
	ll gcd(ll a,ll b){ return b==0?a:gcd(b,a%b); }
	SemiPerfectPower(){
		// 预处理F,G ,w[i]=\mu(i) S[k][i]=\sum_{j \in [1,i]}  F(i,k)=1
		rep(i,1,sqrt(N)) rep(j,1,(N-1)/i/i) F[i*i*j]=i;
		rep(i,2,pow(N,1/3.0)) rep(j,1,(N-1)/i/i/i) G[i*i*i*j]=1;
		rep(i,1,M-1) for(int j=i;j<M;j+=i) Fac[j].pb(i);
		w[1]=1;
		rep(i,2,N-1) {
			if(!notpri[i]) pri[++pc]=i,w[i]=-1;
			for(int j=1;i*pri[j]<N && j<=pc;++j) {
				notpri[i*pri[j]]=1;
				if(i%pri[j]==0) {
					w[i*pri[j]]=0;
					break;
				} 
				w[i*pri[j]]=-w[i];
			}
		}
		rep(i,1,N-1) if(w[i]) {
			S[i].resize(N/i+1);
			rep(j,1,(N-1)/i) S[i][j]=S[i][j-1]+(F[i*j]==1);
		}
	}

	ll Solve2(ll n){
		ll ans=0,UX=pow(n,1/3.0);
        // 防止浮点误差
		if((UX+1)*(UX+1)*(UX+1)<=n) UX++;
		if(UX*UX*UX>n) UX--;
		rep(i,1,UX) if(F[i]==1) {
			ll UY=sqrt(n/i);
            // 防止浮点误差
			if(i*(UY+1)*(UY+1)<=n) UY++;
			if(i*UY*UY>n) UY--;
			ans+=max(0ll,UY-i);
		}
		return ans;
	}

	ll Solve3(ll n){
		ll UX=pow(n,0.25); 
		// 枚举c的上界
		ll ans=0;
		if((UX+1)*(UX+1)*(UX+1)*(UX+1)<=n) UX++;
		rep(x,1,UX) if(!G[x]) {
			ll UY=pow(n/x,1.0/3.0),U=pow(x,1/3.0);
			
			// 防止浮点误差
			if((U+1)*(U+1)*(U+1)<=x) U++;
			if(U*U*U>x) U--;
			if(x*(UY+1)*(UY+1)*(UY+1)<=n) UY++;
			if(x*UY*UY*UY>n) UY--;

			int L=F[x],D=x/L/L;
			for(int G:Fac[D]) {
				for(int g:Fac[G]) {
					if(g*L>U) break;
					rep(f,1,U/g/L) if(gcd(f,D/g)==1) {
						int L=x/G/f/f,R=UY/G/f/f;
						ans+=w[G/g]*(S[G/g][R]-S[G/g][L]);
					}
				}
			}
		}
		return ans;
	}
	ll Solve(ll n) {
		return Solve2(n)+Solve3(n);
	}
	ll count(ll L,ll R) {
		return Solve(R)-Solve(L-1);
	}
};