[Luogu P3704] [BZOJ 4816] [SDOI2017]数字表格

题目描述

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

f[0]=0,f[1]=1f[0]=0,f[1]=1f[0]=0,f[1]=1,

f[n]=f[n−1]+f[n−2],n≥2f[n]=f[n-1]+f[n-2],n\geq 2f[n]=f[n−1]+f[n−2],n≥2

Doris用老师的超级计算机生成了一个n×mn×mn×m的表格，

第iii行第jjj列的格子中的数是f[gcd⁡(i,j)]f[\gcd(i,j)]f[gcd(i,j)]，其中gcd⁡(i,j)\gcd(i,j)gcd(i,j)表示i,ji,ji,j的最大公约数。

Doris的表格中共有n×mn×mn×m个数，她想知道这些数的乘积是多少。

答案对109+710^9+7109+7取模。

输入输出格式

输入格式：

有多组测试数据。

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

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

输出格式：

输出TTT行，第iii行的数是第iii组数据的结果

输入输出样例

输入样例#1：

输出样例#1：

1
6
960

说明

对10%10\%10%的数据，1≤n,m≤1001\leq n,m\leq 1001≤n,m≤100

对30%30\%30%的数据，1≤n,m1\leq n,m1≤n,m

另外存在30%30\%30%的数据，T≤3T\leq 3T≤3

对100%100\%100%的数据，T≤1000,1≤n,m≤106T\leq1000,1\leq n,m\leq 10^6T≤1000,1≤n,m≤106

解题分析

∏i=1n∏j=1mF[gcd(i,j)]=∏d=1nF[d]∑i=1n∑j=1m[gcd(i,j)=d]=∏d=1nF[d]∑p=1⌊nd⌋μ(p)⌊ndp⌋⌊mdp⌋\prod_{i=1}^{n}\prod_{j=1}^mF[gcd(i,j)] \\ =\prod_{d=1}^{n}F[d]^{\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=d]} \\ =\prod_{d=1}^{n}F[d]^{\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\lfloor\frac{n}{dp}\rfloor\lfloor\frac{m}{dp}\rfloor} i=1∏nj=1∏mF[gcd(i,j)]=d=1∏nF[d]∑i=1n∑j=1m[gcd(i,j)=d]=d=1∏nF[d]∑p=1⌊dn⌋μ(p)⌊dpn⌋⌊dpm⌋

套路地设T=dpT=dpT=dp
∏d=1nF[d]∑p=1⌊nd⌋μ(p)⌊ndp⌋⌊mdp⌋=∏T=1n∏d∣TF[d]μ(Td)⌊nT⌋⌊mT⌋=∏T=1n∏d∣T(F[d]μ(Td))⌊nT⌋⌊mT⌋\prod_{d=1}^{n}F[d]^{\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\lfloor\frac{n}{dp}\rfloor\lfloor\frac{m}{dp}\rfloor} \\ =\prod_{T=1}^{n}\prod_{d|T}F[d]^{\mu(\frac{T}{d})\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} \\ =\prod_{T=1}^{n}\prod_{d|T}(F[d]^{\mu(\frac{T}{d})})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} d=1∏nF[d]∑p=1⌊dn⌋μ(p)⌊dpn⌋⌊dpm⌋=T=1∏nd∣T∏F[d]μ(dT)⌊Tn⌋⌊Tm⌋=T=1∏nd∣T∏(F[d]μ(dT))⌊Tn⌋⌊Tm⌋
然后就可以筛出μ(i)\mu(i)μ(i)后O(Nln(N))O(Nln(N))O(Nln(N))预处理， O(TNlog(N))O(T\sqrt N log(N))O(TNlog(N))回答询问了。

代码如下：

#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define MOD 1000000007ll
#define ll long long
#define MX 1000500
template <class T>
IN void in(T &x)
{x = 0; R char c = gc;for (; !isdigit(c); c = gc);for (;  isdigit(c); c = gc)x = (x << 1) + (x << 3) + c - 48;
}
template <class T> IN T max(T a, T b) {return a > b ? a : b;}
template <class T> IN T min(T a, T b) {return a < b ? a : b;}
IN void add(R int ad, int &tar) {tar += ad; if (tar > MOD) tar -= MOD;}
int pcnt;
int miu[MX], pri[MX];
ll sum[MX], F[MX], inv[MX];
bool npr[MX];
IN ll fpow(R ll base, R ll tim)
{ll ret = 1;W (tim){if (tim & 1) ret = ret * base % MOD;base = base * base % MOD, tim >>= 1;}return ret;
}
void get()
{miu[1] = 1; R int tar;for (R int i = 2; i <= 1e6; ++i){if (!npr[i]) pri[++pcnt] = i, miu[i] = -1;for (R int j = 1; j <= pcnt; ++j){tar = i * pri[j];if (tar > 1e6) break;npr[tar] = true;if (!(i % pri[j])) {miu[tar] = 0; break;}miu[tar] = -miu[i];}}sum[0] = F[1] = 1;for (R int i = 2; i <= 1e6; ++i) F[i] = (F[i - 1] + F[i - 2]) % MOD;for (R int i = 1; i <= 1e6; ++i) inv[i] = fpow(F[i], MOD - 2), sum[i] = 1;for (R int i = 1; i <= 1e6; ++i)for (R int j = 1, tar; (tar = j * i) <= 1e6; ++j){if (miu[j] == 1) sum[tar] = sum[tar] * F[i] % MOD;else if (miu[j] == -1) sum[tar] = sum[tar] * inv[i] % MOD;}for (R int i = 2; i <= 1e6; ++i) sum[i] = sum[i - 1] * sum[i] % MOD;
}
int main(void)
{ll ans;int T, n, m; in(T); get();W (T--){in(n), in(m); ans = 1;if (n > m) std::swap(n, m);for (R int lef = 1, rig; lef <= n; lef = rig + 1){rig = std::min(n / (n / lef), m / (m / lef));ans = ans * fpow(sum[rig] * fpow(sum[lef - 1], MOD - 2) % MOD, 1ll * (n / lef) * (m / lef)) % MOD;}printf("%lld\n", ans);}
}