不知道有没有用的科技。

$O(v)-O(1)\gcd$，其中 $v$ 是值域。

预处理

对于一个数 $n$，我们把它分解为 $(x,y,z)$，满足 $xyz=n,x\leq y\leq z$，且 $z\leq\sqrt n$ 或 $z$ 是质数。

下面证明其存在性：

若 $z>\sqrt n$ 且 $z$ 是合数，那么 $z$ 存在一个质因数 $\leq\sqrt n$，可以接着分解。

构造方法：

使用递归构造。

  • 若 $n=1$，构造 $(1,1,1)$；

  • 若 $n$ 是质数，构造 $(1,1,n)$；

  • 若 $n$ 是合数，找到 $n$ 的最小质因子，$p$，找到 $m=n/p$ 的一个分解方法 $(x,y,z)$，那么 $n$ 的合法分解即为 $(px,y,z)$ 的升序排序。

证明：

易得 $x\leq\sqrt[3]m$。

  • 若 $p\leq\sqrt[4]{n}$，则 $px\leq\sqrt[3]{np^2}\leq\sqrt n$。

  • 否则 $p>\sqrt[4]n$

    • 若 $x=1$，显然成立。

    • 否则 $m$ 不是质数，可以得到 $x,y,z\ge p$，那么 $pxyz>\sqrt[4]{n}=n$，矛盾，不存在此情况。

询问

要求 $\gcd(a,b)$，分别枚举 $b$ 的三个因子，求出和 $a$ 的 $\gcd$，然后把 $a$ 除以这个 $\gcd$，三个因子的计算结果乘起来就是答案。

然后具体求法：

  • $x\leq\sqrt n$：预处理。

  • $x>\sqrt n,b\bmod x=0$：结果为 $x$。

  • $x>\sqrt n,b\bmod x\neq0$：结果为 $1$。

代码

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <vector>
#include <cstring>
using namespace std;

namespace solve
{
    const int maxn = 1e6 + 10;
    const int lim = 1000;
    const int mod = 998244353;

    int vis[maxn], pri[maxn], tot;
    int div[maxn][3];
    int table[lim + 1][lim + 1];

    void init(int n = 1e6)
    {
        div[1][1] = div[1][2] = div[1][0] = 1;
        for (int i = 2; i <= n; i++)
        {
            if (!vis[i])
            {
                pri[++tot] = i;
                div[i][0] = 1, div[i][1] = 1, div[i][2] = i;
            }
            for (int j = 1; j <= tot && i * pri[j] <= n; j++)
            {
                int x = i * pri[j];
                vis[x] = 1;
                memcpy(div[x], div[i], sizeof(div[i]));
                div[x][0] *= pri[j];
                sort(div[x], div[x] + 3);
                if (i % pri[j] == 0)
                    break;
            }
        }
        for (int i = 0; i <= lim; i++)
            for (int j = 0; j <= lim; j++)
                table[i][j] = __gcd(i, j);
    }

    int get(int a, int b)
    {
        if (a <= lim && b <= lim)
            return table[a][b];
        int res = 1;
        for (int i = 0; i < 3; i++)
        {
            int x;
            if (div[b][i] <= lim)
                x = table[div[b][i]][a % div[b][i]];
            else if (a % div[b][i] == 0)
                x = div[b][i];
            else
                x = 1;
            res *= x;
            a /= x;
        }
        return res;
    }

    int n;
    int a[maxn], b[maxn];

    void main()
    {
        init();
        cin >> n;
        for (int i = 1; i <= n; i++)
            cin >> a[i];
        for (int i = 1; i <= n; i++)
            cin >> b[i];
        for (int i = 1; i <= n; i++)
        {
            long long res = 0;
            for (int j = 1, val = i; j <= n; j++, val = 1ll * val * i % mod)
                res += 1ll * val * get(a[i], b[j]) % mod;
            cout << res % mod << '\n';
        }
    }
}

int main()
{
    ios::sync_with_stdio(false), cin.tie(0);
    int T = 1;
    // cin >> T;
    while (T--)
        solve::main();
}