zoukankan      html  css  js  c++  java
  • 非递归的扩展欧几里得算法

    $DeclareMathOperator{extgcd}{extgcd}$

    作者按:写这篇随笔是为了解释 tourist 的逆元模板

    template <typename T>
    T inverse(T a, T m) {
        T u = 0, v = 1;
        while (a != 0) {
            T t = m / a;
            m -= t * a; swap(a, m);
            u -= t * v; swap(u, v);
        }
        assert(m == 1);
        return u; // 注意:u 可能为负数
    }
    

    扩展欧几里德算法是求不定方程 $ax + by = gcd(a,b)$ 的一组解 $(x, y)$;若 $a$ 和 $b$ 互素,$x$ 即 $a$ 在模 $b$ 下的逆元。注意到方程 $ax + by = gcd(a, b)$ 和方程 $frac{a}{gcd(a, b)} x + frac{b}{gcd(a, b)} y = 1$ 同解。以下只考虑 $a, b$ 互素的情形。常见的求逆元算法实现如下

    int extgcd(int a, int b, int& x, int& y) {
      int d = a;
      if (b != 0) {
        d = extgcd(b, a % b, y, x);
        y -= (a / b) * x;
      } else {
        x = 1; y = 0;
      }
      return d;
    }
    
    int mod_inverse(int a, int m) {
        int x, int y;
        extgcd(a, m, x, y);
        return (m + x % m) % m;
    }
    

    求 $a$ 在模 $m$ 下的逆元时,对于方程 $ax + my = 1$,我们并不需要把 $x, y$ 求出来,而只需要求出 $x$。相比于 tourist 的实现,上述代码是不够高效的,一则含有冗余计算,二则采用递归实现。下面分析 tourist 的非递归实现之原理。

    一般地,可以通过下述过程求出 $gcd(a, m)$。
    令 $r_0 = m$,$r_1 = a$。
    $r_0 = r_1 q_2 + r_ 2 $
    $r_1 = r_2 q_3 + r_3$
    $r_2 = r_3 q_4 + r_4$
    $vdots$
    $color{blue}{r_{i - 2} = r_{i - 1} q_i + r_i}$
    $r_{i - 1} = r_{i} q_{i + 1} + r_{i + 1}$
    $r_{i} = r_{i + 1} q_{i + 2} + r_{i + 2}$
    $vdots$
    $r_{n - 2} = r_{n - 1} q_n + r_{n}$
    $r_{n} = 0$

    我们有

    $gcd(r_0, r_1) = gcd(r_1, r_2) = dots = gcd(r_{n - 1}, r_{n}) = gcd(r_{n-1}, 0) = r_{n - 1} = 1$

    设 $r_i = s_i a + t_i m$,注意到 $r_i = r_{i - 2} - r_{i - 1} q_i $,故有
    egin{aligned}
    r_i &= s_{i - 2} a + t_{i - 2} m - (s_{i - 1} a + t_{i - 1} m) q_i \
    &= (s_{i - 2} - q_i s_{i - 1}) a + (t_{i - 2} - q_i t_{i - 1}) m
    end{aligned}
    即有递推
    egin{aligned}
    s_i &= s_{i - 2} - q_i s_{i - 1}, \
    t_i &= t_{i - 2} - q_i t_{i - 1}.
    end{aligned}

    我们所感兴趣的是等式 $r_{n - 1} = s_{n - 1} a + t_{n - 1} m = 1$,易见 $s_{n-1}$ 即我们欲求的 $a$ 在模 $m$ 下的逆元。注意到 $s_i$ 的递推式中未出现 $t$,因此可以不管 $t$,只考虑 $s$。

    下面考虑边界条件。
    注意到 $r_0 = s_0 a + t_0 m = m$,$r_1 = s_1 a + t_1 m = a$,迳取 $s_0 = 0, s_1 = 1$。
    循环开始之前,变量 u 存放 $s_0$,变量 v 存放 $s_1$。

    对于 $i ge 0$,第 $i$ 轮循环:

    1. 根据等式 $r_{i} = r_{i + 1} q_{i+2} + r_{i + 2} $,算出 $q_{i+2}$(即变量 t)和 $r_{i+2}$;
    2. 根据递推式 $s_{i+2} = s_i - q_{i+2} s_{i + 1} $,算出 $s_{i+2}$。

    第 $i$ 轮循环结束之后

    • 最新的 $r$,即 $r_{i + 2}$,保存在变量 a 中,$r_{i + 1}$ 存放在 m 中;
    • $s_{i + 1}$ 存放在 u 中,$s_{i + 2}$ 存放在 v 中。

    通过上述分析可以知道当循环结束时,即 a 的值是 $r_n = 0$ 时,m 的值恰是 $r_{n - 1} = gcd(a, m)$。

    题外话

    我们看到两种实现所依据的都是欧几里得算法。作为补充,下面来讨论扩展欧几里得算法算出的 $ax + by = gcd(a,b)$ 的解的大小

    用归纳法可以证明,若 $ab e 0$,则 $|x| le b$ 且 $|y|le a$。

    在 $ b = 0$ 的前一步,即 $a \% b = 0$ 时有 $x = 0$ 且 $y = 1$,结论显然成立。假设调用 extgcd(b, a % b, y', x') 后有 $|x'| le b$ 且 $|y'| le a \% b$ 。在 extgcd(a, b, x, y) 中 $x = x', y = y' - (a / b) x '$,所以有如下不等式成立
    $|x| = |x'| le b$,
    $|y| = |y' - (a / b)x'| le |y'| + (a /b )|x'| le a \% b + (a / b) imes b = a$

    References

    https://brilliant.org/wiki/extended-euclidean-algorithm/#extended-euclidean-algorithm

  • 相关阅读:
    Single Number II
    Pascal's Triangle
    Remove Duplicates from Sorted Array
    Populating Next Right Pointers in Each Node
    Minimum Depth of Binary Tree
    Unique Paths
    Sort Colors
    Swap Nodes in Pairs
    Merge Two Sorted Lists
    Climbing Stairs
  • 原文地址:https://www.cnblogs.com/Patt/p/11638063.html
Copyright © 2011-2022 走看看