1# machineparameters.tcl --
2#   Compute double precision machine parameters.
3#
4# Description
5#   This the Tcl equivalent of the DLAMCH LAPCK function.
6#   In floating point systems, a floating point number is represented
7#   by
8#   x = +/- d1 d2 ... dt basis^e
9#   where digits satisfy
10#   0 <= di <= basis - 1, i = 1, t
11#   with the convention :
12#   - t is the size of the mantissa
13#   - basis is the basis (the "radix")
14#
15# References
16#
17#   "Algorithms to Reveal Properties of Floating-Point Arithmetic"
18#   Michael A. Malcolm
19#   Stanford University
20#   Communications of the ACM
21#   Volume 15 ,  Issue 11  (November 1972)
22#   Pages: 949 - 951
23#
24#   "More on Algorithms that Reveal Properties of Floating
25#   Point Arithmetic Units"
26#   W. Morven Gentleman, University of Waterloo
27#   Scott B. Marovich, Purdue University
28#   Communications of the ACM
29#   Volume 17 ,  Issue 5  (May 1974)
30#   Pages: 276 - 277
31#
32# Example
33#
34#   In the following example, one compute the parameters of a desktop
35#   under Linux with the following Tcl 8.4.19 properties :
36#
37#% parray tcl_platform
38#tcl_platform(byteOrder) = littleEndian
39#tcl_platform(machine)   = i686
40#tcl_platform(os)        = Linux
41#tcl_platform(osVersion) = 2.6.24-19-generic
42#tcl_platform(platform)  = unix
43#tcl_platform(tip,268)   = 1
44#tcl_platform(tip,280)   = 1
45#tcl_platform(user)      = <username>
46#tcl_platform(wordSize)  = 4
47#
48#   The following example creates a machineparameters object,
49#   computes the properties and displays it.
50#
51#     set pp [machineparameters create %AUTO%]
52#     $pp compute
53#     $pp print
54#     $pp destroy
55#
56#   This prints out :
57#
58#     Machine parameters
59#     Epsilon : 1.11022302463e-16
60#     Beta : 2
61#     Rounding : proper
62#     Mantissa : 53
63#     Maximum exponent : 1024
64#     Minimum exponent : -1021
65#     Overflow threshold : 8.98846567431e+307
66#     Underflow threshold : 2.22507385851e-308
67#
68#   That compares well with the results produced by Lapack 3.1.1 :
69#
70#     Epsilon                      =   1.11022302462515654E-016
71#     Safe minimum                 =   2.22507385850720138E-308
72#     Base                         =    2.0000000000000000
73#     Precision                    =   2.22044604925031308E-016
74#     Number of digits in mantissa =    53.000000000000000
75#     Rounding mode                =   1.00000000000000000
76#     Minimum exponent             =   -1021.0000000000000
77#     Underflow threshold          =   2.22507385850720138E-308
78#     Largest exponent             =    1024.0000000000000
79#     Overflow threshold           =   1.79769313486231571E+308
80#     Reciprocal of safe minimum   =   4.49423283715578977E+307
81#
82# Copyright 2008 Michael Baudin
83#
84package require snit
85package provide math::machineparameters 0.1
86
87snit::type machineparameters {
88  # Epsilon is the smallest value so that 1+epsilon>1 is false
89  variable epsilon 0
90  # basis is the basis of the floating-point representation.
91  # basis is usually 2, i.e. binary representation (for example IEEE 754 machines),
92  # but some machines (like HP calculators for example) uses 10, or 16, etc...
93  variable basis 0
94  # The rounding mode used on the machine.
95  # The rounding occurs when more than t digits would be required to
96  # represent the number.
97  # Two modes can be determined with the current system :
98  # "chop" means than only t digits are kept, no matter the value of the number
99  # "proper" means that another rounding mode is used, be it "round to nearest",
100  #   "round up", "round down".
101  variable rounding ""
102  # the size of the mantissa
103  variable mantissa 0
104  # The first non-integer is A = 2^m with m is the
105  # smallest positive integer so that fl(A+1)=A
106  variable firstnoninteger 0
107  # Maximum number of iterations in loops
108  option -maxiteration 10000
109  # Set to 1 to enable verbose logging
110  option -verbose -default 0
111  # The largest positive exponent before overflow occurs
112  variable exponentmax 0
113  # The largest negative exponent before (gradual) underflow occurs
114  variable exponentmin 0
115  # Largest positive value before overflow occurs
116  variable vmax
117  # Largest negative value before (gradual) underflow occurs
118  variable vmin
119  #
120  # compute --
121  #   Computes the machine parameters.
122  #
123  method compute {} {
124    $self log "compute"
125    $self computeepsilon
126    $self computefirstnoninteger
127    $self computebasis
128    $self computerounding
129    $self computemantissa
130    $self computeemax
131    $self computeemin
132    return ""
133  }
134  #
135  # computeepsilon --
136  #   Find epsilon the minimum value for which 1.0 + epsilon > 1.0
137  #
138  method computeepsilon {} {
139    $self log "computeepsilon"
140    set factor 2.
141    set epsilon 0.5
142    for {set i 0} {$i<$options(-maxiteration)} {incr i} {
143      $self log "$i/$options(-maxiteration) : $epsilon"
144      set epsilon [expr {$epsilon / $factor}]
145      set inequality [expr {1.0+$epsilon>1.0}]
146      if {$inequality==0} then {
147        break
148      }
149    }
150    $self log "epsilon : $epsilon (after $i loops)"
151    return ""
152  }
153  #
154  # computefirstnoninteger --
155  #   Compute the first positive non-integer real.
156  #   It is the smallest a such that (a+1)-a is different from 1
157  #
158  method computefirstnoninteger {} {
159    $self log "computefirstnoninteger"
160    set firstnoninteger 2.
161    for {set i 0} {$i < $options(-maxiteration)} {incr i} {
162      $self log "$i/$options(-maxiteration) : $firstnoninteger"
163      set firstnoninteger [expr {2.*$firstnoninteger}]
164      set one [expr {($firstnoninteger+1.)-$firstnoninteger}]
165      if {$one!=1.} then {
166        break
167      }
168    }
169    $self log "Found firstnoninteger : $firstnoninteger"
170    return ""
171  }
172  #
173  # computebasis --
174  #   Compute the basis (basis)
175  #
176  method computebasis {} {
177    $self log "computebasis"
178    #
179    # Compute b where b is the smallest real so that fl(a+b)> a,
180    # where a is the first non integer.
181    # Note :
182    #  With floating point numbers, a+1==a !
183    #  b is denoted by "B" in Malcolm's algorithm
184    #
185    set b 1
186    for {set i 0} {$i < $options(-maxiteration)} {incr i} {
187      $self log "$i/$options(-maxiteration) : $b"
188      set basis [expr {int(($firstnoninteger+$b)-$firstnoninteger)}]
189      if {$basis!=0.} then {
190        break
191      }
192      incr b
193    }
194    $self log "Found basis : $basis"
195    return ""
196  }
197  #
198  # computerounding --
199  #   Compute the rounding mode.
200  # Note:
201  #   This corresponds to DLAMCH implementation (DLAMC1 exactly).
202  #
203  method computerounding {} {
204    $self log "computerounding"
205    # Now determine whether rounding or chopping occurs,  by adding a
206    # bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a (=firstnoninteger).
207    set F [expr {$basis/2.0 - $basis/100.0}]
208    set C [expr {$F + $firstnoninteger}]
209    if {$C==$firstnoninteger} then {
210      set rounding "proper"
211    } else {
212      set rounding "chop"
213    }
214    set F [expr {$basis/2.0 + $basis/100.0}]
215    set C [expr {$F + $firstnoninteger}]
216    if {$rounding=="proper" && $C==$firstnoninteger} then {
217      set rounding "chop"
218    }
219    $self log "Found rounding : $rounding"
220    return ""
221  }
222  #
223  # computemantissa --
224  #   Compute the mantissa size
225  #
226  method computemantissa {} {
227    $self log "computemantissa"
228    set a 1.
229    set mantissa 0
230    for {set i 0} {$i < $options(-maxiteration)} {incr i} {
231      incr mantissa
232      $self log "$i/$options(-maxiteration) : $mantissa"
233      set a [expr {$a * double($basis)}]
234      set one [expr {($a+1)-$a}]
235      if {$one!=1.} then {
236        break
237      }
238    }
239    $self log "Found mantissa : $mantissa"
240    return ""
241  }
242  #
243  # computeemax --
244  #   Compute the maximum exponent before overflow
245  #
246  method computeemax {} {
247    $self log "computeemax"
248    set vmax 1.
249    set exponentmax 1
250    for {set i 0} {$i < $options(-maxiteration)} {incr i} {
251      $self log "Iteration #$i , exponentmax = $exponentmax, vmax = $vmax"
252      incr exponentmax
253      # Condition #1 : no exception is generated
254      set errflag [catch {
255        set new [expr {$vmax * $basis}]
256      }]
257      if {$errflag!=0} then {
258        break
259      }
260      # Condition #2 : one can recover the original number
261      if {$new / $basis != $vmax} then {
262        break
263      }
264      set vmax $new
265    }
266    incr exponentmax -1
267    $self log "Exponent maximum : $exponentmax"
268    $self log "Value maximum : $vmax"
269    return ""
270  }
271  #
272  # computeemin --
273  #   Compute the minimum exponent before underflow
274  #
275  method computeemin {} {
276    $self log "computeemin"
277    set vmin 1.
278    set exponentmin 1
279    for {set i 0} {$i < $options(-maxiteration)} {incr i} {
280      $self log "Iteration #$i , exponentmin = $exponentmin, vmin = $vmin"
281      incr exponentmin -1
282      # Condition #1 : no exception is generated
283      set errflag [catch {
284        set new [expr {$vmin / $basis}]
285      }]
286      if {$errflag!=0} then {
287        break
288      }
289      # Condition #2 : one can recover the original number
290      if {$new * $basis != $vmin} then {
291        break
292      }
293      set vmin $new
294    }
295    incr exponentmin +1
296    # See in DMALCH.f, DLAMC2 relative to IEEE machines.
297    # TODO : what happens on non-IEEE machine ?
298    set exponentmin [expr {$exponentmin - 1 + $mantissa}]
299    set vmin [expr {$vmin * pow($basis,$mantissa-1)}]
300    $self log "Exponent minimum : $exponentmin"
301    $self log "Value minimum : $vmin"
302    return ""
303  }
304  #
305  # log --
306  #   Puts the given message on standard output.
307  #
308  method log {msg} {
309    if {$options(-verbose)==1} then {
310      puts "(mp) $msg"
311    }
312    return ""
313  }
314  #
315  # get --
316  #   Return value for key
317  #
318  method get {key} {
319    $self log "get $key"
320    switch -- $key {
321      -epsilon {
322        set result $epsilon
323      }
324      -rounding {
325        set result $rounding
326      }
327      -basis {
328        set result $basis
329      }
330      -mantissa {
331        set result $mantissa
332      }
333      -exponentmax {
334        set result $exponentmax
335      }
336      -exponentmin {
337        set result $exponentmin
338      }
339      -vmax {
340        set result $vmax
341      }
342      -vmin {
343        set result $vmin
344      }
345      default {
346        error "Unknown key $key"
347      }
348    }
349    return $result
350  }
351  #
352  # print --
353  #   Print machine parameters on standard output
354  #
355  method print {} {
356    set str [$self tostring]
357    puts "$str"
358    return ""
359  }
360  #
361  # tostring --
362  #   Return a report for machine parameters
363  #
364  method tostring {} {
365    set str ""
366    append str "Machine parameters\n"
367    append str  "Epsilon : $epsilon\n"
368    append str "Basis : $basis\n"
369    append str "Rounding : $rounding\n"
370    append str "Mantissa : $mantissa\n"
371    append str "Maximum exponent before overflow : $exponentmax\n"
372    append str "Minimum exponent before underflow : $exponentmin\n"
373    append str "Overflow threshold : $vmax\n"
374    append str "Underflow threshold : $vmin\n"
375    return $str
376  }
377}
378