4

I'm working on an implementation of the Fibonacci sequence in Numpy using the Q-Matrix method. The results are fine up till n = 47. At this point, the matrix_power function is returning incorrect results. Any explanation about why this is happening?

import numpy
def fibonacci(n):
    qmatrix = numpy.matrix([[1, 1], [1, 0]])
    (a,b,c,d) = numpy.linalg.matrix_power(qmatrix,n).flatten().tolist()[0]
    return b
print fibonacci(47) # Outputs -1323752223
Hamman Samuel
  • 2,350
  • 4
  • 30
  • 41
  • 1
    Are you sure? I get `2971215073`. If it is a platform issue, try casting the matrix to another type, like `numpy.linalg.matrix_power(qmatrix,n).astype(numpy.uint64).flatten().tolist()[0]` – alexpeits Jan 27 '17 at 09:20
  • Thanks, indeed it seems like a platform issue. For `n = 47`, casting to `numpy.uint32` fixes the result, but it again gets wrong for larger values. I've not had to do casting before, what's the best approach for getting this working? – Hamman Samuel Jan 27 '17 at 09:27
  • You can read up on [this](https://docs.scipy.org/doc/numpy/user/basics.types.html). Basically a small integer, let's say 4 bit, cannot hold a very large value, and it does weird stuff when you ask to display it, so you cast the array in order to surpass this issue. From what I see, `uint64` is the largest you can get, but I'm no expert. – alexpeits Jan 27 '17 at 09:37

1 Answers1

5

If you are going to be playing around with the Fibonacci numbers, it is probably warranted to sacrifice some speed and use Python's arbitrarily large integers. You can do it by setting your matrix's dtype to object.

You also don't really need to use the np.matrix object, it is almost always better to stick with normal arrays. And you can extract the relevant item without converting your array to a list:

def fibonacci(n):
    qmatrix = numpy.array([[1, 1], [1, 0]], dtype=object)
    return numpy.linalg.matrix_power(qmatrix, n)[0, 1]
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Awesome, this works perfectly. One edit needed in your code is the return should be at index [0,1], according to the Q-matrix theorem. SO won't let me edit because it's just one character :\ – Hamman Samuel Jan 27 '17 at 15:37
  • 1
    Ah, you are right! Fixed it already. – Jaime Jan 27 '17 at 15:39
  • `cumprod` can be used to get a whole range of these `qmatrix` powers (using the `np.matrix` version) https://stackoverflow.com/a/47444215/901925 – hpaulj Nov 22 '17 at 21:54