fast_fibonacci/
lib.rs

1//! # fast-fibonacci
2//!
3//! `fast-fibonacci` uses linear recurrence to quickly find fib(n, mod) in O(log(n)) time.
4//!
5//! Adapted from http://fusharblog.com/solving-linear-recurrence-for-programming-contest/
6
7use ndarray::arr2;
8use ndarray::Array2;
9use num_bigint::BigUint;
10use num::FromPrimitive;
11
12/// Finds the nth fibonacci number with modulo. Runtime O(log(n))
13///
14/// Uses linear recurrence under the covers.
15/// # Examples
16///
17/// ```
18/// assert_eq!(0, fast_fibonacci::fib_with_mod(0, 10));
19/// assert_eq!(1, fast_fibonacci::fib_with_mod(1, 10));
20/// assert_eq!(1, fast_fibonacci::fib_with_mod(2, 10));
21/// assert_eq!(2, fast_fibonacci::fib_with_mod(3, 10));
22/// assert_eq!(546_875, fast_fibonacci::fib_with_mod(1_000_000_000_000_000, 1_000_000));
23/// assert_eq!(875, fast_fibonacci::fib_with_mod(1_000_000_000_000_000, 1_000));
24/// ```
25pub fn fib_with_mod(n: u64, modulo: u64) -> u64 {
26    if n == 0 {
27        return 0;
28    }
29    if n == 1 {
30        return 1;
31    }
32
33    let f = vec![0, 1];
34    let t = arr2(&[
35        [0, 1], 
36        [1, 1]
37    ]);
38    let power_t = matrix_power_with_mod(&t, n, modulo);
39    let mut answer = 0;
40    for i in 0..2 {
41        answer = (answer + (power_t[[0, i]] * f[i])) % modulo;
42    }
43    answer
44}
45
46
47/// BigUint version of fib_with_mod. Finds the nth fibonacci number with modulo. Runtime O(log(n))
48///
49/// Uses linear recurrence under the covers.
50///
51/// # Examples
52/// ```
53/// use num::FromPrimitive;
54///
55/// let ns               = vec![0, 1, 2, 3, 10, 1_000_000_000_000_000, 1_000_000_000_000_000, 1_955_995_342_096_516];
56/// let modulos          = vec![10, 10, 20, 30, 100, 1_000_000, 1_000, u64::MAX];
57/// let expected_results = vec![0, 1, 1, 2, 55, 546_875, 875, 2_886_946_313_980_141_317];
58///
59/// for i in 0..ns.len() {
60///     assert_eq!(
61///         fast_fibonacci::bigfib_with_mod(
62///             &FromPrimitive::from_u64(ns[i]).unwrap(),
63///             &FromPrimitive::from_u64(modulos[i]).unwrap()
64///         ),
65///         FromPrimitive::from_u64(expected_results[i]).unwrap()
66///     );
67/// }
68/// ```
69///
70/// ```
71/// use num_bigint::BigUint;
72///
73/// let big_n: BigUint = BigUint::from_slice(&[100u32, 100, 100, 100, 15129, 12319]);
74/// let big_modulo: BigUint = BigUint::from_slice(&[14u32, 12, 1923876123, 12]);
75/// let expected_result: BigUint = BigUint::from_slice(&[2743227343u32, 920986447, 1158660944, 5]);
76///
77/// assert_eq!(
78///     fast_fibonacci::bigfib_with_mod(&big_n, &big_modulo),
79///     expected_result
80/// );
81/// ```
82pub fn bigfib_with_mod(n: &BigUint, modulo: &BigUint) -> BigUint {
83    let ZERO: BigUint = FromPrimitive::from_u64(0).unwrap();
84    let ONE: BigUint = FromPrimitive::from_u64(1).unwrap();
85    if n == &ZERO || n == &ONE {
86        return n.clone();
87    }
88
89    let f: Vec<BigUint> = vec![ZERO.clone(), ONE.clone()];
90    let t: Array2<BigUint> = arr2(&[
91        [ZERO.clone(), ONE.clone()],
92        [ONE.clone(), ONE.clone()]
93    ]);
94    let power_t = bigfib_matrix_power(&t, n, modulo);
95    let mut answer: BigUint = ZERO.clone();
96    for i in 0..2 {
97        answer = (answer + (&power_t[[0, i]] * &f[i])) % modulo;
98    }
99    return answer;
100}
101
102
103fn bigfib_matrix_power(mat: &Array2<BigUint>, pow: &BigUint, modulo: &BigUint) -> Array2<BigUint> {
104    let ONE: BigUint = FromPrimitive::from_u64(1).unwrap();
105    let TWO: BigUint = FromPrimitive::from_u64(2).unwrap();
106    if pow == &ONE {
107        return mat.clone();
108    }
109    if pow % &TWO == ONE {
110        return bigfib_multiply(
111            &mat, 
112            &bigfib_matrix_power(mat, &(pow - ONE), modulo),
113            modulo
114        );
115    }
116    let x = bigfib_matrix_power(mat, &(pow / TWO), modulo);
117    bigfib_multiply(&x, &x, modulo)
118}
119
120
121fn bigfib_multiply(a: &Array2<BigUint>, b: &Array2<BigUint>, modulo: &BigUint) -> Array2<BigUint> {
122    let ZERO: BigUint = FromPrimitive::from_u64(0).unwrap();
123    let mut return_mat: Array2<BigUint> = arr2(&[
124        [ZERO.clone(), ZERO.clone()],
125        [ZERO.clone(), ZERO.clone()]
126    ]);
127
128    for i in 0..2 {
129        for j in 0..2 {
130            for k in 0..2 {
131                let big_val: BigUint = &return_mat[[i, j]] + (&a[[i, k]] * &b[[k, j]]);
132                return_mat[[i, j]] = big_val % modulo;
133            }
134        }
135    }
136    return_mat
137}
138
139
140fn multiply_with_mod(a: &Array2<u64>, b: &Array2<u64>, modulo: u64) -> Array2<u64> {
141    let mut return_mat: Array2<u64> = Array2::zeros((2, 2));
142
143    let big_mod: BigUint = FromPrimitive::from_u64(modulo).unwrap();
144    for i in 0..2 {
145        for j in 0..2 {
146            for k in 0..2 {
147                let mat_ij: BigUint = FromPrimitive::from_u64(return_mat[[i, j]]).unwrap();
148                let a_ik: BigUint = FromPrimitive::from_u64(a[[i, k]]).unwrap();
149                let b_kj: BigUint = FromPrimitive::from_u64(b[[k, j]]).unwrap();
150
151                let big_val: BigUint = (mat_ij + (
152                    a_ik * b_kj
153                )) % &big_mod;
154
155                return_mat[[i, j]] = small_big_int_to_u64(&big_val);
156            }
157        }
158    }
159    return_mat
160}
161
162
163fn matrix_power_with_mod(mat: &Array2<u64>, pow: u64, modulo: u64) -> Array2<u64> {
164    if pow == 1 {
165        return mat.clone();
166    }
167    if pow % 2 == 1 {
168        return multiply_with_mod(
169            &mat, 
170            &matrix_power_with_mod(
171                mat, 
172                pow - 1, 
173                modulo
174            ), 
175            modulo
176        );
177    }
178    let x = matrix_power_with_mod(mat, pow / 2, modulo);
179    multiply_with_mod(&x, &x, modulo)
180}
181
182
183fn small_big_int_to_u64(big_int: &BigUint) -> u64 {
184    let mut result: u64 = 0;
185
186    let digits = big_int.to_radix_be(10);
187	for i in 0..digits.len() - 1 {
188		result = result + digits[i] as u64;
189		result = result * 10;
190	}
191	result + digits[digits.len() - 1] as u64
192}
193
194#[cfg(test)]
195mod tests {
196    use crate::*;
197    
198    #[test]
199    fn test_first_few() {
200        assert_eq!(fib_with_mod(0, 10), 0);
201        assert_eq!(fib_with_mod(1, 10), 1);
202        assert_eq!(fib_with_mod(2, 10), 1);
203        assert_eq!(fib_with_mod(3, 10), 2);
204        assert_eq!(fib_with_mod(4, 10), 3);
205        assert_eq!(fib_with_mod(5, 10), 5);
206    }
207
208    #[test]
209    fn test_modulo() {
210        assert_eq!(fib_with_mod(100, 1_000_000_000), 261_915_075);
211        assert_eq!(fib_with_mod(100, 1_000_000), 915_075);
212        assert_eq!(fib_with_mod(100, 1_000), 75);
213        assert_eq!(fib_with_mod(100, 10), 5);
214    }
215    
216    #[test]
217    fn test_big() {
218        assert_eq!(fib_with_mod(1_000_000_000_000_000, 1_000_000), 546_875);
219        assert_eq!(fib_with_mod(1_955_995_342_096_516, u64::MAX), 2_886_946_313_980_141_317);
220    }
221
222    #[test]
223    fn test_bigfib() {
224        let ns = vec![0, 1, 2, 3, 10, 1_000_000_000_000_000, 1_955_995_342_096_516];
225        let modulos = vec![10, 10, 20, 30, 100, 1_000_000, u64::MAX];
226        let expected_results = vec![0, 1, 1, 2, 55, 546_875, 2_886_946_313_980_141_317];
227
228        for i in 0..ns.len() {
229            assert_eq!(
230                bigfib_with_mod(
231                    &FromPrimitive::from_u64(ns[i]).unwrap(), 
232                    &FromPrimitive::from_u64(modulos[i]).unwrap()
233                ),
234                FromPrimitive::from_u64(expected_results[i]).unwrap()
235            );
236        }
237    }
238
239    #[test]
240    fn test_large_bigfib() {
241        let n: BigUint = BigUint::from_slice(&[100u32, 100, 100, 100, 15129, 12319]);
242        let modulo: BigUint = BigUint::from_slice(&[14u32, 12, 1923876123, 12]);
243        let expected_result: BigUint = BigUint::from_slice(&[2743227343u32, 920986447, 1158660944, 5]);
244        assert_eq!( 
245            bigfib_with_mod(&n, &modulo),
246            expected_result
247        );
248    }
249}