【51nod 1986 Jason曾不想做的数论题】【数论】

    技术2022-07-10  92

    题意

    给定 n , m n, m n,m,求 ∏ X ∈ [ 1 , m ] n l c m ( X 1 , ⋯   , X n ) gcd ⁡ ( X 1 , ⋯   , X n ) \prod_{X\in [1,m]^n}\mathrm{lcm}(X_1,\cdots,X_n)^{\gcd(X_1,\cdots,X_n)} X[1,m]nlcm(X1,,Xn)gcd(X1,,Xn)

    X ∈ [ 1 , m ] n X\in [1,m]^n X[1,m]n表示 X X X取遍所有长度为 n n n的序列,满足 X i X_i Xi 1 1 1 m m m之间的整数。 n ≤ 1 0 9 , m ≤ 1 0 8 n\le 10^9,m\le 10^8 n109,m108

    分析

    先考虑这个问题的弱化版:求

    S ( m ) = ∏ X ∈ [ 1 , m ] n l c m ( X 1 , ⋯   , X n ) S(m)=\prod_{X\in [1,m]^n}\mathrm{lcm}(X_1,\cdots,X_n) S(m)=X[1,m]nlcm(X1,,Xn)

    第一种算法是

    S ( m ) = ∏ p ∈ P , p k ≤ m p m n − ( m − ⌊ m p k ⌋ ) n S(m)=\prod_{p\in \mathbb{P},p^k\le m}p^{m^n-(m-\lfloor\frac{m}{p^k}\rfloor)^n} S(m)=pP,pkmpmn(mpkm)n

    可以预处理

    f ( n ) = { p n = p k , p ∈ P 1 o t h e r w i s e f(n)=\begin{cases} p & n = p^k,p\in \mathbb{P}\\ 1 & otherwise \end{cases} f(n)={p1n=pk,pPotherwise

    的前缀积,然后通过分块 O ( m log ⁡ P ) O(\sqrt m\log P) O(m logP)求出答案。

    第二种算法是通过min-max反演得到

    l c m ( X 1 , ⋯   , X n ) = ∏ k = 1 n ∏ 1 ≤ y 1 < ⋯ < y k ≤ n gcd ⁡ ( X y 1 , ⋯   , X y k ) ( − 1 ) k − 1 \mathrm{lcm}(X_1,\cdots,X_n)=\prod_{k=1}^n\prod_{1\le y_1<\cdots <y_k\le n}\gcd(X_{y_1},\cdots,X_{y_k})^{(-1)^{k-1}} lcm(X1,,Xn)=k=1n1y1<<ykngcd(Xy1,,Xyk)(1)k1

    如果令

    g ( p ) = ∑ k = 1 n ∑ 1 ≤ y 1 < ⋯ < y k ≤ n , 1 ≤ X y i ≤ p , 1 ≤ X i ≤ m [ gcd ⁡ ( X y 1 , ⋯   , X y k ) = 1 ] ( − 1 ) k − 1 g(p)=\sum_{k=1}^n\sum_{1\le y_1<\cdots <y_k\le n,1\le X_{y_i}\le p,1\le X_i\le m}[\gcd(X_{y_1},\cdots,X_{y_k})=1](-1)^{k-1} g(p)=k=1n1y1<<ykn,1Xyip,1Xim[gcd(Xy1,,Xyk)=1](1)k1

    那么

    S ( m ) = ∏ k = 1 m d g ( ⌊ m d ⌋ ) S(m)=\prod_{k=1}^md^{g(\lfloor\frac{m}{d}\rfloor)} S(m)=k=1mdg(dm)

    注意到 g g g只有 O ( m ) O(\sqrt m) O(m )种取值,只要求出这些 g g g,就能够通过分块 O ( m log ⁡ P ) O(\sqrt m\log P) O(m logP)求出答案。问题在于如何求 g g g。令

    g t ( p ) = ∑ k = 1 n ∑ 1 ≤ y 1 < ⋯ < y k ≤ n , 1 ≤ X y i ≤ p , 1 ≤ X i ≤ m [ gcd ⁡ ( X y 1 , ⋯   , X y k ) = t ] ( − 1 ) k − 1 g_t(p)=\sum_{k=1}^n\sum_{1\le y_1<\cdots <y_k\le n,1\le X_{y_i}\le p,1\le X_i\le m}[\gcd(X_{y_1},\cdots,X_{y_k})=t](-1)^{k-1} gt(p)=k=1n1y1<<ykn,1Xyip,1Xim[gcd(Xy1,,Xyk)=t](1)k1

    A ( p ) = ∑ t = 1 p g t ( p ) = ∑ k = 1 n ( n k ) m n − k p k ( − 1 ) k − 1 A(p)=\sum_{t=1}^pg_t(p)=\sum_{k=1}^n\binom{n}{k}m^{n-k}p^k(-1)^{k-1} A(p)=t=1pgt(p)=k=1n(kn)mnkpk(1)k1

    = m n − ∑ k = 0 n ( n k ) m n − k ( − p ) k = m n − ( m − p ) n =m^n-\sum_{k=0}^n\binom{n}{k}m^{n-k}(-p)^k=m^n-(m-p)^n =mnk=0n(kn)mnk(p)k=mn(mp)n

    又注意到 g t ( p ) = g ( ⌊ p t ⌋ ) g_t(p)=g(\lfloor\frac{p}{t}\rfloor) gt(p)=g(tp),因此

    g ( p ) = g 1 ( p ) = m n − ( m − p ) n − ∑ t = 2 p g ( ⌊ p t ⌋ ) g(p)=g_1(p)=m^n-(m-p)^n-\sum_{t=2}^pg(\lfloor\frac{p}{t}\rfloor) g(p)=g1(p)=mn(mp)nt=2pg(tp)

    用上式计算 g g g,复杂度和杜教筛类似,都是 O ( m 3 4 ) O(m^{\frac{3}{4}}) O(m43)

    回到原问题,类似可设

    G ( p ) = ∏ X ∈ [ 1 , p ] n l c m ( X 1 , ⋯   , X n ) [ gcd ⁡ ( X 1 , ⋯   , X n ) = 1 ] G(p)=\prod_{X\in [1,p]^n}\mathrm{lcm}(X_1,\cdots,X_n)^{[\gcd(X_1,\cdots,X_n)=1]} G(p)=X[1,p]nlcm(X1,,Xn)[gcd(X1,,Xn)=1]

    G 1 ( p ) = ∑ X ∈ [ 1 , p ] n [ gcd ⁡ ( X 1 , ⋯   , X n ) = 1 ] G1(p)=\sum_{X\in [1,p]^n}[\gcd(X_1,\cdots,X_n)=1] G1(p)=X[1,p]n[gcd(X1,,Xn)=1]

    那么

    a n s = ∏ d = 1 m G ( ⌊ m d ⌋ ) d ∗ ( d d ) G 1 ( ⌊ m d ⌋ ) ans=\prod_{d=1}^m G(\lfloor\frac{m}{d}\rfloor)^d*(d^d)^{G1(\lfloor\frac{m}{d}\rfloor)} ans=d=1mG(dm)d(dd)G1(dm)

    如果求出了 G G G G 1 G1 G1以及 d d d^d dd的前缀积,就可以分块来求 a n s ans ans了。求 d d d^d dd的前缀积可以通过

    ∏ i = 1 n i i = ∏ i = 1 n ( i ! ( i − 1 ) ! ) i = ( n ! ) n ∏ i = 1 n − 1 i ! \prod_{i=1}^ni^i=\prod_{i=1}^n\Big(\frac{i!}{(i-1)!}\Big)^i=\frac{(n!)^n}{\prod_{i=1}^{n-1}i!} i=1nii=i=1n((i1)!i!)i=i=1n1i!(n!)n

    因此只需要 O ( m ) O(m) O(m)预处理阶乘的前缀积。

    用求 g g g相同的方式求求 G G G G 1 G1 G1。令

    G 1 t ( p ) = ∑ X ∈ [ 1 , p ] n [ gcd ⁡ ( X 1 , ⋯   , X n ) = t ] G1_t(p)=\sum_{X\in [1,p]^n}[\gcd(X_1,\cdots,X_n)=t] G1t(p)=X[1,p]n[gcd(X1,,Xn)=t]

    那么

    ∑ t = 1 p G 1 t ( p ) = p n ⇒ G 1 ( p ) = p n − ∑ t = 2 p G 1 ( ⌊ p t ⌋ ) \sum_{t=1}^pG1_t(p)=p^n\Rightarrow G1(p)=p^n-\sum_{t=2}^pG1(\lfloor\frac{p}{t}\rfloor) t=1pG1t(p)=pnG1(p)=pnt=2pG1(tp)

    G t ( p ) = ∏ X ∈ [ 1 , p ] n l c m ( X 1 , ⋯   , X n ) gcd ⁡ ( X 1 , ⋯   , X n ) = t G_t(p)=\prod_{X\in [1,p]^n}\mathrm{lcm}(X_1,\cdots,X_n)^{\gcd(X_1,\cdots,X_n)=t} Gt(p)=X[1,p]nlcm(X1,,Xn)gcd(X1,,Xn)=t

    那么

    ∏ t = 1 p G t ( p ) = ∏ X ∈ [ 1 , p ] n l c m ( X 1 , ⋯   , X n ) = S ( p ) \prod_{t=1}^pG_t(p)=\prod_{X\in [1,p]^n}\mathrm{lcm}(X_1,\cdots,X_n)=S(p) t=1pGt(p)=X[1,p]nlcm(X1,,Xn)=S(p)

    得到

    G ( p ) = S ( p ) ∗ ( ∏ t = 2 p G ( ⌊ p t ⌋ ) ∗ t G 1 ( ⌊ p t ⌋ ) ) − 1 G(p)=S(p)*\Big(\prod_{t=2}^pG(\lfloor\frac{p}{t}\rfloor)*t^{G1(\lfloor\frac{p}{t}\rfloor)}\Big)^{-1} G(p)=S(p)(t=2pG(tp)tG1(tp))1

    S ( p ) S(p) S(p)的时候,当 p p p较小时用第一种算法,较大时用第二种算法。

    代码

    #include<bits/stdc++.h> using namespace std; typedef long long LL; const int MOD = 1000000007; const int MOD1 = MOD - 1; const int N = 1e8; const int M = 5e6; const int K = 2e4 + 5; int n, m, jc1[N + 5], jc2[N + 5], po[M + 5], tot, prime[M / 10], inv[M + 5], pro[M + 5], pro_inv[M + 5]; int w[K], w1[K], g[K], G[K], G1[K], sqrtm, h1[K], h2[K]; bool not_prime[M + 5]; int ksm(int x, int y, int mo) { int ans = 1; while (y) { if (y & 1) ans = (LL)ans * x % mo; x = (LL)x * x % mo; y >>= 1; } return ans; } void pre() { po[1] = 1; for (int i = 2; i <= M; i++) { if (!not_prime[i]) po[i] = ksm(i, n, MOD1); for (int j = 1; j <= tot && i * prime[j] <= M; j++) { po[i * prime[j]] = (LL)po[i] * po[prime[j]] % MOD1; if (i % prime[j] == 0) break; } } } void init() { jc1[0] = jc2[0] = inv[0] = inv[1] = 1; for (int i = 1; i <= N; i++) jc1[i] = (LL)jc1[i - 1] * i % MOD, jc2[i] = (LL)jc2[i - 1] * jc1[i] % MOD;; for (int i = 2; i <= M; i++) inv[i] = (LL)(MOD - MOD / i) * inv[MOD % i] % MOD; for (int i = 0; i <= M; i++) pro[i] = pro_inv[i] = 1; for (int i = 2; i <= M; i++) { if (!not_prime[i]) { prime[++tot] = i; for (int j = 1; j <= M / i; j *= i) pro[i * j] = i, pro_inv[i * j] = inv[i]; } pro[i] = (LL)pro[i - 1] * pro[i] % MOD; pro_inv[i] = (LL)pro_inv[i - 1] * pro_inv[i] % MOD; for (int j = 1; j <= tot && i * prime[j] <= M; j++) { not_prime[i * prime[j]] = 1; if (i % prime[j] == 0) break; } } } void pre_calc_g(int m) { int num = 0; for (int i = 1; i <= m; i = m / (m / i) + 1) w1[++num] = m / i; int t = ksm(m, n, MOD1); for (int i = num; i >= 1; i--) { int p = w1[i], id = p <= sqrtm ? h1[p] : h2[m / p]; g[id] = (t + MOD1 - ksm(m - p, n, MOD1)) % MOD1; for (int j = 2, last; j <= p; j = last + 1) { last = p / (p / j); int q = p / j, id1 = q <= sqrtm ? h1[q] : h2[m / q]; (g[id] += MOD1 - (LL)(last - j + 1) * g[id1] % MOD1) %= MOD1; } } } int get_S(int m) { if (m <= M) { int ans = ksm(pro[m], po[m], MOD); for (int i = 1, last; i <= m; i = last + 1) { last = m / (m / i); ans = (LL)ans * ksm((LL)pro_inv[last] * pro[i - 1] % MOD, po[m - m / i], MOD) % MOD; } return ans; } pre_calc_g(m); int ans = 1; for (int i = 1, last; i <= m; i = last + 1) { last = m / (m / i); int p = m / i, id = p <= sqrtm ? h1[p] : h2[m / p]; ans = (LL)ans * ksm((LL)jc1[last] * ksm(jc1[i - 1], MOD - 2, MOD) % MOD, g[id], MOD) % MOD; } return ans; } void pre_calc_G(int num) { for (int i = num; i >= 1; i--) { int p = w[i], id = p <= sqrtm ? h1[p] : h2[m / p]; G1[id] = ksm(p, n, MOD1); G[id] = 1; for (int j = 2, last; j <= p; j = last + 1) { last = p / (p / j); int q = p / j, id1 = q <= sqrtm ? h1[q] : h2[m / q]; (G1[id] += MOD1 - (LL)(last - j + 1) * G1[id1] % MOD1) %= MOD1; G[id] = (LL)G[id] * ksm(G[id1], last - j + 1, MOD) % MOD; G[id] = (LL)G[id] * ksm((LL)jc1[last] * ksm(jc1[j - 1], MOD - 2, MOD) % MOD, G1[id1], MOD) % MOD; } G[id] = (LL)get_S(p) * ksm(G[id], MOD - 2, MOD) % MOD; } } int get_pro(int l, int r) { int ans = (LL)ksm(jc1[r], r, MOD) * ksm(jc2[r - 1], MOD - 2, MOD) % MOD; if (l <= 1) return ans; int w = (LL)ksm(jc1[l - 1], l - 1, MOD) * ksm(jc2[l - 2], MOD - 2, MOD) % MOD; return (LL)ans * ksm(w, MOD - 2, MOD) % MOD; } int solve() { pre(); int num = 0; sqrtm = sqrt(m); for (int i = 1, last; i <= m; i = last + 1) { last = m / (m / i); w[++num] = m / i; if (m / i <= sqrtm) h1[m / i] = num; else h2[last] = num; } pre_calc_G(num); int ans = 1; for (int i = 1, last; i <= m; i = last + 1) { last = m / (m / i); int p = m / i, id = p <= sqrtm ? h1[p] : h2[m / p]; ans = (LL)ans * ksm(G[id], (LL)(last - i + 1) * (last + i) / 2 % MOD1, MOD) % MOD; ans = (LL)ans * ksm(get_pro(i, last), G1[id], MOD) % MOD; } return ans; } int main() { int T; scanf("%d", &T); init(); while (T--) { scanf("%d%d", &n, &m); printf("%d\n", solve()); } return 0; }
    Processed: 0.028, SQL: 9