Luogu6788 「EZEC-3」四月櫻花
2020.8.29更新
被命題人@了，原來的做法已經過不去了……
這題是一道很不錯的數論題，對整除分塊的考察及其時間複雜度分析與優化較為深入。
演算法考察
整除分塊，常見數論函式性質。
演算法分析
我們來推一推柿子。
顯然\(d(y)=\sum_{z|y}1\)，因此\(y^{d(y)}=\prod_{z|y}y\)。
我們把\(y\)拆開：
\[\prod_{z|y}y=\prod_{z|y}z\times\prod_{z|y}\frac{y}{z}=\left(\prod_{z|y}z\right)^2=\prod_{z|y}z^2 \]
因此原式可推導如下：
\[\prod_{x=1}^t\prod_{y|x}\prod_{z|y}\frac{z^2}{(z+1)^2} \]
改為列舉\(y\)
\[\left(\prod_{y=1}^t\prod_{z|y}\frac{z^2}{(z+1)^2}\right)^{\left\lfloor\frac{t}{y}\right\rfloor} \]
改為列舉\(z\)
\[\left(\prod_{z=1}^t\frac{z^2}{(z+1)^2}\right)^{\sum_{y=1}^{\left\lfloor\frac{t}{z}\right\rfloor}\left\lfloor\frac{t}{yz}\right\rfloor} \]
我們設\(f(n)=\sum_{i=1}^n \left\lfloor\frac{n}{i}\right\rfloor\)
\[\left(\prod_{z=1}^t\frac{z^2}{(z+1)^2}\right)^{f\left(\left\lfloor\frac{t}{z}\right\rfloor\right)} \]
整除分塊基礎——求解\(f(n)\)
考慮\(f(n)\)如何求解。（有整除分塊基礎的同學可以跳過）
通過觀察，我們可以發現\(\left\lfloor\frac{n}{i}\right\rfloor\)最多有\(O(\sqrt n)\)種不同的取值，且隨\(i\)增大單調不減。
單調性顯然，取值個數簡單證明一下：
- \(i\le \sqrt n\)時，有\(\sqrt n\)個不同的\(i\)，對應至多\(\sqrt n\)個不同的\(\left\lfloor\frac{n}{i}\right\rfloor\)；
- \(i>\sqrt n\)時，發現\(\left\lfloor\frac{n}{i}\right\rfloor<\sqrt n\)，因此也至多有\(\sqrt n\)個不同的取值。
因此我們考慮列舉這\(\sqrt n\)種不同的取值即可求解，單次求解的時間複雜度為\(O(\sqrt n)\)。
求解\(f(n)\)的程式碼如下：
int get_f(int n){
int ret=0;
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
// l,r表示區間[l,r]內的整數i，floor(n/i)相等（均為 n/l）
// l,r 的值的推導本處不展開
// 可以參考下面相關題目的題解
ret=(ret+1ll*(r-l+1)*(n/l))%mod;
// 區間[l,r]內有(r-l+1)個整數，它們的整除值都是 n/l
}
return ret;
}
相關題目：[CQOI2007]餘數求和
回到原問題，現在我們已經可以用\(O(\sqrt n)\)的時間複雜度求出\(f(n)\)的值。
整除分塊進階——整除分塊巢狀
我們再觀察式子：
\[\left(\prod_{z=1}^t\frac{z^2}{(z+1)^2}\right)^{f\left(\color{red}{\left\lfloor\frac{t}{z}\right\rfloor}\right)} \]
我們發現\(f(n)\)的自變數也出現整除，也就意味著會出現一段連續的\(z\)，使得\(\left\lfloor\frac{t}{z}\right\rfloor\)的值相同，即\(f\left({\left\lfloor\frac{t}{z}\right\rfloor}\right)\)也一定相同。
因此對於指數相同的部分，我們只需要知道底數的部分積即可。
觀察底數部分的部分積。
\[\prod_{z=l}^r\frac{z^2}{(z+1)^2}=\frac{l^2}{(r+1)^2} \]
發現這個值可以配合費馬小定理求逆元快速求出。
類比整除分塊的過程，我們列舉\(O(\sqrt t)\)個不同的\(\left\lfloor\frac{t}{z}\right\rfloor\)的值，指數部分可以用整除分塊求解。
這就是整除分塊套整除分塊。
類比杜教篩，時間複雜度為\(O(t^{\frac{3}{4}}\log t)\)，前面是整除分塊套整除分塊的複雜度，即\(\sum_{i=1}^{\sqrt t}O(\sqrt i)+O(\sqrt{\frac{n}{i}})=O(t^{\frac{3}{4}})\)，後面\(\log t\)是過程中求逆元。
演算法優化
本演算法還有進一步優化的空間。
繼續類比杜教篩。在杜教篩中，我們預處理出部分函式值及其字首和，從而將時間複雜度由\(O(n^\frac{3}{4})\)優化到\(O(n^\frac{2}{3})\)。對於本題，我們能不能也預處理出一部分\(f(n)\)，從而優化時間複雜度呢？
因數與倍數——\(f(n)\)求法優化
我們考慮\(\left\lfloor\frac{a}{b}\right\rfloor\)的實際含義，其可以表示不超過\(a\)的數中有多少個是\(b\)的倍數，即\(\sum_{i=1}^a [b|i]\)。
因此我們把\(f(n)\)改寫如下：
\[f(n)=\sum_{i=1}^n \left\lfloor\frac{n}{i}\right\rfloor \]
\[=\sum_{i=1}^n \sum_{j=1}^n [i|j] \]
\[=\sum_{j=1}^n \sum_{i=1}^n [i|j] \]
\[=\sum_{j=1}^n \sum_{i|j}1 \]
\[=\sum_{j=1}^n \sigma(j) \]
其中\(\sigma(n)\)表示\(n\)有多少個因數。
根據\(\sigma(n)\)的定義，我們不難寫出單次詢問\(O(\sqrt n)\)的演算法求出單個\(\sigma(n)\)，但這個顯然不能滿足我們的需求。
現在我們考慮如何批量求出\(\sigma(n)\)。
注意到對每一個數\(i\)，都會對\(i\)的倍數產生\(1\)的貢獻，因此我們再次列舉倍數，即可批量求出\(\sigma(n)\)，最後求一遍字首和即可批量求出\(f(n)\)。
具體實現如下：
for(int i=1;i<=N;i++){
for(int j=1;i*j<=N;j++){
sig[i*j]++;
}
}
時間複雜度為\(\sum_{i=1}^n\frac{n}{i}=O(n\log n)\)，可用調和級數證明。
上述過程也可以用狄利克雷卷積解釋，即\(1\ *\ 1\)，而求狄利克雷卷積的通用方法正是列舉倍數。
我們也不難看出\(\sigma(n)\)是一個積性函式，並且\(\sigma(p^k)=k+1\)，這意味著我們可以線性篩批量求出\(\sigma(n)\)，時間複雜度為\(O(n)\)。
相關題目：[AHOI2005]約數研究
綜上，我們可以用\(O(n\log n)\)的時間複雜度求出前\(n\)個\(f(n)\)。假設我們預處理前\(M\)個\(f(n)\)的值，那麼總時間複雜度為
\[M\log M+\sqrt t\log t+\sum_{i=1}^{t/M}\sqrt{\frac{t}{i}}\approx(2\sqrt{\frac{t^2}{M}}+M)\log t \]
由不等式相關知識得\(M=t^\frac{1}{3}\)，原式取最小值為\(O(t^{\frac{2}{3}}\log t)\)，實測\(M=500000\)較優。
程式碼實現
- 注意\(2.5\times10^9\)超出了
int範圍，要使用unsigned int； - 注意對指數上的整除分塊求和模數應為\(\varphi(p)=p-1\)；
- 注意常數優化。
#include<bits/stdc++.h>
using namespace std;
#define maxn 1000005
#define ui unsigned
template <typename Tp> void read(Tp &x){
int fh=1;char c=getchar();x=0;
while(c>'9'||c<'0'){if(c=='-'){fh=-1;}c=getchar();}
while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+(c&15);c=getchar();}x*=fh;
}
ui n,mod;
ui ans=1;
ui ksm(ui B,ui P){ui ret=1;while(P){if(P&1)ret=1llu*ret*B%mod;B=1llu*B*B%mod;P>>=1;}return ret;}
ui f[maxn];
ui calc(ui x){//整除分塊求f(n)
if(x<=1000000)return f[x];
ui ret=0;
for(ui l=1,r;l<=x;l=r+1){
r=x/(x/l);
ret=(ret+1ll*(r-l+1)*(x/l))%(mod-1);
}
return ret;
}
void preprocess(int N){//預處理求f(n)
for(int i=1;i<=N;i++){
for(int j=1;i*j<=N;j++){
f[i*j]++;
}
}
for(int i=1;i<=N;i++)f[i]=(f[i-1]+f[i])%mod;
}
signed main(){
read(n);read(mod);
preprocess(1000000);
for(ui l=1,r;l<=n;l=r+1){
r=n/(n/l);//整除分塊巢狀
ui invr=1llu*l*ksm((r+1)%mod,mod-2)%mod;
invr=1llu*invr*invr%mod;
ans=1llu*ans*ksm(invr,calc(n/l))%mod;
}
printf("%u\n",ans);
return 0;
}
後記
本題解中附有一些相關題目，其本身可能難度不是很大，但其思想與本題相同。
對於求解不同的數論問題，列舉因數與列舉倍數各有所長，有時可能要反覆轉化。