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