6

I want to have a numpy array with mpz/mpfr values. Because my code:

import numpy as np
import gmpy2
A=np.ones((5,5));
print A/gmpy2.mpfr(1);

generates:

RuntimeWarning: invalid value encountered in divide
  print A/gmpy2.mpfr(1);
[[1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]
 [1.0 1.0 1.0 1.0 1.0]]

Which as I can understand is the impossibility to convert gmpy mpfr to numpy float64. So how can I get a numpy array with mpfr values in the first place?

Thanks.

Cupitor
  • 11,007
  • 19
  • 65
  • 91

2 Answers2

8

You will need to create your array with dtype=object, and then you can use any python type inside your array. I don't have gmpy2 installed, but the following example should show how it works:

In [3]: a = np.ones((5, 5), dtype=object)

In [5]: import fractions

In [6]: a *= fractions.Fraction(3, 4)

In [7]: a
Out[7]: 
array([[3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4],
       [3/4, 3/4, 3/4, 3/4, 3/4]], dtype=object)

Having a numpy array of dtype=object can be a liitle misleading, because the powerful numpy machinery that makes operations with the standard dtypes super fast, is now taken care of by the default object's python operators, which means that the speed will not be there anymore:

In [12]: b = np.ones((5, 5)) * 0.75

In [13]: %timeit np.sum(a)
1000 loops, best of 3: 1.25 ms per loop

In [14]: %timeit np.sum(b)
10000 loops, best of 3: 23.9 us per loop
Jaime
  • 65,696
  • 17
  • 124
  • 159
  • Then again, `fractions.Fraction` is not an especially fast class. I wonder what the speed delta between native Numpy arrays and an `mpfr` array, seeing as `mpfr` is a relatively low-level C wrapper class. – nneonneo Mar 09 '13 at 07:19
  • 1
    @nneonneo I believe the problem is not that much the speed of the operations, but the fact that there are Python function calls involved in every single one of them, something that doesn't happen with the other numpy dtypes. – Jaime Mar 09 '13 at 07:31
  • Yes, there are Python function calls, but for a class implemented in C, the overhead of these calls might be pretty small. `Fraction` is implemented in pure Python so each call is many bytecode instructions. – nneonneo Mar 09 '13 at 07:33
  • @nneonneo I just repeated the timings with an array of numpy floats and another one of dtype object, so of Python floats, and you do have a point. For a `(5, 5)` array, the performance differences are less than 5%, although numpy comes ahead. But with a `(500, 500)` array numpy is over 25x faster. – Jaime Mar 09 '13 at 14:52
  • If you were using `=` instead of `*=`, you would not need to specify `dtype`. The issue is probably a problem in one of the two libraries. – leewz Jul 04 '15 at 00:37
1

I believe this is a bug in one of the two libraries. I also believe it is fixed.

Input:

import sys
import numpy as np
import gmpy2

print(sys.version)
print(np.__version__)
print(gmpy2.version)

A=np.ones((5,5));
print(A/gmpy2.mpfr(1))

Output:

3.4.2 (v3.4.2:ab2c023a9432, Oct  6 2014, 22:15:05) [MSC v.1600 32 bit (Intel)]
1.9.1
2.0.5
[[mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]
 [mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0') mpfr('1.0')]]

Either Numpy didn't properly say what to do when it encountered an unknown type, or gmpy2 didn't specify how to get divided by something (__rdiv__).

It is not necessary to specify the dtype of an ndarray unless you intend to write over its elements. Operations like multiplication will result in a new ndarray, and Numpy will figure out what dtype to use.

leewz
  • 3,201
  • 1
  • 18
  • 38