数论3：线性同余方程组

欧几里得算法

又称辗转相处法，用以求两个数的最大公约数，引入符号\((a,b)\)表示a、b的最大公约数，\(a|b\)表示a能整除b，\(a\not |\ \ b\)表示a不能整除b。

现要证明\((a,b) = (b,a \mod b)\)
即证\(a = bq + r\)中\((a,b)=(b,r)\)
令\((a,b) = d_1\)，\((b,r) = d_2\)
有\(d_1(a^, - b^,q) = r\)，即\(d_1 | r\)
其已知\(d_1 | b\)，所以\(d_1 | (b,r)\)，即\(d_1 | d_2\)
同理有\(d_2 | d_1\)，所以得\(d_1 = d_2\)

因为\((n, 0) = n\)，所以这也是欧几里得算法的出口。

时间复杂度\(O(\lg \max(a,b))\)，因为每次迭代数字的大小至少减少一半。

int gcdr(int a, int b)
{
    if(b == 0) return a;
    return gcd(b, a%b);
}

int gcd(int a, int b)
{
    while(b) {
        int r = a % b;
        a = b, b = r;
    }
    return a;
}

扩展欧几里得算法

用以求不定方程\(ax + by = (a,b)\)的一组特解，为方便令\((a,b) = d\)，引入符号\(\lfloor \frac a b \rfloor\)表示下取整，讨论其解法。

当\(b=0\)时，特解\(x=1, y=0\)。
当\(b\neq 0\)，令\(r=a \mod b\)考虑方程

\[by_1 + rx_1 = (b, r) \]

因为\((a,b)=(b,r)\)，所以

\[ax + by = by_1 + rx_1 \]

而\(r = a \mod b = a - \lfloor \frac a b \rfloor b\)，所以

\[ax + by = by_1 + (a - \lfloor \frac a b \rfloor b)x_1 \]

整理得

\[ax + by = ax_1 + b(y_1 - \lfloor \frac a b \rfloor x_1) \]

得递归解

\[x = x_1 \qquad y = y_1 - \lfloor \frac a b \rfloor x_1 \]

现一组特解为\(x_0, y_0\)，则通解为

\[x = x_0-bt \quad y = y_0 + at \]

时间复杂度同样为\(O(\lg\max(a,b))\)，因为递归的执行次数的决定因素和欧几里得算法一样。

int exgcd(int a, int b, int &x, int &y)
{
    if(b == 0) { x = 1, y = 0; return a; }
    int d = exgcd(b, a%b, y, x);
    y -= (a / b) * x;
    return d;
}

扩展欧几里得算法求\(ax + by = c\)

可以通过\(ax_0 + by_0 = d\)的特解放大倍数求得。
若\(d \not |\ \ c\)，则无解；
若\(d \ \ \, |\ \ c\)，则无解有以下讨论：
两边同时乘\(\frac c d\)，可有

\[a\frac {x_0 c}{d} + b \frac {y_0c}{d} = c \]

特解为

\[x_1 = \frac {x_0 c}{d} \quad y_1 = \frac {y_0c}{d} \]

通解为

\[x = x_1-\frac b d t \quad y = y_1 + \frac a d t \]

// 有解返回true，否则false
bool solve1(int a, int b, int c, int &x, int &y)
{
    int d = exgcd(a, b, x, y);
    if(c % d != 0) return 0;
    int t = c / d;
    x *= t, y *= t;
    return 1;
}

扩展欧几里得算法求\(ax \equiv b \pmod m\)

\(a \equiv b \pmod m\)，a、b同余代表除m时余数相同，即\(a - \lfloor \frac a m \rfloor m = b - \lfloor \frac b m \rfloor m\)，整理得

\[mk = a-b，即m | (a-b) （k为整数） \]

所以上述同余方程可写为：

\[ax - my = b \]

这就可以用扩展欧几里得算法求\(ax + by = c\)方法求解了。

// 有解返回true，否则false
bool solve2(int a, int b, int m, int &x)
{
    int y;
    return solve1(a, -m, b, x, y);
}

线性同余方程组\(a_ix \equiv b_i \pmod{m_i}\)

我们无法做到一次求解多个方程，那么考虑如何将其化为单一方程求解。

我们的目标是\(x = b \pmod m\)，即b和m。

为方便解决，引入方程\(x \equiv 0 \pmod 1\)，现只需要考虑以下两个方程的合并问题

\[\begin{cases} x &\equiv b_0 & \pmod {m_0} \quad &①\\ a_1x &\equiv b_1 & \pmod {m_1} \quad &② \end{cases} \]

由①得到

\[x=b_0 + m_0u \quad ③ \]

代入②得

\[a_1m_0u \equiv b_1 - a_1b_0 \pmod {m_1} \]

仅当\(d=(a_1m_0, m_1) | (b_1 - a_1b_0)\)时，整个方程组有解。
通过扩展欧几里得算法求\(ax \equiv b \pmod m\)的方法求得解\(u = u_0 \pmod {\frac {m_1}{d}}\)，带回③得，\(x = b_0 + m_0u_0 + \frac {m_0m_1}{d}t\)，即

\[x \equiv b_0 + m_0u_0 \pmod {\frac {m_0m_1}{(a_1m_0, m_1)}} \]

至此合并完成。

int n, A[maxn], B[maxn], M[maxn];

// 有解返回true，否则false
bool solve_equ_set(int &b0, int &m0)
{
    int u0;
    b0 = 0, m0 = 1;

    for(int i = 0; i < n; ++i) {
        int na = A[i] * m0;
        int nb = B[i] - A[i]*b0;
        int d = gcd(na, M[i]);
        if(!solve2(na, nb, M[i], u0)) return 0;
        b0 += m0 * u0;
        m0 *= M[i] / d;
        b0 %= m0;
    }
    return 1;
}