算法1

(快速gcd) $O(n^2)$

一些约定

• N：N为询问的值域。
• Prime：Prime为全体素数集合。

集合

设集合 $\{a_1, a_2, \ldots, a_m\}$ 是n的分解，当且仅当 $a_1 \times a_2 \times \ldots \times a_m = n$。

原理

定理一

内容：可以将值域中的每个 $x$ 分解成 $\{a, b, c\}$，满足 $a, b, c \leq x$ 或 $\in Prime$（定义这种分解为合法分解）。

证明：不妨设 $a \leq b \leq c$。若 $c \notin Prime$ 且 $c > x$，则 $c$ 可分解为 $\{d, e\}$ 且 $d \leq e$，有 $d \leq x$，而 $a \times b = \frac{x}{c} < x$，则有n的分解 $\{d, a \times b, e\}$。若 $e > x$，则可按该规律一直分解直到 $e \in Prime$ 或 $\leq x$。

定理二

内容：对于询问 $\gcd(x, y)$，分别考虑 $a, b, c$ 对答案的贡献，$a$ 对答案的贡献为 $\gcd(a, y)$，再将 $y$ 除以 $\gcd(a, y)$（这个因子已经被算过，不能重复计算），再对 $b, c$ 干同样的事，三个贡献相乘即为 $\gcd(x, y)$。

证明：易得引理若 $r \mid p, r \mid q$ 则 $\gcd(p, q) = r \times \gcd\left(\frac{p}{r}, \frac{q}{r}\right)$。

分别代入：
$$\begin{align*} p_1 &= a \times b \times c, \\ q_1 &= y, \\ r_1 &= \gcd(a, q_1), \\ p_2 &= b \times c, \\ q_2 &= \frac{q_1}{r_1}, \\ r_2 &= \gcd(b, q_2), \\ p_3 &= c, \\ q_3 &= \frac{q_2}{r_2}, \\ r_3 &= \gcd(c, q_3). \end{align*}$$

实现

我们发现实现的难点在于如何在 $O(N)$ 时间内进行分解，询问部分较为容易。

分解

对于 $x \geq 2$，找到 $x$ 的最小质因子 $p$ 以及 $\frac{x}{p}$ 的合法分解 $\{a_0, b_0, c_0\}$ 且 $a_0 \leq b_0 \leq c_0$，$x$ 的一种合法分解即为 $\{a_0 \times p, b_0, c_0\}$ 的升序排序。

考虑证明：

• 当 $x \in Prime$ 时显然成立，分解为 $\{1, 1, x\}$。
• 当 $p \leq \sqrt[4]{x}$ 时将 $a_0 \leq \sqrt[3]{\frac{x}{p}}$ 带入有 $a_0 \times p \leq x$。
• 考虑 $p > \sqrt[4]{x}$ 的情况：
• $a_0 = 1$，显然有 $a_0 \times p = p \leq x$；
• $a_0 \neq 1$，由于 $x$ 不是素数，$\frac{x}{p}$ 的最小质因子 $q$ 即为 $x$ 的第二小质因子，一定 $\geq p$，而 $a_0, b_0, c_0$ 都为 $\frac{x}{p}$ 的非1因子，有 $p \leq q \leq a_0 \leq b_0 \leq c_0$，$p \times a_0 \times b_0 \times c_0 > (\sqrt[4]{x})^4 = x$ 与 $p \times a_0 \times b_0 \times c_0 = x$ 相矛盾，故不存在此情况。

所以只用跑一次线性筛，用最小质因子更新即可，然后预处理出 $\sqrt n \times \sqrt n$ 的 $\gcd$ 数组。

C++ 代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
typedef pair<int,int> PII;
typedef long long LL;
const int N=5010,M=1e6+10,T=1010,mod=998244353;
int pre[T][T];
int a[N],b[N];
int fac[M][3];  // fac为合法分解
int primes[M/10],cnt;
bool st[M];
int n;
void init(int n)
{
    fac[1][0]=fac[1][1]=fac[1][2]=1;
    for(int i=2;i<=n;i++)
    {
        if(!st[i])
        {
            primes[cnt++]=i;
            fac[i][0]=fac[i][1]=1;
            fac[i][2]=i;
        }
        for(int j=0;primes[j]*i<=n;j++)
        {
            int tmp=primes[j]*i;
            st[tmp]=true;
            fac[tmp][0]=primes[j]*fac[i][0];
            fac[tmp][1]=fac[i][1];
            fac[tmp][2]=fac[i][2];
            if(fac[tmp][0]>fac[tmp][1])
                fac[tmp][0]^=fac[tmp][1]^=fac[tmp][0]^=fac[tmp][1];
            if(fac[tmp][1]>fac[tmp][2])
                fac[tmp][1]^=fac[tmp][2]^=fac[tmp][1]^=fac[tmp][2];
            // 对于整数运算a ^= b ^= a ^= b等价于swap(a, b)这里就是手动进行length = 3的排序
            if(i%primes[j]==0)break;
        }
    }
    for(int i=0;i<T;i++)pre[0][i]=pre[i][0]=i;
    for(int i=1;i<T;i++)
        for(int j=1;j<=i;j++)
            pre[i][j]=pre[j][i]=pre[j][i%j];
}
int gcd(int a,int b)
{
    int res=1;
    for(int i=0;i<3;i++)
    {
        int tmp;
        if(fac[a][i]>T)
        {
            if(b%fac[a][i])tmp=1;
            else tmp=fac[a][i];
        }
        else tmp=pre[fac[a][i]][b%fac[a][i]];
        b/=tmp;
        res*=tmp;
    }
    return res;
}
int main()
{
    init(M-1);
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%d",&a[i]);
    for(int i=1;i<=n;i++)scanf("%d",&b[i]);
    for(int i=1;i<=n;i++)
    {
        int pf=1,ans=0;
        for(int j=1;j<=n;j++)
        {
            pf=(LL)pf*i%mod;
            ans=(ans+(LL)pf*gcd(a[i],b[j]))%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}