# Evgenii B. Rudnyi, http://MatrixProgramming.com

from numpy import *
from scipy.linalg import *

A = matrix([
	[1.0,  -0.026625,   0.00070892,   -1.8875E-05    ],
	[0.0,   1.000000,  -0.05325100,    0.00212670    ],
	[1.0,  -0.024859,   0.00061796,   -1.5362E-05    ],
	[0.0,   1.000000,  -0.04971800,    0.00185390    ]
	])

print svdvals(A)

b = array([0.99814 , -0.13086 , 0.98479 , 132.5085 ])

x = solve(A, b)

print norm(A*x - b)

def ndigit(x, n):
	e = 10**(round(math.log10(x)))
	b = round(x/e, n)
	return b*e

y = array([ndigit(i, 5) for i in x])

print x
print y

print norm(A*y - b)

for i in range(5, 16):
	y = array([ndigit(j, i) for j in x])
	print i, norm(A*y - b)