1. 程式人生 > 其它 >[SDOI2017]數字表格 & [MtOI2019]幽靈樂團

[SDOI2017]數字表格 & [MtOI2019]幽靈樂團

莫比烏斯反演

P3704 [SDOI2017]數字表格

首先根據題意寫出答案的表示式

\[\large\prod_{i=1}^n\prod_{j=1}^mf_{\gcd(i,j)} \]

按常規套路改為列舉 \(d=\gcd(i,j)\)
(不妨設 \(n\le m\)

\[\large\prod_{d=1}^n{f_d}^{\sum_{i=1}^n\sum_{j=1}^m~[(i,j)=d]} \]

指數上的式子很熟悉了,單獨拿出來推一下

\[\begin{aligned} \sum_{i=1}^n\sum_{j=1}^m[(i,j)=d] &=\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[(i,j)=1]\\ &=\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}\sum_{t\mid (i,j)}\mu(t)\\ &=\sum_{t=1}^{n/d}\mu(t)\sum_{i=1}^{n/d}\sum_{j=1}^{n/d}[t|i][t|j]\\ &=\sum_{t=1}^{n/d}\mu(t)(n/dt)(m/dt) \end{aligned}\]

代回原式

\[\large\prod_{d=1}^n{f_d}^{\sum_{t=1}^{n/d}\mu(t)(n/dt)(m/dt)} \]

\(k=dt\) ,改為列舉 \(k,d\)

\[\large\prod_{k=1}^n\prod_{d\mid k}{f_d}^{\mu(k/d)(n/k)(m/k)} \]\[\large\prod_{k=1}^n\left(\prod_{d\mid k}{f_d}^{\mu(k/d)}\right)^{(n/k)(m/k)} \]\[\large\prod_{k=1}^n\left(\prod_{d\mid k}{f_{k/d}}^{\mu(d)}\right)^{(n/k)(m/k)} \]

括號裡的東西可以 \(O(n\log n)\)

預處理,然後整除分塊計算就可以了

時間複雜度 \(O(n\log n+T\sqrt n\log n)\)
\(\log\) 是因為用到了快速冪

到這裡已經可以 AC 了
但還是介紹一下 \(O(n\log\log n+T\sqrt n\log n)\) 的做法

可以理解為 dp ,設:

\[\large s_{i,n}=\prod_{d\mid n,d 只含前 i 種質因子}{f_{n/d}}^{\mu(d)} \]\[\large inv_{i,n}={s_{i,n}}^{-1}=\prod_{d\mid n,d 只含前 i 種質因子}{f_{n/d}}^{-\mu(d)} \]

列舉質數 \(p_i\)

\(p_i\nmid n\) 時,顯然 \(s_{i,n}=s_{i-1,n},inv_{i,n}=inv_{i-1,n}\)

\(p_i\mid n\) 時,需要考慮 \(d\) 中含因子 \(p_i\) 的個數

\(p_i\nmid d\) 的情況對 \(s_{i,n}\) 的貢獻即為 \(s_{i-1,n}\)

\({p_i}^2\mid d\)\(\mu(d)=0\) ,因此這時不產生貢獻

\(p_i\mid d\)\({p_i}^2\nmid d\) ,設 \(d=p_id'\) ,則貢獻為

\[\prod_{d'\mid (n/p_i)}{f_{n/p_id'}}^{\mu(p_id')} \]

\(p_i\)\(d'\) 互質,因此 \(\mu(p_id')=\mu(p_i)\mu(d')=-\mu(d')\)

代回之後發現貢獻就是 \(inv_{i-1,n/p_i}\)

\(s_{i,n}=s_{i-1,n}\times inv_{i-1,n/p_i}\)
顯然 \(inv_{i,n}=inv_{i-1,n}\times s_{i-1,n/p_i}\)

綜上,轉移方程為

\[s_{i,n}=\begin{cases}s_{i-1,n}&(p_i\nmid n)\\ s_{i-1,n}\times invs_{i-1,n/p_i} & (p_i\mid n)\end{cases} \]\[inv_{i,n}=\begin{cases}inv_{i-1,n}&(p_i\nmid n)\\ inv_{i-1,n}\times s_{i-1,n/p_i} &(p_i\mid n)\end{cases} \]

初始狀態為 \(s_{0,n}=f_n,inv_{0,n}={f_n}^{-1}\)
但是直接求逆元會導致複雜度升到 \(O(n\log P)\) ,就前功盡棄了
所以還要用一下線性求逆元的方法

#include<stdio.h>
const int N=1000010,P=1e9+7; int T,n,m,ans,t1,t2;
int t,cnt,suf,I,p[100000],s[N],inv[N],v[N],pre[N];
inline int min(int x,int y) { return x<y?x:y; } 
inline int power(int x,int y) {
	int s=1;
	while (y) (y&1)&&(s=1ll*s*x%P),x=1ll*x*x%P,y>>=1;
	return s;
}

int main() {
	s[1]=pre[0]=pre[1]=inv[0]=suf=1;
	for (int i=2; i<N; ++i) { //線性篩、初始化 s
		v[i]||(p[++cnt]=i);
		for (int j=1; j<=cnt&&(t=p[j]*i)<N; ++j) {
			v[t]=1;
			if (i%p[j]==0) break;
		}
		(s[i]=s[i-2]+s[i-1])>=P&&(s[i]-=P);
		pre[i]=1ll*pre[i-1]*s[i]%P;
	}
	I=power(pre[N-1],P-2);
	for (int i=N-1; i; --i) //初始化 inv
		inv[i]=1ll*pre[i-1]*suf%P*I%P,
		suf=1ll*suf*s[i]%P;
	for (int i=1,j; i<=cnt; ++i) // dp
		for (j=(N-1)/p[i],t=j*p[i]; j; --j,t-=p[i])
			s[t]=1ll*s[t]*inv[j]%P,
			inv[t]=1ll*inv[t]*s[j]%P;
	for (int i=2; i<N; ++i) s[i]=1ll*s[i]*s[i-1]%P;
	for (int i=2; i<N; ++i) inv[i]=1ll*inv[i]*inv[i-1]%P;
	//儲存 s 和 inv 的字首積,便於整除分塊計算
	for (scanf("%d",&T); T; --T) {
		scanf("%d%d",&n,&m),ans=1;
		for (int l=1,r,mn=min(n,m); l<=mn; l=r+1) //整除分塊
			r=min(n/(t1=n/l),m/(t2=m/l)),
			t=1ll*s[r]*inv[l-1]%P,
			ans=1ll*ans*power(t,1ll*t1*t2%(P-1))%P;
		printf("%d\n",ans);
	}
	return 0;
}

P5518 [MtOI2019]幽靈樂團 / 莫比烏斯反演基礎練習題

這個題要麻煩得多
全程大力推式子 雖然非常不基礎 但確實是針對性很強的練習題

\[\large\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)} \]

等於

\[\large\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\left(\dfrac{ij}{\gcd(i,j)\gcd(i,k)}\right)^{f(type)} \]

可以轉化為以下兩個子問題:

\[\large M(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Ci^{f(type)} \]\[\large N(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\gcd(i,j)^{f(type)} \]

答案為

\[\large \dfrac{M(A,B,C)\times M(B,A,C)}{N(A,B,C)\times N(A,C,B)} \]

下面按 \(type\) 的取值分為三種情況討論

\(\Large \textbf{1. type=0,f(type)=1}\)

\[\large M(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Ci=\prod_{i=1}^Ai^{BC}=(A!)^{BC} \]

預處理階乘,快速冪回答詢問,複雜度為 \(O(n+T\log n)\)\(n\)\(A,B,C\) 同級, \(T\) 為詢問組數,下同)

\[\large N(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\gcd(i,j)=\left(\prod_{i=1}^A\prod_{j=1}^B\gcd(i,j)\right)^C \]

括號內的式子和上題基本一致 過程不詳細寫了
\(D=\min(A,B)\) ,則

\[\large N(A,B,C)=\prod_{t=1}^D\left(\prod_{d \mid t}d^{\mu(t/d)}\right)^{(A/t)(B/t)C} \]

依然是預處理括號內的式子 然後整除分塊
時間複雜度 \(O(n\log n+T\sqrt n\log n)\)
預處理也可以做到 \(O(n\log \log n)\) ,但對本題而言基本沒啥優化效果

\(\Large \textbf{2. type=1,f(type)=i×j×k}\)

\[\large M(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Ci^{i\times j\times k}=\prod_{i=1}^Ai^{i(\sum_{j=1}^Bj)(\sum_{k=1}^Ck)} \]

\(S(n)=\sum\limits_{i=1}^ni=\dfrac{n(n+1)}2\)

\[\large M(A,B,C)=\left(\prod_{i=1}^Ai^i\right)^{S(B)S(C)} \]

括號內的式子可以 \(O(n\log n)\) 預處理,然後快速冪回答詢問

\[\large N(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\gcd(i,j)^{i\times j\times k}=\left(\prod_{i=1}^A\prod_{j=1}^B\gcd(i,j)^{i\times j}\right)^{S(C)} \]

括號裡的式子拿出來推一下

\[\large \begin{aligned}\prod_{i=1}^A\prod_{j=1}^B\gcd(i,j)^{i\times j} &=\prod_{d=1}^D\prod_{i=1}^A\prod_{j=1}^Bd^{i\times j\times [(i,j)=d]} =\prod_{d=1}^D\prod_{i=1}^{A/d}\prod_{j=1}^{B/d}d^{id\times jd\times [(i,j)=1]} =\prod_{d=1}^Dd^{d^2\sum_{i=1}^{A/d}\sum_{j=1}^{B/d}ij[(i,j)=1]}\\ &=\prod_{d=1}^Dd^{d^2\sum_{i=1}^{A/d}\sum_{j=1}^{B/d}ij\sum_{t\mid (i,j)}~\mu(t)} =\prod_{d=1}^Dd^{d^2\sum_{t=1}^{D/d}\mu(t)\sum_{i=1}^{A/dt}it\sum_{j=1}^{B/dt}jt} =\prod_{d=1}^Dd^{d^2\sum_{t=1}^{D/d}\mu(t)t^2S(A/dt)S(B/dt)}\\ &=\prod_{d=1}^D\prod_{t=1}^{D/d}d^{\mu(t)(dt)^2S(A/dt)S(B/dt)} =\prod_{t'=1}^D\prod_{d\mid t'}d^{\mu(t'/d)t'^2S(A/t')S(B/t')} =\prod_{t'=1}^D\left(\prod_{d\mid t'}d^{\mu(t'/d)}\right)^{t'^2S(A/t')S(B/t')} \end{aligned}\]

推到這裡就可以了,代回 \(N(A,B,C)\) 的表示式

\[\large N(A,B,C)=\prod_{t=1}^D\left(\prod_{d\mid t}d^{\mu(t/d)}\right)^{t^2S(A/t)S(B/t)S(C)} \]

括號裡的式子預處理過了 然後還是整除分塊
時間複雜度 \(O(n\log n+T\sqrt n\log n)\)

\(\Large\textbf{3. type=2,f(type)=gcd(i,j,k)}\)

\[\large M(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C i^{\gcd(i,j,k)}=\prod_{i=1}^Ai^{\sum_{j=1}^B\sum_{k=1}^C\gcd(i,j,k)} \]

指數可以考慮尤拉反演

\[\large \sum_{j=1}^B\sum_{k=1}^C\gcd(i,j,k) =\sum_{j=1}^B\sum_{k=1}^C\sum_{d\mid \gcd(i,j,k)}\varphi(d) =\sum_{d\mid i}\varphi(d)\sum_{j=1}^B\sum_{k=1}^C[d|j][d|k] =\sum_{d\mid i}\varphi(d)(B/d)(C/d)\]

\(E=\min(A,B,C)\)

\[\large M(A,B,C)=\prod_{i=1}^Ai^{\sum_{d\mid i}\varphi(d)(B/d)(C/d)} =\prod_{d=1}^E\left(\prod_{i'=1}^{A/d}i'd\right)^{\varphi(d)(B/d)(C/d)} =\prod_{d=1}^E\left((A/d)!\times d^{A/d}\right)^{\varphi(d)(B/d)(C/d)}\]\[\large M(A,B,C)=\prod_{d=1}^E(A/d)!^{\varphi(d)(B/d)(C/d)}\times \prod_{d=1}^Ed^{\varphi(d)(A/d)(B/d)(C/d)} \]

兩部分都可以整除分塊做

接下來是本題最難處理的式子了

\[\large N(A,B,C)=\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\gcd(i,j)^{\gcd(i,j,k)} \]

改為列舉 \(\gcd(i,j)\)

\[\large N(A,B,C)=\prod_{d=1}^Dd^{\sum_{i=1}^A\sum_{j=1}^B[(i,j)=d]\sum_{k=1}^C\gcd(d,k)} \]

指數還是常規反演

\[\large \sum_{i=1}^A\sum_{j=1}^B[(i,j)=d]=\sum_{i=1}^{A/d}\sum_{j=1}^{B/d}\sum_{t\mid (i,j)}\mu(t)=\sum_{t=1}^{D/d}\mu(t)(A/dt)(B/dt) \]\[\large \sum_{k=1}^C\gcd(d,k)=\sum_{k=1}^C\sum_{t'\mid (d,k)}\varphi(t')=\sum_{t'\mid d}\varphi(t')(C/t') \]\[\large N(A,B,C)=\prod_{d=1}^Dd^{\sum_{t=1}^{D/d}\mu(t)(A/dt)(B/dt)\sum_{t'\mid d}\varphi(t')(C/t')} \]

列舉因數看起來不對勁,換一下列舉方式

\[\large N(A,B,C)=\prod_{t'=1}^E\prod_{d'=1}^{D/t'}(d't')^{\varphi(t')(C/t')\sum_{t=1}^{D/d't'}\mu(t)(A/d't't)(B/d't't)} \]

式子有點醜 所以替換一次字母(

\[\large N(A,B,C)=\prod_{i=1}^E\prod_{j=1}^{D/i}(ij)^{\varphi(i)(C/i)\sum_{k=1}^{D/ij}\mu(k)(A/ijk)(B/ijk)} \]

似乎很難繼續推了
通過看官方題解我們瞭解到一個技巧,把 \(i,j\) 拆開分別算貢獻

\[\large N(A,B,C)=\prod_{i=1}^E\prod_{j=1}^{D/i}i^{\varphi(i)(C/i)\sum_{k=1}^{D/ij}\mu(k)(A/ijk)(B/ijk)}\times \prod_{i=1}^E\prod_{j=1}^{D/i}j^{\varphi(i)(C/i)\sum_{k=1}^{D/ij}\mu(k)(A/ijk)(B/ijk)} \]

對於第一部分,列舉 \(t=jk\)

\[\large\prod_{i=1}^Ei^{\varphi(i)(C/i)\sum_{t=1}^{D/i}~(A/it)(B/it)\sum_{k\mid t}\mu(k)} \]\[\large\prod_{i=1}^Ei^{\varphi(i)(C/i)\sum_{t=1}^{D/i}~(A/it)(B/it)[t=1]} \]

發現只有 \(t=1\) 的情況能對答案產生貢獻

\[\large\prod_{i=1}^Ei^{\varphi(i)(A/i)(B/i)(C/i)} \]

而我們在 \(M(A,B,C)\) 中也算出過相同的式子:

\[\large M(A,B,C)=\prod_{d=1}^E(A/d)!^{\varphi(d)(B/d)(C/d)}\times \prod_{d=1}^Ed^{\varphi(d)(A/d)(B/d)(C/d)} \]

所以直接約掉就好了

最後來處理第二部分

\[\large \prod_{i=1}^E\prod_{j=1}^{D/i}j^{\varphi(i)(C/i)\sum_{k=1}^{D/ij}\mu(k)(A/ijk)(B/ijk)} \]\[\large \prod_{i=1}^E\left(\prod_{j=1}^{D/i}j^{\sum_{k=1}^{D/ij}\mu(k)(A/ijk)(B/ijk)}\right)^{\varphi(i)(C/i)} \]

同樣列舉 \(t=jk\)

\[\large \prod_{i=1}^E\left(\prod_{t=1}^{D/i}\prod_{j\mid t}j^{\mu(t/j)(A/it)(B/it)}\right)^{\varphi(i)(C/i)} \]

簡單整理一下

\[\large \prod_{i=1}^E\left(\prod_{t=1}^{D/i}\left(\prod_{d\mid t}d^{\mu(t/d)}\right)^{(A/it)(B/it)}\right)^{\varphi(i)(C/i)} \]

最內層括號中的式子預處理過了
可以兩層整除分塊做,時間複雜度 \(O(n^{3/4}\log n)\)

還有個 \(O(n^{2/3}\log n)\) 的做法,雖然我沒有寫,而且也證不來複雜度,但還是提一下

之前我們算過

\[\large\prod_{i=1}^A\prod_{j=1}^B\gcd(i,j)=\prod_{t=1}^D\left(\prod_{d\mid T}d^{\mu(t/d)}\right)^{(A/t)(B/t)} \]

帶入到現在的式子裡就是

\[\large\prod_{i=1}^E\left(\prod_{j=1}^{A/i}\prod_{k=1}^{B/i}\gcd(j,k)\right)^{\varphi(i)(C/i)} \]

預處理 \(A,B\le n^{2/3}\) 的所有 \(\prod\limits_{i=1}^A\prod\limits_{j=1}^B\gcd(i,j)\)

\(i\) 需要一層整除分塊
對於括號內的部分,當 \(i\ge n^{1/3}\) 時直接用預處理的值即可,當 \(i<n^{1/3}\) 時需要再套一層整除分塊
時間複雜度證明可類比杜教篩 我都不會證

至此這道題的推式子部分終於結束了

\(O(n\log n+Tn^{3/4}\log n)\) 或者 \(O(n^{4/3}+Tn^{2/3}\log n)\) 的做法都能通過
程式碼其實寫完之後除了編譯器誰都看不懂 但還是象徵性地放一下

#include<stdio.h>
#define M0(A,B) power(fac[A],1ll*B*C%P2)
#define M1(A,B) power(pow[A],1ll*S1(B)*S1(C)%P2) 
const int N=100001; int A,B,C,D,E,T,P,P2,t,cnt,ans,p[10000];
int s[N],inv[N],S[N],INV[N],fac[N],pow[N],phi[N],v[N];
inline int min(int x,int y) { return x<y?x:y; }
inline int S1(int x) { return (1ll*x*(x+1)>>1)%P2; }
inline int power(int x,int y) {
	int ans=1;
	while (y) (y&1)&&(ans=1ll*ans*x%P),x=1ll*x*x%P,y>>=1;
	return ans;
}
void prework() {
	phi[1]=fac[0]=pow[0]=S[0]=INV[0]=1; 
	s[0]=s[1]=inv[0]=inv[1]=1;
	for (int i=2; i<N; ++i) {
		v[i]||(p[++cnt]=i,phi[i]=i-1);
		for (int j=1; j<=cnt&&(t=p[j]*i)<N; ++j) {
			v[t]=1,phi[t]=phi[i]*(p[j]-1);
			if (i%p[j]==0) { phi[t]+=phi[i]; break; }
		}
		(phi[i]+=phi[i-1])>=P2&&(phi[i]-=P2);
	}
	for (int i=2; i<N; ++i) s[i]=i,inv[i]=1ll*(P-P/i)*inv[P%i]%P;
	for (int i=1; i<=cnt; ++i)
		for (int j=(N-1)/p[i],t=j*p[i]; j; --j,t-=p[i])
			s[t]=1ll*s[t]*inv[j]%P,inv[t]=1ll*inv[t]*s[j]%P;
	for (int i=1; i<N; ++i)
		S[i]=1ll*S[i-1]*power(s[i],1ll*i*i%P2)%P,
		INV[i]=1ll*INV[i-1]*power(inv[i],1ll*i*i%P2)%P,
		s[i]=1ll*s[i-1]*s[i]%P,inv[i]=1ll*inv[i-1]*inv[i]%P,
		fac[i]=1ll*fac[i-1]*i%P,pow[i]=1ll*pow[i-1]*power(i,i)%P;
}

int N0(int A,int B,int C) {
	int D=min(A,B),ans=1;
	for (int l=1,r,t1,t2; l<=D; l=r+1)
		r=min(A/(t1=A/l),B/(t2=B/l)),
		ans=1ll*ans*power(1ll*s[r]*inv[l-1]%P,1ll*t1*t2%P2)%P;
	return power(ans,P2-C);
}
int N1(int B,int C) {
	D=min(A,B),ans=1;
	for (int l=1,r,t,t1,t2; l<=D; l=r+1)
		r=min(A/(t1=A/l),B/(t2=B/l)),
		t=1ll*S1(t1)*S1(t2)%P2,
		ans=1ll*ans*power(1ll*S[r]*INV[l-1]%P,t)%P;
	return power(ans,P2-S1(C));
}
int M2(int A,int B) {
	E=min(min(A,B),C),ans=1;
	for (int l=1,r,t,t1,t2,t3; l<=E; l=r+1)
		r=min(min(A/(t1=A/l),B/(t2=B/l)),C/(t3=C/l)),
		t=1ll*t2*t3%P2*(phi[r]-phi[l-1]+P2)%P2,
		ans=1ll*ans*power(fac[t1],t)%P;
	return ans;
}
int N2(int B,int C) {
	E=min(D=min(A,B),C),ans=1;
	for (int l=1,r,t1,t2,t3; l<=E; l=r+1)
		r=min(min(A/(t1=A/l),B/(t2=B/l)),C/(t3=C/l)),
		ans=1ll*ans*N0(t1,t2,1ll*t3*(phi[r]-phi[l-1]+P2)%P2)%P;
	return ans;
}
int main() {
	scanf("%d%d",&T,&P),P2=P-1,prework();
	while (T--)
		scanf("%d%d%d",&A,&B,&C),t=1ll*M0(A,B)*M0(B,A)%P,
		printf("%d ",1ll*t*N0(A,B,C)%P*N0(A,C,B)%P),
		printf("%d ",1ll*M1(A,B)*M1(B,A)%P*N1(B,C)%P*N1(C,B)%P),
		printf("%d\n",1ll*M2(A,B)*M2(B,A)%P*N2(B,C)%P*N2(C,B)%P);
	return 0;
}