题解说明

问题背景：
该程序处理一个与数论和矩阵相关的问题。输入一个数组 $a[1..n]$，需要通过 $\gcd$、欧拉函数、矩阵行列式等运算，计算出一个模 $10^9+7$ 的结果。

核心思路：
$1.$ 快速读入：

• 使用缓冲区优化的 $getC()$ 和 $read()$，加快大数据输入。

$2.$ 欧拉函数与素数筛：

• $init\_p(m)$ 使用线性筛预处理 $1..m$ 的最小质因子 $v[i]$、素数表 $pr[]$ 和欧拉函数 $\varphi[i]$。
• $\varphi$ 的计算遵循积性函数性质。

$3.$ $\gcd$ 表预处理：

• $g[i][j] = \gcd(i,j)$，通过递归公式 $g[i][j] = g[j][i \bmod j]$ 填充。
• 这样后续 $\gcd$ 查询为 $O(1)$。

$4.$ $d[i]$ 的计算：

• $d[i] = \sum_{j=1}^n \gcd(a[i], a[j])$。
• 表示 $a[i]$ 与所有元素的 $\gcd$ 之和。
• 然后对每个 $a[i]$，在 $t[a[i]]$ 中累加 $d[i]$ 的逆元。

$5.$ $s[i]$ 的计算：

• $s[i] = \sum_{j:\, i \mid j} t[j]$。
• 即对每个 $i$，累加所有倍数 $j$ 的 $t[j]$。

$6.$ 矩阵 $b$ 的构造：

• 初始为单位矩阵。
• 对每个 $(i,j)$，计算 $k = \mathrm{lcm}(i,j) = \dfrac{i \cdot j}{\gcd(i,j)}$。
• 若 $k \le m$，则 $b[i][j] \mathrel{-}= \varphi[j] \cdot s[k]$。
• 这样得到一个 $m \times m$ 的矩阵。

$7.$ 答案 $ans$ 的初始值：

• $ans = \prod_{i=1}^{n-1} d[i]$。

$8.$ 高斯消元求行列式：

• 对矩阵 $b$ 做上三角化。
• 每次找到非零主元，若需要交换行则 $ans$ 取负。
• $ans \mathrel{*}=$ 对角线元素。
• 使用模逆元消去下方行。
• 最终 $ans$ 即为行列式乘上之前的积。

复杂度分析：

• 筛法 $O(m)$。
• $\gcd$ 表 $O(m^2)$。
• $d[i]$ 计算 $O(n^2)$。
• 矩阵构造 $O(m^2)$。
• 高斯消元 $O(m^3)$。
  适合 $n,m \le 5000$ 的范围。

实现细节：

• 所有运算取模 $10^9+7$。
• 使用 $qpow$ 计算逆元。
• 高斯消元时注意行交换对行列式符号的影响。
• $\gcd$ 表和 $\varphi$ 的预处理保证后续计算高效。

总结：
该程序结合了数论（欧拉函数、$\gcd$、$\mathrm{lcm}$）、矩阵行列式（高斯消元）、模运算（快速幂、逆元）等知识点，最终输出一个复杂表达式的模值。整体流程是：
输入数组 $\to$ 预处理 $\varphi$ 和 $\gcd$ $\to$ 计算 $d[i], t[], s[]$ $\to$ 构造矩阵 $b$ $\to$ 高斯消元求行列式 $\to$ 输出答案。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// 快速读取字符（缓冲区优化）
inline char getC() {
    static char *p1, *p2, buf[1 << 20];  // 缓冲区及指针
    // 若缓冲区为空，从标准输入读取数据到缓冲区
    return p1 == p2 ? (p2 = (p1 = buf) + fread(buf, 1, 1 << 20, stdin), p1 == p2 ? '\n' : *p1++) : *p1++;
}

// 快速读取整数（基于getC()，处理正负号）
inline ll read() {
    ll f = 1, r = 0;  // f: 符号标记（1为正，-1为负），r: 结果
    char c = getC();

    // 跳过非数字字符，处理负号
    while (!isdigit(c))
        f ^= c == '-', c = getC();  // 若遇到'-'，f取反

    // 读取数字部分
    while (isdigit(c))
        r = r * 10 + c - 48, c = getC();  // 累加数字

    return f ? r : -r;  // 返回带符号的结果
}

#define mk make_pair
#define pii pair<int, int>
#define fi first
#define se second

// 取最小值模板函数
template<class T> inline void ckmin(T &x, T y) {
    if (y < x)
        x = y;
}

// 取最大值模板函数
template<class T> inline void ckmax(T &x, T y) {
    if (x < y)
        x = y;
}

const int N = 5007, inf = 0x3f3f3f3f, mod = 1e9 + 7;

// 模加法：x += y (mod mod)
inline void inc(int &x, int y) {
    x += y;
    if (x >= mod)
        x -= mod;
}

// 模减法：x -= y (mod mod)
inline void dec(int &x, int y) {
    x -= y;
    if (x < 0)
        x += mod;
}

// 模乘法：x = x * y (mod mod)
inline void mul(int &x, int y) {
    x = (ll)x * y % mod;
}

// 模加法：返回(x + y) mod mod
inline int add(int x, int y) {
    return inc(x, y), x;
}

// 模减法：返回(x - y) mod mod
inline int sub(int x, int y) {
    return dec(x, y), x;
}

// 快速幂：返回a^b mod mod
inline int qpow(int a, ll b) {
    int res = 1;  // 结果初始化为1
    for (; b; b >>= 1, mul(a, a))  // 指数二进制分解，底数平方
        if (b & 1)  // 若当前位为1，乘以当前底数
            mul(res, a);
    return res;
}

int s_p, pr[N], v[N], phi[N];  // s_p: 素数个数，pr: 素数表，v: 最小质因子，phi: 欧拉函数

// 初始化素数表和欧拉函数（埃氏筛法）
inline void init_p(int n) {
    phi[1] = 1;  // 1的欧拉函数为1
    for (int i = 2; i <= n; i++) {
        if (!v[i]) {  // 若i为素数
            v[i] = i;  // 最小质因子为自身
            pr[++s_p] = i;  // 加入素数表
            phi[i] = i - 1;  // 素数的欧拉函数为i-1
        }
        // 遍历已找到的素数，更新i*pr[j]的信息
        for (int j = 1; j <= s_p && i * pr[j] <= n; j++) {
            v[i * pr[j]] = pr[j];  // 标记最小质因子
            if (i % pr[j] == 0) {  // 若pr[j]是i的质因子
                phi[i * pr[j]] = phi[i] * pr[j];  // 欧拉函数公式（p^k的倍数）
                break;
            }
            // 否则，欧拉函数为phi[i] * phi[pr[j]]（积性函数）
            phi[i * pr[j]] = phi[i] * phi[pr[j]];
        }
    }
}

int n, m, a[N], b[N][N], d[N], s[N], t[N], g[N][N];  // 全局变量定义

int main() {
#ifdef LOCAL  // 本地调试模式：从文件读写
    freopen("1.in", "r", stdin);
    freopen("1.out", "w", stdout);
#endif
    n = read();  // 读取n

    // 读取数组a，并记录最大值m
    for (int i = 1; i <= n; i++)
        ckmax(m, a[i] = read());

    init_p(m);  // 初始化素数和欧拉函数，上限为m

    // 预处理gcd表：g[i][j] = gcd(i, j)（递归实现）
    for (int i = 0; i <= m; i++)
        for (int j = 0; j <= i; j++)
            if (!j)  // gcd(i, 0) = i
                g[i][j] = g[j][i] = i;
            else  // 欧几里得算法：gcd(i,j) = gcd(j,i%j)
                g[i][j] = g[j][i] = g[j][i % j];

    // 计算d[i] = 所有a[j]与a[i]的gcd之和（j从1到n）
    for (int i = 1; i < n; i++) {
        for (int j = 1; j <= n; j++)
            d[i] += g[a[i]][a[j]];
        // t[x]累加d[i]的模逆（x = a[i]）
        inc(t[a[i]], qpow(d[i], mod - 2));
    }

    // 计算s[i]：所有j是i的倍数的t[j]之和（即s[i] = sum_{j|i} t[j]）
    for (int i = 1; i <= m; i++)
        for (int j = i; j <= m; j += i)
            inc(s[i], t[j]);

    // 初始化矩阵b
    for (int i = 1; i <= m; i++)
        for (int j = 1; j <= m; j++) {
            b[i][j] = (i == j) ? 1 : 0;  // 初始为单位矩阵
            int k = i * j / g[i][j];  // k = lcm(i,j) = i*j/gcd(i,j)
            if (k <= m)  // 若k在范围内，更新b[i][j]
                dec(b[i][j], (ll)phi[j] * s[k] % mod);
        }

    int ans = 1;  // 最终答案，初始为1

    // 计算所有d[i]的乘积（i从1到n-1）
    for (int i = 1; i < n; i++)
        mul(ans, d[i]);

    // 高斯消元：对矩阵b进行上三角化，同时计算行列式
    for (int i = m; i; i--) {
        // 寻找第i列中下方（包括i）非零元素的行p
        int p = i;
        while (p && !b[p][i])
            p--;
        if (!p) {  // 若第i列全为0，行列式为0
            puts("0"), exit(0);
        }
        if (i ^ p) {  // 若p != i，交换两行，行列式符号取反
            swap(b[i], b[p]), ans = mod - ans;
        }

        mul(ans, b[i][i]);  // 累乘主对角线元素
        int iv = qpow(b[i][i], mod - 2);  // 主对角线元素的模逆

        // 收集第i行中非零元素的列索引
        vector<int> vec;
        for (int j = i; j; j--)
            if (b[i][j])
                vec.push_back(j);

        // 消去下方行（i-1到1）的第i列元素
        for (int j = i - 1; j; j--) {
            if (!b[j][i])  // 若当前行第i列已为0，跳过
                continue;
            // 计算消元系数
            int x = (ll)(mod - b[j][i]) * iv % mod;
            // 更新当前行的所有非零列
            for (auto k : vec)
                inc(b[j][k], (ll)x * b[i][k] % mod);
        }
    }

    printf("%d\n", ans);  // 输出结果
    return 0;
}