[SDOI 2017]数字表格

Description

题库链接

记 (f_i) 为 (fibonacci) 数列的第 (i) 项。

求 [prod_{i=1}^nprod_{j=1}^mf_{gcd(i,j)}]

对质数取模，多组询问。

(1leq tleq 1000,1leq n,mleq 10^6)

Solution

[egin{aligned}Rightarrow&prod_{d=1}^{min{n,m}}f(d)^{sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)=1]}\=&prod_{d=1}^{min{n,m}}f(d)^{sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}sumlimits_{kmid gcd(i,j)}mu(k)}\=&prod_{d=1}^{min{n,m}}f(d)^{sumlimits_{k=1}^{minleft{leftlfloorfrac{n}{d} ight floor,leftlfloorfrac{m}{d} ight floor ight}}mu(k)leftlfloorfrac{n}{kd} ight floorleftlfloorfrac{m}{kd} ight floor}end{aligned}]

令 (T=kd) [prod_{T=1}^{min{n,m}}prod_{dmid T}f(d)^{muleft(frac{T}{d} ight)leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floor}]

我们把其中类似于狄利克雷卷积形式的东西记做 (F(T)) [prod_{T=1}^{min{n,m}}F(T)^{leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floor}]

那么可以枚举因子来求 (F(T)) ，显然可以在近似于 (O(nln~n)) 的时限内预处理出来。然后数论分块的复杂度为 (O(t sqrt n)) 。

Code

//It is made by Awson on 2018.2.22
#include <bits/stdc++.h>
#define LL long long
#define dob complex<double>
#define Abs(a) ((a) < 0 ? (-(a)) : (a))
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
#define Swap(a, b) ((a) ^= (b), (b) ^= (a), (a) ^= (b))
#define writeln(x) (write(x), putchar('
'))
#define lowbit(x) ((x)&(-(x)))
using namespace std;
const int N = 1e6;
const int yzh = 1e9+7;
void read(int &x) {
    char ch; bool flag = 0;
    for (ch = getchar(); !isdigit(ch) && ((flag |= (ch == '-')) || 1); ch = getchar());
    for (x = 0; isdigit(ch); x = (x<<1)+(x<<3)+ch-48, ch = getchar());
    x *= 1-2*flag;
}
void print(LL x) {if (x > 9) print(x/10); putchar(x%10+48); }
void write(LL x) {if (x < 0) putchar('-'); print(Abs(x)); }

int n, m, t, mu[N+5], f[N+5], F[N+5];
int isprime[N+5], prime[N+5], tot, inv[N+5];

int quick_pow(int a, int b) {
    int ans = 1;
    while (b) {
    if (b&1) ans = 1ll*ans*a%yzh;
    b >>= 1, a = 1ll*a*a%yzh;
    }
    return ans;
}
void get_pre() {
    inv[1] = inv[2] = f[1] = f[2] = 1; for (int i = 3; i <= N; i++) f[i] = (f[i-1]+f[i-2])%yzh, inv[i] = quick_pow(f[i], yzh-2);
    for (int i = 0; i <= N; i++) isprime[i] = F[i] = 1; isprime[1] = 0, mu[1] = 1;
    for (int i = 2; i <= N; i++) {
    if (isprime[i]) prime[++tot] = i, mu[i] = -1;
    for (int j = 1; j <= tot && i*prime[j] <= N; j++)
        if (i%prime[j] != 0) isprime[i*prime[j]] = 0, mu[i*prime[j]] = -mu[i];
        else {isprime[i*prime[j]] = 0, mu[i*prime[j]] = 0; break; }
    }
    for (int i = 1; i <= N; i++)
    for (int j = 1; j*i <= N; j++)
        if (mu[j] == 1) F[i*j] = 1ll*F[i*j]*f[i]%yzh;
        else if (mu[j] == -1) F[i*j] = 1ll*F[i*j]*inv[i]%yzh;
    for (int i = 2; i <= N; i++) F[i] = 1ll*F[i]*F[i-1]%yzh;
}
void work() {
    get_pre(); read(t);
    while (t--) {
    read(n), read(m); if (n > m) Swap(n, m);
    int ans = 1;
    for (int i = 1, last; i <= n; i = last+1) {
        last = Min(n/(n/i), m/(m/i));
        ans = 1ll*ans*quick_pow(1ll*F[last]*quick_pow(F[i-1], yzh-2)%yzh, 1ll*(n/i)*(m/i)%(yzh-1))%yzh;
    }
    writeln(ans);
    }
}
int main() {
    work(); return 0;
}