SDOI2017 数字表格

#include <cstdio>
#include <cstring>

#include <algorithm>

using namespace std;

typedef long long ll;

const int N = 1000010, P = 1e9L + 7, PC = 500010;

int pc;
int mu[N], fib[N], prod[N];
int p[PC];
bool f[N];

void pre();
void sieve();
int rev(int x);
int mpow(int x, ll k);
void ex_gcd(int a, int b, int& x, int& y);

int main() {
    int t, n, m, ans;
    pre();
    scanf("%d", &t);
    while (t--) {
        scanf("%d%d", &n, &m);
        if (n > m)
            swap(n, m);
        ans = 1;
        for (int l = 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            ans = (ll)ans * mpow((ll)prod[r] * rev(prod[l - 1]) % P, (ll)(n / l) * (m / l)) % P;
        }
        printf("%d\n", ans);
    }
    return 0;
}

int mpow(int x, ll k) {
    int ret = 1;
    for (; k; k >>= 1) {
        if (k & 1)
            ret = (ll)ret * x % P;
        x = (ll)x * x % P;
    }
    return ret;
}

void ex_gcd(int a, int b, int& x, int& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return;
    }
    ex_gcd(b, a % b, y, x);
    y -= a / b * x;
}

int rev(int a) {
    int x, y;
    ex_gcd(a, P, x, y);
    if (x < 0)
        x += P;
    return x;
}

void sieve() {
    mu[1] = 1;
    for (int i = 2; i < N; ++i) {
        if (!f[i]) {
            p[++pc] = i;
            mu[i] = -1;
        }
        for (int j = 1; (ll)i * p[j] < N; ++j) {
            f[i * p[j]] = true;
            mu[i * p[j]] = -mu[i];
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            }
        }
    }
}

void pre() {
    sieve();
    fib[1] = 1;
    for (int i = 2; i < N; ++i) {
        fib[i] = fib[i - 1] + fib[i - 2];
        if (fib[i] >= P)
            fib[i] -= P;
    }
    fill(prod, prod + N, 1);
    for (int i = 1; i < N; ++i)
        for (int j = 1; i * j < N; ++j)
            switch (mu[j]) {
            case 1:
                prod[i * j] = (ll)prod[i * j] * fib[i] % P;
                break;
            case -1:
                prod[i * j] = (ll)prod[i * j] * rev(fib[i]) % P;
                break;
            }
    for (int i = 2; i < N; ++i)
        prod[i] = (ll)prod[i] * prod[i - 1] % P;
}