Python mpmath not arbitrary precision?

Question

I'm trying to continue on my previous question in which I'm trying to calculate Fibonacci numbers using Benet's algorithm. To work with arbitrary precision I found mpmath. However the implementation seems to fail above certain value. For instance the 99th value gives:

218922995834555891712

This should be (ref):

218922995834555169026

Here is my code:

from mpmath import *

def Phi():
    return (1 + sqrt(5)) / 2

def phi():
    return (1 - sqrt(5)) / 2

def F(n):
    return (power(Phi(), n) - power(phi(), n)) / sqrt(5)

start = 99
end = 100

for x in range(start, end):
    print(x, int(F(x)))

score 4 · Answer 1 · answered Aug 16 '16 at 16:32

mpmath does do arbitrary precision math, and it does do it accurately to any precision (as described above) if you are using the arbitrary precision math module and not the default behavior.

mpmath has more than one module which determines the accuracy and speed of the results (to be chosen depending on what you need), and by default it uses Python floats, which is what I believe you saw above.

If you call mpmath's fib( ) having set mp.dps high enough, you will get the correct answer as stated above.

>>> from mpmath import mp
>>> mp.dps = 25
>>> mp.nprint( mp.fib( 99 ), 25 )
218922995834555169026.0
>>> mp.nprint( mpmath.fib( 99 ), 25 )
218922995834555169026.0

Whereas, if you don't use the mp module, you will only get results as accurate as a Python double.

>>> import mpmath
>>> mpmath.dps = 25
>>> mpmath.nprint( mpmath.fib( 99 ), 25
218922995834555170816.0

score 3 · Accepted Answer · answered Sep 07 '14 at 21:19

3

mpmath provides arbitrary precision (as set in mpmath.mp.dps), but still inaccuate calculation. For example, mpmath.sqrt(5) is not accurate, so any calculation based on that will also be inaccurate.

To get an accurate result for sqrt(5), you have to use a library which supports abstract calculation, e.g. http://sympy.org/ .

To get an accurate result for Fibonacci numbers, probably the simplest way is using an algorithm which does only integer arithmetics. For example:

def fib(n):
  if n < 0:
    raise ValueError

  def fib_rec(n):
    if n == 0:
      return 0, 1
    else:
      a, b = fib_rec(n >> 1)
      c = a * ((b << 1) - a)
      d = b * b + a * a
      if n & 1:
        return d, c + d
      else:
        return c, d

  return fib_rec(n)[0]

answered Sep 07 '14 at 21:19

pts

80,836
20
110
183

Ah I understand where the problem is going. Is there any way to calculate the amount of precision needed for a `sqrt`? (I need to calculate it that way) – Ropstah Sep 07 '14 at 21:27
Your algorithm obviously works, and for this part of the problem I'm gonna use that. – Ropstah Sep 07 '14 at 21:28
It's possible to calculate the required precision, but such calculations are usually much more complicated to devise than the actual `sqrt(5)` and Fibonacci calculation. – pts Sep 07 '14 at 21:35
You may want to run a benchmark whether `power` (with the required precision) or `fib_rec` is faster. My expectation is that they are both about O(log(n)). – pts Sep 07 '14 at 21:36
As a matter of fact it seems to work like this! However I do need to measure execution time and start optimizing. – Ropstah Sep 07 '14 at 21:54
I can start playing with the precision indeed which allows me to 'tune' the performance of the code. Using a precision of `n / 4` seems to give correct results. `n / 10` starts to lose precision. I've been running for 14 hours, but calculating `n > 1*10^9` seems impossible using `sympy`. Might `fib_rec` be faster or does your feeling say otherwise? – Ropstah Sep 08 '14 at 19:01
I think `fib_rec` will be the fastest because it's the most direct, without any abstraction or overhead. But measure it for yourself. – pts Sep 08 '14 at 22:09

score 2 · Answer 3 · answered Aug 22 '16 at 01:51

Actually mpmath's default precision is 15 which I think is not enough if you want to get the result of up to 21-digit precision.

One thing you can do is set the precision to be a higher value and use mpmath's defined arithmetic functions for addition, subtraction, etc.

    from mpmath import mp
    mp.dps = 50
    sqrt5 = mp.sqrt(5)
    def Phi():
        return 0.5*mp.fadd(1, sqrt5)

    def phi():
        return 0.5*mp.fsub(1, sqrt5)

    def F(n):
        return mp.fdiv(mp.power(Phi(), n) - mp.power(phi(), n), sqrt5)

    print int(F(99))

This will give you

    218922995834555169026L

Python mpmath not arbitrary precision?

3 Answers3

Linked