问题描述

给定一个 $n \times m$ 的表格，其中格子 $(i,j)$ 的值为 $f[\gcd(i,j)]$，其中 $f$ 是斐波那契数列（$f[0]=0, f[1]=1, f[n]=f[n-1]+f[n-2]$）。要求计算表格中所有数的乘积，对 $10^9+7$ 取模。

解题思路

问题转化

设 $P(n,m) = \prod_{i=1}^{n} \prod_{j=1}^{m} f[\gcd(i,j)]$，我们需要计算 $P(n,m) \mod (10^9+7)$。

核心思想是按 $\gcd(i,j)$ 的值进行分组计算：

• 当 $\gcd(i,j) = d$ 时，可令 $i = d \cdot x, j = d \cdot y$，则 $\gcd(x,y) = 1$
• 满足条件的 $(x,y)$ 对数量为 $\sum_{x=1}^{n/d} \sum_{y=1}^{m/d} [\gcd(x,y) = 1]$

数学推导

利用莫比乌斯反演，上述计数可转化为：

$$\sum_{k=1}^{\min(n/d,m/d)} \mu(k) \cdot \left\lfloor \frac{n}{dk} \right\rfloor \cdot \left\lfloor \frac{m}{dk} \right\rfloor $$

因此乘积可表示为：

$$P(n,m) = \prod_{d=1}^{\min(n,m)} f[d]^{\sum_{k=1}^{\min(n/d,m/d)} \mu(k) \cdot \left\lfloor \frac{n}{dk} \right\rfloor \cdot \left\lfloor \frac{m}{dk} \right\rfloor} $$

通过换元 $T = dk$，进一步化简为：

$$P(n,m) = \prod_{T=1}^{\min(n,m)} \left( \prod_{g|T} f[g]^{\mu(T/g)} \right)^{\left\lfloor \frac{n}{T} \right\rfloor \cdot \left\lfloor \frac{m}{T} \right\rfloor} $$

算法实现

1. 预处理：

   • 线性筛计算莫比乌斯函数 $\mu$
   • 预处理斐波那契数列 $f$ 模 $10^9+7$
   • 计算 $s[T] = \prod_{g|T} f[g]^{\mu(T/g)}$

2. 查询计算：

   • 利用数论分块优化乘积计算
   • 对相同的 $\left\lfloor \frac{n}{T} \right\rfloor \cdot \left\lfloor \frac{m}{T} \right\rfloor$ 区间统一处理
   • 使用前缀积和快速幂提高效率

代码实现

import sys
MOD = 10**9 + 7

def main():
    input = sys.stdin.read().split()
    ptr = 0
    n = int(input[ptr])
    ptr += 1
    m = int(input[ptr])
    ptr += 1
    max_nm = max(n, m)
    
    # 线性筛求莫比乌斯函数
    max_mu = max_nm
    mu = [1] * (max_mu + 1)
    is_prime = [True] * (max_mu + 1)
    prime = []
    for i in range(2, max_mu + 1):
        if is_prime[i]:
            prime.append(i)
            mu[i] = -1
        for p in prime:
            if i * p > max_mu:
                break
            is_prime[i * p] = False
            if i % p == 0:
                mu[i * p] = 0
                break
            else:
                mu[i * p] = -mu[i]
    
    # 预处理斐波那契数列
    f = [0] * (max_mu + 1)
    if max_mu >= 1:
        f[1] = 1
    for i in range(2, max_mu + 1):
        f[i] = (f[i-1] + f[i-2]) % MOD
    
    # 计算f的逆元
    inv_f = [1] * (max_mu + 1)
    for i in range(1, max_mu + 1):
        inv_f[i] = pow(f[i], MOD-2, MOD)
    
    # 计算s[T] = product_{g|T} f[g]^μ(T/g)
    s = [1] * (max_mu + 1)
    for g in range(1, max_mu + 1):
        for j in range(1, max_mu // g + 1):
            T = g * j
            if mu[j] == 1:
                s[T] = s[T] * f[g] % MOD
            elif mu[j] == -1:
                s[T] = s[T] * inv_f[g] % MOD
    
    # 计算前缀积
    prefix = [1] * (max_mu + 1)
    for i in range(1, max_mu + 1):
        prefix[i] = prefix[i-1] * s[i] % MOD
    
    # 数论分块计算结果
    res = 1
    min_nm = min(n, m)
    T = 1
    while T <= min_nm:
        nt = n // T
        mt = m // T
        if nt == 0 or mt == 0:
            break
        
        # 找到最大的r使得n//r == nt且m//r == mt
        r1 = n // nt if nt != 0 else min_nm
        r2 = m // mt if mt != 0 else min_nm
        r = min(r1, r2, min_nm)
        
        # 计算区间[T, r]的乘积
        product = prefix[r] * pow(prefix[T-1], MOD-2, MOD) % MOD
        
        # 计算指数
        exponent = nt * mt % (MOD-1)  # 费马小定理优化
        res = res * pow(product, exponent, MOD) % MOD
        
        T = r + 1
    
    print(res)

if __name__ == "__main__":
    main()

算法复杂度分析

• 预处理时间：$O(n \log n)$，主要用于计算 $s[T]$
• 查询时间：$O(\sqrt{n} + \sqrt{m})$，得益于数论分块优化
• 空间复杂度：$O(n + m)$，用于存储各种数组

该算法能够高效处理 $n, m \leq 10^6$ 的情况，通过数学转化和优化技巧，将原本难以直接计算的双重乘积问题转化为可高效求解的形式。