LOJ #6222. 幂数 ！（加强版）
题目链接
题意
给定整数 \(n(1\le n\le 10^{25})\)，求 \(n\) 以内 Powerful Number 的个数，以及它们的和。
题解
Part 1
如果 \(x\) 是一个 Powerful Number，那么它一定可以表示成 \(a^2b^3\) 的形式。
我们限制 \(b\) 不含（大于 \(1\) 的）平方因子，那么表示方案唯一。
证明是简单的，考虑 \(x\) 每一项质因子 \(p^k(k\ge 2)\) 如何分配即可。
如果 \(k\) 是偶数，将 \(p^k\) 全部分配给 \(a^2\)；
如果 \(k\) 是奇数，将 \(p^{k-3}\) 分配给 \(a^2\)，将 \(p^3\) 分配给 \(b^3\)。
显然 \(b\) 的每个质因子的幂次最多为 \(1\)，并且这是唯一的分配方案。证毕。
至此我们可以将所求转化为易于计算的式子：
其中 \(\mu\) 为莫比乌斯函数。
Part 2
先来看 \(cnt\) 怎么算。
我们引入一个重要的优化求和的 trick：
取一对实数 \(x,y\)，使得 \(x^2y^3=n\)。
若 \(i^2j^3\le n\)，则 \(i\le x\) 和 \(j\le y\) 中至少有一个成立。
对 \(i\le x\) 和 \(j\le y\) 的情况分别求和，再减去公共部分。式子如下：
记 \(\displaystyle{\rm S}\mu^2(n)=\sum_{i=1}^n\mu^2(i)\)，那么式子可以写成：
预处理之后可以 \(O(x+y)\) 计算。
然而预处理的复杂度是 \(O(\sqrt[3]n+y)\)，不能接受。
再考虑一个用于计算 \({\rm S}\mu^2\) 的恒等式：
证明：
实际上，\(\displaystyle\mu^2(i)=\sum_{j^2~|i}\mu(j)\)。
设 \(i\) 的最大平方因子为 \(t\)，则 \(\displaystyle\sum_{j^2~|i}\mu(j)=\sum_{j|t}\mu(j)=[t=1]=\mu^2(i)\)。
用上式代换 \(\mu^2(i)\)，化简：
证毕。
（本质上就是对平方因子进行容斥）
如果直接用这个式子计算，那么求 \({\rm S}\mu^2(n)\) 的复杂度是 \(O(\sqrt n)\)。事实上可以继续优化：
设 \(m=\sqrt[3]n\)，套用前面的 trick：
在 \(O(\sqrt n)\) 的预处理之后，我们可以 \(O(\sqrt[3]n)\) 求出 \({\rm S}\mu^2(n)\)。
回忆一下 \(cnt\) 的计算式：
现在我们的复杂度变为 \(\displaystyle O\left(y+\sum_{i=1}^x\sqrt[9]{n/i^2}\right)\)。
积分一下，也就是 \(O(y+\sqrt[9]{nx^7})\)。其中 \(x^2y^3=n\)。
为了平衡复杂度，可以通过均值不等式或者求导来寻找最值。
易得 \(x=O(n^{2/13})\) 时复杂度最优，为 \(O(n^{3/13})\)。
Part 2
\(sum\) 的处理方法与 \(cnt\) 基本一致。我们快速过一遍式子。
记 \({\rm id}^k(n)=n^k,f(n)=\mu^2(n)n^3\)。
其中 \({\rm Sid}^2(n)=\dfrac{n(n+1)(2n+1)}6\) 可以直接计算。
对于 \({\rm S}f(n)\)，我们有：
设 \(g(n)=\mu(j)j^6\)：
其中 \({\rm Sid}^3(n)=\left(\dfrac{n(n+1)}2\right)^2\) 同样可以直接计算。
\({\rm S}g\) 直接预处理即可，复杂度为 \(O(\sqrt[6]n)\)。
Part 3
到这里就做完了。
时间复杂度 \(O(n^{3/13})\)，空间复杂度 \(O(n^{3/13})\)。
代码实现如下。需要注意开方的精度问题。
#include<bits/stdc++.h>
#define rep(i,l,r) for (int i=l; i<=r; ++i)
typedef long long ll;
typedef __int128 Lll;
const int N=6e5+5,M=15005;
int t1; ll cnt,t2; Lll sum;
std::bitset<N>mu2;
int mu[N],smu[N]; Lll g[M],sg[M];
inline Lll read() {
Lll x=0; char ch=getchar();
while (!isdigit(ch)) ch=getchar();
while (isdigit(ch)) x=x*10+(ch^48),ch=getchar();
return x;
}
void write(Lll x) {
if (x>9) write(x/10);
putchar(x%10|48);
}
void prework(const int n,const int m) {
mu[1]=1;
rep(i,2,n) if (!mu[i])
for (int j=i; j<=n; j+=i) mu[j]=(mu[j]<0?1:-1);
for (int i=2,t; (t=i*i)<=n; ++i) if (mu[t])
for (int j=t; j<=n; j+=t) mu[j]=0;
rep(i,1,n) {
smu[i]=smu[i-1];
if (mu[i]) mu2[i]=1,smu[i]+=mu[i];
}
rep(i,1,m) {
sg[i]=sg[i-1];
if (mu[i]) sg[i]+=(g[i]=mu[i]*i*Lll(i)*i*i*i*i);
}
}
inline Lll S2(ll x) { return x*Lll(x+1)*(x<<1|1)/6; }
inline Lll S3(ll x) { return x=x*(x+1)>>1,Lll(x)*x; }
inline ll Sqrt(Lll x) { ll t=sqrtl(x); return t-(Lll(t)*t>x); }
inline int Cbrt(Lll x) { int t=cbrtl(x); return t-(Lll(t)*t*t>x); }
//精确开方，在 sqrtl,cbrtl 基础上微调
int Smu2(const int n) {
const int m=cbrt(n); int s=0,t;
rep(i,1,m) t=n/i,s+=t/i*mu[i]+smu[(int)sqrt(t)];
return s-m*smu[m];
}
Lll Sf(const int n) {
const int m=cbrt(n); Lll s=0; int t;
rep(i,1,m) t=n/i,s+=S3(t/i)*g[i]+sg[(int)sqrt(t)]*i*i*i;
return s-S3(m)*sg[m];
}
int main() {
const Lll n=read();
const int x=pow(n,2./13),y=cbrtl(n/(x*x));
prework(y+100,pow(n,1./6)+100);
rep(i,1,x)
t1=Cbrt(n/(i*i)),
cnt+=Smu2(t1),sum+=Sf(t1)*i*i;
rep(i,1,y) if (mu2[i])
t2=Sqrt(n/(ll(i)*i*i)),
cnt+=t2,sum+=S2(t2)*i*i*i;
printf("%lld\n",cnt-ll(x)*Smu2(y));
write(sum-S2(x)*Sf(y));
return 0;
}