BZOJ4816 数字表格

4816: [Sdoi2017]数字表格

Time Limit: 50 Sec Memory Limit: 128 MB

Description

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么

f[0]=0

f[1]=1

f[n]=f[n-1]+f[n-2],n>=2

Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

Input

有多组测试数据。

第一个一个数T，表示数据组数。

接下来T行，每行两个数n,m

T<=1000,1<=n,m<=10^6

Output

输出T行，第i行的数是第i组数据的结果

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

题解

这道题很好地延续了SDOI的优良传统，考了一道莫比乌斯反演以供娱乐。

由于我们一眼发现了这是一道莫比乌斯反演水题，正如做所有的莫比乌斯反演一样，我们把要求的式子先写出来并推导

\begin{eqnarray*}
ans & = & \prod_{i}^{n}\prod _{j}^{m}f( \gcd(i,j) ) \\
& = & \prod_{k}^{n}f(k) ^ {\sum_{i}^{\frac{n}{k}}\sum_{j}^{\frac{m}{k}}[\gcd(i,j)=1]} \\
& = & \prod_{k}^{n}f(k) ^ {\sum_{i}^{\frac{n}{k}}\sum_{j}^{\frac{m}{k}}\sum_{x\mid{i}\&x\mid{j}}~~\mu(x)} \\
& = & \prod_{k}^{n}\prod_{x}^{\frac{n}{k}}(f(k) ^ {\mu(x)})^{\frac{n}{kx}\frac{m}{kx}}
\end{eqnarray*}

为了能够把$f(k) ^ {\mu(x)}$提出来，显然，我们可以设$T=kx$，$g(T)=\prod_{k\mid{T}}f(k) ^ {\mu(\frac{T}{k})}$

化简得到$ans = \prod_{T}^{n}g(T)^{\frac{n}{T}\frac{m}{T}}$

求出$f$以及$f$的逆元，线性筛求$\mu$，Dirichlet卷积求出$g$，然后计算$g$的前缀积$g'$以及$g'$的逆元

查询使用大众喜闻乐见的分块，至此，我们切掉了这道水题。

代码

#include<bits/stdc++.h>
using namespace std;
template <class _T> inline void read(_T &_x) {int _t; bool flag = false;while ((_t = getchar()) != '-' && (_t < '0' || _t > '9')) ;if (_t == '-') _t = getchar(), flag = true; _x = _t - '0';while ((_t = getchar()) >= '0' && _t <= '9') _x = _x * 10 + _t - '0';if (flag) _x = -_x;
}
typedef long long LL;
const int maxn = 1000010;
const int mod = 1000000007;
int f[maxn], f_re[maxn], mu[maxn], g[maxn], g_re[maxn];
bool vis[maxn]; int prime[maxn / 10], pcnt;
#define trans(x) ((int)((LL)x % mod))
#define reg register
inline void update(int &a, reg int b) {a = trans(a * b);if (a < 0) a += mod;
}
inline int qpower(int a, reg LL b) {int ret = 1;while (b) {if (b & 1) update(ret, a);update(a, a), b >>= 1;}return ret;
}
inline int calc(reg int a, reg int b) {if (b == 0) return 1;return b < 0 ? f_re[a] : f[a];
}
void init() {f[0] = 0, f[1] = f_re[1] = g[1] = mu[1] = 1;reg int i, j;for (i = 2; i < maxn; ++i) {f[i] = f[i - 1] + f[i - 2];if (f[i] >= mod) f[i] -= mod;f_re[i] = qpower(f[i], mod - 2);}for (i = 2; i < maxn; ++i) {g[i] = 1;if (!vis[i]) {prime[++pcnt] = i;mu[i] = -1;}for (j = 1; j <= pcnt && prime[j] * i < maxn; ++j) {vis[i * prime[j]] = true;if (i % prime[j] == 0) {mu[i * prime[j]] = 0;break;}mu[i * prime[j]] = -mu[i];}}for (i = 1; i * i < maxn; ++i) {update(g[i * i], calc(i, mu[i]));for (j = i + 1; i * j < maxn; ++j)update(g[i * j], trans(calc(i, mu[j]) * calc(j, mu[i])));}g[0] = g_re[0] = 1;for (i = 1; i < maxn; ++i) {update(g[i], g[i - 1]);g_re[i] = qpower(g[i], mod - 2);}
}
inline int query(reg int a, reg int b) {if (a > b) {int t = a; a = b, b = t; }int ret = 1;for (reg int i = 1, j, x, y; i <= a; i = j + 1) {x = a / i, y = b / i;j = min(a / x, b / y);update(ret, qpower(trans(g[j] * g_re[i - 1]), (LL)x * y));}return ret;
}
int main() {//freopen(".in", "r", stdin);//freopen(".out", "w", stdout);
    init();int T, a, b; read(T);while (T--) {read(a), read(b);printf("%d\n", query(a, b));}return 0;
}

View Code

转载于:https://www.cnblogs.com/akhpl/p/6730864.html