模算術

模算術(modular arithmetic)又可稱為同餘算術,是一種特別的整數算術系統。
這個系統要先指定一個整數 m 做為模數(modulus),而系統裡的數值就會像是在 0m1 的圈圈裡面一樣,如果互相運算後的數值要超出了這個範圍,就會繞回範圍內對應的數值。

日常生活中,表示時刻周而復始的時候就很常應用模算術。
比如說 22 點過了 8 個小時之後會說是 6 點,而不會說是 30 點(好吧,事實上有人這麼說)。
一天 24 小時,24 就是這裡的模數,而我們講「幾點」的時候只在乎它距離半夜的 0 點過了多少小時,這種觀點下的 306 會是一樣的意思。

同餘

同餘關係

如同上面一天 24 小時的例子,在這樣周而復始的系統裡面,一個數值不管中間繞了多少圈,我們在乎的只有它比起點多了多少。
換句話說,我們只在乎一個數值除以模數之後的餘數。
所以在談論模算術之前,我們可以先認識同餘關係(congruence),顧名思義就是「餘數相同的關係」。

對於整數 abm,如果 m 整除 ab,就說(在模 m 之下) ab 同餘,記為 ab(modm), 讀做「a 同餘 bm」。

上面的定義並沒有直接提及「餘數」,但意思其實是相同的。
根據除法原理,任何整數都可以寫成 qm+r 的形式,其中 q 是除以模數 m 之後的商,而餘數 r 則滿足不等式 0r<m
如果兩個整數 ab 除以模數 m 時有相同的餘數 r,就代表各別有整數 qaqb 使得 {a=qam+r,b=qbm+r, 那麼 m 就會整除 (qaqb)m=ab,反之亦然。

按照同餘符號的定義,之前的例子就會表達成「306(mod24)」。
這裡表示同餘的符號「」寫得像是等號「=」一樣,就是暗示我們同餘關係是種等價關係:當我們只在乎餘數的時候,同餘的兩數就「視為是相同的」。

不過,互相同餘的整數也是有無限多個的,比如在模 24 之下跟 6 同餘的就有 30542667018 等等。
我們往往會需要一個代表元素來概括這些等價的分身,就像是我們講 6 點而不是 30 點一樣。
在模數 m 的系統中,整數 a 的代表元素 r 通常會取與 a 同餘的數之中滿足 0r<m 的那一個,其實也就是取 a 除以 m 之後的餘數,數學上有時候也會用 amodm 來表示這樣的 r

同餘關係何不直接用餘數定義?

我們討論有理數的時候並不在乎分子分母的數值組合,而只在乎它們的比值,所以說 2/4=3/6
要說兩組數字 (a,b)(c,d) 的比值 a/bc/d 是相等的,有兩種效果一樣的方式:一種是「交叉相乘相等 ad=cb」、一種是「兩者的最簡分數相等」。
雖然平常我們都是以最簡分數當代表來去思考這些比值,但實際要寫數學定義或定理的時候,依賴最簡分數反而不太好用。
比如要說 2/4=3/6,當我們是用「交叉相乘相等」作為比值相等的定義,直接說明 2×6=3×4 就足夠;但如果我們依賴最簡分數,則往往需要多幾個「化成最簡分數」的操作,如果要完整地敘述則通常還需要談到公因數:

gcd(2,4)=2,可得 2/4 的最簡分數是 1/2
並且由 gcd(3,6)=3,可得 3/6 的最簡分數也是 1/2
因為 2/43/6 的最簡分數相等,所以 2/4=3/6

而「同餘關係」與「餘數」的關聯就類似於「比值相等」與「最簡分數」的關聯。
如果我們討論模算術的時候都依賴餘數,就必須先用除法原理講清楚何謂餘數才能進一步討論,並且最後也必須回到餘數才能做結論。
相比之下,雖然將同餘定義為「兩數的差是模數的倍數」可能不太直覺,但論述起來其實相對容易許多。

另一方面,同餘關係的觀點更注重在「視為相同」這件事上。
或許,刻意避免直觀卻過於細節的操作,能讓我們更容易理解它之後衍伸出的各種性質與定理。

程式語言中的 mod

程式語言通常也有提供取餘數的運算子或函式,相當於前面提到的二元運算 amodm
我們可以透過比較兩數的餘數是否相等來判斷兩數是否同餘,也能透過判斷餘數是否為 0 來判斷因倍數關係。

C++ 範例
30 % 24  // 6
54 % 24 == 26670 % 24  // true

bool isEven(int n) {
    return n % 2 == 0;
}
Python 範例
30 % 24  # 6
54 % 24 == 26670 % 24  # True

def is_even(n):
    return n % 2 == 0
Haskell 範例
30 `mod` 24  -- 6
54 `mod` 24 == 26670 `mod` 24  -- True

isEven :: Int -> Bool
isEven n = n `mod` 2 == 0

模算術的加、減、乘

加法、減法、乘法

模算術系統不單單只是「拿餘數表示數值」的系統,既然被稱作算術系統,就代表它在運算方面也是有規則的。

比如 6 點經過 72 小時之後會回到 6 點,因為 72 小時是 24 小時的倍數、是 3 天整。
又比如相對於 9 點,在 30 小時之前的時刻是 3 點,而 30 小時之後的時刻則是 15 點,因為 30 小時是 1 天又 6 小時。

這告訴我們,當我們只在乎運算結果的餘數時,我們不必等到所有加法或減法運算結束之後才取餘數,而是能直接取那些加數或減數的餘數來做運算。
在實際應用時,雖然「先取餘數做運算,再取結果的餘數」可能比起「全部算完再取餘數」多了好幾次取餘數的步驟,但卻會因為減少運算中途的大數字而省下了整體運算所花費的時間與空間。

這樣的好事不僅僅發生在加減法上,乘法也能直接拿餘數來做運算。
而事實上這些性質並不侷限在餘數,只要是互相同餘的數都可以自由替換,這也進一步體現了「視為相同」的觀點。
於是我們有下面的性質,並且在之後附上乘法性質的證明。

如果 a0a1(modm)b0b1(modm),那麼就有
a0+b0a1+b1(modm),a0b0a1b1(modm),a0×b0a1×b1(modm).

如果 a0a1(modm)b0b1(modm),根據定義,代表 a1a0b1b0 都是 m 的倍數。
並且因為
a1b1a0b0=a1b1a1b0+a1b0a0b0=a1(b1b0)+b0(a1a0),
而兩個 m 倍數的整係數組合依然是 m 倍數,於是就得到 a0b0a1b1(modm)

我們也不難驗證說普通整數運算的某些常見性質在模算術系統中也是存在的,比如加法與乘法的交換律與結合律,或乘法對加減法的分配律。

減法與相反數

普通整數運算時,我們會說 ab=a+(b),也就是說「減去某數就是加上某數的相反數」。
在模算術系統中也有一樣的規則,比如說 9 點同時是 17 點之前 8 小時的時刻,也是 17 點之後 16 小時的時刻,用同餘式來表達就是 917817+16(mod24), 這裡的 16 就相當於 8 的相反數,並且我們的確有 816(mod24)

雖然模算術系統所謂的相反數不像是普通整數加法那樣,互相位在數線上的正負兩邊,但它們的性質是一樣的:兩者相加之後會互相消滅而變成加法單位元素。
這裡的「加法單位元素」就是指跟任何其他數相加都不會讓那個其他數有所增減的數,在普通整數加法裡面是 0,而模算術裡的則是與 0 同餘的數,換句話說就是模數的倍數。

如果 a+b0(modm),那麼對任何整數 x 都有 xax+b(modm).

模算術的除法

除法與倒數

不同於加、減、乘,我們會發現在模算術系統裡面,除法並不能跟普通運算直接做對應。
比如說我們知道 7×53511(mod24),但 11 卻不能直接除以 5 得到 7,我們似乎就不能「不在乎一個數值繞了多少圈」。
另外,就算是能夠直接整除的情況,似乎也不能像其他運算那樣自由替換。
比如我們知道 14020(mod24),並且有 140÷10=1420÷10=2,但是 142(mod24)

會發生這些問題都是因為我們直接拿普通運算的方式來做除法,但像這樣直接對應是錯誤的。
為了要弄清楚在模算術的世界裡該怎麼做除法,我們應該從另一個角度看待除法:「除以某數就是乘以某數的倒數」。

普通有理數運算裡面,非零整數 a 的倒數是有理數 1/a,但重點並不是分數,而是互為倒數的兩個數值所滿足的運算性質:兩者相乘之後會互相抵消而變成乘法單位元素。
這裡的「乘法單位元素」就是指跟任何其他數相乘都不會讓那個其他數有所改變的數,也就是 1

上面這段的關於倒數的敘述跟前一節談論相反數的敘述是差不多的,它們都是指運算上相對應的反元素(inverse),互為反元素的數在運算後能夠互相抵消變成不做事的單位元素(identity)。
所以,減法與除法並不是獨立的存在,它們分別是代表能夠抵消加法跟乘法的逆運算:
x+aa=x+a+(a)=x+0=x,x×a÷a=x×a×(1/a)=x×1=x.

既然如此,要在模算術系統裡面做除法,我們也要用相同的運算性質去得到它。
也就是說,當我們要在模算術系統中除以某數,就是要乘以它所謂的倒數。
不難檢查,模算術世界裡的乘法單位元素也是 1(以及與 1 同餘的其他數)。

對於整數 abm,如果 a×b1(modm),就稱(在模 m 之下) ba 的倒數。
我們習慣用 1 次方代表倒數,記為 a1b(modm).

回到這節一開始 11 除以 5 的問題。
我們從 5×5251(mod24) 知道 5 就是自己的倒數,所以說 11÷511×5557(mod24), 於是我們真的從 11 除以 5 得到 7 了。

倒數存在的條件

但我們會發現除以 10 的問題並沒有被解決。
在模 24 之下,我們找不到 10 的倒數,我們甚至會發現如果 10 真的有倒數 r 就會出問題。
例如
1212×112×10×r12×2×5×r0×5×r0(mod24)
120(mod24),顯然這告訴我們並不是所有數都能有倒數。
但要注意,這並不代表在模 24 之下 10 沒有倍數,像是 164×10(mod24) 就是 10 的倍數,只不過在目前的模算術系統下我們沒有辦法從 16 還原回 4 而已。

那麼,倒數存在的條件又是什麼呢?
從上面的例子看到,我們從 10 裡面分解出因數 2 之後,能讓 212 乘出 0,而這就是問題的關鍵。

對於整數 am,如果 am 有公因數 d1,那麼 a 就沒有模 m 之下的倒數。

假設 a 有倒數 r
因為 dam 的公因數,所以存在整數 bn 使得有 a=bdm=nd
那麼就有
nn×1n×a×rn×d×b×r0×b×r0(modm),
但因為 d1,所以 0<n<m,也就是說 n0(modm),有矛盾。
所以說一開始的假設「a 有倒數 r」是錯誤的。

在數學上,如果兩個不同餘 0 的數相乘能同餘 0,那麼這兩個數就被稱為「零因子(zero divisor)」,意思是 0 的因數。
而除以零因子與除以 0 在算術系統上是同樣的罪過。
上面的性質告訴我們如果一個數是零因子的倍數,那麼它就沒有倒數。

那麼反過來,如果一個數不是零因子的倍數,或換句話說它與模數互質,它就會有倒數嗎?
答案是肯定的,下面提供一個用到 Bézout’s 引理的證明,但我們不會去證明那個引理。

對於整數 abc,二元一次方程式 ax+by=c 有整數解 (x,y),若且唯若 cgcd(a,b) 的倍數。

對於整數 ama 有模 m 之下的倒數,若且唯若 am 互質。

我們已經在前面的性質推導過,am 不互質時 a 沒有倒數。

現在假設 am 互質,也就是說 gcd(a,m)=1
根據 Bézout’s 引理,我們知道二元一次方程式 ax+my=1 有整數解,假設為 (r,q)
那麼就有 1=ar+mqar(modm).
也就是說,a 有倒數 a1r(modm)

額外補充:除以不與模數互質的數

前面我們知道了做除法的前提是除數必須要與模數互質。
但如果不互質又會怎麼樣?

我們回去看在模 24 之下的情況。
之前說過透過乘法我們能知道 4×1016(mod24),但沒有除法可以讓我們從 16 還原回 4
如果我們考慮 10 的所有倍數甚至會發現 16×1016(mod24)
顯然,如果想要從 16 除以 10 還原回某個數,就會遇到 416 兩種不同的答案。
這個情況可以非常糟糕,比如考慮 18 的倍數,會發現
122×186×1810×1814×1818×1822×18(mod24),
代表如果想算 12 除以 18 竟然會有 6 種不同的答案!

我們可以稍微深入研究一下。
假設兩個整數 am 不互質,也就是有 gcd(a,m)=d1,那麼就存在整數 bn 使得有 a=bdm=nd
接著就會發現,對於任意的整數 kx 都有
(x+k×n)×ax×a+k×n×ax×a+k×n×d×bx×a+k×0×bx×a(modm).
在模 m 之下,上面的 k 可以代入 0d1 得到不同的值,代表能與 a 乘出同個結果的共有 d 種組合。
所以說,如果除數不與模數互質,那麼除法就一定不會有唯一解。

除非,我們把這些 x+kn 都「視為是相同的」。
如果我們放寬條件,不再固守模數 m,而是換到更小的模數 n 系統裡,我們會獲得除法的唯一解,但同時我們也失去了分辨 x+kn 的能力。
下面的命題是我們勉強能做到的事情。

對於整數 abcm,如果 cab 的公因數,並且 ab(modm),則 acbc(modmgcd(c,m)).

求倒數

前一節我們提到只要除數與模數互質我們就有倒數能做除法,但我們還沒提到如何找出對應的倒數。
我們當然可以窮舉並測試算術系統裡所有的數,但這種暴力做法在模數大的時候顯然會曠日廢時。
那我們有沒有更有效率的方法呢?

下面兩個小節我們分別提供兩種不同的算法,並且都附上對應的程式語言實作範例。

Python 3.8 之後可以直接使用 pow(a, -1, m) 算出 a1(modm)

擴充歐幾里得算法

擴充歐幾里得算法(extended Euclidean algorithm)是歐幾里得算法的擴充。
歐幾里得算法就是俗稱的輾轉相除法,可能可以稱為西方數學最早的演算法,它被使用來求得給定兩個整數 ab 的最大公因數 gcd(a,b)
而當我們把演算法步驟倒推回來時,我們能進一步得到最大公因數的整係數組合 gcd(a,b)=μa+νb,這組整係數 (μ,ν) 恰好就是 Bézout’s 引理提過的整數解。
擴充算法就是建立在這件事上,在算出 gcd(a,b) 的同時算出整係數 μν

有了擴充歐幾里得算法,那就如同之前的證明一樣,如果要求 a1(modm),我們只要對 am 套用算法並得到 a 的係數即可。
下面提供程式語言的實作範例,其中 extEuclid(a,b)=(d,μ,ν) 滿足 μa+νb=d=gcd(a,b)

C++ 範例
std::tuple<int, int, int> extEuclid(int a, int b) {
    if (b == 0)
        return {a, 1, 0};
    auto [d, m, n] = extEuclid(b, a % b);
    return {d, n, m - (a / b) * n};
}

int reciprocal(int a, int m) {
    return std::get<1>(extEuclid(a, m));
}
Python 範例
def ext_euclid(a, b):
    if b == 0:
        return (a, 1, 0)
    q, r = divmod(a, b)
    d, m, n = ext_euclid(b, r)
    return (d, n, m - q * n)

def reciprocal(a, m):
    return ext_euclid(a, m)[1]
Haskell 範例
extEuclid :: Int -> Int -> (Int, Int, Int)
extEuclid a 0 = (a, 1, 0)
extEuclid a b = (d, n, m - n * q)
  where (q, r)    = a `divMod` b
        (d, m, n) = extEuclid b r

reciprocal :: Int -> Int -> Int
reciprocal a m = r
  where (_, r, _) = extEuclid a m

費馬歐拉定理

數論裡面還有費馬歐拉定理,它告訴我們任意整數 a 只要跟模數 m 互質,那麼 a 升到某個次方一定會同餘 1,這個次方可以只跟 m 相關。
而費馬小定理是當 m 剛好是質數時的特例。

下面的 φ 是數論裡的歐拉函數,φ(m) 代表小於 m 且與 m 互質的正整數個數。
比如小於 24 且與 24 互質的正整數恰好是 15711131719238 個,所以 φ(24)=8
我們也不難看出對於任何質數 p,都有等式 φ(p)=p1

對於整數 a 與質數 p,如果 a 不是 p 的倍數,則 ap11(modp).

對於整數 am,如果 am 互質,則 aφ(m)1(modm).

根據定理,我們能看到 a×aφ(m)11(modm),於是
a1aφ(m)1(modm);a1ap2(modp).

實務上通常會利用像是反覆平方法(repeated squaring,或叫做快速冪)的求次方法,並配合前面提過的「先取餘數做運算,再取結果的餘數」來提高計算次方的效率。

雖然 φ(m)1 次方並不是達到 1 所需要的最小次方數,所需的次方數往往可以更小。
例如 φ(24)1=7,但在模 24 之下,除了 1 以外,與 24 互質的其他 7 個數都只要平方就會同餘 1,而不必算到 7 次方。
不過,若因此要去分析各個數值分別要取幾次方就不切實際了,所以實務上的算法都還是直接取 φ(m)1 次方。
這也讓這種算法的運算次數能不依賴於 a 的值,並且當固定模數 m 時,我們可以事先算出 φ(m) 的值。

下面提供程式語言的實作範例,這裡我們只提供模數是質數的狀況,也就是利用到費馬小定理求倒數。

C++ 範例
int modpow(int b, int e, int p) {
    b %= p;
    int result = 1;
    while (e) {
        if (e & 1)
            result = result * b % p;
        b = b * b % p;
        e >>= 1;
    }
    return result;
}

int reciprocal(int a, int p) {
    return modpow(a, p - 2, p);
}
Python 範例
def reciprocal(a, p):
    return pow(a, p - 2, p)
Haskell 範例
modpow :: Int -> Int -> Int -> Int
modpow _ 0 _ = 1
modpow b e p | odd e     = b * r `mod` p
             | otherwise = r
  where r = modpow (b * b `mod` p) (e `div` 2) p

reciprocal :: Int -> Int -> Int
reciprocal a p = modpow a (p - 2) p