LOJ #2000. 「SDOI2017」数字表格

zrzring 2020-07-23 PM 95℃ 0条

这个数据范围一看就是要对每个fib[i]计算贡献,那么

$$ ans = \prod_{k = 1}^{min(n ,m)} fib[k]^{\sum_{i = 1}^{n}\sum_{j = 1}^{m}[gcd(i, j)==k]} $$

出现了非常套路的莫反,我们对fib[k]的指数变换一下

$$ \sum_{i = 1}^{n}\sum_{j = 1}^{m}[gcd(i, j)==k] = \sum_{d = 1}^{\lfloor\frac{n}{k}\rfloor}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\mu(d) $$

于是计算

$$ ans = \prod_{k = 1}^{min(n ,m)} fib[k]^{\sum_{d = 1}^{\lfloor\frac{n}{k}\rfloor}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\mu(d)} $$

但是这个指数预处理是$O(N\sqrt{N})$的,不能接受

$$ \prod_{k = 1}^{min(n ,m)} fib[k]^{\sum_{d = 1}^{\lfloor\frac{n}{k}\rfloor}\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor\mu(d)} = \prod_{k = 1}^{min(n, m)}\left(\prod_{d|k}fib[d]^{\mu(\lfloor\frac{k}{d}\rfloor)}\right)^{\lfloor\frac{n}{k}\rfloor\lfloor\frac{m}{k}\rfloor} $$

这样一来括号内的式子可以做到$O(NlogN)$预处理,单次计算$O(\sqrt{N})$的复杂度,可以通过此题

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <queue>
#define inf 1e9
using namespace std;
void file() {
    freopen("read.in", "r", stdin);
    freopen("write.out", "w", stdout);
}
const int N = 1e6 + 10, dsq = 1e9 + 7;
inline int read() {
    int sym = 0, res = 0; char ch = getchar();
    while (!isdigit(ch)) sym |= (ch == '-'), ch = getchar();
    while (isdigit(ch)) res = (res << 3) + (res << 1) + (ch ^ 48), ch = getchar();
    return sym ? -res : res;
}
int n, m, T, cnt, mu[N], flag[N], prime[N], invg[N], g[N], invf[N], f[N];
long long qpow(long long a, long long x) {
    long long res = 1;
    while (x) {if (x & 1) res = res * a % dsq; a = a * a % dsq; x >>= 1;}
    return res;
}
void init() {
    mu[1] = 1;
    for (int i = 2; i <= N - 10; i++) {
        if (!flag[i]) {
            prime[++cnt] = i; mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= N - 10; j++) {
            flag[i * prime[j]] = 1;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0; break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    f[1] = invf[1] = 1;
    for (int i = 2; i <= N - 10; i++) f[i] = (f[i - 1] + f[i - 2]) % dsq, invf[i] = qpow(f[i], dsq - 2);
    for (int i = 0; i <= N - 10; i++) g[i] = 1, invg[i] = 1;
    for (int i = 2; i <= N - 10; i++) {
        for (int j = i, k = 1; j <= N - 10; j += i, k++) {
            if (mu[k] == 1) {
                g[j] = 1ll * g[j] * f[i] % dsq;
                invg[j] = 1ll * invg[j] * invf[i] % dsq;
            }
            if (mu[k] == -1) {
                g[j] = 1ll * g[j] * invf[i] % dsq;
                invg[j] = 1ll * invg[j] * f[i] % dsq;
            }
        }
    }
    for (int i = 2; i <= N - 10; i++) {
        g[i] = 1ll * g[i] * g[i - 1] % dsq;
        invg[i] = 1ll * invg[i] * invg[i - 1] % dsq;
    }
}
long long solve(int n, int m) {
    long long ans = 1;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans = ans * qpow(1ll * g[r] * invg[l - 1] % dsq, 1ll * (n / l) * (m / l) % (dsq - 1)) % dsq;
    }
    return ans;
}
int main() {
    T = read(); init();
    while (T--) {
        n = read(); m = read(); if (n > m) swap(n, m); printf("%lld\n", solve(n, m));
    }
    return 0;
}

评论


您还没有登录呦


登录 注册
Title - Artist
0:00