Precission issue converting SymPy expressions to NumPy longdoubles

user_1069658 0
Viewed 44

Upon converting SymPy expressions to NumPy longdoubles (code snippet below), I noticed that the last three digits in the NumPy results are typically rubbish. I don't understand why this happens, as I'm explicitly indicating that the NumPy results should be of dtype=np.longdouble. What am I doing wrong, and is there a (straightforward) way to fix it?

There doesn't seem to be a single mention of SymPy in NumPy's documentation, so I'm not sure what exactly NumPy's array constructor does when given SymPy input.

    from sympy import *
    import numpy as np

    print(np.finfo(np.longdouble))

    symMat = Matrix([[sqrt(5), sqrt(7)]])
    print(symMat)
    print(N(symMat, 19))

    np.set_printoptions(precision = 18)

    numMat = np.array(symMat, dtype=np.longdouble)
    print("", repr(numMat))
    numMat = np.array(N(symMat, 19), dtype=np.longdouble)
    print("", repr(numMat))

This produces the following output:

    Machine parameters for float128
    ---------------------------------------------------------------
    precision =  18   resolution = 1e-18
    machep =    -63   eps =        1.084202172485504434e-19
    negep =     -64   epsneg =     5.42101086242752217e-20
    minexp = -16382   tiny =       3.3621031431120935063e-4932
    maxexp =  16384   max =        1.189731495357231765e+4932
    nexp =       15   min =        -max
    ---------------------------------------------------------------

    Matrix([[sqrt(5), sqrt(7)]])
    Matrix([[2.236067977499789696, 2.645751311064590591]])
     array([[2.236067977499789805, 2.645751311064590716]], dtype=float128)
     array([[2.236067977499789805, 2.645751311064590716]], dtype=float128)
sympy numpy python
1 Answers

I'm pretty sure np.longdouble just converts the inputs via __float__:

In [4]: np.array([float(N(sqrt(5), 19))], dtype=np.longdouble)
Out[4]: array([2.236067977499789805], dtype=float128)

Python doesn't generally have a good way of interoperating between different non-integer types: everything usually goes through float. That being said numpy doesn't recognise that sympy's Rational follows the numbers ABC and so a better conversion would be possible in that case. For SymPy's Float there is no good way because the Real ABC doesn't give any way to interoperate floating point types apart from converting them to float:

https://docs.python.org/3/library/numbers.html#numbers.Real https://bugs.python.org/issue43602#msg391183

You can do this manually by turning the Float into Rational and extracting the numerators and denominators as integers which np.longdouble can handle:

In [28]: M = Matrix([[sqrt(5), sqrt(7)]])

In [29]: Mfr = M.applyfunc(lambda e: Rational(N(e, 19)))

In [30]: Mfr
Out[30]: 
⎡10312043428088987147  24402748658945006457⎤
⎢────────────────────  ────────────────────⎥
⎣4611686018427387904   9223372036854775808 ⎦

In [31]: Mnum = Mfr.applyfunc(lambda e: e.numerator)

In [32]: Mden = Mfr.applyfunc(lambda e: e.denominator)

In [33]: Mnum
Out[33]: [10312043428088987147  24402748658945006457]

In [34]: Mden
Out[34]: [4611686018427387904  9223372036854775808]

In [35]: to_np = lambda A: np.array([[int(e) for e in row] for row in A.tolist()], dtype=np.longdouble)

In [36]: a = to_np(Mnum) / to_np(Mden)

In [37]: a
Out[37]: array([[2.236067977499789696, 2.64575131106459059 ]], dtype=float128)
user_94509910
Related