Modular multiplicative inverse
Keywords: 模反元素、費馬小定理、輾轉相除法
前言
排列組合中，一定很常遇到類似下面的式子：
total = (n!) / (x1! * x2! * x3! * …)
不論是分子，或是分母的每一個階乘，可能都是非常大的數字，若今天存在題目問說，給一個prime value P，想知道total mod p為多少時，該如何找到答案？
這問題看似很簡單，但如果寫成程式就意味著，所有中間過程的數值，都要受到32-bit integer或64bit integer的限制，為了避免溢位，如果運算是加減乘，都可以輕易的在計算的過程中，不斷重複取mod。但除法卻沒有這種特性，該怎麼辦？
這時候就可以輪到模反元素登場。
模反元素
Modular multiplicative inverse
根據維基百科: 模反元素的定義如下：
其中有解的充分必要條件為a和n必須互質。
擴展歐幾里得算法
// 複製用Python程式碼
def ext_euclid(a, b):
s0, s1 = 1, 0
t0, t1 = 0, 1
r0, r1 = a, b while(r1 != 0):
q = r0 // r1
r0, r1 = r1, r0 - q * r1
s0, s1 = s1, s0 - q * s1
t0, t1 = t1, t0 - q * t1
return s0, t0, r0
假設a, b為兩個整數，有gcd(a, b)=m，該算法可得另外兩個整數s和t且
s * a + t * b = m
若a, b互質，s即為整數a對b的模反元素。
其證明的方式很暴力，如下圖，將輾轉相除法(歐幾里得算法)的過程展開，統整出s, t, r的疊代歷程。
費馬小定理
若取模的n是一個質數，則(a^-1) mod n，也就是a的模反元素可使用費馬小定理求得。
a^(p-1) =1 (mod p)
根據上式兩側同乘a^(-1)，可以再擴展為下式：
a^(p-2) = a^(-1) (mod p)
也就是若p是質數，求解a^(-1)等價求解a^(p-2)：
# (a^(p-2)) % ppower(a, p-2, p) // C++, <algorithm>pow(a, p-2, p) # Python, default
結論
a^(-1) mod p不論是透過輾轉相除法或是費馬小定理，都可以高效的得解。
輾轉相除法被證明疊代次數不會超過較小數值的10進位長度的5倍。因此複雜度為log10(min(a, p))。
power可以做對半的切分，因此是log2(p)，雖然複雜度略高一些，但勝在lib是現成的，實務上較為簡單。
最後我們回顧最開始的問題total % p，將分母的部分全部改為模反元素(MMI)，得到：
total % p =
((n!) * MMI(x1!) * MMI(x2!) * MMI(x3!) * …) % p
