TopCoder - 12584 SRM 582 Div 1 SemiPerfectPower (莫比乌斯反演)
TopCoder - 12584 SRM 582 Div 1 SemiPerfectPower (莫比乌斯反演)
题目大意：
给定\(L,R\)，求\([L,R]\)中能够表示为\(a\cdot b^c(1\leq a<b,c>1)\)的数(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)\)
由于\(a<b\)，所以显然\(n>a^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\)并不多，可以直接枚举
#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);
}
};
