$nCk \mod m$ の組み合わせの求め方を3つ紹介します。
間違っている、もっといい方法があるなどありましたらコメントで知らせていただけると嬉しいです。
[求め方1] もっとも一般的な方法
フェルマーの小定理より
a^p \equiv a \mod m
よって
a^{-1} \equiv a^{p-2} \mod m
// 冪乗余計算アルゴリズム
// a^b mod m を計算します。
fn pow_mod(a: i64, b: i64, m: i64) -> i64 {
if b == 0 {
1
} else if b % 2 == 0 {
let d = pow_mod(a, b/2, m);
(d * d) % m
} else {
(a * pow_mod(a, b-1, m)) % m
}
}
// nCc mod m を計算します。
fn combination(mut n: i64, c: i64, m: i64) -> i64 {
let mut upe = 1;
let mut dow = 1;
for i in 1..c + 1 {
upe = upe * n % m;
dow = dow * i % m;
n -= 1;
}
upe * pow_mod(dow, m-2, m) % m
}
[求め方2] バイナリ法
[求め方1]より早いです。
バイナリ法はRSA暗号で冪乗余を計算する際に使用されたりします。
combination関数は上と変わりません。
trait MOD<T> {
// べき乗余計算アルゴリズム
// a^n mod m を計算する
// バイナリ法
fn pow_mod(a: T, n: T, m: T) -> T;
}
impl MOD<i64> for i64 {
fn pow_mod(mut a: i64, mut n: i64, m: i64) -> i64 {
let mut res: i64 = 1;
while n > 0 {
if n & 1 == 1 {
res = res * a % m;
}
a = a * a % m;
n >>= 1;
}
res
}
}
impl MOD<usize> for usize {
fn pow_mod(mut a: usize, mut n: usize, m: usize) -> usize {
let mut res: usize = 1;
while n > 0 {
if n & 1 == 1 {
res = res * a % m;
}
a = a * a % m;
n >>= 1;
}
res
}
}
[求め方3] 拡張ユークリッド互除法
今回紹介した3つの方法の中でもっとも速いです。
別記事に詳しく説明しました。
// 拡張ユークリッド互除法
// 入力: 整数 u, v (u > v > 0)
// 返り値: e := v^{-1} mod u を満たす e (=t)
fn extended_euclidean(u: i64, v: i64) -> i64 {
let mut r0 = u;
let mut r1 = v;
let mut s0 = 1;
let mut s1 = 0;
let mut t0 = 0;
let mut t1 = 1;
while r1 != 0 {
let q = r0 / r1;
let r = r0 - q * r1;
let s = s0 - q * s1;
let t = t0 - q * t1;
r0 = r1;
s0 = s1;
t0 = t1;
r1 = r;
s1 = s;
t1 = t;
}
if t0 < 0 {
t0 + u
} else {
t0
}
}
// nCc mod m を計算します。
fn combination(mut n: i64, c: i64, m: i64) -> i64 {
let mut upe = 1;
let mut dow = 1;
for i in 1..c + 1 {
upe = upe * n % m;
dow = dow * i % m;
n -= 1;
}
upe * extended_euclidean(m, dow) % m
}
[PS] パスカルの三角形を使った方法
過去AtCoderでの出題
など