题解：P15378 二次根式 / sosqrt
bcdmwSjy
·
·
题解
首先，我们发现 F^2\times G=N!，这非常容易得出，所以我们只需要求出 F 和 G 中的一个即可，这里我们先求出 F。
对于每个正整数 a，考虑所有 i=a^2\times b\le n，且 b 无平方因子的数，这样的 i 的个数即为 c(a)，则 F=\prod\limits_{a=1}^{\sqrt N}a^{c(a)}。
考虑如何计算不超过 $n$ 的无平方因子的数的个数，容斥一下得结果为 $n-\sum\limits_{p_1\in P}\lfloor\frac{n}{p_1^2}\rfloor+\sum\limits_{p_1,p_2 \in P}\lfloor\frac{n}{p_1^2p_2^2}\rfloor-\dots$，由莫比乌斯函数的定义可以得到原式的值为 $\sum\limits_{d=1}^{\sqrt n}\mu(d)\lfloor\frac{n}{d^2}\rfloor$。
这样，我们就可以快速计算 $F$ 了，时间复杂度为 $O(\sqrt N\log N)$，常数非常小。算出 $F$ 后，$G=\dfrac{N!}{F^2}$，$N!$ 可以分段打表计算，但是赛时我以为洛谷打不了几个数，常数太大过不了，就~~不小心~~写了个快速阶乘算法上去。
代码：
```cpp
constexpr int mod=1000000007;
inline ll qpow(ll a,ll b){
ll res=1;
for (;b;b>>=1,a=a*a%mod) if (b&1) res=res*a%mod;
return res;
}
int mu[40001];
bitset<40001> np;
vector<int> pr;
void init(int n){
np[0]=np[1]=true;
mu[1]=1;
for (int i=2;i<=n;i++){
if (not np[i]){
pr.push_back(i);
mu[i]=-1;
}
for (auto j:pr){
int s=i*j;
if (s>n) break;
np[s]=true;
if (i%j==0){
mu[s]=0;
break;
}
mu[s]=-mu[i];
}
}
}
inline int sqfree(int n){
if (n<=0) return 0;
int res=0;
for (int d=1;d*d<=n;d++){
if (mu[d]) res+=mu[d]*(n/(d*d));
}
return res;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(nullptr);
cout.tie(nullptr);
init(40000);
factorial::init(mod); // 这里可以换成你喜欢的阶乘算法
int T;
cin>>T;
while (T--){
int n;
cin>>n;
ll ans=1;
for (int i=1;i*i<=n;i++){
int m=n/(i*i);
ans=ans*qpow(i,sqfree(m))%mod;
}
cout<<ans<<" "<<factorial::calc(n)*qpow(qpow(ans,mod-2),2)%mod<<"\n";
}
return 0;
}
```
都写到这里了，那就再说一下快速阶乘算法吧，其实这个算法没有那么难以理解。
我们先分块，设块长为 $B$，再设 $f_m(x)=\prod\limits_{i=x+1}^{x+m}i$，只要计算出 $f_B(0),f_B(B),f_B(2B),\dots,f_B(B^2)$，把它们乘起来，再暴力计算剩下的部分即可。
我们发现我们只需要 $f$ 的一些点值，不需要维护 $f$ 的系数。假设我们已经求出了 $f_m$，考虑如何得到 $f_{2m}$，这也非常简单，$f_{2m}(x)=f_m(x)\times f_m(x+m)$（这里的乘法是直接把系数相乘），只需要进行[连续点值平移](https://www.luogu.com.cn/problem/P5667)就好了。
但是我们现在维护的点值并不连续，但是连个点之间差是一定的，只需要设 $p_m(x)=f_m(Bx)$，这样 $p_m$ 的点值就是连续的了。
在实现中，我取了 $B=2^{\lceil\log_2(\sqrt{0.5mod})\rceil}$（有威尔逊定理优化，只需要处理到 $\frac12mod$），这样就没有了从 $p_m$ 到 $p_{m+1}$ 的过程，在不影响复杂度的情况下使代码更简洁。
我们发现求出的多项式在每次计算阶乘时是相同的，可以先预处理一次，然后后续直接调用，这样我们就得到了一个 $O(\sqrt n\log n)$ 预处理 $O(\sqrt n)$ 查询的阶乘算法。
给个代码吧：
```cpp
namespace factorial{
MPoly f;
int B;
void init(const int mod){
if (mod==f.mod) return;
f.setmod(mod);
f.resize(2);
B=1<<int(ceil(0.5*log2(mod*0.5)));
int r=qpow(B,mod-2,f.fm);
f.a[0]=1,f.a[1]=B+1;
for (int l=2;l<=B;l<<=1){
f.extend(pointValueTranslation(f,(l>>1)+1));
f.mulcoefficient(pointValueTranslation(f,r*ll(l>>1)%f.fm));
f.resize(l+1);
}
}
int calc(int n){
if (n>f.mod) return 0;
if (n>f.mod-1-n){
int res=qpow(calc(f.mod-1-n),f.mod-2,f.fm);
return res?((f.mod-n)&1?f.mod-res:res):0;
}else{
ll res=1;
for (int i=0;i<n/B;i++) res=res*f.a[i]%f.fm;
for (int i=n/B*B+1;i<=n;i++) res=res*i%f.fm;
return res;
}
}
}
```