|
| 1 | +## Fibonacci Numbers |
| 2 | +## =========================================================================== |
| 3 | +## Fibonacci numbers are numbers that are part of the Fibonacci sequence. |
| 4 | +## The Fibonacci sequence generally starts with 0 and 1, every next number |
| 5 | +## is the sum of the two preceding numbers: |
| 6 | +## - F(0) = 0; |
| 7 | +## - F(1) = 1; |
| 8 | +## - F(n) = F(n-1) + F(n-2); |
| 9 | +## |
| 10 | +## References: |
| 11 | +## - https://en.wikipedia.org/wiki/Fibonacci_sequence |
| 12 | +## - https://oeis.org/A000045 |
| 13 | +{.push raises: [].} |
| 14 | + |
| 15 | +runnableExamples: |
| 16 | + |
| 17 | + doAssert fibonacciSeqIterative(5) == @[Natural(0), 1, 1, 2, 3] |
| 18 | + |
| 19 | + doAssert fibonacciClosure(4) == Natural(3) |
| 20 | + |
| 21 | + for f in fibonacciIterator(1.Natural..4.Natural): |
| 22 | + echo f # 1, 1, 2, 3 |
| 23 | + |
| 24 | +import std/math |
| 25 | + |
| 26 | +type |
| 27 | + Matrix2x2 = array[2, array[2, int]] |
| 28 | + |
| 29 | +const |
| 30 | + IdentityMatrix: Matrix2x2 = [[1, 0], |
| 31 | + [0, 1]] |
| 32 | + InitialFibonacciMatrix: Matrix2x2 = [[1, 1], |
| 33 | + [1, 0]] |
| 34 | + |
| 35 | + |
| 36 | +## Recursive Implementation |
| 37 | +## --------------------------------------------------------------------------- |
| 38 | + |
| 39 | +func fibonacciRecursive*(nth: Natural): Natural = |
| 40 | + ## Calculates the n-th fibonacci number in a recursive manner. |
| 41 | + ## Recursive algorithm is extremely slow for bigger numbers of n. |
| 42 | + if nth <= 1: return nth |
| 43 | + fibonacciRecursive(nth-2) + fibonacciRecursive(nth-1) |
| 44 | + |
| 45 | + |
| 46 | +## List of first terms |
| 47 | +## --------------------------------------------------------------------------- |
| 48 | + |
| 49 | +func fibonacciSeqIterative*(n: Positive): seq[Natural] = |
| 50 | + ## Generates a list of n first fibonacci numbers in iterative manner. |
| 51 | + result = newSeq[Natural](n) |
| 52 | + result[0] = 0 |
| 53 | + if n > 1: |
| 54 | + result[1] = 1 |
| 55 | + for i in 2..<n: |
| 56 | + result[i] = result[i-1] + result[i-2] |
| 57 | + |
| 58 | + |
| 59 | +## Closed-form approximation of N-th term |
| 60 | +## --------------------------------------------------------------------------- |
| 61 | + |
| 62 | +func fibonacciClosedFormApproximation*(nth: Natural): Natural = |
| 63 | + ## Approximates the n-th fibonacci number in constant time O(1) |
| 64 | + ## with use of a closed-form expression also known as Binet's formula. |
| 65 | + ## Will be incorrect for large numbers of n, due to rounding error. |
| 66 | + const Sqrt5 = sqrt(5'f) |
| 67 | + const Phi = (Sqrt5 + 1) / 2 # golden ratio |
| 68 | + let powPhi = pow(Phi, nth.float) |
| 69 | + Natural(powPhi/Sqrt5 + 0.5) |
| 70 | + |
| 71 | + |
| 72 | +## Closure and Iterator Implementations |
| 73 | +## --------------------------------------------------------------------------- |
| 74 | + |
| 75 | +func makeFibClosure(): proc(): Natural = |
| 76 | + ## Closure constructor. Returns procedure which can be called to get |
| 77 | + ## the next fibonacci number, starting with F(2). |
| 78 | + var prev = 0 |
| 79 | + var current = 1 |
| 80 | + proc(): Natural = |
| 81 | + swap(prev, current) |
| 82 | + current += prev |
| 83 | + result = current |
| 84 | + |
| 85 | +proc fibonacciClosure*(nth: Natural): Natural = |
| 86 | + ## Calculates the n-th fibonacci number with use of a closure. |
| 87 | + if nth <= 1: return nth |
| 88 | + let fib = makeFibClosure() |
| 89 | + for _ in 2..<nth: |
| 90 | + discard fib() |
| 91 | + fib() |
| 92 | + |
| 93 | +proc fibonacciSeqClosure*(n: Positive): seq[Natural] = |
| 94 | + ## Generates a list of n first fibonacci numbers with use of a closure. |
| 95 | + result = newSeq[Natural](n) |
| 96 | + result[0] = 0 |
| 97 | + if n > 1: |
| 98 | + result[1] = 1 |
| 99 | + let fib = makeFibClosure() |
| 100 | + for i in 2..<n: |
| 101 | + result[i] = fib() |
| 102 | + |
| 103 | +iterator fibonacciIterator*(s: HSlice[Natural, Natural]): Natural = |
| 104 | + ## Nim iterator. |
| 105 | + ## Returns fibonacci numbers from F(s.a) to F(s.b). |
| 106 | + var prev = 0 |
| 107 | + var current = 1 |
| 108 | + for i in 0..s.b: |
| 109 | + if i <= 1: |
| 110 | + if i >= s.a: yield i |
| 111 | + else: |
| 112 | + swap(prev, current) |
| 113 | + current += prev |
| 114 | + if i >= s.a: yield current |
| 115 | + |
| 116 | + |
| 117 | +## An asymptotic faster matrix algorithm |
| 118 | +## --------------------------------------------------------------------------- |
| 119 | + |
| 120 | +func `*`(m1, m2: Matrix2x2): Matrix2x2 = |
| 121 | + let |
| 122 | + a = m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0] |
| 123 | + b = m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1] |
| 124 | + z = m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0] |
| 125 | + y = m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1] |
| 126 | + |
| 127 | + [[a, b], [z, y]] |
| 128 | + |
| 129 | +func pow(matrix: Matrix2x2, n: Natural): Matrix2x2 = |
| 130 | + ## Fast binary matrix exponentiation (divide-and-conquer approach) |
| 131 | + var matrix = matrix |
| 132 | + var n = n |
| 133 | + |
| 134 | + result = IdentityMatrix |
| 135 | + while n != 0: |
| 136 | + if n mod 2 == 1: |
| 137 | + result = result * matrix |
| 138 | + matrix = matrix * matrix |
| 139 | + n = n div 2 |
| 140 | + |
| 141 | +func fibonacciMatrix*(nth: Natural): Natural = |
| 142 | + ## Calculates the n-th fibonacci number with use of matrix arithmetic. |
| 143 | + if nth <= 1: return nth |
| 144 | + var matrix = InitialFibonacciMatrix |
| 145 | + matrix.pow(nth - 1)[0][0] |
| 146 | + |
| 147 | + |
| 148 | +when isMainModule: |
| 149 | + import std/unittest |
| 150 | + import std/sequtils |
| 151 | + import std/sugar |
| 152 | + |
| 153 | + const |
| 154 | + GeneralMatrix = [[1, 2], [3, 4]] |
| 155 | + ExpectedMatrixPow4 = [[199, 290], [435, 634]] |
| 156 | + ExpectedMatrixPow5 = [[1069, 1558], [2337, 3406]] |
| 157 | + |
| 158 | + LowerNth: Natural = 0 |
| 159 | + UpperNth: Natural = 31 |
| 160 | + OverflowNth: Natural = 93 |
| 161 | + |
| 162 | + Count = 32 |
| 163 | + OverflowCount = 94 |
| 164 | + |
| 165 | + HighSlice = Natural(32)..Natural(40) |
| 166 | + |
| 167 | + Expected = @[Natural(0), 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, |
| 168 | + 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, |
| 169 | + 75025, 121393, 196418, 317811, 514229, 832040, 1346269] |
| 170 | + ## F0 .. F31 |
| 171 | + ExpectedHigh = @[Natural(2178309), 3524578, 5702887, 9227465, 14930352, |
| 172 | + 24157817, 39088169, 63245986, 102334155] |
| 173 | + ## F32 .. F40 |
| 174 | + |
| 175 | + template checkFib(list: openArray[Natural]) = |
| 176 | + check list == Expected |
| 177 | + |
| 178 | + template checkFib(calcTerm: proc(n: Natural): Natural, |
| 179 | + range = LowerNth..UpperNth) = |
| 180 | + let list = collect(for i in range: calcTerm(i)) |
| 181 | + check list == Expected |
| 182 | + |
| 183 | + template checkFibOverflow(code: typed) = |
| 184 | + expect OverflowDefect: |
| 185 | + discard code |
| 186 | + |
| 187 | + suite "Matrix exponentiation": |
| 188 | + test "Matrix to the power of 0": |
| 189 | + check pow(GeneralMatrix, 0) == IdentityMatrix |
| 190 | + test "Matrix to the power of 1": |
| 191 | + check pow(GeneralMatrix, 1) == GeneralMatrix |
| 192 | + test "Matrix to the power of 4": |
| 193 | + check pow(GeneralMatrix, 4) == ExpectedMatrixPow4 |
| 194 | + test "Matrix to the power of 5": |
| 195 | + check pow(GeneralMatrix, 5) == ExpectedMatrixPow5 |
| 196 | + |
| 197 | + suite "Fibonacci Numbers": |
| 198 | + test "F0..F31 - Recursive Version": |
| 199 | + checkFib(fibonacciRecursive) |
| 200 | + test "F0..F31 - Closure Version": |
| 201 | + checkFib(fibonacciClosure) |
| 202 | + test "F0..F31 - Matrix Version": |
| 203 | + checkFib(fibonacciMatrix) |
| 204 | + test "F0..F31 - Iterative Sequence Version": |
| 205 | + checkFib(fibonacciSeqIterative(Count)) |
| 206 | + test "F0..F31 - Closure Sequence Version": |
| 207 | + checkFib(fibonacciSeqClosure(Count)) |
| 208 | + test "F0..F31 - Closed-form Approximation": |
| 209 | + checkFib(fibonacciClosedFormApproximation) |
| 210 | + test "F0..F31 - Nim Iterator": |
| 211 | + checkFib(fibonacciIterator(LowerNth..UpperNth).toSeq) |
| 212 | + |
| 213 | + test "Closed-form approximation fails when nth >= 32": |
| 214 | + let list = collect(for i in HighSlice: fibonacciClosedFormApproximation(i)) |
| 215 | + check list != ExpectedHigh |
| 216 | + |
| 217 | + #test "Recursive procedure overflows when nth >= 93": # too slow at this point |
| 218 | + # checkFibOverflow(fibonacciRecursive(OverflowNth)) |
| 219 | + test "Closure procedure overflows when nth >= 93": |
| 220 | + checkFibOverflow(fibonacciClosure(OverflowNth)) |
| 221 | + test "Matrix procedure overflows when nth >= 93": |
| 222 | + checkFibOverflow(fibonacciMatrix(OverflowNth)) |
| 223 | + test "Iterative Sequence function overflows when n >= 94": |
| 224 | + checkFibOverflow(fibonacciSeqIterative(OverflowCount)) |
| 225 | + test "Closure Sequence procedure overflows when n >= 94": |
| 226 | + checkFibOverflow(fibonacciSeqClosure(OverflowCount)) |
| 227 | + test "Nim Iterator overflows when one or both slice indexes >= 93": |
| 228 | + checkFibOverflow(fibonacciIterator(LowerNth..OverflowNth).toSeq()) |
0 commit comments