Compute matrix inverse with decimal object

There's a related questions Matrix inverse with Decimal type NumPy 2015 a while ago which did not have a definite answer. There's a second question from me Is there a way for python to perform a matrix inversion at 500 decimal precision where hpaulj provided some updated suggestions.

Basically decimal is a standard python library capable of computing arbitrary precession value at arbitrary orders. It can be operated by most of the Numpy function such as poly to evaluate the polynomials

np.polyval([Decimal(1),Decimal(2)], Decimal(3.1) )

Decimal('5.100000000')

It can also be cased to a numpy array or being initiated as a dtype object array(Are Decimal 'dtypes' available in NumPy? 2011).

np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])

array([[Decimal('1'), Decimal('2')],
       [Decimal('3'), Decimal('4')]], dtype=object)

matrix_m=np.zeros((2,2) ,dtype=np.dtype)
for ix in range(0,2):
    for iy in range(0,2):
        matrix_m[ix,iy]=Decimal(ix)+Decimal(iy);

array([[Decimal('0'), Decimal('1')],
       [Decimal('1'), Decimal('2')]], dtype=object)

Some array operation from numpy also worked when Decimal was the element,

np.exp( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]]) )

array([[Decimal('2.718281828'), Decimal('7.389056099')],
       [Decimal('20.08553692'), Decimal('54.59815003')]], dtype=object)

np.sqrt( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]]) )

array([[Decimal('1'), Decimal('1.414213562')],
       [Decimal('1.732050808'), Decimal('2')]], dtype=object)

and, at single element, the numpy calculation agreed with decimal's native function

np.exp(Decimal(1))==Decimal(1).exp()
True

The useful constant was also provided

def pi():
    """Compute Pi to the current precision.
    #https://docs.python.org/3/library/decimal.html

    >>> print(pi())
    3.141592653589793238462643383

    """
    getcontext().prec += 2  # extra digits for intermediate steps
    three = Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n+na, na+8
        d, da = d+da, da+32
        t = (t * n) / d
        s += t
    getcontext().prec -= 2
    return +s               # unary plus applies the new precision

However, it turned out that both the determinate and the inverse of the matrix in the numpy

np.linalg.det(np.array([[Decimal(1),Decimal(2)],[Decimal(1),Decimal(3)]]))

File <__array_function__ internals>:180, in det(*args, **kwargs)
File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:2154, in det(a)
   2152 t, result_t = _commonType(a)
   2153 signature = 'D->D' if isComplexType(t) else 'd->d'
-> 2154 r = _umath_linalg.det(a, signature=signature)
   2155 r = r.astype(result_t, copy=False)
   2156 return r

np.linalg.inv(np.array([[Decimal(1),Decimal(2)],[Decimal(1),Decimal(3)]]))

File <__array_function__ internals>:180, in inv(*args, **kwargs)
File ~\anaconda3\lib\site-packages\numpy\linalg\linalg.py:552, in inv(a)
    550 signature = 'D->D' if isComplexType(t) else 'd->d'
    551 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 552 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    553 return wrap(ainv.astype(result_t, copy=False))

returned the same error

UFuncTypeError: Cannot cast ufunc 'inv' input from dtype('O') to dtype('float64') with casting rule 'same_kind'

Which is not what was intended. It should just calculate the object according to the arithmetic and decimal should be able to compute the value itself. hpaulj's post provided an alternative method to cast the decimal object to mfp object of mpmath package

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) 

matrix(
[['1.0', '2.0'],
 ['3.0', '4.0']])

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) [0,0]

mpf('1.0')

and then perform the inverse in mpmath package.

mp.matrix( np.array([[Decimal(1),Decimal(2)],[Decimal(3),Decimal(4)]])) **(-1)

matrix(
[['-2.0', '1.0'],
 ['1.5', '-0.5']])

This could work, however, it lost the nice functionally of decimal package and involved large amount casting elements from mpmath to numpy and decimal objects. The mpf() object's computational speed is also significantly slower than the calculation Decimal() object's.

Is there an easy way to write or improve the code from the numpy package directly so that a np.inverse() could be used on decimal array? Is there any way to compute the matrix inverse with decimal object?