Automatically round arithmetic operations to eight decimals

Question

I am doing some numerical analysis exercise where I need calculate solution of linear system using a specific algorithm. My answer differs from the answer of the book by some decimal places which I believe is due to rounding errors. Is there a way where I can automatically set arithmetic to round eight decimal places after each arithmetic operation? The following is my python code.

import numpy as np
A1 = [4, -1, 0, 0, -1, 4, -1, 0,\
     0, -1, 4, -1, 0, 0, -1, 4]
A1 = np.array(A1).reshape([4,4])
I = -np.identity(4)
O = np.zeros([4,4])

A = np.block([[A1, I, O, O],
             [I, A1, I, O],
             [O, I, A1, I],
             [O, O, I, A1]])

b = np.array([1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6])

def conj_solve(A, b, pre=False):
    n = len(A)
    
    C = np.identity(n)        
    
    if pre == True:
        for i in range(n):
            C[i, i] = np.sqrt(A[i, i])
            
    Ci = np.linalg.inv(C)
    Ct = np.transpose(Ci)
    
    x = np.zeros(n)
    r = b - np.matmul(A, x)
    w = np.matmul(Ci, r)
    v = np.matmul(Ct, w)
    alpha = np.dot(w, w)
    
    for i in range(MAX_ITER):
        if np.linalg.norm(v, np.infty) < TOL:            
            print(i+1, "steps")            
            print(x)
            print(r)
            return
        u = np.matmul(A, v)
        t = alpha/np.dot(v, u)
        x = x + t*v
        r = r - t*u
        w = np.matmul(Ci, r)
        beta = np.dot(w, w)
        
        if np.abs(beta) < TOL:
            if np.linalg.norm(r, np.infty) < TOL:
                print(i+1, "steps")
                print(x)
                print(r)
                return
        s = beta/alpha
        v = np.matmul(Ct, w) + s*v
        alpha = beta
    print("Max iteration exceeded")
    return x

MAX_ITER = 1000
TOL = 0.05

sol = conj_solve(A, b, pre=True)

Using this, I get 2.55516527 as first element of array which should be 2.55613420.

OR, is there a language/program where I can specify the precision of arithmetic?

Answer 1

Precision/rounding during the calculation is unlikely to be the issue.

To test this I ran the calculation with precisions that bracket the precision you are aiming for: once with np.float64, and once with np.float32. Here is a table of the printed results, their approximate decimal precision, and the result of the calculation (ie, the first printed array value).

numpy type       decimal places         result
-------------------------------------------------
np.float64               15              2.55516527
np.float32                6              2.5551653

Given that these are so much in agreement, I doubt an intermediate precision of 8 decimal places is going to give an answer that's not between these two results (ie, 2.55613420 that's off in the 4th digit).

---

This isn't part isn't part of my answer, but is a comment on using mpmath. The questioner suggested it in the comments, and it was my first thought too, so I ran a quick test to see if it behaved how I expected with low precision calculations. It didn't, so I abandoned it (but I'm not an expert with it).

Here's my test function, basically multiplying 1/N by N and 1/N repeatedly to emphasise the error in 1/N.

def precision_test(dps=100, N=19, t=mpmath.mpf):
    with mpmath.workdps(dps):  
        x = t(1)/t(N)
        print(x)
        y = x
        for i in range(10000):
            y *= x
            y *= N
        print(y)

This works as expected with, eg, np.float32:

precision_test(dps=2, N=3, t=np.float32)
# 0.33333334
# 0.3334327041164994

Note that the error has propagated into more significant digits, as expected.

But with mpmath, I could never get that to happen (testing with a range of dps and a various prime N values):

precision_test(dps=2, N=3)
# 0.33
# 0.33

Because of this test, I decided mpmath is not going to give normal results for low precision calculations.

TL;DR:
mpmath didn't behave how I expected at low precision so I abandoned it.