Path: csiph.com!usenet.pasdenom.info!weretis.net!feeder4.news.weretis.net!rt.uk.eu.org!newsfeed.xs4all.nl!newsfeed1.news.xs4all.nl!xs4all!post.news.xs4all.nl!not-for-mail
Date: Sat, 10 May 2014 08:48:17 -0700
From: Gary Herron <gary.herron@islandtraining.com>
User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:24.0) Gecko/20100101 Thunderbird/24.5.0
MIME-Version: 1.0
To: python-list@python.org
Subject: Re: How can this assert() ever trigger?
References: <536e44c1$0$27147$e4fe514c@dreader35.news.xs4all.nl>
In-Reply-To: <536e44c1$0$27147$e4fe514c@dreader35.news.xs4all.nl>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
Precedence: list
Newsgroups: comp.lang.python
Message-ID: <mailman.9852.1399737403.18130.python-list@python.org>
Lines: 80
NNTP-Posting-Host: 2001:888:2000:d::a6
Xref: csiph.com comp.lang.python:71252

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