How to set the precision for a float in Numpy to 52 decimals?

Question

I need to solve a differential equation numerically, however I have to use a number which is precise up to 52 decimals and I only managed to get 16 decimals with np.rounddouble.

I really want to use numpy for this problem, since I was thought so in class, however if it only works in Maple (which my professor used) I would have to write him about this problem.

Any better method to get more than 16 decimals is also appreciated.

Edit: I tried to use pythons decimals as proposed by you, however after about 16 digits, my number seized to be accurate.

imgur.com/a/ZnjRF2M

Does this have to do with conversion to float between steps?

Last edit: YES it was a conversion to float since I didnt put apostrophs before and after the number!

Answer 1

You can use python's Decimal. Decimal arithmetic is implemented in software so it doesn't have the speed of a double, but then you'll have this same problem in any language. Decimal numbers are generated from some context where you can set precision to any number. There is a default context that works for module level decimals. This code sets the default, builds and array and then does a simple operation on the array.

>>> import numpy as np
>>> from decimal import Decimal
>>> import decimal
>>> decimal.getcontext().prec = 52
>>> arr = np.array([Decimal(i) for i in range(10)])
>>> arr
array([Decimal('0'), Decimal('1'), Decimal('2'), Decimal('3'),
       Decimal('4'), Decimal('5'), Decimal('6'), Decimal('7'),
       Decimal('8'), Decimal('9')], dtype=object)
>>> arr + 3
array([Decimal('3'), Decimal('4'), Decimal('5'), Decimal('6'),
       Decimal('7'), Decimal('8'), Decimal('9'), Decimal('10'),
       Decimal('11'), Decimal('12')], dtype=object)
>>> 

Answer 2

Let's compare some ways of working with decimal.Decimal:

In [217]: import decimal
In [218]: Decimal = decimal.Decimal
In [219]: alist = [Decimal(i) for i in range(10)]
In [220]: alist
Out[220]: 
[Decimal('0'),
 Decimal('1'),
 Decimal('2'),
 Decimal('3'),
 Decimal('4'),
 Decimal('5'),
 Decimal('6'),
 Decimal('7'),
 Decimal('8'),
 Decimal('9')]
In [221]: arr = np.array(alist)
In [222]: arr
Out[222]: 
array([Decimal('0'), Decimal('1'), Decimal('2'), Decimal('3'),
       Decimal('4'), Decimal('5'), Decimal('6'), Decimal('7'),
       Decimal('8'), Decimal('9')], dtype=object)

Some math works fine with object dtype array. Basically it delegates the action to the objects. If they implement the right method, great. If they don't, it will fail (np.exp(arr) will fail for nearly every kind of object, even numbers and numpy arrays, because no one implements a x.exp() method.)

In [223]: arr+10
Out[223]: 
array([Decimal('10'), Decimal('11'), Decimal('12'), Decimal('13'),
       Decimal('14'), Decimal('15'), Decimal('16'), Decimal('17'),
       Decimal('18'), Decimal('19')], dtype=object)
In [224]: timeit arr+10
6.21 µs ± 71.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

But compare that with a list comprehension:

In [225]: [x+10 for x in alist]
Out[225]: 
[Decimal('10'),
 Decimal('11'),
 Decimal('12'),
 Decimal('13'),
 Decimal('14'),
 Decimal('15'),
 Decimal('16'),
 Decimal('17'),
 Decimal('18'),
 Decimal('19')]
In [226]: timeit [x+10 for x in alist]
2.23 µs ± 76.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

The object array may be easier to code, but it isn't faster. The usual scaling disclaimers apply, but I don't think it will make a difference.

Now compare those with regular numeric numpy:

In [227]: y = np.arange(10)
In [228]: y+10
Out[228]: array([10, 11, 12, 13, 14, 15, 16, 17, 18, 19])
In [229]: timeit y+10
1.54 µs ± 23.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

---

I have to take one comment back. decimal does implement a exp method.

In [250]: alist[5].exp()
Out[250]: Decimal('148.4131591025766034211155800')
In [251]: np.exp(arr)
Out[251]: 
array([Decimal('1'), Decimal('2.718281828459045235360287471'),
       Decimal('7.389056098930650227230427461'),
       Decimal('20.08553692318766774092852965'), 
        ...], dtype=object)
In [252]: timeit np.exp(arr)
180 µs ± 6.33 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [254]: timeit [i.exp() for i in alist]
165 µs ± 621 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Again compare that with the pure numpy:

In [255]: timeit np.exp(np.arange(10))
4.22 µs ± 107 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

===

Another option is mpmath. The symbolic package sympy makes a lot of use of mpmath. It seems to be slower than decimal, but it has a lot more mathematical functions. So for differential equation work it may be better, with or without sympy. Beware of trying to use scipy ode solvers with either of these extensions. scipy is an extension of numpy, so it's usefulness with sympy, mpmath (and presumably decimal) will be hit-or-miss.

In [268]: mlist = [mpmath.mpf(i) for i in range(10)]
In [269]: ten=mpmath.mpf(10)
In [270]: timeit [x+ten for x in mlist]
21.6 µs ± 1.28 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [271]: timeit [mpmath.exp(i) for i in mlist]
101 µs ± 281 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Answer 3

I don't think you can.

"NumPy does not provide a dtype with more precision than C’s long double". https://numpy.org/doc/stable/user/basics.types.html

How much precision does a long double give?

"depends on the hardware and on the development environment: specifically, x86 machines provide hardware floating-point with 80-bit precision, and while most C compilers provide this as their long double type,"...

But you want 56 decimals, and 10 to the power of 56 is much more than 2 to the power of 80.