[求め方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
}