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
2 3
4 5
6 7
Sample Output
1
6
960
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 < '' || _t > '')) ;
if (_t == '-') _t = getchar(), flag = true; _x = _t - '';
while ((_t = getchar()) >= '' && _t <= '') _x = _x * + _t - '';
if (flag) _x = -_x;
}
typedef long long LL;
const int maxn = ;
const int mod = ;
int f[maxn], f_re[maxn], mu[maxn], g[maxn], g_re[maxn];
bool vis[maxn]; int prime[maxn / ], 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 < ) a += mod;
}
inline int qpower(int a, reg LL b) {
int ret = ;
while (b) {
if (b & ) update(ret, a);
update(a, a), b >>= ;
}
return ret;
}
inline int calc(reg int a, reg int b) {
if (b == ) return ;
return b < ? f_re[a] : f[a];
}
void init() {
f[] = , f[] = f_re[] = g[] = mu[] = ;
reg int i, j;
for (i = ; i < maxn; ++i) {
f[i] = f[i - ] + f[i - ];
if (f[i] >= mod) f[i] -= mod;
f_re[i] = qpower(f[i], mod - );
}
for (i = ; i < maxn; ++i) {
g[i] = ;
if (!vis[i]) {
prime[++pcnt] = i;
mu[i] = -;
}
for (j = ; j <= pcnt && prime[j] * i < maxn; ++j) {
vis[i * prime[j]] = true;
if (i % prime[j] == ) {
mu[i * prime[j]] = ;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
for (i = ; i * i < maxn; ++i) {
update(g[i * i], calc(i, mu[i]));
for (j = i + ; i * j < maxn; ++j)
update(g[i * j], trans(calc(i, mu[j]) * calc(j, mu[i])));
}
g[] = g_re[] = ;
for (i = ; i < maxn; ++i) {
update(g[i], g[i - ]);
g_re[i] = qpower(g[i], mod - );
}
}
inline int query(reg int a, reg int b) {
if (a > b) {int t = a; a = b, b = t; }
int ret = ;
for (reg int i = , j, x, y; i <= a; i = j + ) {
x = a / i, y = b / i;
j = min(a / x, b / y);
update(ret, qpower(trans(g[j] * g_re[i - ]), (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 ;
}