On 05/10/2014 08:24 AM, Albert van der Horst wrote:
I have the following code for calculating the determinant of a matrix. It works inasfar that it gives the same result as an octave program on a same matrix./ ---------------------------------------------------------------- def determinant( mat ): ''' Return the determinant of the n by n matrix mat i row j column Destroys mat ! ''' #print "getting determinat of", mat n=len(mat) nom = 1. if n == 1: return mat[0][0] lastr = mat.pop() jx=-1 for j in xrange(n): if lastr[j]: jx=j break if jx==-1: return 0. result = lastr[jx] assert(result<>0.) # Make column jx zero by subtracting a multiple of the last row. for i in xrange(n-1): pivot = mat[i][jx] if 0. == pivot: continue assert(result<>0.) nom *= result # Compenstate for multiplying a row. for j in xrange(n): mat[i][j] *= result for j in xrange(n): mat[i][j] -= pivot*lastr[j] # Remove colunm jx for i in xrange(n-1): x= mat[i].pop(jx) assert( x==0 ) if (n-1+jx)%2<>0: result = -result det = determinant( mat ) assert(nom<>0.) return result*det/nom /----------------------------------------- Now on some matrices the assert triggers, meaning that nom is zero. How can that ever happen? mon start out as 1. and gets multiplied with a number that is asserted to be not zero.
Easily due to *underflow* precision trouble. Your "result" may never be zero, but it can be very small. Take the product of many of such tiny values, and the result can be less then the smallest value representable by a float, at which point it becomes zero.
To see this clearly, try this Python code: >>> a = 1.0 >>> while a > 0: ... a = a*1.0e-50 ... print(a) ... 1e-50 1e-100 1e-150 1e-200 1e-250 1e-300 0.0 Gary Herron
Any hints appreciated. Groetjes Albert
-- https://mail.python.org/mailman/listinfo/python-list
