P5221 Product 另解
推式子的方向不一样，最后呈现的结果不一样。
为了方便，\(\log\) 均表示以原根 \(g\) 为底的离散对数。
\[\begin{aligned}
\prod_{i=1}^n\prod_{j=1}^n\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)} &= {\prod_{i=1}^n\prod_{j=1}^n (i \times j)} \times {\prod_{i=1}^n\prod_{j=1}^n \gcd (i,j)^{-2}}
\end{aligned}
\]
拆开分别考虑
前一项直接处理。
\[\begin{aligned}
\prod_{i=1}^n\prod_{j=1}^n (i \times j) &= \prod _{i=1} ^ n i^n \prod _{j=1} ^n j \\ &= (n!)^n \prod _{i=1} ^ n i ^ n \\ &= (n!)^{2n}
\end{aligned}
\]
后一项看着 \(\prod\) 就烦，考虑取离散对数。（此处省略指数上的 \(-2\)）
\[\begin{aligned}
\log \prod _{i=1} ^{n} \prod _{j=1} ^n \gcd (i,j) &= \sum _{i=1} ^n\sum_{j=1}^n \log \gcd (i,j) \\&= \sum _{d=1} ^ n \log d \sum _{i=1} ^n \sum _{j=1} ^n [\gcd (i,j)=1] \\&= \sum _{d=1} ^n \log d \sum _{i=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} \sum _{j=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} [\gcd (i,j)=1] \\&= \sum _{d=1}^n \log d \sum _{i=1} ^ {\left \lfloor \frac{n}{d} \right \rfloor} \sum _{j=1} ^ {\left \lfloor \frac{n}{d} \right \rfloor} \sum _{l \mid \gcd (i,j) } \mu (l) \\ &= \sum _{d=1} ^n \log d \sum _{i=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} \sum _{j=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} \sum _{l=1} ^n \mu (l) [l \mid i][l \mid j] \\&= \sum _{d=1} ^n \log d \sum _{l=1} ^n \mu(l) \sum _{i=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} [l \mid i] \sum _{j=1} ^{\left \lfloor \frac{n}{d} \right \rfloor} [l \mid j] \\&= \sum _{d=1} ^n \log d \sum _{l=1} ^ n \mu (l) {\left \lfloor \frac{n}{dl} \right \rfloor}^2 \\&= \sum_{T=1} ^n {\left \lfloor \frac{n}{T} \right \rfloor}^2 \sum _{d \mid T } \log d \times \mu (\frac{T}{d})
\end{aligned}
\]
这个式子直接计算枚举 \(d\) 是 \(O(n \ln n)\) 的，前提是每个 \(d\) 只枚举到有用的。但是一般枚举一个数的因子是 \(O(\sqrt n)\) 的，凭空增加复杂度。于是考虑质因数分解优化。发扬人类智慧。
表面上计算每个数的质因数分解是 \(\sum \limits _{i=1} ^n \sqrt i \approx 6\times 10^8\) 的，实际上用筛了质数之后优化，再使用这段程序测算后发现只有 \(1 \times 10^6\)。
测试代码
pri =[]
def breakdown(x:int) -> int:
ans = 0
for i in pri:
ans +=1
if(x%i==0):
while x%i==0:
ans+=1
x//=i
if(i*i>x): break
if x>1:ans+=1
return ans
n=int(1000000)
r=0
not_prime = [False] * n
def pre(n):
for i in range(2, n + 1):
if not not_prime[i]:
pri.append(i)
for pri_j in pri:
if i * pri_j > n:
break
not_prime[i * pri_j] = True
if i % pri_j == 0:
break
for i in range(n):
r+=breakdown(i)
print(r)
至于求 \([1,n]\) 离散对数，考虑 BSGS 共用哈希表，设往哈希表中插入 \(B\) 个数，只有质数是需要查询哈希表，总共有 \(\frac{nP}{B\ln n}\) 次查询。当 \(B = \sqrt \frac{nP}{\ln n}\) 时有最优复杂度 \(O(\sqrt \frac{nP}{\ln n}) \approx 2 \times 10 ^6\)。
然后我们有了 \(1 \sim n\) 的质因数分解，只需要 DFS 遍历所有因子即可。
复杂度比较玄学，注意常数优化。
AC 代码
#include <bits/stdc++.h>
#include <bits/extc++.h>
using namespace std;
const int mod=104857601,G=3,phi=mod-1;
const int N=1e6+10;
const int B=1200000;
const int S=B-1,step=100867123;//step = G^{B-1}
bool __st__;
int pr[N/10];//质数表
bitset<N> vis;
short mu[N];
int lg[N];//以 G 为底的离散对数
int n;
namespace HST{
const int HSTMOD=2524939;
struct FastMod{
typedef unsigned long long ull;
typedef __uint128_t lll;
ull b,m;
void init(ull _b) {this->b=_b,m=ull((lll(1)<<64)/_b);}
ull operator()(ull a){
if(b==1) return 0;
ull q=(ull)((lll(m)*a)>>64),r=a-q*b;
return r>=b?r-b:r;
}
}moder;
struct node{
int key,val,nxt;
}pool[B+1];
int pcnt;
struct hash_table{
signed char s;
int head[HSTMOD];
void insert(int x,int w){
int t=moder(x);
for(int i=head[t];i;i=pool[i].nxt){
if(pool[i].key==x){
pool[i].val=w;
return ;
}
}
pcnt++;
pool[pcnt].key=x;
pool[pcnt].val=w;
pool[pcnt].nxt=head[t];
head[t]=pcnt;
}
void init(){
moder.init(HSTMOD);
}
int qry(int x){
assert(x>=0);
int t=moder(x);
for(int i=head[t];i;i=pool[i].nxt){
if(pool[i].key==x){
return pool[i].val;
}
}
return -1;
}
};
}
int fpow(int x,int p){
int res=1;
while(p){
if(p&1) res=1ll*res*x%mod;
x=1ll*x*x%mod;
p>>=1;
}
return res;
}
HST::hash_table hs;//G ^ i -> i
int BSGS(int b){
int y=fpow(b,mod-2);
for(int i=0;;i++){
if(hs.qry(y)!=-1){
return i*S-hs.qry(y);
}
y=1ll*y*step%mod;
}
}
void euler(){
hs.init();
int x=1;
{//初始化哈希表
for(int i=1;i<=B;i++){
x=x*G%mod;
hs.insert(x,i);
}
}
int cnt=0;
for(int i=2;i<N;i++){
if(!vis.test(i)){
mu[i]=-1;
lg[i]=BSGS(i);
pr[++cnt]=i;
}
for(int j=1;pr[j]*i<N;j++){
vis.set(pr[j]*i);
lg[pr[j]*i]=lg[pr[j]]+lg[i];
if(i%pr[j]==0){
mu[pr[j]*i]=0;
break;
}else{
mu[pr[j]*i]=mu[pr[j]]*mu[i];
}
}
}
mu[1]=1;
}
int ds[17][2];
int siz;
void breakdown(int x){
siz=0;
if(x==1){
ds[siz][0]=1;
ds[siz][1]=0;
siz++;
return;
}
for(int i=1;pr[i]*pr[i]<=x;i++){
if(x%pr[i]==0){
int cnt=0;
while(x%pr[i]==0){
cnt++;
x/=pr[i];
}
ds[siz][0]=pr[i];
ds[siz][1]=cnt;
siz++;
}
}
if(x>1){
ds[siz][0]=x;
ds[siz][1]=1;
siz++;
}
}
int moder(int t){
return t<0?t+phi:t;
}
int dfs(int id,int d,int T){
if(id==siz){
return 1ll*lg[d]*moder(mu[T/d])%phi;
}
int res=0;
int base=1;
for(int i=0;i<=ds[id][1];i++){
res+=dfs(id+1,d*base,T);
if(res>=phi) res-=phi;
base*=ds[id][0];
}
return res;
}
int alpha(int T){
breakdown(T);
return 1ll*dfs(0,1,T)*(n/T)%phi*(n/T)%phi;
}
bool __ed__;
int main(){
cerr<<(&__ed__-&__st__)/1024.0/1024.0<<endl;
euler();
cin>>n;
int res=1;
for(int i=1;i<=n;i++){
res=1ll*res*i%mod;
}
res=fpow(res,2*n);
int ww=0;
for(int i=1;i<=n;i++){
ww+=alpha(i);
if(ww>=phi)ww-=phi;
}
res=1ll*res*fpow(fpow(G,ww),mod-3)%mod;
printf("%d\n",res);
cerr<<clock()<<endl;
return 0;
}
最后以 457ms 的最大时间，30.10MB 的最大空间险通过。
浙公网安备 33010602011771号