高斯整数及其应用

高斯整数

高斯整数定义：形如\(a+b\cdot i\)的复数被称为高斯整数，其中\(a,b \subseteq \mathbb{Z}\)，高斯整数的全体记作\(\mathbb{Z}[i]\)
四则运算：高斯整数的四则运算规则同复数的四则运算规则。
封闭性：设\(\alpha=a+b\cdot i,\beta=c+d\cdot i\)是高斯整数，则\(\alpha+\beta,\alpha-\beta,\alpha\cdot \beta\)均是高斯整数。
范数：对于高斯整数\(\alpha=a+b\cdot i\)，记\(N(\alpha)=a^{2}+b^{2}\)为\(\alpha\)的范数。显然有\(\alpha\cdot \overline{\alpha}=N(\alpha),\quad N(\alpha\cdot \beta)=N(\alpha)\cdot N(\beta)\)。
整除性：高斯整数\(\alpha\)整除高斯整数\(\beta\)当且仅当\(\exists \gamma \subseteq \mathbb{Z}[i],\alpha \cdot \gamma=\beta\)，记作\(\alpha|\beta\)
单位与相伴：若\(\epsilon|1\)，则称高斯整数\(\epsilon\)是单位，若\(\epsilon\)是单位，则称\(\epsilon\cdot \alpha\)是\(\alpha\)的一个相伴。高斯整数的单位有且仅有\(1,-1,i,-i\)。
高斯素数定义：若非零的高斯整数\(\pi\)不是单位，且只能够被单位和他的相伴整除，则称\(\pi\)为高斯素数。
定理1：若高斯整数\(\pi\)，满足\(N(\pi)=p\)，其中\(p\)是有理素数，则\(\pi,\overline{\pi}\)是高斯素数。
证明：

\[\begin{aligned} &假设\pi=\alpha\cdot \beta\Rightarrow N(\pi)=N(\alpha)\cdot N(\beta)=p,那么N(\alpha)=1或N(\beta)=1\\ &这说明\alpha和\beta中有一个是单位,\pi不能分解成两个非单位的高斯素数之积,他是高斯素数 \end{aligned} \]

最大公约数和唯一分解

最大公约数：设\(\alpha,\beta \subseteq \mathbb{Z}[i]\)则称\(\alpha,\beta\)的最大公约数\(\gamma\)为满足以下两个条件的高斯整数：\((i)\gamma|\alpha,\gamma|\beta\quad(ii)\delta|\alpha \wedge \delta|\beta\Rightarrow\delta|\gamma\)
带余除法：设\(\alpha,\beta \subseteq \mathbb{Z}[i]\)且\(\beta\neq 0\)，则存在高斯整数\(\gamma,\rho\)满足\(\alpha=\beta\cdot \gamma+\rho\)且\(0\leq N(\rho)<N(\beta)\)我们称\(\gamma\)为商，\(\rho\)为余数。
证明：

\[\begin{aligned} &假设\alpha=a+b\cdot i,\beta=c+d\cdot i,\frac{a+b\cdot i}{c+d\cdot i}=x+y\cdot i则z+y\cdot i为高斯整数当且仅当\beta|\alpha\\ &令s=\left[ x+\frac{1}{2} \right] ,t=\left[ y+\frac{1}{2} \right]\Rightarrow x+y\cdot i=(s+f)+(t+g)\cdot i，其中f,g \subseteq \mathbb{R}且\left\vert f \right\vert \leq \frac{1}{2} \left\vert g \right\vert \leq \frac{1}{2}\\ &现在令\gamma=s+t\cdot i, \rho=\alpha-\beta\cdot \gamma\\ &那么N(\rho)=N(\beta\cdot (x+y\cdot i-\gamma))=N(\beta)\cdot N(f+g\cdot i)\leq \frac{1}{2}\cdot N(\beta) \end{aligned} \]

定理2：设\(\alpha,\beta \subseteq \mathbb{Z}[i]\)，\((i)\exists \gamma \subseteq \mathbb{Z}[i],满足\gamma是\alpha,\beta的最大公约数\quad (ii)\exists \nu,\mu \subseteq \mathbb{Z}[i],\mu\cdot \alpha+\nu\cdot \beta=\gamma\)。
证明：

\[\begin{aligned} &设集合S=\left\{ N(\mu\cdot \alpha+\nu\cdot \beta)|\mu,\nu \subseteq \mathbb{Z}[i] \right\},那么集合中一定存在最小元\\ &令\gamma=\mu_0\cdot \alpha+\nu_0\cdot \beta,假设N(\gamma)是S中的最小元,我们证明\gamma就是\alpha,\beta的最大公约数\\ &假设\delta|\alpha,\delta|\beta\Rightarrow\alpha=\delta\cdot \omicron,\beta=\delta\cdot \lambda\Rightarrow\gamma=\mu_0\cdot \omicron\cdot \delta+\nu_0\cdot \lambda\cdot \delta\Rightarrow\delta|\gamma\\ &假设\tau=\mu_1\cdot \alpha+\nu_1\cdot \beta，将\tau对\gamma做带余除法，\tau=\omega\cdot \gamma+\xi,N(\xi)\leq N(\gamma)\\ &容易知道\xi也是形如\mu\cdot \alpha+\nu\cdot \beta的高斯整数，因此\xi只能是0，否则将与N(\gamma)是集合S中最小元素矛盾。\\ &因此有\gamma|\tau\Rightarrow \gamma|\alpha,\gamma|\beta \end{aligned} \]

欧几里得算法：令\(\rho_0=\alpha,\rho_1=\beta\)为非零高斯整数，若连续使用高斯整数的带余除法，可以得到\(\rho_j=\rho_{j+1}\cdot \gamma_{j+1}+\rho_{j+2}\)其中\(N(\rho_{j+2})<N(\rho_{j+1}),j=0,1,2, \ldots ,n-2\)并且\(\rho_{n+1}=0\)。则最后一个非零余数\(\rho_n\)就是\(\alpha,\beta\)的最大公约数。
引理1：若\(\pi\)是高斯素数，\(\alpha,\beta\)是高斯整数，且有\(\pi|\alpha\cdot \beta\)，那么有\(\pi|\alpha\)或者\(\pi|\beta\)。
证明：

\[\begin{aligned} &假设\pi\not|\beta，由于\pi是高斯素数，\pi和\beta的最大公约数只能是单位\\ &所以由定理2存在\mu,\nu满足\mu\cdot \beta+\nu\cdot \pi=1\Rightarrow \mu\cdot \beta\cdot \alpha+\nu\cdot \pi\cdot \alpha=\alpha\Rightarrow \pi|\alpha\\ &同理，\pi\not|\alpha\Rightarrow \pi|\beta，引理得证。\\ &这个引理还可以推广到多个数的情形：高斯素数\pi|\alpha_1\cdot \alpha_2\cdots\alpha_n\Rightarrow \exists j \subseteq [1,n],\pi|\alpha_j \end{aligned} \]

唯一分解定理：设\(\gamma\)是非零非单位的高斯整数，则\((i)\gamma能分解为一些高斯素数的乘积\quad (ii)该分解在某种意义上是唯一的，即若\gamma=\rho_1\cdot \rho_2\cdots\rho_s=\pi_1\cdot \pi_2\cdots\pi_t\Rightarrow s=t,并且对这些项重新排列，可以使得\rho_i与\pi_i相伴,i=1,2, \cdots ,s\)
证明：

\[\begin{aligned} &当N(\gamma)=2时，由定理1可知，\gamma是高斯素数\\ &现在假设N(\gamma)>2,且任意范数小于N(\gamma)的高斯整数都能表示为一些高斯素数的乘积\\ &若\gamma是高斯素数，则他显然可以表示为高斯素数的乘积，若他不是高斯素数\\ &则\gamma=\theta\cdot \mu,1<N(\theta)<N(\gamma),1<N(\mu)<N(\gamma),由归纳假设，\theta,\mu都能分解成一些高斯素数之积，故\gamma也能分解成一些高斯素数之积。\\ &我们同样使用第二数学归纳法证明(ii)\\ &若\gamma是高斯素数，则他的分解显然唯一，若其不是高斯素数\\ &假设\gamma=\rho_1\cdot \rho_2\cdots\rho_s=\pi_1\cdot \pi_2\cdots\pi_t\Rightarrow \rho_1|\pi_1\cdot \pi_2\cdots\pi_t\Rightarrow\exists 1\leq j\leq t,\rho_1|\pi_j\\ &我们对\pi_1\cdots \pi_t重新排列，使得\rho_1|\pi_1\Rightarrow\rho_1和\pi_1相伴。然后两边同时除以\pi_1对剩余部分使用归纳假设。\\ &这里我们可以看出唯一分解成立的关键就是引理1的成立。 \end{aligned} \]

高斯整数与平方和

求出所有的正整数\(a,b\)满足\(a^{2}+b^{2}=n\)

定理3：形如\(4k+1\)的有理质数\(p\)可以分解为两个有理整数的平方和。
证明：

\[\begin{aligned} &-1为4k+1形素数的二次剩余，于是存在有理整数t满足t^2\equiv -1(mod p)\\ &p|(t-i)\cdot (t+i)，若p是高斯素数，则由引理1，p|t-i或p|t+i，而能整除p的高斯整数都形如p\cdot (a+b\cdot i)，t+i,t-i显然不满足这种形式。\\ &于是p不是高斯素数，那么存在\alpha,\beta \subseteq \mathbb{Z}[i],N(\alpha)>1,N(\beta)>1,\alpha\cdot \beta=p\Rightarrow N(p)=p^2=N(\alpha)\cdot N(\beta)\Rightarrow N(\alpha)=N(\beta)=p\\ &从而p可以写成两个有理整数的平方和。\\ &假设p=a^2+b^2则p=(a+b\cdot i)\cdot (a-b\cdot i)，由定理1，a+b\cdot i与a-b\cdot i均是高斯素数。\\ &考虑如何求出a和b，如果能找到t，满足t^2\equiv-1\left( \mod p\right)，那么只需要求出\operatorname{gcd}(t+i,p)就能得到a,b\\ &而有一半的数都是p的二次剩余，因此我们可以随机，找到p的一个二次非剩余k，则k^{\frac{p-1}{2}}\equiv -1(\mod p)\\ &取t=k^{\frac{p-1}{4}}即可。 \end{aligned} \]

定理4：形如\(4k+3\)的有理素数\(p\)是高斯素数。
证明：

\[\begin{aligned} &首先对于任意正整数a，a^2\equiv 0或1(\mod 4)\Rightarrow p不能表示为两个有理整数的平方和。\\ &假设p不是高斯素数，则存在\alpha,\beta \subseteq \mathbb{Z}[i],N(\alpha)>1,N(\beta)>1,\alpha\cdot \beta=p\Rightarrow N(p)=p^2=N(\alpha)\cdot N(\beta)\Rightarrow N(\alpha)=N(\beta)=p\\ &从而p可以写成两个有理整数的平方和，这产生了矛盾\\ &于是p是高斯素数。 \end{aligned} \]

平方和：

\[\begin{aligned} &假设n是正整数，且有如下分解形式n=2^w p_1^{a_1}p_2^{a_2}\cdots p_s^{a_s}q_1^{b_1}q_2^{b_2}\cdots q_t^{b_t}，其中p_i是形如4k+1的有理素数，q_i是形如4k+3的有理素数。\\ &接下来将n分解为高斯素数之积n=i^w(1-i)^{2w}(c_1+d_1i)^{a_1}(c_1-d_1i)^{a_1}\cdots ( c_s+d_si)^{a_s}(c_s-d_si)^{a_s}q_1^{b_1}q_2^{b_2}\cdots q_t^{b_t}\\ &将n分解为a^2+b^2即将n分解为(a+bi)\cdot (a-bi)\\ &设a-bi=\epsilon\cdot (1-i)^e(c_1+d_1i)^{f_1}(c_1-d_1i)^{g_1}\cdots (c_s+d_si)^{f_s}(c_s-d_si)^{g_s}q_1^{h_1}q_2^{h_2}\cdots q_t^{h_t}\\ &则a+bi=\overline{\epsilon}\cdot (1+i)^e(c_1-d_1i)^{f_1}(c_1+d_1i)^{g_1}\cdots (c_s-d_si)^{f_s}(c_s+d_si)^{g_s}q_1^{h_1}q_2^{h_2}\cdots q_t^{h_t}\\ &(a-bi)\cdot (a+bi)=2^e(c_1+d_1i)^{f_1+g_1}(c_1-d_1i)^{f_1+g_1}\cdots (c_s+d_si)^{f_s+g_s}(c_s-d_si)^{f_s+g_s}q_1^{2h_1}q_2^{2h_2}\cdots q_t^{2h_t}\\ &因此需要满足e=w,f_i+g_i=a_i(i=1,2, \cdots ,s),2h_i=b_i(i=1,2, \cdots,t),这样我们就得到了要求的a,b\\ &我们可以看到，正整数n要想分解成两个正整数的平方和，充要条件是n的形如4k+3的质因子的次数均为偶数。 \end{aligned} \]

一些例子

将\(n\times m\)的矩形放入坐标系中，使得矩形的四个角坐标均为整点，问所有不同的放置方式下矩形包含的完整\(1\times 1\)正方形数量之和。两个放置方式相同当且仅当他们之间可以仅通过平移重合。\(n,m\leq 10^{18}\)2022牛客暑期多校训练营7 D.The Pool

\[\begin{aligned} &将矩形ABCD中顶点A放置到原点，假设AB=n,AD=m\\ &考虑B的可行位置，即找到整数a,b满足a^2+b^2=n^2，同时需要判断一下对应的D点是否在整点上。\\ &假设B坐标(a,b)，D坐标(c,d)\\ &s_1=以a,b为边长的直角三角形中1\times 1方格的数量\\ &s_2=以c,d为边长的直角三角形中1\times 1方格的数量\\ &则ans=(s_1+s_2)\times 2+(a-c)\times (d-b)\\ &s_1=\sum\limits_{i=0}^{a-1}\left\lfloor \frac{b\cdot x}{a} \right\rfloor可以方便的用类欧几里得算法求解。 \end{aligned} \]

code:

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iostream>
#define int long long

using namespace std;

const int mo = 998244353;
const int d[12] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};

struct gauss {
    __int128 a, b;
    gauss() {}
    gauss(int o, int p) {
        a = o;
        b = p;
    }
    bool operator==(const gauss& o) const { return (a == o.a) && (b == o.b); }
    bool operator<(const gauss& o) const {
        return (a < o.a) || ((a == o.a) && (b < o.b));
    }
    gauss operator-(const gauss& o) const { return gauss(a - o.a, b - o.b); }
    gauss operator*(const gauss& o) const {
        return gauss(a * o.a - b * o.b, a * o.b + b * o.a);
    }
    gauss operator%(const gauss& o) const {
        int t1 = floor(((long double)a / o.b + (long double)b / o.a) /
                           ((long double)o.a / o.b + (long double)o.b / o.a) +
                       0.5);
        int t2 = floor(((long double)b / o.b - (long double)a / o.a) /
                           ((long double)o.a / o.b + (long double)o.b / o.a) +
                       0.5);
        return gauss(a, b) - o * gauss(t1, t2);
    }
};

int n, m, bo[1000005], prime[200005], s0 = 0, cnt = 0, div1[1005][2],
                                      div3[1005][3], s1 = 0, s3 = 0, s2 = 0;
gauss divid[1005], possible[2000005], st[1005][1005];

__int128 abss(__int128 o) { return (o > 0) ? o : (-o); }

int gcd(int o, int p) { return (!p) ? o : gcd(p, o % p); }

int mi(int o, int p, int mo) {
    int yu = 1;
    while (p) {
        if (p & 1) yu = (__int128)yu * o % mo;
        o = (__int128)o * o % mo;
        p >>= 1;
    }
    return yu;
}

bool test(int o, int p) {
    if (mi(o, p - 1, p) != 1) return 0;
    int yu = p - 1;
    while (!(yu & 1)) {
        yu >>= 1;
        if (mi(o, yu, p) == (p - 1)) return 1;
    }
    return mi(o, yu, p) == 1;
}

bool Miller_Rabin(int o) {
    if (o <= 1) return 1;
    for (int i = 0; i < 12; i++)
        if (o == d[i]) return 1;
    for (int i = 0; i < 12; i++)
        if (!test(d[i], o)) return 0;
    return 1;
}

int Pollard_Rho(int o) {
    int x1 = 1ll * rand() * rand() * rand() * rand() % o + 1, c = rand() + 1,
        y = x1;
    for (int lim = 1;; lim = min(lim, lim << 1)) {
        int cnt = 1;
        for (int i = 1; i <= lim; i++) {
            x1 = ((__int128)x1 * x1 + c) % o;
            y = ((__int128)y * y + c) % o;
            y = ((__int128)y * y + c) % o;
            cnt = (__int128)cnt * abs(x1 - y) % o;
        }
        int yu = gcd(cnt, o);
        if (yu > 1) return yu;
    }
    return o;
}

gauss gauss_mi(gauss o, int p) {
    gauss yu = gauss(1, 0);
    while (p) {
        if (p & 1) yu = yu * o;
        o = o * o;
        p >>= 1;
    }
    return yu;
}

gauss gauss_gcd(gauss o, gauss p) {
    if (p == gauss(0, 0)) return o;
    return gauss_gcd(p, o % p);
}

gauss split(int o) {
    int yu = 1ll * rand() * rand() * rand() * rand() % o + 1;
    while (mi(yu, (o - 1) >> 1, o) != (o - 1)) {
        yu = 1ll * rand() * rand() * rand() * rand() % o + 1;
    }
    yu = mi(yu, (o - 1) >> 2, o);
    gauss t = gauss_gcd(gauss(yu, -1), gauss(o, 0));
    int a1 = abss(t.a), a2 = abss(t.b);
    if (a1 > a2) swap(a1, a2);
    return gauss(a1, a2);
}

void insert(gauss o) { possible[++s0] = gauss(abss(o.a), abss(o.b)); }

void dfs(int o, gauss p) {
    if (o > s1) {
        gauss yu = p;
        insert(yu);
        insert(yu * gauss(0, 1));
        return;
    }
    gauss yu = gauss(1, 0);
    st[o][0] = gauss(1, 0);
    for (int i = 1; i <= div1[o][1]; i++)
        st[o][i] = st[o][i - 1] * gauss(divid[o].a, -divid[o].b);
    for (int i = 0; i <= div1[o][1]; i++) {
        dfs(o + 1, p * st[o][div1[o][1] - i] * yu);
        yu = yu * divid[o];
    }
}

int Mo(int o) { return (o >= mo) ? (o - mo) : o; }

int calc(__int128 a, __int128 b, __int128 c, __int128 n) {
    if (a == 0) return ((b / c) % mo) * ((n + 1) % mo) % mo;
    if ((a >= c) || (b >= c)) {
        int s1 = ((n * (n + 1) / 2) % mo * ((a / c) % mo)) % mo;
        int s2 = ((n + 1) % mo) * ((b / c) % mo) % mo;
        return Mo(Mo(s1 + s2) + calc(a % c, b % c, c, n));
    }
    int s1 = ((n % mo) * ((a * n + b) / c) % mo -
              calc(c, c - b - 1, a, (a * n + b) / c - 1)) %
             mo;
    return Mo(s1 + mo);
}

signed main() {
    srand((unsigned)time(NULL));
    for (int i = 2; i <= 100000; i++) {
        if (!bo[i]) prime[++cnt] = i;
        for (int j = 1; (j <= cnt) && (i * prime[j] <= 100000); j++) {
            bo[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }
    int T;
    cin >> T;
    while (T--) {
        scanf("%lld%lld", &n, &m);
        int gg = gcd(n, m);
        s1 = s2 = s3 = s0 = 0;
        for (int i = 1; i <= cnt; i++)
            if (gg % prime[i] == 0) {
                if (prime[i] % 4 == 1) {
                    s1++;
                    div1[s1][0] = prime[i];
                    div1[s1][1] = 0;
                    while (gg % prime[i] == 0) div1[s1][1]++, gg /= prime[i];
                } else if (prime[i] % 4 == 3) {
                    s3++;
                    div3[s3][0] = prime[i];
                    div3[s3][1] = 0;
                    div3[s3][2] = 1;
                    while (gg % prime[i] == 0) {
                        div3[s3][1]++;
                        gg /= prime[i];
                        div3[s3][2] *= prime[i];
                    }
                } else
                    while (gg % 2 == 0) gg >>= 1, s2++;
            }
        while (!Miller_Rabin(gg)) {
            int yu = Pollard_Rho(gg);
            gg /= yu;
            if (yu > gg) swap(yu, gg);
            if (yu > 1) {
                if (yu % 4 == 1) {
                    s1++;
                    div1[s1][0] = yu;
                    div1[s1][1] = 1;
                    while (gg % yu == 0) div1[s1][1]++, gg /= yu;
                } else {
                    s3++;
                    div3[s3][0] = yu;
                    div3[s3][1] = div3[s3][2] = 1;
                    while (gg % yu == 0) {
                        div3[s3][1]++;
                        gg /= yu;
                        div3[s3][2] *= yu;
                    }
                }
            }
        }
        if (gg > 1) {
            if (gg % 4 == 1) {
                s1++;
                div1[s1][0] = gg;
                div1[s1][1] = 1;
            } else {
                s3++;
                div3[s3][0] = gg;
                div3[s3][1] = 1;
                div3[s3][2] = gg;
            }
        }
        s2 <<= 1;
        for (int i = 1; i <= s1; i++) div1[i][1] <<= 1;
        for (int i = 1; i <= s3; i++) div3[i][1] <<= 1;
        for (int i = 1; i <= s1; i++) divid[i] = split(div1[i][0]);
        gauss yu = gauss_mi(gauss(1, -1), s2);
        for (int i = 1; i <= s3; i++) yu = yu * gauss(div3[i][2], 0);
        dfs(1, yu);
        for (int i = 1; i <= s0; i++) {
            possible[i].a = abss(possible[i].a);
            possible[i].b = abss(possible[i].b);
        }
        for (int i = 1; i <= s0; i++)
            if (possible[i].a > possible[i].b)
                swap(possible[i].a, possible[i].b);
        sort(possible + 1, possible + s0 + 1);
        possible[0] = gauss(0, 0);
        int ans = 0;
        gg = gcd(n, m);
        for (int i = 1; i <= s0; i++) {
            if (possible[i] == possible[i - 1]) continue;
            if (!possible[i].a) continue;
            if (!possible[i].b) continue;
            int a = possible[i].a * (n / gg), b = possible[i].b * (n / gg);
            int c = possible[i].b * (m / gg), d = possible[i].a * (m / gg);
            int yu = calc(b, 0, a, a - 1) * 2 % mo;
            yu = (yu + calc(c, 0, d, d - 1) * 2) % mo;
            yu = (yu + (__int128)(c - a) * (b - d) % mo) % mo;
            yu = (yu + mo) % mo;
            if (possible[i].a != possible[i].b) yu = yu * 2 % mo;
            ans = (ans + yu) % mo;
        }
        ans = (ans + (__int128)n * m % mo) % mo;
        if (n != m) ans = ans * 2 % mo;
        printf("%lld\n", ans);
    }
    return 0;
}