Updated on 2020.8.6
巨幅更新，對積性函數和狄利克雷卷積部分進行重構。
新增對一類特殊求和式的講解。
Updated on 2020.9.9
添加了幾個例題。
前置知識
小碎骨
艾佛森括號 \([P] = \begin{cases} 1 &\text{If P is true}\\ 0 &\text{Otherwise} \end{cases}\)
此處 \(P\) 是一個可真可假的命題。
引理1
\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]
證明
\[\dfrac{a}{b} = \left\lfloor{\dfrac{a}{b}}\right\rfloor + r(0\le r < 1) \]
\[\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor\dfrac{a}{b}\times\dfrac{1}{c}\right\rfloor = \left\lfloor{\dfrac{1}{c}\times\left({\left\lfloor{\dfrac{a}{b} + r}\right\rfloor}\right)}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c} + \dfrac{r}{c}}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor{\dfrac{a}{b}}\right\rfloor}{c}}\right\rfloor \]
數論分塊
內容獨立了出來，詳細內容見 數論分塊 - Luckyblock
對於一類含有\(\left\lfloor\frac{n}{i}\right\rfloor\)的求和式 (\(n\) 為常數)，由於\(\left\lfloor\frac{n}{i}\right\rfloor\)單調不增，故存在多個區間\([l,r]\), 使得\(\left\lfloor\frac{n}{i}\right\rfloor = \left\lfloor\frac{n}{j}\right\rfloor(i,j\in [l,r])\) 成立。
對於任意一個\(i\)，最大的滿足上式的 \(j=\left\lfloor{\dfrac{n}{\left\lfloor{\dfrac{n}{i}}\right\rfloor}}\right\rfloor\)
積性函數
定義
若 \(\gcd(x,y) = 1\) 且 \(f(xy)=f(x)f(y)\)，則 \(f(n)\) 為積性函數。
性質
若 \(f(x)\)，\(g(x)\)均為積性函數，則以下函數也為積性函數：
\[\begin{aligned} & h(x) = f(x^p)\\ & h(x) = f^p(x)\\ & h(x) = f(x)g(x)\\ & h(x) = \sum_{d\mid x} f(d)g\left(\dfrac{x}{d}\right) \end{aligned}\]
常見積性函數
- 單位函數 \(e(n) = [n = 1]\)。
- 冪函數 \(\operatorname{Id}_{k}(n) = n^k\)， \(\operatorname{id}_1(n)\) 通常簡記為\(\operatorname{id}(n)\)。
- 常數函數 \(1(n) = 1\)。
- 因數個數 \(\operatorname{d}(n) = \sum\limits_{d\mid n} 1\)。
- 除數函數 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。
\(k=1\) 時為因數和函數，通常簡記為 \(\sigma(n)\)。
\(k=0\) 時為因數個數函數 \(\sigma_{0}(n)\)。 - 歐拉函數 \(\varphi(n) = \sum\limits_{i=1}^{n} [\gcd(i,n) = 1]\)。
- 莫比烏斯函數 \(\mu(n) = \begin{cases}1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{為}\ n\ \text{的本質不同質因子個數} \end{cases}\)
不是很懂上面寫的什么玩意？
不用深究，有個印象繼續往下看就好。
莫比烏斯函數
定義
\(\mu\) 為莫比烏斯函數，定義為
\[\mu(n) = \begin{cases} 1 &n=1\\0 &n\ \text{含有平方因子}\\(-1)^k &k\text{為}\ n\ \text{的本質不同質因子個數} \end{cases}\]
解釋
令 \(n = \prod\limits_{i=1}^{k} p_{i}^{c_i}\)，其中\(p_i\)為質因子，\(c_i\ge 1\)。
- \(n=1\)時，\(\mu (n) = 1\)。
- \(n\not ={1}\)時 ，
- \(\exist i\in [1,k], c_i > 1\) 時，\(\mu (n) = 0\)。
當某質因子出現次數大於\(1\)時，\(\mu (n) = 0\)。 - \(\forall i\in [1,k], c_i = 1\) 時，\(\mu (n) = (-1)^k\)。
當每個質因子只出現一次時，即\(n = \prod\limits_{i=1}^{k}p_i\)，\(\{p_i\}\)中元素唯一。
\(\mu (n) = (-1)^k\)，此處\(k\)為質因子的種類數。
- \(\exist i\in [1,k], c_i > 1\) 時，\(\mu (n) = 0\)。
性質
莫比烏斯函數是積性函數，且具有以下性質
\[\large \sum_{d\mid n} \mu (d) = [n=1] \]
證明，設 \(n = \prod\limits_{i=1}^{k}{p_i^{c_i}}, n' = \prod\limits_{i=1}^{k}{p_i}\)。
- 根據莫比烏斯函數定義，則有：\(\sum\limits_{d\mid n}{\mu(d)} = \sum\limits_{d\mid n'}{\mu(d)}\)。
- 若 \(n'\) 的某因子 \(d\)，有 \(\mu (d) = (-1)^i\)，則它由 \(i\) 個 本質不同的質因子組成。
由於質因子總數為 \(k\)，滿足上式的因子數為 \(C_{k}^{i}\)。 - 對於原求和式，轉為枚舉 \(\mu(d)\) 的值。
則 \(\sum\limits_{d\mid n'}{\mu(d)} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i} = \sum\limits_{i=0}^{k}{C_{k}^{i} \times (-1)^i\times 1^{k-i}}\)
根據二項式定理，上式 \(= (1+(-1))^k\)，
易知該式在 \(k=0\)，即 \(n=0\) 時為 \(1\)，否則為 \(0\)。
反演常用結論
一個反演常用結論:
\[[\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)} \]
證明 1：
設 \(n = \gcd(i,j)\)，則右\(= \sum\limits_{d\mid n} {\mu (d)} = [n = 1] = [\gcd(i,j) = 1]=\) 左。
證明 2：
暴力展開：\([\gcd(i,j) = 1] = e(\gcd(i,j)) = \sum\limits_{d\mid \gcd(i,j)}\mu(d)\)
線性篩求莫比烏斯函數
\(\mu\) 為積性函數，因此可以線性篩莫比烏斯函數。
int pnum, mu[kMaxn], p[kMaxn];
bool vis[kMaxn];
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
}
狄利克雷(Dirichlet)卷積
建議閱讀 算法學習筆記(35): 狄利克雷卷積 By: Pecco。
定義兩個數論函數 \(f,g\) 的狄利克雷卷積為
\[\large(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]
性質
- 顯然滿足 交換律，結合律，分配律。
- \(f \ast g = g \ast f\)
- \((f \ast g) \ast h = f\ast (g\ast h)\)
- \(f\ast (g+h) = f\ast g + f\ast h\)
- \(e\) 為狄利克雷卷積的單位元，有\((f\ast e)(n) = f(n)\)。
- 若 \(f,g\) 為積性函數，則 \(f\ast g\) 為積性函數。
關於單位元 \(e\)
有：
\[e = \mu \ast 1=\sum\limits_{d\mid n} \mu (d) \]
證明
\[\begin{aligned} (f\ast e)(n) = \sum_{d\mid n} f(d)e(\dfrac{n}{d}) = \sum_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] \end{aligned}\]
- 對於\([\dfrac{n}{d} = 1]\)，當且僅當 \(\dfrac{n}{d}=1\)，即 \(d=n\) 時為 \(1\)，否則為\(0\)。
- 則當 \(d=n\) 時，\(f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)。
當 \(d\not ={n}\) 時，\(f(d)\left[\dfrac{n}{d} = 1\right] = 0\)。
綜上，\((f\ast e)(n) = \sum\limits_{d\mid n} f(d)\left[\dfrac{n}{d} = 1\right] = f(n)\)，滿足單位元的性質。
則 \(e = \mu \ast 1\) 成立。
除數函數與冪函數
冪函數 \(\operatorname{Id}_{k}(n) = n^k\)。
除數函數 \(\sigma_{k}(n) = \sum\limits_{d\mid n} d^k\)。
顯然有：
\[(\operatorname{Id}_k\ast 1)(n) = \sum_{d\mid n} \operatorname{Id_k}(d) = \sum_{d\mid n} d^k = \sigma_k(n) \]
當 \(k=0\) 時，\(\operatorname{Id_0}=1\)，\(\sigma_0\) 為因數個數函數，有：
\[(1\ast 1)(n) = \sum_{d\mid n}1 = \sigma_0(n) \]
歐拉函數與恆等函數
\[\begin{aligned} \varphi \ast 1 =& \operatorname{Id}\\ \varphi =& \operatorname{Id}\ast \mu \end{aligned}\]
對於一式，當 \(n=p^m\) 時（\(p\) 為質數），有：
\[(\varphi \ast 1)(p^m) = \sum_{d\mid n}\varphi(d) = \varphi(1) +\sum_{i=1}^{m}\varphi(p^i) = 1 +\sum_{i=1}^{m}(p^i-p^{i-1}) = p^m \]
\(p^i\) 的因子有 \(p^{i-1}\) 個，為 \(1\sim p^{i-1}\)，故 \(\varphi(p^i) = p^i-p^{i-1}\)。
又 \((\varphi \ast 1)(n)\) 為積性函數，則對於任意正整數 \(n\)，有：
\[(\varphi \ast 1)(n) = (\varphi \ast 1)\left(\prod p^m\right) = \prod\left(\varphi \ast 1\right)(p^m) = \prod p^m = n \]
得證。
對於 2 式，在 1 式基礎上兩側同時 \(\ast \mu\) 即得。
左側變為 \(\varphi \ast 1\ast \mu = \varphi \ast e = \varphi\)。
計算
\[(f\ast g) (n) = \sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right) \]
考慮枚舉 \(n\) 的因子，將貢獻累加即可。
顯然可以使用埃氏篩篩出所有前綴狄利克雷卷積，復雜度 \(O(nk\log n)\)，其中 \(k\) 是計算函數值的復雜度。
莫比烏斯反演
反演是個啥？反演。
公式
設\(f(n),g(n)\)為兩個數論函數。
如果有
\[\large f(n) = (g\ast 1)(n) = \sum\limits_{d\mid n}{g(d)} \]
那么有
\[\large g(n) = (\mu \ast f)(n)=\sum\limits_{d\mid n} {\mu(d)f\left(\dfrac{n}{d}\right)} \]
證明
方法一：運用卷積。
原問題為：已知 \(f = g\ast 1\)，證明 \(g = f\ast \mu\)。
易知如下轉化：\(f\ast \mu = g\ast 1 \ast \mu \Longrightarrow f\ast \mu = g\ast e = g\)。
方法二：對原式進行數論變換。
-
用 \(\sum\limits_{d\mid n}g(d)\) 替換\(f\left(\dfrac{n}{d}\right)\)。
\[\sum_{d\mid n}{\mu(d)\sum_{k\mid \frac{n}{d}}g(k)} \] -
變換求和順序。
\[\sum_{k\mid n}g(k)\sum_{d\mid \frac{n}{k}}{\mu(d)} \] -
依據 \(\sum\limits_{d\mid n}{\mu(d)} = [n=1]\)，僅當 \(\dfrac{n}{k} = 1\) 時，\(\sum\limits_{d\mid \frac{n}{k}}{\mu(d)} = 1\)，否則為 \(0\)。
此時\(k=n\)，故上式等價於 \(\sum\limits_{k\mid n} {[n=k]\cdot g(k)} = g(n)\)。
舉例
從 狄利克雷(Dirichlet)卷積 部分可以知道一些積性函數的關系。
嘗試對它們進行反演：
\[e = \mu \ast 1\iff\mu = \mu\ast e = \mu \]
\[\sigma_k = \operatorname{Id}_k\ast 1\iff \operatorname{Id}_k = \mu \ast \sigma_k \]
\[\operatorname{Id}=\varphi \ast 1\iff \varphi = \operatorname{Id}\ast \mu \]
關於一類求和式
\[\sum_{i=1}^{n}\sum_{j=1}^{m}f(\gcd(i,j)) \]
一般套路
考慮構造出函數 \(g\)，滿足下列形式：
\[f(n) = g\ast 1 = \sum_{d\mid n}g(d) \]
則求和式變為：
\[\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}g(d) \]
考慮算術基本定理，發現若滿足 \(d\mid \gcd (i,j)\)，則 \(d\mid i\) 且 \(d\mid j\) 成立。
考慮 \(g(d)\) 在何時做出貢獻，調整枚舉順序：
\[\sum_{d=1}g(d)\sum_{i=1}^n[d\mid i]\sum_{j=1}^m[d\mid j] \]
\(\sum\limits_{i=1}^{n}[d\mid i]\) 等價於 \(1\sim n\) 中 \(d\) 的倍數的個數，則上式等價於：
\[\sum_{d=1}g(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \]
數論分塊求解即可。
例 1
\[\sum_{i=1}^{n}\sum_{j=1}^{m}\gcd(i,j) \]
發現此題的 \(f\) 等價於 \(\operatorname{Id}\)，則上式等價於：
\[\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{Id}(\gcd(i,j))\\ =& \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid \gcd(i,j)}\varphi(d)\\ =& \sum_{d=1}\varphi(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \end{aligned}\]
例 2
\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}e(\gcd(i,j))\\ =& \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum_{d\mid \gcd(i,j)}\mu (d)\\ =& \sum_{d=1}\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor \end{aligned}\]
感悟
卷點什么東西，把 \(g\) 卷出來。
\(g\) 不一定是特殊意義的函數。
例題
[HAOI2011] Problem b
\(n\) 組詢問，每次給定參數 \(a,b,c,d,k\)，求：
\[\sum\limits_{x=a}^{b}\sum\limits_{y=c}^{d}[\gcd(x,y) = k] \]\(1\le n,k,a,b,c,d\le 5\times 10^4\)，\(a\le b,c\le d\)。
令 \(f(n,m) = \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = k]\)。
根據容斥原理，則原式等價於 \(f(b,d) - f(a-1,d) - f(b,d-1) + f(a-1,d-1)\)。
\(f\) 變成了上述一類求和式的形式，考慮化簡 \(f\)。
易知原式等價於
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\gcd(i,j) = 1] \]
代入反演常用結論 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，上式化為
\[\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d\mid \gcd(i,j)}{\mu(d)} \]
變換求和順序，先枚舉\(d\mid \gcd(i,j)\)，可得
\[\sum_{d=1}\mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}[d\mid i]\sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[d\mid j] \]
對於上式后半的解釋：當\(d\mid i\) 且 \(d\mid j\) 時，\(d\mid \gcd(i,j)\)。
易知\(1\sim \left\lfloor\dfrac{n}{k}\right\rfloor\) 中 \(d\) 的倍數有 \(\left\lfloor\dfrac{\left\lfloor\dfrac{n}{k}\right\rfloor}{d}\right\rfloor = \left\lfloor\dfrac{n}{kd}\right\rfloor\) 個（由引理 1 可知），原式變為
\[\sum_{d=1}\mu(d)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \]
預處理 \(\mu\) 后，顯然可以數論分塊求解，復雜度\(O(n + T\sqrt{n})\)。
代碼
//知識點：莫比烏斯反演
/*
//By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <algorithm>
#define ll long long
const int MARX = 6e4 + 10;
//=============================================================
int N, a, b, c, d, k;
int cnt, Mobius[MARX], Prime[MARX], Sum[MARX];
bool vis[MARX];
//=============================================================
inline int read()
{
int f = 1, w = 0; char ch = getchar();
for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMobius()
{
Mobius[1] = 1;
int MAX = MARX - 10;
for(int i = 2; i <= MAX; i ++)
{
if(! vis[i]) Mobius[i] = - 1, Prime[++ cnt] = i;
for(int j = 1; j <= cnt && i * Prime[j] <= MAX; j ++)
{
vis[i * Prime[j]] = true;
if(i % Prime[j] == 0) break;
Mobius[i * Prime[j]] = - Mobius[i];
}
}
for(int i = 1; i <= MAX; i ++) Sum[i] = Sum[i - 1] + Mobius[i];
}
ll Calc(int x, int y)
{
ll ans = 0ll; int max = std :: min(x, y);
for(int l = 1, r; l <= max; l = r + 1)
r = std :: min(x / (x / l), y / (y / l)),
ans += (1ll * x / (1ll * l * k)) * (1ll * y / (1ll * l * k)) * (Sum[r] - Sum[l - 1]);
return ans;
}
ll Solve()
{
a = read(), b = read(), c = read(), d = read(), k = read();
return Calc(b, d) - Calc(b, c - 1) - Calc(a - 1, d) + Calc(a - 1, c - 1);
}
//=============================================================
int main()
{
N = read(); GetMobius();
while(N --) printf("%lld\n", Solve());
return 0;
}
[國家集訓隊]Crash的數字表格
給定 \(n,m\) , 求：
\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]\(1\le n,m\le 10^7\)。
易知原式等價於：
\[\sum_{i=1}^{n}\sum_{j=1}^{m}\dfrac{ij}{\gcd (i,j)} \]
考慮枚舉 \(\gcd(i,j)\)，設枚舉量為 \(d\)。
則 \(d=\gcd(i,j)\) 的充要條件是滿足 \(d|i, d|j\) 且 \(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\)，則原式等價於：
\[\sum_{i=1}^n\sum_{j=1}^m \sum_{d=1} \dfrac{ij}{d}[d|i][d|j][\gcd(\frac{i}{d}, \frac{j}{d})=1] \]
先枚舉 \(d\)，則原式等價於：
\[\sum_{d=1}\sum_{i=1}^{n}[d\mid i]\sum_{j=1}^m [d\mid j][\gcd(\dfrac{i}{d}, \dfrac{j}{d}=1)] \dfrac{ij}{d} \]
這個 \(d\) 很煩人，把 \(i,j\) 中的 \(d\) 提出來，變為枚舉 \(\frac{i}{d}\)，\(\frac{j}{d}\)。
消去 \(d\mid i\)，\(d\mid j\) 的限定條件，則原式等價於：
\[\begin{aligned} &\sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]\dfrac{id\times jd}{d}\\ =& \sum_{d=1}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ijd\\ =& \sum_{d=1}d \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]
單獨考慮后半部分，設 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)。
發現 \(f(x,y)\) 的形式與 [HAOI2011] Problem b 中的式子類似，代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波：
\[\begin{aligned} f(x,y) =& \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\\ =& \sum_{i=1}^{x} \sum_{j=1}^{y}\sum_{d\mid \gcd(i,j)}\mu(d) ij\\ =& \sum_{d=1}\mu(d)\sum_{i=1}^{x}[d\mid i] \sum_{j=1}^{y}[d\mid j] ij\\ =& \sum_{d=1}\mu(d) d^2\sum_{i=1}^{\left\lfloor \frac{x}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{y}{d}\right\rfloor}ij \end{aligned}\]
前半部分 \(\sum\limits_{d=1}\mu(d) d^2\)，可以考慮篩出 \(\mu(d)\) 后求前綴和。
后半部分是等差數列乘等差數列的形式，設 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)，\(g_{p,q}\) 可以通過下式 \(O(1)\) 計算：
\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]
則對於 \(f(x,y)\)，有：
\[f(x,y) = \sum_{d=1}\mu(d) d^2\cdot g(\left\lfloor \frac{x}{d}\right\rfloor, \left\lfloor\frac{y}{d}\right\rfloor) \]
數論分塊求解即可。
再看回原式，原式等價於：
\[\sum_{d=1}d\cdot f(\left\lfloor\frac{n}{d}\right\rfloor, \left\lfloor\frac{m}{d}\right\rfloor) \]
又是一個可以數論分塊求解的形式。
線性篩預處理后 數論分塊套數論分塊，復雜度 \(O(n + m)\)，瓶頸是線性篩。
一些注意的點
處理時會出現求平方的運算，需要特別注意取模問題，ll 都會爆，被坑慘了。
在預處理前綴和的這個地方：
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //僅令 mu + kMod
注意先對平方取模，在乘上 \(\mu\)，否則就會爆掉。
以及可以僅令 \(mu + mod\)。
以及這個地方：
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
平方計算，注意隨時取模。
代碼
//知識點：莫比烏斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
sum[1] = 1ll;
for (int i = 1; i <= lim_; ++ i) {
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //僅令 mu + kMod
}
}
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
int lim = std :: min(n_, m_), ret = 0;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
}
return ret;
}
int Sum(ll l, ll r) {
return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() {
int n = read(), m = read();
int lim = std :: min(n, m), ans = 0;
Euler(lim);
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
}
printf("%d", ans);
return 0;
}
/*
7718820 8445880
*/
SP5971 LCMSUM - LCM Sum
\(T\) 次詢問，每次詢問給定 \(n\)，求：
\[\sum_{i=1}^{n}\operatorname{lcm}(i,n) \]\(1<T\le 3\times 10^5\)，\(1\le n\le 10^6\)。
法一：無腦暴力
先拆 \(\operatorname{lcm}\)，原式等價於：
\[n\sum_{i=1}^{n}\dfrac{i}{\gcd(i,n)} \]
套路的枚舉 \(\gcd(i,n)\)，調換枚舉順序，原式等價於：
\[\begin{aligned} &n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd(\dfrac{i}{d},\dfrac{n}{d}) = 1]\dfrac{i}{d} \end{aligned}\]
把 \(i,n\) 中的 \(d\) 提出來，變為枚舉 \(\frac{i}{d}\)，消去整除的條件，原式等價於：
\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \]
代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等價於：
\[n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} i \sum_{t|\gcd(i,\frac{n}{d})}\mu (t) \]
值得注意的是 \(t\) 的上界為 \(\frac{n}{d}\)，\(dt\le n\)。
調換枚舉順序，先枚舉 \(t\)，原式等價於：
\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t) \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} [t|i] i \]
套路地消去整除的條件，把 \(i\) 中的 \(t\) 提出來，原式等價於：
\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t \sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor} i \]
對於最后的一個求和項，設 \(g(x) = \sum\limits_{i=1}^{x}i = \frac{x(x+1)}{2}\)，顯然可以 \(O(1)\) 求解，原式等價於：
\[n\sum_{d|n}\sum_{t|\frac{n}{d}} \mu(t)t\cdot g(\left\lfloor\dfrac{n}{dt}\right\rfloor) \]
考慮枚舉 \(T = dt\)，顯然 \(T\le n\)。
\(\mu(t)t\) 與 \(d\) 無關，可以直接考慮枚舉 \(t|T\)，原式等價於：
\[n\sum_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)\sum_{t|T}\mu(t)t \]
前半塊是一個數論分塊的形式，可以 \(O(\sqrt{n})\) 求解。
考慮后半塊，設 \(f(n)=\sum\limits_{d|n}\mu(d)d\)，發現它是一個積性函數，可以線性篩篩出，具體地：
\[f(n)= \begin{cases} 1-n &n\in \mathrm{primes} \\ f(\frac{x}{p}) &p^2\mid n\\ f(\frac{x}{p})f(p) &p^2\nmid n \end{cases} \]
其中 \(p\) 為 \(n\) 的最小質因子。
此時已經可以線性篩 + 數論分塊求解，復雜度 \(O(n+T\sqrt{n})\)，比較菜雞，時限 500ms 過不了。
考慮篩出 \(f\) 后再用埃氏篩預處理 \(\sum\limits_{T=1}^{n}g(\left\lfloor\dfrac{n}{T}\right\rfloor)f(T)\)，輸出時乘上 \(n\)，復雜度變為 \(O(n\log^2 n + n)\)。
法二：
同樣先拆 \(\operatorname{lcm}\)，枚舉 \(\gcd(i,n)\)，調換枚舉順序，原式等價於：
\[\begin{aligned} &n\sum_{i=1}^{n}\dfrac{i}{\gcd(i,n)}\\ =& n\sum_{i=1}\sum_{d|i} [\gcd(i,n) = d]\dfrac{i}{d}\\ =& n\sum_{d|n}\sum_{i=1}^{n}[d|i][\gcd({i},{n}) = d]\dfrac{i}{d} \end{aligned}\]
把 \(i,n\) 中的 \(d\) 提出來，變為枚舉 \(\frac{i}{d}\)，消去整除的條件，原式等價於：
\[\begin{aligned} &n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(id,n) = d]i\\ =& n\sum_{d|n}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,\dfrac{n}{d}) = 1]i \end{aligned}\]
調整枚舉對象，上式等價於：
\[n\sum_{d|n}\sum_{i=1}^{d}[\gcd(i,d) = 1]i \]
考慮 \(\sum\limits_{i=1}^{d}[\gcd(i,d) = 1]i\) 的實際意義，表示 \([1,d]\) 中與 \(d\) 互質的數的和。
當 \(d>1\) 時，與 \(d\) 互質的數總是成對存在，即若 \(\gcd(i,d)=1\) 成立，則 \(\gcd(d-i,d)=1\) 成立。
每對這樣的數的和為 \(d\)，共有 \(\frac{\varphi(d)}{2}\) 對這樣的數。
則原式等價於：
\[n\sum_{d|n}\dfrac{\varphi(d)d}{2} \]
可以直接預處理答案。
預處理時先線性篩出 \(\varphi\)，再埃氏篩枚舉 \(i\) 的倍數，令它們的答案加上 \(\frac{\varphi(i)i}{2}\)，最后輸出時乘上 \(n\)。
復雜度 \(O(n\log^2 n + T)\)。
法二代碼
//知識點：莫比烏斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 1e6;
//=============================================================
ll phi[kMaxn + 10], ans[kMaxn + 10];
int pnum, p[kMaxn + 10];
bool flag[kMaxn + 10];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void GetMax(int &fir, int sec) {
if (sec > fir) fir = sec;
}
void GetMin(int &fir, int sec) {
if (sec < fir) fir = sec;
}
void GetPrime() {
phi[1] = 1, flag[1] = true; //注意初始化
for (int i = 2; i <= kMaxn; ++ i) {
if (! flag[i]) {
p[++ pnum] = i;
phi[i] = i - 1ll;
}
for (int j = 1; j <= pnum && i * p[j] <= kMaxn; ++ j) {
flag[i * p[j]] = true;
if (i % p[j]) {
phi[i * p[j]] = phi[i] * phi[p[j]];
} else {
phi[i * p[j]] = phi[i] * p[j];
break;
}
}
}
for (int i = 1; i <= kMaxn; ++ i) {
for (int j = 1; i * j <= kMaxn; ++ j) {
ans[i * j] += (i == 1 ? 1 : 1ll * phi[i] * i / 2);
}
}
}
//=============================================================
int main() {
GetPrime();
int T = read();
while (T --) {
int n = read();
printf("%lld\n", 1ll * ans[n] * n);
}
return 0;
}
[SDOI2015]約數個數和
\(T\) 次詢問，每次詢問給定 \(n,m\)。
定義 >\(\operatorname{d}(i)\) 為 \(i\) 的約數個數，求：\[\sum_{i=1}^{n}\sum_{j=1}^m\operatorname{d}(ij) \]\(1<T,n\le 5\times 10^4\)。
一個結論：
\[\operatorname{d}(ij) = \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1] \]
證明
先考慮 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 的情況，有：
\[\operatorname{d}(p^{a+b})=\sum_{x|p^a}\sum_{y|p^b}[\gcd(x,y)=1] \]
對於等式左側，\(p^{a+b}\) 的約數個數為 \(a+b+1\)。
對於等式右側，在保證 \(\gcd(x,y)=1\) 成立的情況下，有貢獻的數對 \((x,y)\) 只能是下列三種形式：
- \(x>0,y-0\)，\(x\) 有 \(a\) 種取值方法。
- \(x=0,y>0\)，\(y\) 有 \(b\) 種取值方法。
- \(x=0,y=0\)。
則等式右側貢獻次數為 \(a+b+1\) 次，等於 \(p^{a+b}\) 的約數個數。
則當 \(i = p^a\)，\(j=p^b(p\in \mathrm{primes})\) 時等式成立。
又不同質因子間相互獨立，上述情況可拓展到一般的情況。
對 \(\operatorname{d}(i,j)\) 進行化簡，代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等價於：
\[\begin{aligned} \operatorname{d}(ij) =& \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\ =& \sum_{x|i}\sum_{y|j}\sum\limits_{d\mid \gcd(x,y)} {\mu (d)} \end{aligned}\]
調換枚舉順序，先枚舉 \(d\)，原式等價於：
\[\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|i}[d|x]\sum_{y|j}[d|y] \]
把各項中的 \(d\) 提出來，消去整除的條件，原式等價於：
\[\begin{aligned} &\sum_{d=1}[d|i][d|j]{\mu (d)}\sum_{x|\frac{i}{d}}\sum_{y|\frac{j}{d}}1\\ =& \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \end{aligned}\]
將 \(\operatorname{d}(ij)\) 代回原式，原式等價於：
\[\sum_{i=1}^{n}\sum_{j=1}^m \sum_{d=1}[d|i][d|j]{\mu (d)}\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]
調換枚舉順序，先枚舉 \(d\)，原式等價於：
\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{n}[d|i]\sum_{j=1}^m [d|j]\cdot \operatorname{d}(\dfrac{i}{d})\operatorname{d}(\dfrac{j}{d}) \]
把 \(i,j\) 中的 \(d\) 提出來，變為枚舉 \(\frac{i}{d}, \frac{j}{d}\)，消去整除的條件，原式等價於：
\[\sum_{d=1}{\mu (d)}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\operatorname{d}(i)\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\operatorname{d}(j) \]
考慮預處理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\)，則原式等價於：
\[\sum_{d=1}{\mu (d)}S\left({\left\lfloor\frac{n}{d}\right\rfloor}\right)S\left({\left\lfloor\frac{m}{d}\right\rfloor}\right) \]
線性篩預處理 \(\mu,\operatorname{d}\)，數論分塊求解即可，復雜度 \(O(n+T\sqrt{n})\)。
代碼
//知識點：莫比烏斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 為最小質因子的次數
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true;
mu[1] = d[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
num[i] = 1;
d[i] = 2;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
num[i * p[j]] = num[i] + 1;
d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
break;
}
mu[i * p[j]] = - mu[i];
num[i * p[j]] = 1;
d[i * p[j]] = 2ll * d[i]; //
}
}
for (int i = 1; i <= lim_; ++ i) {
summu[i] = summu[i - 1] + mu[i];
sumd[i] = sumd[i - 1] + d[i];
}
}
//=============================================================
int main() {
Euler(kMaxn - 10);
int T = read();
while (T --) {
int n = read(), m = read(), lim = std :: min(n, m);
ll ans = 0ll;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
}
printf("%lld\n", ans);
}
return 0;
}
/*
in
1
32 43
*/
/*
out
15420
*/
P3768 簡單的數學題
給定參數 \(n\)、\(p\)，求：
\[\left(\sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot \gcd(i,j)\right)\bmod p \]\(n\leq10^{10}\)，\(5\times10^8\leq p\leq1.1\times10^9\)，\(p\in \mathrm{primes}\)。
時限 4s。
無腦套路暴力。
考慮先枚舉 \(\gcd(i,j)\)，原式等價於：
\[\begin{aligned} &\sum_{d=1}d\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=d]ij\\ =& \sum_{d=1}d\sum_{i=1}^{n}[d|i]\sum_{j=1}^{n}[d|j][\gcd(\dfrac{i}{d},\dfrac{j}{d})=1]ij \end{aligned}\]
提出 \(d\)，變為枚舉 \(\frac{i}{d}\) 與 \(\frac{j}{d}\)，消去整除的條件，原式等價於：
\[\begin{aligned} &\sum_{d=1}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]id\cdot jd\\ =& \sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[\gcd(i,j)=1]ij \end{aligned}\]
代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\)，原式等價於：
\[\sum_{d=1} d^3\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}ij\sum_{t|\gcd(i,j)}\mu (t) \]
值得注意的是 \(t\) 的上界為 \(\frac{n}{d}\)，\(dt\le n\)。
調換枚舉順序，先枚舉 \(t\)，原式等價於：
\[\sum_{d=1} d^3\sum_{t=1}\mu(t)\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|i] \sum_{j=1}^{\left\lfloor\frac{n}{d}\right\rfloor}[t|j]ij \]
和上面一樣，提出 \(t\)，套路地消去整除的條件，原式等價於：
\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\sum_{i=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{dt}\right\rfloor}ij \]
發現后面兩個求和是等差數列乘等差數列的形式。
設 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\)，\(g_{p,q}\) 可以通過下式 \(O(1)\) 計算：
\[g(p,q) = \sum_{i=1}^{p} i \sum_{j=1}^{q}j= \dfrac{(1 + p)\times p}{2}\times \dfrac{(1+q)\times q}{2} \]
代入原式，原式等價於：
\[\sum_{d=1} d^3\sum_{t=1}\mu(t)t^2\cdot g\left({\left\lfloor\frac{n}{dt}\right\rfloor},{\left\lfloor\frac{n}{dt}\right\rfloor}\right) \]
考慮枚舉 \(T = dt\)，顯然 \(T\le n\)。
再考慮枚舉 \(d|T\)，即可得到 \(t = \frac{T}{d}\)，原式等價於：
\[\begin{aligned} &\sum_{T=1}^{n}g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d^3\mu{\left(\dfrac{T}{d}\right)}\left(\dfrac{T}{d}\right)^2\\ =& \sum_{T=1}^{n}T^2 g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right)\sum_{d|T}d\cdot \mu{\left(\dfrac{T}{d}\right)} \end{aligned}\]
對於后面這一坨，用 \(\sum\limits_{d|T}d\cdot \mu{\left(\frac{T}{d}\right)} = \operatorname{Id} \ast \mu(T)= \varphi(T)\) 反演，則原式等價於：
\[\sum_{T=1}^{n}T^2 \varphi(T) \cdot g\left({\left\lfloor\frac{n}{T}\right\rfloor},{\left\lfloor\frac{n}{T}\right\rfloor}\right) \]
后半塊可以數論分塊，考慮前半塊。
發現前半段即為 \(\operatorname{Id}^2(T)\varphi(T)\)，又是前綴和形式，考慮杜教篩。
有：
\[f(n) = \operatorname{Id}^2\varphi(n) \]
考慮找到一個函數 \(g\)，構造函數 \(h = f\ast g\) 使其便於求值，有：
\[h(n) = \sum_{d|n} d^2\varphi(d)\cdot g\left(\dfrac{n}{d}\right) \]
看到同時存在 \(d\) 和 \(\frac{n}{d}\)，考慮把 \(d^2\) 消去。
令 \(g = \operatorname{Id}^2\)，有：
\[\begin{aligned} h(n) =& \sum_{d|n} d^2\varphi(d)\cdot \left(\dfrac{n}{d}\right)^2\\ =& n^2\sum_{d|n} \varphi(d)\\ =& n^2 \cdot \varphi \ast 1(n)\\ \end{aligned}\]
又 \(\varphi \ast 1 = \operatorname{Id}\)，則有：
\[h(n) = n^3 \]
找到了合適的 \(g\)，套杜教篩的公式。
\[\begin{aligned} g(1)S(n) &= \sum_{i=1}^{n}h(i) - \sum_{i=2}^{n}g(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ S(n) &= \sum_{i=1}^{n} i^3 - \sum_{i=2}^{n} i^2\cdot S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right) \end{aligned}\]
前一項是自然數的立方和，有 \(\sum\limits_{i=1}^{n} i^3 = (\frac{n(n+1)}{2})^2\)。證明詳見：自然數前n項平方和、立方和公式及證明 - 百度文庫。
后一項直接等差數列求和 + 數論分塊求解即可。
「SDOI2017」數字表格
設 \(f_{i}\) 表示斐波那契數列的第 \(i\) 項。
\(T\) 組數據，每次給定參數 \(n,m\)，求：\[\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)} \pmod {10^9 + 7} \]\(1\le T\le 10^3\)，\(1\le n,m\le 10^6\)。
5S，256MB。
以下欽定 \(n\ge m\)。
大力化式子，先套路地枚舉 \(\gcd(i,j)\)，用初中知識把兩個 \(\prod\) 化到指數位置，原式等於：
\[\large\begin{aligned} &\prod_{d = 1}^{n}\prod_{i=1}^{n}\prod_{j=1}^{m}f_{d}[\gcd(i,j)=d]\\ =&\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=d]\right)} \end{aligned}\]
分母套路一波，有：
\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \end{aligned}\]
代回原式，原式等於：
\[\large \begin{aligned} &\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor\right)}\\ =& \prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor} \end{aligned}\]
考慮再暴力拆一波，原式等於：
\[\large \begin{aligned} &\prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\\ =& \prod_{d = 1}^{n}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f_{d}^{\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\right) \end{aligned}\]
做不動了，但發現變量僅有 \(k,d,kd\)，考慮更換枚舉對象改為枚舉 \(t = kd\) 與 \(d\)，則原式等於：
\[\large\prod_{t=1}^{n}\left(\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}}\right)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]
枚舉對象變成了約數形式。從后面的式子推前面的式子是比較顯然的，可以認為這種枚舉 \(t=kd\) 的形式是一種逆向操作。
設：
\[\large g(t)=\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}} \]
\(g(t)\) 可以用類似埃氏篩的方法 \(O(n\log ^2 n)\) 地預處理出來。再把 \(g\) 代回原式，原式等於：
\[\large\prod_{t=1}^{n}g(t)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]
可以考慮預處理 \(g(t)\) 的前綴積，數論分塊枚舉指數求解即可。
總時間復雜度 \(O(n\log ^2 n + T\sqrt n)\)，輕微卡常可以過。
//知識點：莫比烏斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
void Euler() {
vis[1] = true, mu[1] = 1;
for (int i = 2; i <= kN; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
}
for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
}
void Prepare() {
g[1] = g[2] = 1;
f[1] = f[2] = 1;
invf[1] = invf[2] = 1;
for (int i = 3; i <= kN; ++ i) {
g[i] = 1;
f[i] = (f[i - 1] + f[i - 2]) % mod;
invf[i] = QPow(f[i], mod - 2);
}
Euler();
for (int d = 1; d <= kN; ++ d) {
for (int j = 1; d * j <= kN; ++ j) {
if (mu[j] == 1) {
g[d * j] = g[d * j] * f[d] % mod;
} else if (mu[j] == -1) {
g[d * j] = g[d * j] * invf[d] % mod;
}
}
}
invp[0] = prod[0] = 1;
for (int i = 1; i <= kN; ++ i) {
prod[i] = prod[i - 1] * g[i] % mod;
invp[i] = QPow(prod[i], mod - 2);
}
}
//=============================================================
int main() {
Prepare();
int T = read();
while (T -- ){
n = read(), m = read(), ans = 1;
if (n < m) std::swap(n, m);
for (LL l = 1, r = 1; l <= m; l = r + 1) {
r = std::min(n / (n / l), m / (m / l));
ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
}
printf("%lld\n", ans);
}
return 0;
}
P5518 [MtOI2019]幽靈樂團 / 莫比烏斯反演基礎練習題
給定參數 \(p\)，有 \(T\) 組數據，每次給定參數 \(A,B,C\)，求：
\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(type)} \]其中 \(f(type)\) 的取值如下：
\[f(type) = \begin{cases} 1 &type = 0\\ i\times j\times k &type = 1\\ \gcd(i,j,k) &type = 2 \end{cases}\]\(1\le A,B,C\le 10^5\)，\(10^7\le p\le 1.05\times 10^9\)，\(p\in \mathbb{P}\)，\(T=70\)。
2.5S，128MB。
先化下原式，原式等於：
\[\prod_{i=1}^{A}\prod_{j=1}^{B}\prod_{k=1}^{C}\left(\dfrac{i\times j }{\gcd(i,j)\times \gcd(i,k)}\right)^{f(type)} \]
發現每一項僅與兩個變量有關，設：
\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{f(type)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{f(type)} \end{aligned}\]
發現 \(\prod\) 可以隨意交換，則原式等價於：
\[\dfrac{f_1(a,b,c)\times f_1(b,a,c)}{f_2(a,b,c)\times f_2(a,c,b)} \]
考慮在 \(type\) 取值不同時，如何快速求得 \(f_1\) 與 \(f_2\)。
一共有 \(6\) 個需要推導的式子，不妨就叫它們 \(1\sim 6\) 面了（
type = 0
\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j) \end{aligned}\]
對於 1 面，顯然有：
\[\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i = \left(\prod_{i=1}^{a}i\right)^{b\times c} \]
預處理階乘 + 快速冪即可，單次計算時間復雜度 \(O(\log n)\)。
再考慮 2 面，套路地枚舉 \(\gcd\)，顯然有：
\[\large \begin{aligned} &\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\\ =& \left(\prod_{d=1} d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\right)}\right)^c \end{aligned}\]
指數是個套路，可以看這里：P3455 [POI2007]ZAP-Queries。於是有：
\[\begin{aligned} &\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{a}{kd}\right\rfloor\left\lfloor\dfrac{b}{kd}\right\rfloor \end{aligned}\]
代回原式，略做處理，則原式等於：
\[\large \begin{aligned} &\left(\prod_{d=1} d^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c\\ =& \left(\prod_{d=1} \left(d^{\sum\limits_{k=1}\mu(k)}\right)^{\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor}\right)^c\\ =& \prod_{d=1} \left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(\mu(k)\left\lfloor\frac{a}{kd}\right\rfloor\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)^c \end{aligned}\]
像「SDOI2017」數字表格 一樣，考慮枚舉 \(t=kd\) 和 \(d\)，則原式等於：
\[\large \prod_{t=1}^{n}\left(\left(\prod_{d|t} d^{\mu{\left(\frac{t}{d}\right)}}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]
設：
\[\large g_0(t) = \prod_{d|t}d^{\mu\left(\frac{t}{d}\right)} \]
線性篩預處理 \(\mu\) 后，\(g_0(t)\) 可以用埃氏篩預處理，時間復雜度 \(O(n\log n)\)。再代回原式，原式等於：
\[\large \prod_{t=1}^{a}\left(g_0(t)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\right)^c \]
預處理 \(g_0(t)\) 的前綴積和前綴積的逆元，復雜度 \(O(n\log n)\)。
數論分塊 + 快速冪計算即可，單次時間復雜度 \(O(\sqrt n\log n)\)。
type = 1
\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{i\times j\times k} \end{aligned}\]
考慮 3 面，把 \(\prod k\) 扔到指數位置，有：
\[\large \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{i\times j\times k} = \prod_{i=1}^{a}\prod_{j=1}^{b}i^{\left(i\times j\times \sum\limits_{k = 1}^{c} k\right)} \]
再把 \(\prod j\) 也扔到指數位置，引入 \(\operatorname{sum}(n) = \sum_{i=1}^{n} i = \frac{n(n+1)}{2}\)，原式等於：
\[\left(\prod_{i=1}^{a}i^i\right)^{\operatorname{sum}(b)\times \operatorname{sum}(c)} \]
預處理 \(i^i\) 的前綴積，復雜度 \(O(n\log n)\)。
指數可以 \(O(1)\) 算出，然后快速冪，單次時間復雜度 \(O(\log n)\)。
根據費馬小定理，指數需要對 \(p - 1\) 取模。注意 \(p-1\) 不是質數，計算 \(\operatorname{sum}\) 時不能用逆元，但乘不爆 LL，直接算就行。
再考慮 4 面，發現 \(k\) 與 \(\gcd\) 無關，則同樣把 \(\prod k\) 扔到指數位置，則有：
\[\large \begin{aligned} &\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)^{\operatorname{sum}(c)} \end{aligned}\]
套路地枚舉 \(\gcd\)，原式等於：
\[\large \left(\prod_{d=1}d^{\left(\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\right)}\right)^{\operatorname{sum}(c)} \]
大力化簡指數，有：
\[\large \begin{aligned} &\sum\limits_{i=1}^a \sum\limits_{j=1}^b i\times j[\gcd(i,j)=d]\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} i\times j[\gcd(i,j)=1\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{t|\gcd(i,j)}\mu(t)\\ =& d^2 \sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j\sum\limits_{k|\gcd(i,j)}\mu(k)\\ =& d^2 \sum\limits_{k=1}\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor} i[k|i] \sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor} j[k|j]\\ =& d^2 \sum\limits_{k=1}k^2\mu(k)\sum\limits_{i=1}^{\left\lfloor\frac{a}{kd}\right\rfloor} i\sum\limits_{j=1}^{\left\lfloor\frac{b}{kd}\right\rfloor} j\\ =& d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\\ \end{aligned}\]
指數化不動了，代回原式，原式等於：
\[\large \left(\prod_{d=1}d^{\left(d^2\sum\limits_{k=1}k^2\mu(k)\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)} \]
同 2 面的情況，先展開一下，再枚舉 \(t=kd\) 和 \(d\)，原式等於：
\[\large \begin{aligned} &\left(\prod_{d=1}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}d^{\left(d^2 k^2\mu(k)\right)}\right)^{\left(\operatorname{sum}{\left(\left\lfloor\frac{a}{kd}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{kd}\right\rfloor\right)}\right)}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(d^2\left(\frac{t}{d}\right)^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)}\\ =& \prod_{t=1}\left(\left(\prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)}\right)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \end{aligned}\]
與二面相同，設：
\[\large g_1(t) = \prod_{d|t}d^{\left(t^2\mu\left(\frac{t}{d}\right)\right)} \]
\(g_1(t)\) 可以用埃氏篩套快速冪篩出，時間復雜度 \(O(n\log^2 n)\)。再代回原式，原式等於：
\[\prod_{t=1}\left(g_1(t)^{\operatorname{sum}{\left(\left\lfloor\frac{a}{t}\right\rfloor\right)} \operatorname{sum}{\left(\left\lfloor\frac{b}{t}\right\rfloor\right)}}\right)^{\operatorname{sum}(c)} \]
同樣預處理 \(g_1(t)\) 的前綴積及其逆元，時間復雜度 \(O(n\log n)\)。
整除分塊 + 快速冪即可，單次時間復雜度 \(O(\sqrt n\log n)\)。
注意指數的取模。
type = 2
\[\begin{aligned} f_1(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ f_2(a,b,c) &= \prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)} \end{aligned}\]
考慮 5 面，手段同上，大力反演化簡一波，再調換枚舉對象，則有：
\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} i^{\gcd(i,j,k)}\\ =&\prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{b}\sum\limits_{k=1}^{c}[\gcd(i,j,k)=d]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{a}i^{\left(\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(\frac{i}{d},j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[\gcd(i,j,k)=1]\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}\sum\limits_{x|\gcd(i,j,k)}{\mu(x)}\right)}\\ =& \prod_{d=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\sum\limits_{x=1}\mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}(id)^{\left(d\times \mu(x)[x|i]\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}[x|j]\sum\limits_{k=1}^{\left\lfloor\frac{c}{d}\right\rfloor}[x|k]\right)}\\ =& \prod_{d=1}\prod_{x=1}\prod\limits_{i=1}^{\left\lfloor\frac{a}{xd}\right\rfloor}(ixd)^{\left(d\times \mu(x){\left\lfloor\frac{b}{xd}\right\rfloor}{\left\lfloor\frac{c}{xd}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}(it)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\prod_{d|T}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\prod_{i=1}^{\left\lfloor\frac{a}{t}\right\rfloor}i\right)^{d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ \end{aligned}\]
引入 \(\operatorname{fac}(n) = \prod_{i=1}^{n} i\)，再根據枚舉對象調整一下指數，原式等於：
\[\large \begin{aligned} &\prod_{t = 1}\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\left(d\times \mu\left(\frac{t}{d}\right){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}\right)}\\ =& \prod_{t = 1}\left(\prod_{d|t}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\sum\limits_{d|t}d\times \mu\left(\frac{t}{d}\right)}\right)^{{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}} \end{aligned}\]
指數中出現了一個經典的狄利克雷卷積的形式，對其進行反演。
將 \((\operatorname{Id}\ast \mu) (n)= \varphi (n)\) 代入原式，原式等於：
\[\large \begin{aligned} &\prod_{t = 1}\left(t^{\left\lfloor\frac{a}{t}\right\rfloor}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\\ =& \prod_{t = 1}\left(t^{\varphi(t)\left\lfloor\frac{a}{t}\right\rfloor{\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\times \operatorname{fac}\left(\left\lfloor\frac{a}{t}\right\rfloor\right)^{\varphi(t){\left\lfloor\frac{b}{t}\right\rfloor}{\left\lfloor\frac{c}{t}\right\rfloor}}\right) \end{aligned}\]
預處理 \(t^{\varphi(t)}\) 的前綴積及逆元，階乘的前綴積及階乘逆元，\(\pmod {p-1}\) 下的 \(\varphi(t)\) 的前綴和（指數
），時間復雜度 \(O(n\log n)\)。
同樣整除分塊 + 快速冪即可，單次時間復雜度 \(O(\sqrt n\log n)\)。
然后是最掉 sans 的 6 面。有 \(\gcd(i,j,k) = \gcd(\gcd(i,j), k)\)，考慮先枚舉 \(\gcd(i,j)\)，然后套路化式子，則有：
\[\large \begin{aligned} &\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} \gcd(i,j)^{\gcd(i,j,k)}\\ =& \prod_{d=1}\prod_{i=1}^{a}\prod_{j=1}^{b}\prod_{k=1}^{c} [\gcd(i,j)=d] d^{\gcd(d,k)}\\ =& \prod_{d=1} \left(d^{\left(\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b} [\gcd(i,j)=d]\right)}\right)^{\sum\limits_{k=1}^{c}\gcd(d,k)} \end{aligned}\]
先考慮最外面的指數，這也是個套路，可以參考 一個例子。用 \(\operatorname{Id} = \varphi \ast 1\) 反演，顯然有：
\[\large \begin{aligned} &\sum\limits_{k=1}^{c}\gcd(d,k)\\ =& \sum\limits_{k=1}^{c}\sum_{x|\gcd(d,k)}\varphi(x)\\ =& \sum_{x=1}\varphi(x)[x|d]\sum_{k=1}^{c}[x|k]\\ =& \sum_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor \end{aligned}\]
再考慮里面的指數，發現這式子在 2 面已經推了一遍了，於是直接拿過來用，有：
\[\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}[\gcd(i,j) = d]=\sum_{y=1}\mu(y)\left\lfloor\dfrac{a}{yd}\right\rfloor\left\lfloor\dfrac{b}{yd}\right\rfloor \]
將化簡后的兩個指數代入原式，原式等於：
\[\large \begin{aligned} &\prod_{d=1} \left(d^{\left(\sum\limits_{y=1}\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]
與 2、4 面同樣套路地，考慮枚舉 \(t=yd\) 和 \(d\)，再略作調整，原式等於：
\[\large \begin{aligned} &\prod_{d=1} \left(\prod\limits_{y=1}d^{\left(\mu(y)\left\lfloor\frac{a}{yd}\right\rfloor\left\lfloor\frac{b}{yd}\right\rfloor\right)}\right)^{\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor}\\ =& \prod_{t=1}\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\\ =& \prod_{t=1}\left(\prod_{d|t} d^{\left(\mu\left(\frac{t}{d}\right)\sum\limits_{x|d}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{t=1}\left(\prod_{d|t} \prod_{x|d}d^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]
發現要同時枚舉 \(d\) 和 \(x\)，化不動了。
從題解里學到一個比較神的技巧，考慮把 \(d\) 拆成 \(x\) 和 \(\frac{d}{x}\) 分別計算貢獻再相乘，即分別計算下兩式：
\[\large \begin{aligned} &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}x^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ &\prod_{t=1}\left(\prod_{d|t} \prod_{x|d}{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor} \end{aligned}\]
先考慮 \(x\) 的情況，首先把枚舉 \(x\) 調整到最外層。設 \(\operatorname{lim}=\max(a,b,c)\)，則原式等於：
\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{x}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t} {x}^{\left(\mu\left(\frac{tx}{dx}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\prod_{d|t} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]
把 \(\prod {t}\) 挪到指數位置，原式等於：
\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\sum\limits_{d|t}\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\mu\left(\frac{t}{d}\right)\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor\sum\limits_{d|t}\mu\left(\frac{t}{d}\right)\right)} \end{aligned}\]
指數中又出現了一個經典的狄利克雷卷積的形式，對其進行反演。
將 \((\mu \ast 1) (n)= \epsilon (n)=[n=1]\) 代入原式，原式等於：
\[\large \begin{aligned} &\prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\sum\limits_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor[t=1]\right)}\\ =& \prod_{x=1} {x}^{\left(\varphi(x)\left\lfloor\frac{a}{x}\right\rfloor\left\lfloor\frac{b}{x}\right\rfloor\left\lfloor\frac{c}{x}\right\rfloor\right)} \end{aligned}\]
得到了一個非常優美的式子，而且發現這個式子是 5 面最終答案的一部分。同 5 面的做法，直接整除分塊即可。
再考慮 \(\frac{d}{x}\) 的情況，同上先把枚舉 \(x\) 放到最外層，並調整一下指數，則原式等於：
\[\large \begin{aligned} &\prod_{x=1} \prod_{t=1}^{\operatorname{lim}}[x|t]\left(\prod_{d|t} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{t}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{t}\right\rfloor\left\lfloor\frac{b}{t}\right\rfloor}\\ =& \prod_{x=1} \prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\left(\mu\left(\frac{tx}{d}\right)\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\\ =& \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|tx} [x|d]{\left(\frac{d}{x}\right)}^{\mu\left(\frac{tx}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \end{aligned}\]
考慮枚舉 \(dx\)，替換原來的 \(d\)，注意一下這里的倍數關系。原式等於：
\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}\left(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\right)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]
發現最內層的式子 \(\prod_{d|t}d^{\mu\left(\frac{t}{d}\right)}\)，即為二面處理過的 \(g_0(t)\)。直接代入，原式等於：
\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{a}{tx}\right\rfloor\left\lfloor\frac{b}{tx}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]
一個小結論，證明可以看 這里：
\[\forall a,b,c\in \mathbb{Z},\left\lfloor\dfrac{a}{bc}\right\rfloor = \left\lfloor{\dfrac{\left\lfloor\dfrac{a}{b}\right\rfloor}{c}}\right\rfloor \]
則原式等於：
\[\large \prod_{x=1} \left(\prod_{t=1}^{\left\lfloor\frac{\operatorname{lim}}{x}\right\rfloor}g_0(t)^{\left\lfloor\frac{\left\lfloor\frac{a}{x}\right\rfloor}{t}\right\rfloor\left\lfloor\frac{\left\lfloor\frac{b}{x}\right\rfloor}{t}\right\rfloor}\right)^{\varphi(x)\left\lfloor\frac{c}{x}\right\rfloor} \]
於是可以先對外層整除分塊，再對內層整除分塊。
然后就做完了，哈哈哈。
一些實現上的小技巧：
- 逆元能預處理就預處理。
- 注意對指數取模，模數為 \(p-1\)。
//知識點：莫比烏斯反演
/*
By:Luckyblock
用了比較清晰易懂的變量名，閱讀體驗應該會比較好。
vsc 的自動補全真是太好啦！
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
using std::min;
using std::max;
#define LL long long
const int Lim = 1e5;
const int kN = 1e5 + 10;
//=============================================================
LL A, B, C, mod, ans;
int T, p_num, p[kN];
bool vis[kN];
LL mu[kN], phi[kN], fac[kN], g[2][kN];
LL sumphi[kN], prodt_phi[kN], prodi_i[kN], prodg[2][kN];
LL inv[kN], inv_fac[kN], inv_prodt_phi[kN], inv_prodg[2][kN];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) {
w = (w << 3) + (w << 1) + (ch ^ '0');
}
return f * w;
}
void Chkmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
x_ %= mod;
y_ %= mod - 1;
LL ret = 1;
for (; y_; y_ >>= 1ll) {
if (y_ & 1) ret = ret * x_ % mod;
x_ = x_ * x_ % mod;
}
return ret;
}
LL Inv(LL x_) {
return QPow(x_, mod - 2);
}
LL Sum(LL n_) {
return ((n_ * (n_ + 1ll)) / 2ll) % (mod - 1);
}
void Euler() {
vis[1] = true, mu[1] = phi[1] = 1; //初值
for (int i = 2; i <= Lim; ++ i) {
if (! vis[i]) {
p[++ p_num] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for (int j = 1; j <= p_num && i * p[j] <= Lim; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
phi[i * p[j]] = phi[i] * p[j];
break;
}
mu[i * p[j]] = -mu[i];
phi[i * p[j]] = phi[i] * (p[j] - 1);
}
}
}
void Prepare() {
Euler();
inv[1] = fac[0] = prodt_phi[0] = prodi_i[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
g[0][i] = g[1][i] = 1;
fac[i] = 1ll * fac[i - 1] * i % mod;
sumphi[i] = (sumphi[i - 1] + phi[i]) % (mod - 1);
prodi_i[i] = prodi_i[i - 1] * QPow(i, i) % mod;
if (i > 1) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
prodt_phi[i] = prodt_phi[i - 1] * QPow(i, phi[i]) % mod;
inv_prodt_phi[i] = Inv(prodt_phi[i]);
}
for (int d = 1; d <= Lim; ++ d) {
for (int j = 1; d * j <= Lim; ++ j) {
int t = d * j;
if (mu[j] == 1) {
g[0][t] = g[0][t] * d % mod;
g[1][t] = g[1][t] * QPow(1ll * d, 1ll * t * t) % mod;
} else if (mu[j] == -1) {
g[0][t] = g[0][t] * inv[d] % mod;
g[1][t] = g[1][t] * Inv(QPow(1ll * d, 1ll * t * t)) % mod;
}
}
}
inv_prodg[0][0] = prodg[0][0] = 1;
inv_prodg[1][0] = prodg[1][0] = 1;
inv_prodt_phi[0] = 1;
for (int i = 1; i <= Lim; ++ i) {
for (int j = 0; j <= 1; ++ j) {
prodg[j][i] = prodg[j][i - 1] * g[j][i] % mod;
inv_prodg[j][i] = Inv(prodg[j][i]);
}
}
}
LL f1(LL a_, LL b_, LL c_, int type) {
if (! type) return QPow(fac[a_], b_ * c_);
if (type == 1) return QPow(prodi_i[a_], Sum(b_) * Sum(c_));
LL ret = 1, lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(fac[a_ / l], (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
return ret;
}
LL f2_2(LL a_, LL b_) {
LL ret = 1;
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
ret = ret * QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l)) % mod;
}
return ret;
}
LL f2(LL a_, LL b_, LL c_, int type) {
LL ret = 1;
if (! type) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[0][r] * inv_prodg[0][l - 1], 1ll * (a_ / l) * (b_ / l));
ret = (ret * QPow(val, c_)) % mod;
}
} else if (type == 1) {
for (LL l = 1, r = 1; l <= min(a_, b_); l = r + 1) {
r = min(a_ / (a_ / l), b_ / (b_ / l));
LL val = QPow(prodg[1][r] * inv_prodg[1][l - 1], Sum(a_ / l) * Sum(b_ / l));
ret = (ret * QPow(val, Sum(c_))) % mod;
}
} else {
LL lim = min(min(a_, b_), c_);
for (LL l = 1, r = 1; l <= lim; l = r + 1) {
r = min(min(a_ / (a_ / l), b_ / (b_ / l)), c_ / (c_ / l));
ret = ret * QPow(f2_2(a_ / l, b_ / l), (sumphi[r] - sumphi[l - 1] + mod - 1) % (mod - 1) * (c_ / l)) % mod;
ret = ret * QPow(prodt_phi[r] * inv_prodt_phi[l - 1], (a_ / l) * (b_ / l) % (mod - 1) * (c_ / l)) % mod;
}
}
return ret;
}
//=============================================================
int main() {
T = read(), mod = read();
Prepare();
while (T -- ) {
A = read(), B = read(), C = read();
for (int i = 0; i <= 2; ++ i) {
ans = f1(A, B, C, i) * f1(B, A, C, i) % mod;
ans = ans * Inv(f2(A, B, C, i)) % mod * Inv(f2(A, C, B, i)) % mod;
printf("%lld ", ans);
}
printf("\n");
}
return 0;
}
寫在最后
參考資料：
Oi-Wiki-莫比烏斯反演
算法學習筆記(35): 狄利克雷卷積 By: Pecco
題解 SP5971 【LCMSUM - LCM Sum】 - BJpers2 的博客
題解 SP5971 【LCMSUM - LCM Sum】 - Venus 的博客
題解 P3327 【[SDOI2015]約數個數和】 - suncongbo 的博客
把 Oi-Wiki 上的內容進行了復制 整理擴充。
我是個沒有感情的復讀機(大霧)