[BZOJ4816]数字表格

来源：互联网发布：js如何给radio赋值编辑：程序博客网时间：2024/05/18 02:02

题目链接：BZOJ4816

题目大意
求∏1≤i≤n,1≤j≤mF(gcd(i,j))，答案对1e9+7取模，其中F为fibonacci数列，F(0)=0,F(1)=1。

分析
1. 枚举gcd，转化为求∏nd=1F(d)h(d)，其中h(d)=∑ni=1∑mj=1[gcd(i,j)==d]。
2. 而求h(d)作为莫比乌斯反演的经典问题，我们知道结论，于是变成了∏nd=1F(d)(∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋)。
3. 设T=dk，同时外移T，就变成了∏nT=1(∏d∣TF(d)μ(Td))⌊nT⌋⌊mT⌋。
4. 设g(T)=∏d∣TF(d)μ(Td)，那么预处理出g(x)的前缀积以及前缀积的逆元即可N−−√logN求出每次询问。
5. 那么如何求g(x)呢？首先初始化g(x)=1，我们枚举倍数d，对所有d∣T，g(T)=g(T)∗F(d)μ(Td)；因此，我们在求出g之前，需要先筛出μ，求出F，以及F的逆元。
6. 然而因为模数过大，不能O(P)预处理逆元，于是我们需要另一种方法求数列的逆元：对于数列a1,a2,a3...an，我们设Si为前i项积；那么我们可以O(logP)求出Sn的逆元invs[n]=(a1a2a3...an)−1，那么an的逆元inva[n]=Sn−1⋅invs[n]=a−1n，而且invs[n−1]=invs[n]⋅an，这样就可以递推出数列元素的逆元以及元素前缀积的逆元了。
7. 总结：预处理筛出μ，求出F以及其前缀积，然后求逆元；枚举因数求出g；然后求前缀积，求逆元；最后就可以O(TN−−√logN)过掉了。虽然步骤多，但是代码不长。

上代码

#include <cstdio>#include <cstring>#include <iostream>#include <algorithm>using namespace std;typedef long long LL;const int N = 1e6 + 10;const int MOD = 1e9 + 7;int n, m;inline int read() {    register char ch = getchar();    register int ans = 0, neg = 1;    for (; !isdigit(ch); ch = getchar())        if (ch == '-') neg = -1;    for (; isdigit(ch); ch = getchar())        ans = ans * 10 + ch - '0';    return ans * neg;}inline int pls(int a, int b) {    LL ans = (LL)a + b;    return ans >= (LL)MOD ? ans - MOD : ans;}inline int mul(int a, int b) {    return (LL)a * b % MOD;}inline int power(int a, int b) {    int ans = 1;    while (b) {        if (b & 1) ans = mul(ans, a);        b >>= 1, a = mul(a, a);    }    return ans;}bool notPrime[N];int prime[N], mu[N];int fi[N], pf[N], invfi[N], invpf[N];int gi[N], pg[N], invpg[N];void preCalc(int a) {    int cnt = 0, mut = 1;    mu[1] = 1, fi[0] = 0, fi[1] = pf[1] = 1;    for (int i = 2; i <= a; ++i) {        if (!notPrime[i]) prime[++cnt] = i, mu[i] = -1;        for (int j = 1, k; j <= cnt && (k = i * prime[j]) <= a; ++j) {            notPrime[k] = true;            if (i % prime[j] == 0) {                mu[k] = 0; break;            }            mu[k] = -mu[i];        }        pf[i] = mul(pf[i - 1], fi[i] = pls(fi[i - 1], fi[i - 2]));    }    invpf[a] = power(pf[a], MOD - 2);    for (int i = a; i >= 1; --i) {        gi[i] = 1;        invpf[i - 1] = mul(invpf[i], fi[i]);        invfi[i] = mul(pf[i - 1], invpf[i]);    }    for (int i = 2; i <= a; ++i)        for (int j = 1, k; (k = i * j) <= a; ++j)            if (mu[j] == 1) gi[k] = mul(gi[k], fi[i]);            else if (mu[j] == -1) gi[k] = mul(gi[k], invfi[i]);    pg[1] = 1;    for (int i = 2; i <= a; ++i)        pg[i] = mul(pg[i - 1], gi[i]);    invpg[a] = power(pg[a], MOD - 2);    for (int i = a; i >= 1; --i)        invpg[i - 1] = mul(invpg[i], gi[i]);}int main() {    int T = read();    preCalc(1e6);    while (T--) {        n = read(), m = read();        if (n > m) swap(n, m);        int ans = 1;        for (int i = 1, nxt; i <= n; i = ++nxt) {            nxt = min(n / (n / i), m / (m / i));            ans = mul(ans, power(mul(pg[nxt], invpg[i - 1]), (LL)(n / i) * (m / i) % (MOD - 1)));        }        printf("%d\n", ans);    }    return 0;}

以上

0 0