machineparameters

MB (11-2008) Machine parameters can be computed with the fortran routine DLAMCH, included in LAPACK. That routine computes some properties of the floating point arithmetic, such as the number of bits in the mantissa, or the epsilon number. It is very simple to do the same in Tcl and that is what I did in the machineparameters SNIT-based class available in the Tclrep project :

http://tclrep.cvs.sourceforge.net/viewvc/tclrep/modules/machineparameters/

One motivation for such a component is to be able to test the convergence of numerical algorithms, with respect to the precision available for floating-point numbers.

The following example creates a machineparameters object, computes the properties and displays it.

     set pp [machineparameters create %AUTO%]
     $pp configure -verbose 1
     $pp compute
     $pp print
     $pp destroy

This prints out :

     Machine parameters
     Epsilon : 1.11022302463e-16
     Beta : 2
     Rounding : chop
     Mantissa : 53
     Maximum exponent : 1023
     Minimum exponent : -1074
     Overflow threshold : 8.98846567431e+307
     Underflow threshold : 4.94065645841e-324

For example, the following is the algorithm which allows to compute epsilon, the largest double such that 1+e>1 is not true.

      set factor 2.
      set epsilon 0.5
      for {set i 0} {$i<$options(-maxiteration)} {incr i} {
        set epsilon [expr {$epsilon / $factor}]
        set inequality [expr {1.0+$epsilon>1.0}]
        if {$inequality==0} then {
          break
        }
      }

The following algorithm is explained in detail in "Algorithms to Reveal Properties of Floating-Point Arithmetic", Michael A. Malcolm and in "More on Algorithms that Reveal Properties of Floating Point Arithmetic Units", W. Morven Gentleman, Scott B. Marovich. It computes the last positive integer n which can be exactly represented with a floating point number, that is, such that (n+1)-n is not equal to 1. This is because all the digits in the mantissa have been used so that the integer which is to be represented requires a digit in the exponent.

       set firstnoninteger 2.
       for {set i 0} {$i < $options(-maxiteration)} {incr i} {
         set firstnoninteger [expr {2.*$firstnoninteger}]
         set one [expr {($firstnoninteger+1.)-$firstnoninteger}]
         if {$one!=1.} then {
           break
         }
       }

In addition to the Tcl implementation, more details are available in the two original articles or in the DLAMCH routine itself.

AM (18 november 2008) This could be a useful addition to Tcllib. Oh, by the way: determining "epsilon" (the smallest value such that 1+epsilon can be distinguished from 1, so a bit the opposite from the above algorithm) by a naive algorithm is challenging with languages like C or Fortran because of the extra precision programs written in these languages tend to use on many modern CPUs.

MB 2008-11-20 Thank you for your comments. I made further inquires in the DLAMCH routine and fixed some differences between the package and the DLAMCH results. Now the package better mimics the results of DLAMCH. Here are updated parameters on my machine. Notice the difference in the minimum exponent before underflow (from -1073 to -1021) and the underflow threshold (from 4.94065645841e-324 to 2.22507385851e-308).

     Epsilon : 1.11022302463e-016
     Basis : 2
     Rounding : chop
     Mantissa : 53
     Maximum exponent before overflow : 1024
     Minimum exponent before underflow : -1021
     Overflow threshold : 8.98846567431e+307
     Underflow threshold : 2.22507385851e-308

It is easy to experiment what happens when such extreme numbers are processed.

     % package require machineparameters
     0.1
     % set pp [machineparameters create %AUTO%]
     ::machineparameters1
     % $pp compute

Now get the maximal value before overflow and try to multiply it. It fails, as expected.

     % set vmax [$pp get -vmax]
     8.98846567431e+307
     % expr {$vmax * 2}
     floating-point value too large to represent

The minimum value before underflow does not have the same property.

     % set vmin [$pp get -vmin]
    2.22507385851e-308

The old value 4.94065645841e-324 had the desired property.

   % set a 4.94065645841e-324
   4.94065645841e-324
   % expr $a / 2
   0.0

I changed the value to reflect the algorithm in DLAMCH. The power is first computed as -1073, then updated to -1073 - 1 + 53, where 53 is the size of the mantissa.

One can also check that epsilon has the property that 1 + epsilon is not greater than 1.

     % set eps [$pp get -epsilon]
     1.11022302463e-016
     % expr {1+$eps>1}
     0

One can also compute epsilon directly with the function pow, by using the fact that the mantissa uses 53 bits :

     % expr pow(2,-53)
     1.11022302463e-016

This is not necessary, but it is cleaner to do some cleanup.

     % $pp destroy