1######################################################################## 2# BigFloat for Tcl 3# Copyright (C) 2003-2005 ARNOLD Stephane 4# It is published with the terms of tcllib's BSD-style license. 5# See the file named license.terms. 6######################################################################## 7 8package require Tcl 8.5 9 10# this line helps when I want to source this file again and again 11catch {namespace delete ::math::bigfloat} 12 13# private namespace 14# this software works only with Tcl v8.4 and higher 15# it is using the package math::bignum 16namespace eval ::math::bigfloat { 17 # cached constants 18 # ln(2) with arbitrary precision 19 variable Log2 20 # Pi with arb. precision 21 variable Pi 22 variable _pi0 23} 24 25 26 27 28################################################################################ 29# procedures that handle floating-point numbers 30# these procedures are sorted by name (after eventually removing the underscores) 31# 32# BigFloats are internally represented as a list : 33# {"F" Mantissa Exponent Delta} where "F" is a character which determins 34# the datatype, Mantissa and Delta are two big integers and Exponent another integer. 35# 36# The BigFloat value equals to (Mantissa +/- Delta)*2^Exponent 37# So the internal representation is binary, but trying to get as close as possible to 38# the decimal one when converted to a string. 39# When calling [fromstr], the Delta parameter is set to the value of 1 at the position 40# of the last decimal digit. 41# Example : 1.50 belongs to [1.49,1.51], but internally Delta may not equal to 1. 42# Because of the binary representation, it is between 1 and 1+(2^-15). 43# 44# So Mantissa and Delta are not limited in size, but in practice Delta is kept under 45# 2^32 by the 'normalize' procedure, to avoid a never-ended growth of memory used. 46# Indeed, when you perform some computations, the Delta parameter (which represent 47# the uncertainty on the value of the Mantissa) may increase. 48# Exponent, as an integer, is limited to 32 bits, and this limit seems fair. 49# The exponent is indeed involved in logarithmic computations, so it may be 50# a mistake to give it a too large value. 51 52# Retrieving the parameters of a BigFloat is often done with that command : 53# foreach {dummy int exp delta} $bigfloat {break} 54# (dummy is not used, it is just used to get the "F" marker). 55# The isInt, isFloat, checkNumber and checkFloat procedures are used 56# to check data types 57# 58# Taylor development are often used to compute the analysis functions (like exp(),log()...) 59# To learn how it is done in practice, take a look at ::math::bigfloat::_asin 60# While doing computation on Mantissas, we do not care about the last digit, 61# because if we compute correctly Deltas, the digits that remain will be exact. 62################################################################################ 63 64 65################################################################################ 66# returns the absolute value 67################################################################################ 68proc ::math::bigfloat::abs {number} { 69 checkNumber $number 70 if {[isInt $number]} { 71 # set sign to positive for a BigInt 72 return [expr {abs($number)}] 73 } 74 # set sign to positive for a BigFloat into the Mantissa (index 1) 75 lset number 1 [expr {abs([lindex $number 1])}] 76 return $number 77} 78 79 80################################################################################ 81# arccosinus of a BigFloat 82################################################################################ 83proc ::math::bigfloat::acos {x} { 84 # handy proc for checking datatype 85 checkFloat $x 86 foreach {dummy entier exp delta} $x {break} 87 set precision [expr {($exp<0)?(-$exp):1}] 88 # acos(0.0)=Pi/2 89 # 26/07/2005 : changed precision from decimal to binary 90 # with the second parameter of pi command 91 set piOverTwo [floatRShift [pi $precision 1]] 92 if {[iszero $x]} { 93 # $x is too close to zero -> acos(0)=PI/2 94 return $piOverTwo 95 } 96 # acos(-x)= Pi/2 + asin(x) 97 if {$entier<0} { 98 return [add $piOverTwo [asin [abs $x]]] 99 } 100 # we always use _asin to compute the result 101 # but as it is a Taylor development, the value given to [_asin] 102 # has to be a bit smaller than 1 ; by using that trick : acos(x)=asin(sqrt(1-x^2)) 103 # we can limit the entry of the Taylor development below 1/sqrt(2) 104 if {[compare $x [fromstr 0.7071]]>0} { 105 # x > sqrt(2)/2 : trying to make _asin converge quickly 106 # creating 0 and 1 with the same precision as the entry 107 set fzero [list F 0 -$precision 1] 108 # 1.000 with $precision zeros 109 set fone [list F [expr {1<<$precision}] -$precision 1] 110 # when $x is close to 1 (acos(1.0)=0.0) 111 if {[equal $fone $x]} { 112 return $fzero 113 } 114 if {[compare $fone $x]<0} { 115 # the behavior assumed because acos(x) is not defined 116 # when |x|>1 117 error "acos on a number greater than 1" 118 } 119 # acos(x) = asin(sqrt(1 - x^2)) 120 # since 1 - cos(x)^2 = sin(x)^2 121 # x> sqrt(2)/2 so x^2 > 1/2 so 1-x^2<1/2 122 set x [sqrt [sub $fone [mul $x $x]]] 123 # the parameter named x is smaller than sqrt(2)/2 124 return [_asin $x] 125 } 126 # acos(x) = Pi/2 - asin(x) 127 # x<sqrt(2)/2 here too 128 return [sub $piOverTwo [_asin $x]] 129} 130 131 132################################################################################ 133# returns A + B 134################################################################################ 135proc ::math::bigfloat::add {a b} { 136 checkNumber $a 137 checkNumber $b 138 if {[isInt $a]} { 139 if {[isInt $b]} { 140 # intAdd adds two BigInts 141 return [incr a $b] 142 } 143 # adds the BigInt a to the BigFloat b 144 return [addInt2Float $b $a] 145 } 146 if {[isInt $b]} { 147 # ... and vice-versa 148 return [addInt2Float $a $b] 149 } 150 # retrieving parameters from A and B 151 foreach {dummy integerA expA deltaA} $a {break} 152 foreach {dummy integerB expB deltaB} $b {break} 153 if {$expA<$expB} { 154 foreach {dummy integerA expA deltaA} $b {break} 155 foreach {dummy integerB expB deltaB} $a {break} 156 } 157 # when we add two numbers which have different digit numbers (after the dot) 158 # for example : 1.0 and 0.00001 159 # We promote the one with the less number of digits (1.0) to the same level as 160 # the other : so 1.00000. 161 # that is why we shift left the number which has the greater exponent 162 # But we do not forget the Delta parameter, which is lshift'ed too. 163 if {$expA>$expB} { 164 set diff [expr {$expA-$expB}] 165 set integerA [expr {$integerA<<$diff}] 166 set deltaA [expr {$deltaA<<$diff}] 167 incr integerA $integerB 168 incr deltaA $deltaB 169 return [normalize [list F $integerA $expB $deltaA]] 170 } elseif {$expA==$expB} { 171 # nothing to shift left 172 return [normalize [list F [incr integerA $integerB] $expA [incr deltaA $deltaB]]] 173 } else { 174 error "internal error" 175 } 176} 177 178################################################################################ 179# returns the sum A(BigFloat) + B(BigInt) 180# the greatest advantage of this method is that the uncertainty 181# of the result remains unchanged, in respect to the entry's uncertainty (deltaA) 182################################################################################ 183proc ::math::bigfloat::addInt2Float {a b} { 184 # type checking 185 checkFloat $a 186 if {![isInt $b]} { 187 error "second argument is not an integer" 188 } 189 # retrieving data from $a 190 foreach {dummy integerA expA deltaA} $a {break} 191 # to add an int to a BigFloat,... 192 if {$expA>0} { 193 # we have to put the integer integerA 194 # to the level of zero exponent : 1e8 --> 100000000e0 195 set shift $expA 196 set integerA [expr {($integerA<<$shift)+$b}] 197 set deltaA [expr {$deltaA<<$shift}] 198 # we have to normalize, because we have shifted the mantissa 199 # and the uncertainty left 200 return [normalize [list F $integerA 0 $deltaA]] 201 } elseif {$expA==0} { 202 # integerA is already at integer level : float=(integerA)e0 203 return [normalize [list F [incr integerA $b] \ 204 0 $deltaA]] 205 } else { 206 # here we have something like 234e-2 + 3 207 # we have to shift the integer left by the exponent |$expA| 208 incr integerA [expr {$b<<(-$expA)}] 209 return [normalize [list F $integerA $expA $deltaA]] 210 } 211} 212 213 214################################################################################ 215# arcsinus of a BigFloat 216################################################################################ 217proc ::math::bigfloat::asin {x} { 218 # type checking 219 checkFloat $x 220 foreach {dummy entier exp delta} $x {break} 221 if {$exp>-1} { 222 error "not enough precision on input (asin)" 223 } 224 set precision [expr {-$exp}] 225 # when x=0, return 0 at the same precision as the input was 226 if {[iszero $x]} { 227 return [list F 0 -$precision 1] 228 } 229 # asin(-x)=-asin(x) 230 if {$entier<0} { 231 return [opp [asin [abs $x]]] 232 } 233 # 26/07/2005 : changed precision from decimal to binary 234 set piOverTwo [floatRShift [pi $precision 1]] 235 # now a little trick : asin(x)=Pi/2-asin(sqrt(1-x^2)) 236 # so we can limit the entry of the Taylor development 237 # to 1/sqrt(2)~0.7071 238 # the comparison is : if x>0.7071 then ... 239 if {[compare $x [fromstr 0.7071]]>0} { 240 set fone [list F [expr {1<<$precision}] -$precision 1] 241 # asin(1)=Pi/2 (with the same precision as the entry has) 242 if {[equal $fone $x]} { 243 return $piOverTwo 244 } 245 if {[compare $x $fone]>0} { 246 error "asin on a number greater than 1" 247 } 248 # asin(x)=Pi/2-asin(sqrt(1-x^2)) 249 set x [sqrt [sub $fone [mul $x $x]]] 250 return [sub $piOverTwo [_asin $x]] 251 } 252 return [normalize [_asin $x]] 253} 254 255################################################################################ 256# _asin : arcsinus of numbers between 0 and +1 257################################################################################ 258proc ::math::bigfloat::_asin {x} { 259 # Taylor development 260 # asin(x)=x + 1/2 x^3/3 + 3/2.4 x^5/5 + 3.5/2.4.6 x^7/7 + ... 261 # into this iterative form : 262 # asin(x)=x * (1 + 1/2 * x^2 * (1/3 + 3/4 *x^2 * (... 263 # ...* (1/(2n-1) + (2n-1)/2n * x^2 / (2n+1))...))) 264 # we show how is really computed the development : 265 # we don't need to set a var with x^n or a product of integers 266 # all we need is : x^2, 2n-1, 2n, 2n+1 and a few variables 267 foreach {dummy mantissa exp delta} $x {break} 268 set precision [expr {-$exp}] 269 if {$precision+1<[bits $mantissa]} { 270 error "sinus greater than 1" 271 } 272 # precision is the number of after-dot digits 273 set result $mantissa 274 set delta_final $delta 275 # resultat is the final result, and delta_final 276 # will contain the uncertainty of the result 277 # square is the square of the mantissa 278 set square [expr {$mantissa*$mantissa>>$precision}] 279 # dt is the uncertainty of Mantissa 280 set dt [expr {$mantissa*$delta>>($precision-1)}] 281 incr dt 282 set num 1 283 # two will be used into the loop 284 set i 3 285 set denom 2 286 # the nth factor equals : $num/$denom* $mantissa/$i 287 set delta [expr {$delta*$square + $dt*($delta+$mantissa)}] 288 set delta [expr {($delta*$num)/ $denom >>$precision}] 289 incr delta 290 # we do not multiply the Mantissa by $num right now because it is 1 ! 291 # but we have Mantissa=$x 292 # and we want Mantissa*$x^2 * $num / $denom / $i 293 set mantissa [expr {($mantissa*$square>>$precision)/$denom}] 294 # do not forget the modified Taylor development : 295 # asin(x)=x * (1 + 1/2*x^2*(1/3 + 3/4*x^2*(...*(1/(2n-1) + (2n-1)/2n*x^2/(2n+1))...))) 296 # all we need is : x^2, 2n-1, 2n, 2n+1 and a few variables 297 # $num=2n-1 $denom=2n $square=x^2 and $i=2n+1 298 set mantissa_temp [expr {$mantissa/$i}] 299 set delta_temp [expr {1+$delta/$i}] 300 # when the Mantissa increment is smaller than the Delta increment, 301 # we would not get much precision by continuing the development 302 while {$mantissa_temp!=0} { 303 # Mantissa = Mantissa * $num/$denom * $square 304 # Add Mantissa/$i, which is stored in $mantissa_temp, to the result 305 incr result $mantissa_temp 306 incr delta_final $delta_temp 307 # here we have $two instead of [fromstr 2] (optimization) 308 # num=num+2,i=i+2,denom=denom+2 309 # because num=2n-1 denom=2n and i=2n+1 310 incr num 2 311 incr i 2 312 incr denom 2 313 # computes precisly the future Delta parameter 314 set delta [expr {$delta*$square+$dt*($delta+$mantissa)}] 315 set delta [expr {($delta*$num)/$denom>>$precision}] 316 incr delta 317 set mantissa [expr {$mantissa*$square>>$precision}] 318 set mantissa [expr {($mantissa*$num)/$denom}] 319 set mantissa_temp [expr {$mantissa/$i}] 320 set delta_temp [expr {1+$delta/$i}] 321 } 322 return [normalize [list F $result $exp $delta_final]] 323} 324 325################################################################################ 326# arctangent : returns atan(x) 327################################################################################ 328proc ::math::bigfloat::atan {x} { 329 checkFloat $x 330 foreach {dummy mantissa exp delta} $x {break} 331 if {$exp>=0} { 332 error "not enough precision to compute atan" 333 } 334 set precision [expr {-$exp}] 335 # atan(0)=0 336 if {[iszero $x]} { 337 return [list F 0 -$precision $delta] 338 } 339 # atan(-x)=-atan(x) 340 if {$mantissa<0} { 341 return [opp [atan [abs $x]]] 342 } 343 # now x is strictly positive 344 # at this moment, we are trying to limit |x| to a fair acceptable number 345 # to ensure that Taylor development will converge quickly 346 set float1 [list F [expr {1<<$precision}] -$precision 1] 347 if {[compare $float1 $x]<0} { 348 # compare x to 2.4142 349 if {[compare $x [fromstr 2.4142]]<0} { 350 # atan(x)=Pi/4 + atan((x-1)/(x+1)) 351 # as 1<x<2.4142 : (x-1)/(x+1)=1-2/(x+1) belongs to 352 # the range : ]0,1-2/3.414[ 353 # that equals ]0,0.414[ 354 set pi_sur_quatre [floatRShift [pi $precision 1] 2] 355 return [add $pi_sur_quatre [atan \ 356 [div [sub $x $float1] [add $x $float1]]]] 357 } 358 # atan(x)=Pi/2-atan(1/x) 359 # 1/x < 1/2.414 so the argument is lower than 0.414 360 set pi_over_two [floatRShift [pi $precision 1]] 361 return [sub $pi_over_two [atan [div $float1 $x]]] 362 } 363 if {[compare $x [fromstr 0.4142]]>0} { 364 # atan(x)=Pi/4 + atan((x-1)/(x+1)) 365 # x>0.420 so (x-1)/(x+1)=1 - 2/(x+1) > 1-2/1.414 366 # > -0.414 367 # x<1 so (x-1)/(x+1)<0 368 set pi_sur_quatre [floatRShift [pi $precision 1] 2] 369 return [add $pi_sur_quatre [atan \ 370 [div [sub $x $float1] [add $x $float1]]]] 371 } 372 # precision increment : to have less uncertainty 373 # we add a little more precision so that the result would be more accurate 374 # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) 375 # when we have n steps in Taylor development : the nth term is : 376 # x^(2n-1)/(2n-1) 377 # and the loss of precision is of 2n (n sums and n divisions) 378 # this command is called with x<sqrt(2)-1 379 # if we add an increment to the precision, say n: 380 # (sqrt(2)-1)^(2n-1)/(2n-1) has to be lower than 2^(-precision-n-1) 381 # (2n-1)*log(sqrt(2)-1)-log(2n-1)<-(precision+n+1)*log(2) 382 # 2n(log(sqrt(2)-1)-log(sqrt(2)))<-(precision-1)*log(2)+log(2n-1)+log(sqrt(2)-1) 383 # 2n*log(1-1/sqrt(2))<-(precision-1)*log(2)+log(2n-1)+log(2)/2 384 # 2n/sqrt(2)>(precision-3/2)*log(2)-log(2n-1) 385 # hence log(2n-1)<2n-1 386 # n*sqrt(2)>(precision-1.5)*log(2)+1-2n 387 # n*(sqrt(2)+2)>(precision-1.5)*log(2)+1 388 set n [expr {int((log(2)*($precision-1.5)+1)/(sqrt(2)+2)+1)}] 389 incr precision $n 390 set mantissa [expr {$mantissa<<$n}] 391 set delta [expr {$delta<<$n}] 392 # end of adding precision increment 393 # now computing Taylor development : 394 # atan(x)=x - x^3/3 + x^5/5 - x^7/7 ... + (-1)^n*x^(2n+1)/(2n+1) 395 # atan(x)=x * (1 - x^2 * (1/3 - x^2 * (1/5 - x^2 * (...*(1/(2n-1) - x^2 / (2n+1))...)))) 396 # what do we need to compute this ? 397 # x^2 ($square), 2n+1 ($divider), $result, the nth term of the development ($t) 398 # and the nth term multiplied by 2n+1 ($temp) 399 # then we do this (with care keeping as much precision as possible): 400 # while ($t <>0) : 401 # $result=$result+$t 402 # $temp=$temp * $square 403 # $divider = $divider+2 404 # $t=$temp/$divider 405 # end-while 406 set result $mantissa 407 set delta_end $delta 408 # we store the square of the integer (mantissa) 409 # Delta of Mantissa^2 = Delta * 2 = Delta << 1 410 set delta_square [expr {$delta<<1}] 411 set square [expr {$mantissa*$mantissa>>$precision}] 412 # the (2n+1) divider 413 set divider 3 414 # computing precisely the uncertainty 415 set delta [expr {1+($delta_square*$mantissa+$delta*$square>>$precision)}] 416 # temp contains (-1)^n*x^(2n+1) 417 set temp [expr {-$mantissa*$square>>$precision}] 418 set t [expr {$temp/$divider}] 419 set dt [expr {1+$delta/$divider}] 420 while {$t!=0} { 421 incr result $t 422 incr delta_end $dt 423 incr divider 2 424 set delta [expr {1+($delta_square*abs($temp)+$delta*($delta_square+$square)>>$precision)}] 425 set temp [expr {-$temp*$square>>$precision}] 426 set t [expr {$temp/$divider}] 427 set dt [expr {1+$delta/$divider}] 428 } 429 # we have to normalize because the uncertainty might be greater than 2**16 430 # moreover it is the most often case 431 return [normalize [list F $result [expr {$exp-$n}] $delta_end]] 432} 433 434 435################################################################################ 436# compute atan(1/integer) at a given precision 437# this proc is only used to compute Pi 438# it is using the same Taylor development as [atan] 439################################################################################ 440proc ::math::bigfloat::_atanfract {integer precision} { 441 # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) 442 # when we have n steps in Taylor development : the nth term is : 443 # 1/denom^(2n+1)/(2n+1) 444 # and the loss of precision is of 2n (n sums and n divisions) 445 # this command is called with integer>=5 446 # 447 # We do not want to compute the Delta parameter, so we just 448 # can increment precision (with lshift) in order for the result to be precise. 449 # Remember : we compute atan2(1,$integer) with $precision bits 450 # $integer has no Delta parameter as it is a BigInt, of course, so 451 # theorically we could compute *any* number of digits. 452 # 453 # if we add an increment to the precision, say n: 454 # (1/5)^(2n-1)/(2n-1) has to be lower than (1/2)^(precision+n-1) 455 # Calculus : 456 # log(left term) < log(right term) 457 # log(1/left term) > log(1/right term) 458 # (2n-1)*log(5)+log(2n-1)>(precision+n-1)*log(2) 459 # n(2log(5)-log(2))>(precision-1)*log(2)-log(2n-1)+log(5) 460 # -log(2n-1)>-(2n-1) 461 # n(2log(5)-log(2)+2)>(precision-1)*log(2)+1+log(5) 462 set n [expr {int((($precision-1)*log(2)+1+log(5))/(2*log(5)-log(2)+2)+1)}] 463 incr precision $n 464 # first term of the development : 1/integer 465 set a [expr {(1<<$precision)/$integer}] 466 # 's' will contain the result 467 set s $a 468 # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) 469 # equals x (1 - x^2 * (1/3 + x^2 * (... * (1/(2n-3) + (-1)^(n+1) * x^2 / (2n-1))...))) 470 # all we need to store is : 2n-1 ($denom), x^(2n+1) and x^2 ($square) and two results : 471 # - the nth term => $u 472 # - the nth term * (2n-1) => $t 473 # + of course, the result $s 474 set square [expr {$integer*$integer}] 475 set denom 3 476 # $t is (-1)^n*x^(2n+1) 477 set t [expr {-$a/$square}] 478 set u [expr {$t/$denom}] 479 # we break the loop when the current term of the development is null 480 while {$u!=0} { 481 incr s $u 482 # denominator= (2n+1) 483 incr denom 2 484 # div $t by x^2 485 set t [expr {-$t/$square}] 486 set u [expr {$t/$denom}] 487 } 488 # go back to the initial precision 489 return [expr {$s>>$n}] 490} 491 492# 493# bits : computes the number of bits of an integer, approx. 494# 495proc ::math::bigfloat::bits {int} { 496 set l [string length [set int [expr {abs($int)}]]] 497 # int<10**l -> log_2(int)=l*log_2(10) 498 set l [expr {int($l*log(10)/log(2))+1}] 499 if {$int>>$l!=0} { 500 error "bad result: $l bits" 501 } 502 while {($int>>($l-1))==0} { 503 incr l -1 504 } 505 return $l 506} 507 508################################################################################ 509# returns the integer part of a BigFloat, as a BigInt 510# the result is the same one you would have 511# if you had called [expr {ceil($x)}] 512################################################################################ 513proc ::math::bigfloat::ceil {number} { 514 checkFloat $number 515 set number [normalize $number] 516 if {[iszero $number]} { 517 return 0 518 } 519 foreach {dummy integer exp delta} $number {break} 520 if {$exp>=0} { 521 error "not enough precision to perform rounding (ceil)" 522 } 523 # saving the sign ... 524 set sign [expr {$integer<0}] 525 set integer [expr {abs($integer)}] 526 # integer part 527 set try [expr {$integer>>(-$exp)}] 528 if {$sign} { 529 return [opp $try] 530 } 531 # fractional part 532 if {($try<<(-$exp))!=$integer} { 533 return [incr try] 534 } 535 return $try 536} 537 538 539################################################################################ 540# checks each variable to be a BigFloat 541# arguments : each argument is the name of a variable to be checked 542################################################################################ 543proc ::math::bigfloat::checkFloat {number} { 544 if {![isFloat $number]} { 545 error "BigFloat expected" 546 } 547} 548 549################################################################################ 550# checks if each number is either a BigFloat or a BigInt 551# arguments : each argument is the name of a variable to be checked 552################################################################################ 553proc ::math::bigfloat::checkNumber {x} { 554 if {![isFloat $x] && ![isInt $x]} { 555 error "input is not an integer, nor a BigFloat" 556 } 557} 558 559 560################################################################################ 561# returns 0 if A and B are equal, else returns 1 or -1 562# accordingly to the sign of (A - B) 563################################################################################ 564proc ::math::bigfloat::compare {a b} { 565 if {[isInt $a] && [isInt $b]} { 566 set diff [expr {$a-$b}] 567 if {$diff>0} {return 1} elseif {$diff<0} {return -1} 568 return 0 569 } 570 checkFloat $a 571 checkFloat $b 572 if {[equal $a $b]} {return 0} 573 if {[lindex [sub $a $b] 1]<0} {return -1} 574 return 1 575} 576 577 578 579 580################################################################################ 581# gets cos(x) 582# throws an error if there is not enough precision on the input 583################################################################################ 584proc ::math::bigfloat::cos {x} { 585 checkFloat $x 586 foreach {dummy integer exp delta} $x {break} 587 if {$exp>-2} { 588 error "not enough precision on floating-point number" 589 } 590 set precision [expr {-$exp}] 591 # cos(2kPi+x)=cos(x) 592 foreach {n integer} [divPiQuarter $integer $precision] {break} 593 # now integer>=0 and <Pi/2 594 set d [expr {$n%4}] 595 # add trigonometric circle turns number to delta 596 incr delta [expr {abs($n)}] 597 set signe 0 598 # cos(Pi-x)=-cos(x) 599 # cos(-x)=cos(x) 600 # cos(Pi/2-x)=sin(x) 601 switch -- $d { 602 1 {set signe 1;set l [_sin2 $integer $precision $delta]} 603 2 {set signe 1;set l [_cos2 $integer $precision $delta]} 604 0 {set l [_cos2 $integer $precision $delta]} 605 3 {set l [_sin2 $integer $precision $delta]} 606 default {error "internal error"} 607 } 608 # precision -> exp (multiplied by -1) 609 #idebug break 610 lset l 1 [expr {-([lindex $l 1])}] 611 # set the sign 612 if {$signe} { 613 lset l 0 [expr {-[lindex $l 0]}] 614 } 615 #idebug break 616 return [normalize [linsert $l 0 F]] 617} 618 619################################################################################ 620# compute cos(x) where 0<=x<Pi/2 621# returns : a list formed with : 622# 1. the mantissa 623# 2. the precision (opposite of the exponent) 624# 3. the uncertainty (doubt range) 625################################################################################ 626proc ::math::bigfloat::_cos2 {x precision delta} { 627 # precision bits after the dot 628 set pi [_pi $precision] 629 set pis2 [expr {$pi>>1}] 630 set pis4 [expr {$pis2>>1}] 631 if {$x>=$pis4} { 632 # cos(Pi/2-x)=sin(x) 633 set x [expr {$pis2-$x}] 634 incr delta 635 return [_sin $x $precision $delta] 636 } 637 #idebug break 638 return [_cos $x $precision $delta] 639} 640 641################################################################################ 642# compute cos(x) where 0<=x<Pi/4 643# returns : a list formed with : 644# 1. the mantissa 645# 2. the precision (opposite of the exponent) 646# 3. the uncertainty (doubt range) 647################################################################################ 648proc ::math::bigfloat::_cos {x precision delta} { 649 set float1 [expr {1<<$precision}] 650 # Taylor development follows : 651 # cos(x)=1-x^2/2 + x^4/4! ... + (-1)^(2n)*x^(2n)/2n! 652 # cos(x)= 1 - x^2/1.2 * (1 - x^2/3.4 * (... * (1 - x^2/(2n.(2n-1))...)) 653 # variables : $s (the Mantissa of the result) 654 # $denom1 & $denom2 (2n-1 & 2n) 655 # $x as the square of what is named x in 'cos(x)' 656 set s $float1 657 # 'd' is the uncertainty on x^2 658 set d [expr {$x*($delta<<1)}] 659 set d [expr {1+($d>>$precision)}] 660 # x=x^2 (because in this Taylor development, there are only even powers of x) 661 set x [expr {$x*$x>>$precision}] 662 set denom1 1 663 set denom2 2 664 set t [expr {-($x>>1)}] 665 set dt $d 666 while {$t!=0} { 667 incr s $t 668 incr delta $dt 669 incr denom1 2 670 incr denom2 2 671 set dt [expr {$x*$dt+($t+$dt)*$d>>$precision}] 672 incr dt 673 set t [expr {$x*$t>>$precision}] 674 set t [expr {-$t/($denom1*$denom2)}] 675 } 676 return [list $s $precision $delta] 677} 678 679################################################################################ 680# cotangent : the trivial algorithm is used 681################################################################################ 682proc ::math::bigfloat::cotan {x} { 683 return [::math::bigfloat::div [::math::bigfloat::cos $x] [::math::bigfloat::sin $x]] 684} 685 686################################################################################ 687# converts angles from degrees to radians 688# deg/180=rad/Pi 689################################################################################ 690proc ::math::bigfloat::deg2rad {x} { 691 checkFloat $x 692 set xLen [expr {-[lindex $x 2]}] 693 if {$xLen<3} { 694 error "number too loose to convert to radians" 695 } 696 set pi [pi $xLen 1] 697 return [div [mul $x $pi] 180] 698} 699 700 701 702################################################################################ 703# private proc to get : x modulo Pi/2 704# and the quotient (x divided by Pi/2) 705# used by cos , sin & others 706################################################################################ 707proc ::math::bigfloat::divPiQuarter {integer precision} { 708 incr precision 2 709 set integer [expr {$integer<<1}] 710 #idebug break 711 set P [_pi $precision] 712 # modulo 2Pi 713 set integer [expr {$integer%$P}] 714 # end modulo 2Pi 715 # 2Pi>>1 = Pi of course! 716 set P [expr {$P>>1}] 717 set n [expr {$integer/$P}] 718 set integer [expr {$integer%$P}] 719 # now divide by Pi/2 720 # multiply n by 2 721 set n [expr {$n<<1}] 722 # pi/2=Pi>>1 723 set P [expr {$P>>1}] 724 return [list [incr n [expr {$integer/$P}]] [expr {($integer%$P)>>1}]] 725} 726 727 728################################################################################ 729# divide A by B and returns the result 730# throw error : divide by zero 731################################################################################ 732proc ::math::bigfloat::div {a b} { 733 checkNumber $a 734 checkNumber $b 735 # dispatch to an appropriate procedure 736 if {[isInt $a]} { 737 if {[isInt $b]} { 738 return [expr {$a/$b}] 739 } 740 error "trying to divide an integer by a BigFloat" 741 } 742 if {[isInt $b]} {return [divFloatByInt $a $b]} 743 foreach {dummy integerA expA deltaA} $a {break} 744 foreach {dummy integerB expB deltaB} $b {break} 745 # computes the limits of the doubt (or uncertainty) interval 746 set BMin [expr {$integerB-$deltaB}] 747 set BMax [expr {$integerB+$deltaB}] 748 if {$BMin>$BMax} { 749 # swap BMin and BMax 750 set temp $BMin 751 set BMin $BMax 752 set BMax $temp 753 } 754 # multiply by zero gives zero 755 if {$integerA==0} { 756 # why not return any number or the integer 0 ? 757 # because there is an exponent that might be different between two BigFloats 758 # 0.00 --> exp = -2, 0.000000 -> exp = -6 759 return $a 760 } 761 # test of the division by zero 762 if {$BMin*$BMax<0 || $BMin==0 || $BMax==0} { 763 error "divide by zero" 764 } 765 # shift A because we need accuracy 766 set l [bits $integerB] 767 set integerA [expr {$integerA<<$l}] 768 set deltaA [expr {$deltaA<<$l}] 769 set exp [expr {$expA-$l-$expB}] 770 # relative uncertainties (dX/X) are added 771 # to give the relative uncertainty of the result 772 # i.e. 3% on A + 2% on B --> 5% on the quotient 773 # d(A/B)/(A/B)=dA/A + dB/B 774 # Q=A/B 775 # dQ=dA/B + dB*A/B*B 776 # dQ is "delta" 777 set delta [expr {($deltaB*abs($integerA))/abs($integerB)}] 778 set delta [expr {([incr delta]+$deltaA)/abs($integerB)}] 779 set quotient [expr {$integerA/$integerB}] 780 if {$integerB*$integerA<0} { 781 incr quotient -1 782 } 783 return [normalize [list F $quotient $exp [incr delta]]] 784} 785 786 787 788 789################################################################################ 790# divide a BigFloat A by a BigInt B 791# throw error : divide by zero 792################################################################################ 793proc ::math::bigfloat::divFloatByInt {a b} { 794 # type check 795 checkFloat $a 796 if {![isInt $b]} { 797 error "second argument is not an integer" 798 } 799 foreach {dummy integer exp delta} $a {break} 800 # zero divider test 801 if {$b==0} { 802 error "divide by zero" 803 } 804 # shift left for accuracy ; see other comments in [div] procedure 805 set l [bits $b] 806 set integer [expr {$integer<<$l}] 807 set delta [expr {$delta<<$l}] 808 incr exp -$l 809 set integer [expr {$integer/$b}] 810 # the uncertainty is always evaluated to the ceil value 811 # and as an absolute value 812 set delta [expr {$delta/abs($b)+1}] 813 return [normalize [list F $integer $exp $delta]] 814} 815 816 817 818 819 820################################################################################ 821# returns 1 if A and B are equal, 0 otherwise 822# IN : a, b (BigFloats) 823################################################################################ 824proc ::math::bigfloat::equal {a b} { 825 if {[isInt $a] && [isInt $b]} { 826 return [expr {$a==$b}] 827 } 828 # now a & b should only be BigFloats 829 checkFloat $a 830 checkFloat $b 831 foreach {dummy aint aexp adelta} $a {break} 832 foreach {dummy bint bexp bdelta} $b {break} 833 # set all Mantissas and Deltas to the same level (exponent) 834 # with lshift 835 set diff [expr {$aexp-$bexp}] 836 if {$diff<0} { 837 set diff [expr {-$diff}] 838 set bint [expr {$bint<<$diff}] 839 set bdelta [expr {$bdelta<<$diff}] 840 } elseif {$diff>0} { 841 set aint [expr {$aint<<$diff}] 842 set adelta [expr {$adelta<<$diff}] 843 } 844 # compute limits of the number's doubt range 845 set asupInt [expr {$aint+$adelta}] 846 set ainfInt [expr {$aint-$adelta}] 847 set bsupInt [expr {$bint+$bdelta}] 848 set binfInt [expr {$bint-$bdelta}] 849 # A & B are equal 850 # if their doubt ranges overlap themselves 851 if {$bint==$aint} { 852 return 1 853 } 854 if {$bint>$aint} { 855 set r [expr {$asupInt>=$binfInt}] 856 } else { 857 set r [expr {$bsupInt>=$ainfInt}] 858 } 859 return $r 860} 861 862################################################################################ 863# returns exp(X) where X is a BigFloat 864################################################################################ 865proc ::math::bigfloat::exp {x} { 866 checkFloat $x 867 foreach {dummy integer exp delta} $x {break} 868 if {$exp>=0} { 869 # shift till exp<0 with respect to the internal representation 870 # of the number 871 incr exp 872 set integer [expr {$integer<<$exp}] 873 set delta [expr {$delta<<$exp}] 874 set exp -1 875 } 876 # add 8 bits of precision for safety 877 set precision [expr {8-$exp}] 878 set integer [expr {$integer<<8}] 879 set delta [expr {$delta<<8}] 880 set Log2 [_log2 $precision] 881 set new_exp [expr {$integer/$Log2}] 882 set integer [expr {$integer%$Log2}] 883 # $new_exp = integer part of x/log(2) 884 # $integer = remainder 885 # exp(K.log(2)+r)=2^K.exp(r) 886 # so we just have to compute exp(r), r is small so 887 # the Taylor development will converge quickly 888 incr delta $new_exp 889 foreach {integer delta} [_exp $integer $precision $delta] {break} 890 set delta [expr {$delta>>8}] 891 incr precision -8 892 # multiply by 2^K , and take care of the sign 893 # example : X=-6.log(2)+0.01 894 # exp(X)=exp(0.01)*2^-6 895 # if {abs($new_exp)>>30!=0} { 896 # error "floating-point overflow due to exp" 897 # } 898 set exp [expr {$new_exp-$precision}] 899 incr delta 900 return [normalize [list F [expr {$integer>>8}] $exp $delta]] 901} 902 903 904################################################################################ 905# private procedure to compute exponentials 906# using Taylor development of exp(x) : 907# exp(x)=1+ x + x^2/2 + x^3/3! +...+x^n/n! 908# input : integer (the mantissa) 909# precision (the number of decimals) 910# delta (the doubt limit, or uncertainty) 911# returns a list : 1. the mantissa of the result 912# 2. the doubt limit, or uncertainty 913################################################################################ 914proc ::math::bigfloat::_exp {integer precision delta} { 915 if {$integer==0} { 916 # exp(0)=1 917 return [list [expr {1<<$precision}] $delta] 918 } 919 set s [expr {(1<<$precision)+$integer}] 920 set d [expr {1+$delta/2}] 921 incr delta $delta 922 # dt = uncertainty on x^2 923 set dt [expr {1+($d*$integer>>$precision)}] 924 # t= x^2/2 = x^2>>1 925 set t [expr {$integer*$integer>>$precision+1}] 926 set denom 2 927 while {$t!=0} { 928 # the sum is called 's' 929 incr s $t 930 incr delta $dt 931 # we do not have to keep trace of the factorial, we just iterate divisions 932 incr denom 933 # add delta 934 set d [expr {1+$d/$denom}] 935 incr dt $d 936 # get x^n from x^(n-1) 937 set t [expr {($integer*$t>>$precision)/$denom}] 938 } 939 return [list $s $delta] 940} 941################################################################################ 942# divide a BigFloat by 2 power 'n' 943################################################################################ 944proc ::math::bigfloat::floatRShift {float {n 1}} { 945 return [lset float 2 [expr {[lindex $float 2]-$n}]] 946} 947 948 949 950################################################################################ 951# procedure floor : identical to [expr floor($x)] in functionality 952# arguments : number IN (a BigFloat) 953# returns : the floor value as a BigInt 954################################################################################ 955proc ::math::bigfloat::floor {number} { 956 checkFloat $number 957 if {[iszero $number]} { 958 # returns the BigInt 0 959 return 0 960 } 961 foreach {dummy integer exp delta} $number {break} 962 if {$exp>=0} { 963 error "not enough precision to perform rounding (floor)" 964 } 965 # floor(n.xxxx)=n when n is positive 966 if {$integer>0} {return [expr {$integer>>(-$exp)}]} 967 set integer [expr {abs($integer)}] 968 # integer part 969 set try [expr {$integer>>(-$exp)}] 970 # floor(-n.xxxx)=-(n+1) when xxxx!=0 971 if {$try<<(-$exp)!=$integer} { 972 incr try 973 } 974 return [expr {-$try}] 975} 976 977 978################################################################################ 979# returns a list formed by an integer and an exponent 980# x = (A +/- C) * 10 power B 981# return [list "F" A B C] (where F is the BigFloat tag) 982# A and C are BigInts, B is a raw integer 983# return also a BigInt when there is neither a dot, nor a 'e' exponent 984# 985# arguments : -base base integer 986# or integer 987# or float 988# or float trailingZeros 989################################################################################ 990proc ::math::bigfloat::fromstr {number {addzeros 0}} { 991 if {$addzeros<0} { 992 error "second argument has to be a positive integer" 993 } 994 # eliminate the sign problem 995 # added on 05/08/2005 996 # setting '$signe' to the sign of the number 997 set number [string trimleft $number +] 998 if {[string index $number 0]=="-"} { 999 set signe 1 1000 set string [string range $number 1 end] 1001 } else { 1002 set signe 0 1003 set string $number 1004 } 1005 # integer case (not a floating-point number) 1006 if {[string is digit $string]} { 1007 if {$addzeros!=0} { 1008 error "second argument not allowed with an integer" 1009 } 1010 # we have completed converting an integer to a BigInt 1011 # please note that most math::bigfloat procs accept BigInts as arguments 1012 return $number 1013 } 1014 # floating-point number : check for an exponent 1015 # scientific notation 1016 set tab [split $string e] 1017 if {[llength $tab]>2} { 1018 # there are more than one 'e' letter in the number 1019 error "syntax error in number : $string" 1020 } 1021 if {[llength $tab]==2} { 1022 set exp [lindex $tab 1] 1023 # now exp can look like +099 so you need to handle octal numbers 1024 # too bad... 1025 # find the sign (if any?) 1026 regexp {^[\+\-]?} $exp expsign 1027 # trim the number with left-side 0's 1028 set found [string length $expsign] 1029 set exp $expsign[string trimleft [string range $exp $found end] 0] 1030 set mantissa [lindex $tab 0] 1031 } else { 1032 set exp 0 1033 set mantissa [lindex $tab 0] 1034 } 1035 # a floating-point number may have a dot 1036 set tab [split [string trimleft $mantissa 0] .] 1037 if {[llength $tab]>2} {error "syntax error in number : $string"} 1038 if {[llength $tab]==2} { 1039 set mantissa [join $tab ""] 1040 # increment by the number of decimals (after the dot) 1041 incr exp -[string length [lindex $tab 1]] 1042 } 1043 # this is necessary to ensure we can call fromstr (recursively) with 1044 # the mantissa ($number) 1045 if {![string is digit $mantissa]} { 1046 error "$number is not a number" 1047 } 1048 # take account of trailing zeros 1049 incr exp -$addzeros 1050 # multiply $number by 10^$trailingZeros 1051 append mantissa [string repeat 0 $addzeros] 1052 # add the sign 1053 # here we avoid octal numbers by trimming the leading zeros! 1054 # 2005-10-28 S.ARNOLD 1055 if {$signe} {set mantissa [expr {-[string trimleft $mantissa 0]}]} 1056 # the F tags a BigFloat 1057 # a BigInt is like any other integer since Tcl 8.5, 1058 # because expr now supports arbitrary length integers 1059 return [_fromstr $mantissa $exp] 1060} 1061 1062################################################################################ 1063# private procedure to transform decimal floats into binary ones 1064# IN : 1065# - number : a BigInt representing the Mantissa 1066# - exp : the decimal exponent (a simple integer) 1067# OUT : 1068# $number * 10^$exp, as the internal binary representation of a BigFloat 1069################################################################################ 1070proc ::math::bigfloat::_fromstr {number exp} { 1071 set number [string trimleft $number 0] 1072 if {$number==""} { 1073 return [list F 0 [expr {int($exp*log(10)/log(2))-15}] [expr {1<<15}]] 1074 } 1075 if {$exp==0} { 1076 return [list F $number 0 1] 1077 } 1078 if {$exp>0} { 1079 # mul by 10^exp, then normalize 1080 set power [expr {10**$exp}] 1081 set number [expr {$number*$power}] 1082 return [normalize [list F $number 0 $power]] 1083 } 1084 # now exp is negative or null 1085 # the closest power of 2 to the 'exp'th power of ten, but greater than it 1086 # 10**$exp<2**$binaryExp 1087 # $binaryExp>$exp*log(10)/log(2) 1088 set binaryExp [expr {int(-$exp*log(10)/log(2))+1+16}] 1089 # then compute n * 2^binaryExp / 10^(-exp) 1090 # (exp is negative) 1091 # equals n * 2^(binaryExp+exp) / 5^(-exp) 1092 set diff [expr {$binaryExp+$exp}] 1093 if {$diff<0} { 1094 error "internal error" 1095 } 1096 set power [expr {5**(-$exp)}] 1097 set number [expr {($number<<$diff)/$power}] 1098 set delta [expr {(1<<$diff)/$power}] 1099 return [normalize [list F $number [expr {-$binaryExp}] [incr delta]]] 1100} 1101 1102 1103################################################################################ 1104# fromdouble : 1105# like fromstr, but for a double scalar value 1106# arguments : 1107# double - the number to convert to a BigFloat 1108# exp (optional) - the total number of digits 1109################################################################################ 1110proc ::math::bigfloat::fromdouble {double {exp {}}} { 1111 set mantissa [lindex [split $double e] 0] 1112 # line added by SArnold on 05/08/2005 1113 set mantissa [string trimleft [string map {+ "" - ""} $mantissa] 0] 1114 set precision [string length [string map {. ""} $mantissa]] 1115 if { $exp != {} && [incr exp]>$precision } { 1116 return [fromstr $double [expr {$exp-$precision}]] 1117 } else { 1118 # tests have failed : not enough precision or no exp argument 1119 return [fromstr $double] 1120 } 1121} 1122 1123 1124################################################################################ 1125# converts a BigInt into a BigFloat with a given decimal precision 1126################################################################################ 1127proc ::math::bigfloat::int2float {int {decimals 1}} { 1128 # it seems like we need some kind of type handling 1129 # very odd in this Tcl world :-( 1130 if {![isInt $int]} { 1131 error "first argument is not an integer" 1132 } 1133 if {$decimals<1} { 1134 error "non-positive decimals number" 1135 } 1136 # the lowest number of decimals is 1, because 1137 # [tostr [fromstr 10.0]] returns 10. 1138 # (we lose 1 digit when converting back to string) 1139 set int [expr {$int*10**$decimals}] 1140 return [_fromstr $int [expr {-$decimals}]] 1141} 1142 1143 1144 1145################################################################################ 1146# multiplies 'leftop' by 'rightop' and rshift the result by 'shift' 1147################################################################################ 1148proc ::math::bigfloat::intMulShift {leftop rightop shift} { 1149 return [::math::bignum::rshift [::math::bignum::mul $leftop $rightop] $shift] 1150} 1151 1152################################################################################ 1153# returns 1 if x is a BigFloat, 0 elsewhere 1154################################################################################ 1155proc ::math::bigfloat::isFloat {x} { 1156 # a BigFloat is a list of : "F" mantissa exponent delta 1157 if {[llength $x]!=4} { 1158 return 0 1159 } 1160 # the marker is the letter "F" 1161 if {[string equal [lindex $x 0] F]} { 1162 return 1 1163 } 1164 return 0 1165} 1166 1167################################################################################ 1168# checks that n is a BigInt (a number create by math::bignum::fromstr) 1169################################################################################ 1170proc ::math::bigfloat::isInt {n} { 1171 if {[llength $n]>1} { 1172 return 0 1173 } 1174 # if {[string is digit $n]} { 1175 # return 1 1176 # } 1177 return 1 1178} 1179 1180 1181 1182################################################################################ 1183# returns 1 if x is null, 0 otherwise 1184################################################################################ 1185proc ::math::bigfloat::iszero {x} { 1186 if {[isInt $x]} { 1187 return [expr {$x==0}] 1188 } 1189 checkFloat $x 1190 # now we do some interval rounding : if a number's interval englobs 0, 1191 # it is considered to be equal to zero 1192 foreach {dummy integer exp delta} $x {break} 1193 if {$delta>=abs($integer)} {return 1} 1194 return 0 1195} 1196 1197 1198################################################################################ 1199# compute log(X) 1200################################################################################ 1201proc ::math::bigfloat::log {x} { 1202 checkFloat $x 1203 foreach {dummy integer exp delta} $x {break} 1204 if {$integer<=0} { 1205 error "zero logarithm error" 1206 } 1207 if {[iszero $x]} { 1208 error "number equals zero" 1209 } 1210 set precision [bits $integer] 1211 # uncertainty of the logarithm 1212 set delta [_logOnePlusEpsilon $delta $integer $precision] 1213 incr delta 1214 # we got : x = 1xxxxxx (binary number with 'precision' bits) * 2^exp 1215 # we need : x = 0.1xxxxxx(binary) *2^(exp+precision) 1216 incr exp $precision 1217 foreach {integer deltaIncr} [_log $integer] {break} 1218 incr delta $deltaIncr 1219 # log(a * 2^exp)= log(a) + exp*log(2) 1220 # result = log(x) + exp*log(2) 1221 # as x<1 log(x)<0 but 'integer' (result of '_log') is the absolute value 1222 # that is why we substract $integer to log(2)*$exp 1223 set integer [expr {[_log2 $precision]*$exp-$integer}] 1224 incr delta [expr {abs($exp)}] 1225 return [normalize [list F $integer -$precision $delta]] 1226} 1227 1228 1229################################################################################ 1230# compute log(1-epsNum/epsDenom)=log(1-'epsilon') 1231# Taylor development gives -x -x^2/2 -x^3/3 -x^4/4 ... 1232# used by 'log' command because log(x+/-epsilon)=log(x)+log(1+/-(epsilon/x)) 1233# so the uncertainty equals abs(log(1-epsilon/x)) 1234# ================================================ 1235# arguments : 1236# epsNum IN (the numerator of epsilon) 1237# epsDenom IN (the denominator of epsilon) 1238# precision IN (the number of bits after the dot) 1239# 1240# 'epsilon' = epsNum*2^-precision/epsDenom 1241################################################################################ 1242proc ::math::bigfloat::_logOnePlusEpsilon {epsNum epsDenom precision} { 1243 if {$epsNum>=$epsDenom} { 1244 error "number is null" 1245 } 1246 set s [expr {($epsNum<<$precision)/$epsDenom}] 1247 set divider 2 1248 set t [expr {$s*$epsNum/$epsDenom}] 1249 set u [expr {$t/$divider}] 1250 # when u (the current term of the development) is zero, we have reached our goal 1251 # it has converged 1252 while {$u!=0} { 1253 incr s $u 1254 # divider = order of the term = 'n' 1255 incr divider 1256 # t = (epsilon)^n 1257 set t [expr {$t*$epsNum/$epsDenom}] 1258 # u = t/n = (epsilon)^n/n and is the nth term of the Taylor development 1259 set u [expr {$t/$divider}] 1260 } 1261 return $s 1262} 1263 1264 1265################################################################################ 1266# compute log(0.xxxxxxxx) : log(1-epsilon)=-eps-eps^2/2-eps^3/3...-eps^n/n 1267################################################################################ 1268proc ::math::bigfloat::_log {integer} { 1269 # the uncertainty is nbSteps with nbSteps<=nbBits 1270 # take nbSteps=nbBits (the worse case) and log(nbBits+increment)=increment 1271 set precision [bits $integer] 1272 set n [expr {int(log($precision+2*log($precision)))}] 1273 set integer [expr {$integer<<$n}] 1274 incr precision $n 1275 set delta 3 1276 # 1-epsilon=integer 1277 set integer [expr {(1<<$precision)-$integer}] 1278 set s $integer 1279 # t=x^2 1280 set t [expr {$integer*$integer>>$precision}] 1281 set denom 2 1282 # u=x^2/2 (second term) 1283 set u [expr {$t/$denom}] 1284 while {$u!=0} { 1285 # while the current term is not zero, it has not converged 1286 incr s $u 1287 incr delta 1288 # t=x^n 1289 set t [expr {$t*$integer>>$precision}] 1290 # denom = n (the order of the current development term) 1291 # u = x^n/n (the nth term of Taylor development) 1292 set u [expr {$t/[incr denom]}] 1293 } 1294 # shift right to restore the precision 1295 set delta 1296 return [list [expr {$s>>$n}] [expr {($delta>>$n)+1}]] 1297} 1298 1299################################################################################ 1300# computes log(num/denom) with 'precision' bits 1301# used to compute some analysis constants with a given accuracy 1302# you might not call this procedure directly : it assumes 'num/denom'>4/5 1303# and 'num/denom'<1 1304################################################################################ 1305proc ::math::bigfloat::__log {num denom precision} { 1306 # Please Note : we here need a precision increment, in order to 1307 # keep accuracy at $precision digits. If we just hold $precision digits, 1308 # each number being precise at the last digit +/- 1, 1309 # we would lose accuracy because small uncertainties add to themselves. 1310 # Example : 0.0001 + 0.0010 = 0.0011 +/- 0.0002 1311 # This is quite the same reason that made tcl_precision defaults to 12 : 1312 # internally, doubles are computed with 17 digits, but to keep precision 1313 # we need to limit our results to 12. 1314 # The solution : given a precision target, increment precision with a 1315 # computed value so that all digits of he result are exacts. 1316 # 1317 # p is the precision 1318 # pk is the precision increment 1319 # 2 power pk is also the maximum number of iterations 1320 # for a number close to 1 but lower than 1, 1321 # (denom-num)/denum is (in our case) lower than 1/5 1322 # so the maximum nb of iterations is for: 1323 # 1/5*(1+1/5*(1/2+1/5*(1/3+1/5*(...)))) 1324 # the last term is 1/n*(1/5)^n 1325 # for the last term to be lower than 2^(-p-pk) 1326 # the number of iterations has to be 1327 # 2^(-pk).(1/5)^(2^pk) < 2^(-p-pk) 1328 # log(1/5).2^pk < -p 1329 # 2^pk > p/log(5) 1330 # pk > log(2)*log(p/log(5)) 1331 # now set the variable n to the precision increment i.e. pk 1332 set n [expr {int(log(2)*log($precision/log(5)))+1}] 1333 incr precision $n 1334 # log(num/denom)=log(1-(denom-num)/denom) 1335 # log(1+x) = x + x^2/2 + x^3/3 + ... + x^n/n 1336 # = x(1 + x(1/2 + x(1/3 + x(...+ x(1/(n-1) + x/n)...)))) 1337 set num [expr {$denom-$num}] 1338 # $s holds the result 1339 set s [expr {($num<<$precision)/$denom}] 1340 # $t holds x^n 1341 set t [expr {$s*$num/$denom}] 1342 set d 2 1343 # $u holds x^n/n 1344 set u [expr {$t/$d}] 1345 while {$u!=0} { 1346 incr s $u 1347 # get x^n * x 1348 set t [expr {$t*$num/$denom}] 1349 # get n+1 1350 incr d 1351 # then : $u = x^(n+1)/(n+1) 1352 set u [expr {$t/$d}] 1353 } 1354 # see head of the proc : we return the value with its target precision 1355 return [expr {$s>>$n}] 1356} 1357 1358################################################################################ 1359# computes log(2) with 'precision' bits and caches it into a namespace variable 1360################################################################################ 1361proc ::math::bigfloat::__logbis {precision} { 1362 set increment [expr {int(log($precision)/log(2)+1)}] 1363 incr precision $increment 1364 # ln(2)=3*ln(1-4/5)+ln(1-125/128) 1365 set a [__log 125 128 $precision] 1366 set b [__log 4 5 $precision] 1367 set r [expr {$b*3+$a}] 1368 set ::math::bigfloat::Log2 [expr {$r>>$increment}] 1369 # formerly (when BigFloats were stored in ten radix) we had to compute log(10) 1370 # ln(10)=10.ln(1-4/5)+3*ln(1-125/128) 1371} 1372 1373 1374################################################################################ 1375# retrieves log(2) with 'precision' bits ; the result is cached 1376################################################################################ 1377proc ::math::bigfloat::_log2 {precision} { 1378 variable Log2 1379 if {![info exists Log2]} { 1380 __logbis $precision 1381 } else { 1382 # the constant is cached and computed again when more precision is needed 1383 set l [bits $Log2] 1384 if {$precision>$l} { 1385 __logbis $precision 1386 } 1387 } 1388 # return log(2) with 'precision' bits even when the cached value has more bits 1389 return [_round $Log2 $precision] 1390} 1391 1392 1393################################################################################ 1394# returns A modulo B (like with fmod() math function) 1395################################################################################ 1396proc ::math::bigfloat::mod {a b} { 1397 checkNumber $a 1398 checkNumber $b 1399 if {[isInt $a] && [isInt $b]} {return [expr {$a%$b}]} 1400 if {[isInt $a]} {error "trying to divide an integer by a BigFloat"} 1401 set quotient [div $a $b] 1402 # examples : fmod(3,2)=1 quotient=1.5 1403 # fmod(1,2)=1 quotient=0.5 1404 # quotient>0 and b>0 : get floor(quotient) 1405 # fmod(-3,-2)=-1 quotient=1.5 1406 # fmod(-1,-2)=-1 quotient=0.5 1407 # quotient>0 and b<0 : get floor(quotient) 1408 # fmod(-3,2)=-1 quotient=-1.5 1409 # fmod(-1,2)=-1 quotient=-0.5 1410 # quotient<0 and b>0 : get ceil(quotient) 1411 # fmod(3,-2)=1 quotient=-1.5 1412 # fmod(1,-2)=1 quotient=-0.5 1413 # quotient<0 and b<0 : get ceil(quotient) 1414 if {[sign $quotient]} { 1415 set quotient [ceil $quotient] 1416 } else { 1417 set quotient [floor $quotient] 1418 } 1419 return [sub $a [mul $quotient $b]] 1420} 1421 1422################################################################################ 1423# returns A times B 1424################################################################################ 1425proc ::math::bigfloat::mul {a b} { 1426 checkNumber $a 1427 checkNumber $b 1428 # dispatch the command to appropriate commands regarding types (BigInt & BigFloat) 1429 if {[isInt $a]} { 1430 if {[isInt $b]} { 1431 return [expr {$a*$b}] 1432 } 1433 return [mulFloatByInt $b $a] 1434 } 1435 if {[isInt $b]} {return [mulFloatByInt $a $b]} 1436 # now we are sure that 'a' and 'b' are BigFloats 1437 foreach {dummy integerA expA deltaA} $a {break} 1438 foreach {dummy integerB expB deltaB} $b {break} 1439 # 2^expA * 2^expB = 2^(expA+expB) 1440 set exp [expr {$expA+$expB}] 1441 # mantissas are multiplied 1442 set integer [expr {$integerA*$integerB}] 1443 # compute precisely the uncertainty 1444 set delta [expr {$deltaA*(abs($integerB)+$deltaB)+abs($integerA)*$deltaB+1}] 1445 # we have to normalize because 'delta' may be too big 1446 return [normalize [list F $integer $exp $delta]] 1447} 1448 1449################################################################################ 1450# returns A times B, where B is a positive integer 1451################################################################################ 1452proc ::math::bigfloat::mulFloatByInt {a b} { 1453 checkFloat $a 1454 foreach {dummy integer exp delta} $a {break} 1455 if {$b==0} { 1456 return [list F 0 $exp $delta] 1457 } 1458 # Mantissa and Delta are simply multplied by $b 1459 set integer [expr {$integer*$b}] 1460 set delta [expr {$delta*$b}] 1461 # We normalize because Delta could have seriously increased 1462 return [normalize [list F $integer $exp $delta]] 1463} 1464 1465################################################################################ 1466# normalizes a number : Delta (accuracy of the BigFloat) 1467# has to be limited, because the memory use increase 1468# quickly when we do some computations, as the Mantissa and Delta 1469# increase together 1470# The solution : limit the size of Delta to 16 bits 1471################################################################################ 1472proc ::math::bigfloat::normalize {number} { 1473 checkFloat $number 1474 foreach {dummy integer exp delta} $number {break} 1475 set l [bits $delta] 1476 if {$l>16} { 1477 incr l -16 1478 # $l holds the supplementary size (in bits) 1479 # now we can shift right by $l bits 1480 # always round upper the Delta 1481 set delta [expr {$delta>>$l}] 1482 incr delta 1483 set integer [expr {$integer>>$l}] 1484 incr exp $l 1485 } 1486 return [list F $integer $exp $delta] 1487} 1488 1489 1490 1491################################################################################ 1492# returns -A (the opposite) 1493################################################################################ 1494proc ::math::bigfloat::opp {a} { 1495 checkNumber $a 1496 if {[iszero $a]} { 1497 return $a 1498 } 1499 if {[isInt $a]} { 1500 return [expr {-$a}] 1501 } 1502 # recursive call 1503 lset a 1 [expr {-[lindex $a 1]}] 1504 return $a 1505} 1506 1507################################################################################ 1508# gets Pi with precision bits 1509# after the dot (after you call [tostr] on the result) 1510################################################################################ 1511proc ::math::bigfloat::pi {precision {binary 0}} { 1512 if {![isInt $precision]} { 1513 error "'$precision' expected to be an integer" 1514 } 1515 if {!$binary} { 1516 # convert decimal digit length into bit length 1517 set precision [expr {int(ceil($precision*log(10)/log(2)))}] 1518 } 1519 return [list F [_pi $precision] -$precision 1] 1520} 1521 1522# 1523# Procedure that resets the stored cached Pi constant 1524# 1525proc ::math::bigfloat::reset {} { 1526 variable _pi0 1527 if {[info exists _pi0]} {unset _pi0} 1528} 1529 1530proc ::math::bigfloat::_pi {precision} { 1531 # the constant Pi begins with 3.xxx 1532 # so we need 2 digits to store the digit '3' 1533 # and then we will have precision+2 bits in the mantissa 1534 variable _pi0 1535 if {![info exists _pi0]} { 1536 set _pi0 [__pi $precision] 1537 } 1538 set lenPiGlobal [bits $_pi0] 1539 if {$lenPiGlobal<$precision} { 1540 set _pi0 [__pi $precision] 1541 } 1542 return [expr {$_pi0 >> [bits $_pi0]-2-$precision}] 1543} 1544 1545################################################################################ 1546# computes an integer representing Pi in binary radix, with precision bits 1547################################################################################ 1548proc ::math::bigfloat::__pi {precision} { 1549 set safetyLimit 8 1550 # for safety and for the better precision, we do so ... 1551 incr precision $safetyLimit 1552 # formula found in the Math litterature (on Wikipedia 1553 # Pi/4 = 44.atan(1/57) + 7.atan(1/239) - 12.atan(1/682) + 24.atan(1/12943) 1554 set a [expr {[_atanfract 57 $precision]*44}] 1555 incr a [expr {[_atanfract 239 $precision]*7}] 1556 set a [expr {$a - [_atanfract 682 $precision]*12}] 1557 incr a [expr {[_atanfract 12943 $precision]*24}] 1558 return [expr {$a>>$safetyLimit-2}] 1559} 1560 1561################################################################################ 1562# shift right an integer until it haves $precision bits 1563# round at the same time 1564################################################################################ 1565proc ::math::bigfloat::_round {integer precision} { 1566 set shift [expr {[bits $integer]-$precision}] 1567 if {$shift==0} { 1568 return $integer 1569 } 1570 # $result holds the shifted integer 1571 set result [expr {$integer>>$shift}] 1572 # $shift-1 is the bit just rights the last bit of the result 1573 # Example : integer=1000010 shift=2 1574 # => result=10000 and the tested bit is '1' 1575 if {$integer & (1<<($shift-1))} { 1576 # we round to the upper limit 1577 return [incr result] 1578 } 1579 return $result 1580} 1581 1582################################################################################ 1583# returns A power B, where B is a positive integer 1584################################################################################ 1585proc ::math::bigfloat::pow {a b} { 1586 checkNumber $a 1587 if {$b<0} { 1588 error "pow : exponent is not a positive integer" 1589 } 1590 # case where it is obvious that we should use the appropriate command 1591 # from math::bignum (added 5th March 2005) 1592 if {[isInt $a]} { 1593 return [expr {$a**$b}] 1594 } 1595 # algorithm : exponent=$b = Sum(i=0..n) b(i)2^i 1596 # $a^$b = $a^( b(0) + 2b(1) + 4b(2) + ... + 2^n*b(n) ) 1597 # we have $a^(x+y)=$a^x * $a^y 1598 # then $a^$b = Product(i=0...n) $a^(2^i*b(i)) 1599 # b(i) is boolean so $a^(2^i*b(i))= 1 when b(i)=0 and = $a^(2^i) when b(i)=1 1600 # then $a^$b = Product(i=0...n and b(i)=1) $a^(2^i) and 1 when $b=0 1601 if {$b==0} {return 1} 1602 # $res holds the result 1603 set res 1 1604 while {1} { 1605 # at the beginning i=0 1606 # $remainder is b(i) 1607 set remainder [expr {$b&1}] 1608 # $b 'rshift'ed by 1 bit : i=i+1 1609 # so next time we will test bit b(i+1) 1610 set b [expr {$b>>1}] 1611 # if b(i)=1 1612 if {$remainder} { 1613 # mul the result by $a^(2^i) 1614 # if i=0 we multiply by $a^(2^0)=$a^1=$a 1615 set res [mul $res $a] 1616 } 1617 # no more bits at '1' in $b : $res is the result 1618 if {$b==0} { 1619 return [normalize $res] 1620 } 1621 # i=i+1 : $a^(2^(i+1)) = square of $a^(2^i) 1622 set a [mul $a $a] 1623 } 1624} 1625 1626################################################################################ 1627# converts angles for radians to degrees 1628################################################################################ 1629proc ::math::bigfloat::rad2deg {x} { 1630 checkFloat $x 1631 set xLen [expr {-[lindex $x 2]}] 1632 if {$xLen<3} { 1633 error "number too loose to convert to degrees" 1634 } 1635 # $rad/Pi=$deg/180 1636 # so result in deg = $radians*180/Pi 1637 return [div [mul $x 180] [pi $xLen 1]] 1638} 1639 1640################################################################################ 1641# retourne la partie enti�re (ou 0) du nombre "number" 1642################################################################################ 1643proc ::math::bigfloat::round {number} { 1644 checkFloat $number 1645 #set number [normalize $number] 1646 # fetching integers (or BigInts) from the internal representation 1647 foreach {dummy integer exp delta} $number {break} 1648 if {$integer==0} { 1649 return 0 1650 } 1651 if {$exp>=0} { 1652 error "not enough precision to round (in round)" 1653 } 1654 set exp [expr {-$exp}] 1655 # saving the sign, ... 1656 set sign [expr {$integer<0}] 1657 set integer [expr {abs($integer)}] 1658 # integer part of the number 1659 set try [expr {$integer>>$exp}] 1660 # first bit after the dot 1661 set way [expr {$integer>>($exp-1)&1}] 1662 # delta is shifted so it gives the integer part of 2*delta 1663 set delta [expr {$delta>>($exp-1)}] 1664 # when delta is too big to compute rounded value ( 1665 if {$delta!=0} { 1666 error "not enough precision to round (in round)" 1667 } 1668 if {$way} { 1669 incr try 1670 } 1671 # ... restore the sign now 1672 if {$sign} {return [expr {-$try}]} 1673 return $try 1674} 1675 1676################################################################################ 1677# round and divide by 10^n 1678################################################################################ 1679proc ::math::bigfloat::roundshift {integer n} { 1680 # $exp= 10^$n 1681 incr n -1 1682 set exp [expr {10**$n}] 1683 set toround [expr {$integer/$exp}] 1684 if {$toround%10>=5} { 1685 return [expr {$toround/10+1}] 1686 } 1687 return [expr {$toround/10}] 1688} 1689 1690################################################################################ 1691# gets the sign of either a bignum, or a BitFloat 1692# we keep the bignum convention : 0 for positive, 1 for negative 1693################################################################################ 1694proc ::math::bigfloat::sign {n} { 1695 if {[isInt $n]} { 1696 return [expr {$n<0}] 1697 } 1698 checkFloat $n 1699 # sign of 0=0 1700 if {[iszero $n]} {return 0} 1701 # the sign of the Mantissa, which is a BigInt 1702 return [sign [lindex $n 1]] 1703} 1704 1705 1706################################################################################ 1707# gets sin(x) 1708################################################################################ 1709proc ::math::bigfloat::sin {x} { 1710 checkFloat $x 1711 foreach {dummy integer exp delta} $x {break} 1712 if {$exp>-2} { 1713 error "sin : not enough precision" 1714 } 1715 set precision [expr {-$exp}] 1716 # sin(2kPi+x)=sin(x) 1717 # $integer is now the modulo of the division of the mantissa by Pi/4 1718 # and $n is the quotient 1719 foreach {n integer} [divPiQuarter $integer $precision] {break} 1720 incr delta $n 1721 set d [expr {$n%4}] 1722 # now integer>=0 1723 # x = $n*Pi/4 + $integer and $n belongs to [0,3] 1724 # sin(2Pi-x)=-sin(x) 1725 # sin(Pi-x)=sin(x) 1726 # sin(Pi/2+x)=cos(x) 1727 set sign 0 1728 switch -- $d { 1729 0 {set l [_sin2 $integer $precision $delta]} 1730 1 {set l [_cos2 $integer $precision $delta]} 1731 2 {set sign 1;set l [_sin2 $integer $precision $delta]} 1732 3 {set sign 1;set l [_cos2 $integer $precision $delta]} 1733 default {error "internal error"} 1734 } 1735 # $l is a list : {Mantissa Precision Delta} 1736 # precision --> the opposite of the exponent 1737 # 1.000 = 1000*10^-3 so exponent=-3 and precision=3 digits 1738 lset l 1 [expr {-([lindex $l 1])}] 1739 # the sign depends on the switch statement below 1740 #::math::bignum::setsign integer $sign 1741 if {$sign} { 1742 lset l 0 [expr {-[lindex $l 0]}] 1743 } 1744 # we insert the Bigfloat tag (F) and normalize the final result 1745 return [normalize [linsert $l 0 F]] 1746} 1747 1748proc ::math::bigfloat::_sin2 {x precision delta} { 1749 set pi [_pi $precision] 1750 # shift right by 1 = divide by 2 1751 # shift right by 2 = divide by 4 1752 set pis2 [expr {$pi>>1}] 1753 set pis4 [expr {$pis2>>1}] 1754 if {$x>=$pis4} { 1755 # sin(Pi/2-x)=cos(x) 1756 incr delta 1757 set x [expr {$pis2-$x}] 1758 return [_cos $x $precision $delta] 1759 } 1760 return [_sin $x $precision $delta] 1761} 1762 1763################################################################################ 1764# sin(x) with 'x' lower than Pi/4 and positive 1765# 'x' is the Mantissa - 'delta' is Delta 1766# 'precision' is the opposite of the exponent 1767################################################################################ 1768proc ::math::bigfloat::_sin {x precision delta} { 1769 # $s holds the result 1770 set s $x 1771 # sin(x) = x - x^3/3! + x^5/5! - ... + (-1)^n*x^(2n+1)/(2n+1)! 1772 # = x * (1 - x^2/(2*3) * (1 - x^2/(4*5) * (...* (1 - x^2/(2n*(2n+1)) )...))) 1773 # The second expression allows us to compute the less we can 1774 1775 # $double holds the uncertainty (Delta) of x^2 : 2*(Mantissa*Delta) + Delta^2 1776 # (Mantissa+Delta)^2=Mantissa^2 + 2*Mantissa*Delta + Delta^2 1777 set double [expr {$x*$delta>>$precision-1}] 1778 incr double [expr {1+$delta*$delta>>$precision}] 1779 # $x holds the Mantissa of x^2 1780 set x [expr {$x*$x>>$precision}] 1781 set dt [expr {$x*$delta+$double*($s+$delta)>>$precision}] 1782 incr dt 1783 # $t holds $s * -(x^2) / (2n*(2n+1)) 1784 # mul by x^2 1785 set t [expr {$s*$x>>$precision}] 1786 set denom2 2 1787 set denom3 3 1788 # mul by -1 (opp) and divide by 2*3 1789 set t [expr {-$t/($denom2*$denom3)}] 1790 while {$t!=0} { 1791 incr s $t 1792 incr delta $dt 1793 # incr n => 2n --> 2n+2 and 2n+1 --> 2n+3 1794 incr denom2 2 1795 incr denom3 2 1796 # $dt is the Delta corresponding to $t 1797 # $double "" "" "" "" $x (x^2) 1798 # ($t+$dt) * ($x+$double) = $t*$x + ($dt*$x + $t*$double) + $dt*$double 1799 # Mantissa^ ^--------Delta-------------------^ 1800 set dt [expr {$x*$dt+($t+$dt)*$double>>$precision}] 1801 set t [expr {$t*$x>>$precision}] 1802 # removed 2005/08/31 by sarnold75 1803 #set dt [::math::bignum::add $dt $double] 1804 set denom [expr {$denom2*$denom3}] 1805 # now computing : div by -2n(2n+1) 1806 set dt [expr {1+$dt/$denom}] 1807 set t [expr {-$t/$denom}] 1808 } 1809 return [list $s $precision $delta] 1810} 1811 1812 1813################################################################################ 1814# procedure for extracting the square root of a BigFloat 1815################################################################################ 1816proc ::math::bigfloat::sqrt {x} { 1817 checkFloat $x 1818 foreach {dummy integer exp delta} $x {break} 1819 # if x=0, return 0 1820 if {[iszero $x]} { 1821 # return zero, taking care of its precision ($exp) 1822 return [list F 0 $exp $delta] 1823 } 1824 # we cannot get sqrt(x) if x<0 1825 if {[lindex $integer 0]<0} { 1826 error "negative sqrt input" 1827 } 1828 # (1+epsilon)^p = 1 + epsilon*(p-1) + epsilon^2*(p-1)*(p-2)/2! + ... 1829 # + epsilon^n*(p-1)*...*(p-n)/n! 1830 # sqrt(1 + epsilon) = (1 + epsilon)^(1/2) 1831 # = 1 - epsilon/2 - epsilon^2*3/(4*2!) - ... 1832 # - epsilon^n*(3*5*..*(2n-1))/(2^n*n!) 1833 # sqrt(1 - epsilon) = 1 + Sum(i=1..infinity) epsilon^i*(3*5*...*(2i-1))/(i!*2^i) 1834 # sqrt(n +/- delta)=sqrt(n) * sqrt(1 +/- delta/n) 1835 # so the uncertainty on sqrt(n +/- delta) equals sqrt(n) * (sqrt(1 - delta/n) - 1) 1836 # sqrt(1+eps) < sqrt(1-eps) because their logarithm compare as : 1837 # -ln(2)(1+eps) < -ln(2)(1-eps) 1838 # finally : 1839 # Delta = sqrt(n) * Sum(i=1..infinity) (delta/n)^i*(3*5*...*(2i-1))/(i!*2^i) 1840 # here we compute the second term of the product by _sqrtOnePlusEpsilon 1841 set delta [_sqrtOnePlusEpsilon $delta $integer] 1842 set intLen [bits $integer] 1843 # removed 2005/08/31 by sarnold75, readded 2005/08/31 1844 set precision $intLen 1845 # intLen + exp = number of bits before the dot 1846 #set precision [expr {-$exp}] 1847 # square root extraction 1848 set integer [expr {$integer<<$intLen}] 1849 incr exp -$intLen 1850 incr intLen $intLen 1851 # there is an exponent 2^$exp : when $exp is odd, we would need to compute sqrt(2) 1852 # so we decrement $exp, in order to get it even, and we do not need sqrt(2) anymore ! 1853 if {$exp&1} { 1854 incr exp -1 1855 set integer [expr {$integer<<1}] 1856 incr intLen 1857 incr precision 1858 } 1859 # using a low-level (taken from math::bignum) root extraction procedure 1860 # using binary operators 1861 set integer [_sqrt $integer] 1862 # delta has to be multiplied by the square root 1863 set delta [expr {$delta*$integer>>$precision}] 1864 # round to the ceiling the uncertainty (worst precision, the fastest to compute) 1865 incr delta 1866 # we are sure that $exp is even, see above 1867 return [normalize [list F $integer [expr {$exp/2}] $delta]] 1868} 1869 1870 1871 1872################################################################################ 1873# compute abs(sqrt(1-delta/integer)-1) 1874# the returned value is a relative uncertainty 1875################################################################################ 1876proc ::math::bigfloat::_sqrtOnePlusEpsilon {delta integer} { 1877 # sqrt(1-x) - 1 = x/2 + x^2*3/(2^2*2!) + x^3*3*5/(2^3*3!) + ... 1878 # = x/2 * (1 + x*3/(2*2) * ( 1 + x*5/(2*3) * 1879 # (...* (1 + x*(2n-1)/(2n) ) )...))) 1880 set l [bits $integer] 1881 # to compute delta/integer we have to shift left to keep the same precision level 1882 # we have a better accuracy computing (delta << lg(integer))/integer 1883 # than computing (delta/integer) << lg(integer) 1884 set x [expr {($delta<<$l)/$integer}] 1885 # denom holds 2n 1886 set denom 4 1887 # x/2 1888 set result [expr {$x>>1}] 1889 # x^2*3/(2!*2^2) 1890 # numerator holds 2n-1 1891 set numerator 3 1892 set temp [expr {($result*$delta*$numerator)/($integer*$denom)}] 1893 incr temp 1894 while {$temp!=0} { 1895 incr result $temp 1896 incr numerator 2 1897 incr denom 2 1898 # n = n+1 ==> num=num+2 denom=denom+2 1899 # num=2n+1 denom=2n+2 1900 set temp [expr {($temp*$delta*$numerator)/($integer*$denom)}] 1901 } 1902 return $result 1903} 1904 1905# 1906# Computes the square root of an integer 1907# Returns an integer 1908# 1909proc ::math::bigfloat::_sqrt {n} { 1910 set i [expr {(([bits $n]-1)/2)+1}] 1911 set b [expr {$i*2}] ; # Bit to set to get 2^i*2^i 1912 1913 set r 0 ; # guess 1914 set x 0 ; # guess^2 1915 set s 0 ; # guess^2 backup 1916 set t 0 ; # intermediate result 1917 for {} {$i >= 0} {incr i -1; incr b -2} { 1918 set x [expr {$s+($t|(1<<$b))}] 1919 if {abs($x)<= abs($n)} { 1920 set s $x 1921 set r [expr {$r|(1<<$i)}] 1922 set t [expr {$t|(1<<$b+1)}] 1923 } 1924 set t [expr {$t>>1}] 1925 } 1926 return $r 1927} 1928 1929################################################################################ 1930# substracts B to A 1931################################################################################ 1932proc ::math::bigfloat::sub {a b} { 1933 checkNumber $a 1934 checkNumber $b 1935 if {[isInt $a] && [isInt $b]} { 1936 # the math::bignum::sub proc is designed to work with BigInts 1937 return [expr {$a-$b}] 1938 } 1939 return [add $a [opp $b]] 1940} 1941 1942################################################################################ 1943# tangent (trivial algorithm) 1944################################################################################ 1945proc ::math::bigfloat::tan {x} { 1946 return [::math::bigfloat::div [::math::bigfloat::sin $x] [::math::bigfloat::cos $x]] 1947} 1948 1949################################################################################ 1950# returns a power of ten 1951################################################################################ 1952proc ::math::bigfloat::tenPow {n} { 1953 return [expr {10**$n}] 1954} 1955 1956 1957################################################################################ 1958# converts a BigInt to a double (basic floating-point type) 1959# with respect to the global variable 'tcl_precision' 1960################################################################################ 1961proc ::math::bigfloat::todouble {x} { 1962 global tcl_precision 1963 set precision $tcl_precision 1964 if {$precision==0} { 1965 # this is a cheat, I must admit, for Tcl 8.5 1966 set precision 16 1967 } 1968 checkFloat $x 1969 # get the string repr of x without the '+' sign 1970 # please note: here we call math::bigfloat::tostr 1971 set result [string trimleft [tostr $x] +] 1972 set minus "" 1973 if {[string index $result 0]=="-"} { 1974 set minus - 1975 set result [string range $result 1 end] 1976 } 1977 1978 set l [split $result e] 1979 set exp 0 1980 if {[llength $l]==2} { 1981 # exp : x=Mantissa*2^Exp 1982 set exp [lindex $l 1] 1983 } 1984 # caution with octal numbers : we have to remove heading zeros 1985 # but count them as digits 1986 regexp {^0*} $result zeros 1987 incr exp -[string length $zeros] 1988 # Mantissa = integerPart.fractionalPart 1989 set l [split [lindex $l 0] .] 1990 set integerPart [lindex $l 0] 1991 set integerLen [string length $integerPart] 1992 set fractionalPart [lindex $l 1] 1993 # The number of digits in Mantissa, excluding the dot and the leading zeros, of course 1994 set integer [string trimleft $integerPart$fractionalPart 0] 1995 if {$integer eq ""} { 1996 set integer 0 1997 } 1998 set len [string length $integer] 1999 # Now Mantissa is stored in $integer 2000 if {$len>$precision} { 2001 set lenDiff [expr {$len-$precision}] 2002 # true when the number begins with a zero 2003 set zeroHead 0 2004 if {[string index $integer 0]==0} { 2005 incr lenDiff -1 2006 set zeroHead 1 2007 } 2008 set integer [roundshift $integer $lenDiff] 2009 if {$zeroHead} { 2010 set integer 0$integer 2011 } 2012 set len [string length $integer] 2013 if {$len<$integerLen} { 2014 set exp [expr {$integerLen-$len}] 2015 # restore the true length 2016 set integerLen $len 2017 } 2018 } 2019 # number = 'sign'*'integer'*10^'exp' 2020 if {$exp==0} { 2021 # no scientific notation 2022 set exp "" 2023 } else { 2024 # scientific notation 2025 set exp e$exp 2026 } 2027 # place the dot just before the index $integerLen in the Mantissa 2028 set result [string range $integer 0 [expr {$integerLen-1}]] 2029 append result .[string range $integer $integerLen end] 2030 # join the Mantissa with the sign before and the exponent after 2031 return $minus$result$exp 2032} 2033 2034################################################################################ 2035# converts a number stored as a list to a string in which all digits are true 2036################################################################################ 2037proc ::math::bigfloat::tostr {args} { 2038 if {[llength $args]==2} { 2039 if {![string equal [lindex $args 0] -nosci]} {error "unknown option: should be -nosci"} 2040 set nosci yes 2041 set number [lindex $args 1] 2042 } else { 2043 if {[llength $args]!=1} {error "syntax error: should be tostr ?-nosci? number"} 2044 set nosci no 2045 set number [lindex $args 0] 2046 } 2047 if {[isInt $number]} { 2048 return $number 2049 } 2050 checkFloat $number 2051 foreach {dummy integer exp delta} $number {break} 2052 if {[iszero $number]} { 2053 # we do matter how much precision $number has : 2054 # it can be 0.0000000 or 0.0, the result is not the same zero 2055 #return 0 2056 } 2057 if {$exp>0} { 2058 # the power of ten the closest but greater than 2^$exp 2059 # if it was lower than the power of 2, we would have more precision 2060 # than existing in the number 2061 set newExp [expr {int(ceil($exp*log(2)/log(10)))}] 2062 # 'integer' <- 'integer' * 2^exp / 10^newExp 2063 # equals 'integer' * 2^(exp-newExp) / 5^newExp 2064 set binExp [expr {$exp-$newExp}] 2065 if {$binExp<0} { 2066 # it cannot happen 2067 error "internal error" 2068 } 2069 # 5^newExp 2070 set fivePower [expr {5**$newExp}] 2071 # 'lshift'ing $integer by $binExp bits is like multiplying it by 2^$binExp 2072 # but much, much faster 2073 set integer [expr {($integer<<$binExp)/$fivePower}] 2074 # $integer is the Mantissa - Delta should follow the same operations 2075 set delta [expr {($delta<<$binExp)/$fivePower}] 2076 set exp $newExp 2077 } elseif {$exp<0} { 2078 # the power of ten the closest but lower than 2^$exp 2079 # same remark about the precision 2080 set newExp [expr {int(floor(-$exp*log(2)/log(10)))}] 2081 # 'integer' <- 'integer' * 10^newExp / 2^(-exp) 2082 # equals 'integer' * 5^(newExp) / 2^(-exp-newExp) 2083 set binShift [expr {-$exp-$newExp}] 2084 set fivePower [expr {5**$newExp}] 2085 # rshifting is like dividing by 2^$binShift, but faster as we said above about lshift 2086 set integer [expr {$integer*$fivePower>>$binShift}] 2087 set delta [expr {$delta*$fivePower>>$binShift}] 2088 set exp -$newExp 2089 } 2090 # saving the sign, to restore it into the result 2091 set result [expr {abs($integer)}] 2092 set sign [expr {$integer<0}] 2093 # rounded 'integer' +/- 'delta' 2094 set up [expr {$result+$delta}] 2095 set down [expr {$result-$delta}] 2096 if {($up<0 && $down>0)||($up>0 && $down<0)} { 2097 # $up>0 and $down<0 or vice-versa : then the number is considered equal to zero 2098 set isZero yes 2099 # delta <= 2**n (n = bits(delta)) 2100 # 2**n <= 10**exp , then 2101 # exp >= n.log(2)/log(10) 2102 # delta <= 10**(n.log(2)/log(10)) 2103 incr exp [expr {int(ceil([bits $delta]*log(2)/log(10)))}] 2104 set result 0 2105 } else { 2106 # iterate until the convergence of the rounding 2107 # we incr $shift until $up and $down are rounded to the same number 2108 # at each pass we lose one digit of precision, so necessarly it will success 2109 for {set shift 1} { 2110 [roundshift $up $shift]!=[roundshift $down $shift] 2111 } { 2112 incr shift 2113 } {} 2114 incr exp $shift 2115 set result [roundshift $up $shift] 2116 set isZero no 2117 } 2118 set l [string length $result] 2119 # now formatting the number the most nicely for having a clear reading 2120 # would'nt we allow a number being constantly displayed 2121 # as : 0.2947497845e+012 , would we ? 2122 if {$nosci} { 2123 if {$exp >= 0} { 2124 append result [string repeat 0 $exp]. 2125 } elseif {$l + $exp > 0} { 2126 set result [string range $result 0 end-[expr {-$exp}]].[string range $result end-[expr {-1-$exp}] end] 2127 } else { 2128 set result 0.[string repeat 0 [expr {-$exp-$l}]]$result 2129 } 2130 } else { 2131 if {$exp>0} { 2132 # we display 423*10^6 as : 4.23e+8 2133 # Length of mantissa : $l 2134 # Increment exp by $l-1 because the first digit is placed before the dot, 2135 # the other ($l-1) digits following the dot. 2136 incr exp [incr l -1] 2137 set result [string index $result 0].[string range $result 1 end] 2138 append result "e+$exp" 2139 } elseif {$exp==0} { 2140 # it must have a dot to be a floating-point number (syntaxically speaking) 2141 append result . 2142 } else { 2143 set exp [expr {-$exp}] 2144 if {$exp < $l} { 2145 # we can display the number nicely as xxxx.yyyy* 2146 # the problem of the sign is solved finally at the bottom of the proc 2147 set n [string range $result 0 end-$exp] 2148 incr exp -1 2149 append n .[string range $result end-$exp end] 2150 set result $n 2151 } elseif {$l==$exp} { 2152 # we avoid to use the scientific notation 2153 # because it is harder to read 2154 set result "0.$result" 2155 } else { 2156 # ... but here there is no choice, we should not represent a number 2157 # with more than one leading zero 2158 set result [string index $result 0].[string range $result 1 end]e-[expr {$exp-$l+1}] 2159 } 2160 } 2161 } 2162 # restore the sign : we only put a minus on numbers that are different from zero 2163 if {$sign==1 && !$isZero} {set result "-$result"} 2164 return $result 2165} 2166 2167################################################################################ 2168# PART IV 2169# HYPERBOLIC FUNCTIONS 2170################################################################################ 2171 2172################################################################################ 2173# hyperbolic cosinus 2174################################################################################ 2175proc ::math::bigfloat::cosh {x} { 2176 # cosh(x) = (exp(x)+exp(-x))/2 2177 # dividing by 2 is done faster by 'rshift'ing 2178 return [floatRShift [add [exp $x] [exp [opp $x]]] 1] 2179} 2180 2181################################################################################ 2182# hyperbolic sinus 2183################################################################################ 2184proc ::math::bigfloat::sinh {x} { 2185 # sinh(x) = (exp(x)-exp(-x))/2 2186 # dividing by 2 is done faster by 'rshift'ing 2187 return [floatRShift [sub [exp $x] [exp [opp $x]]] 1] 2188} 2189 2190################################################################################ 2191# hyperbolic tangent 2192################################################################################ 2193proc ::math::bigfloat::tanh {x} { 2194 set up [exp $x] 2195 set down [exp [opp $x]] 2196 # tanh(x)=sinh(x)/cosh(x)= (exp(x)-exp(-x))/2/ [(exp(x)+exp(-x))/2] 2197 # =(exp(x)-exp(-x))/(exp(x)+exp(-x)) 2198 # =($up-$down)/($up+$down) 2199 return [div [sub $up $down] [add $up $down]] 2200} 2201 2202# exporting public interface 2203namespace eval ::math::bigfloat { 2204 foreach function { 2205 add mul sub div mod pow 2206 iszero compare equal 2207 fromstr tostr fromdouble todouble 2208 int2float isInt isFloat 2209 exp log sqrt round ceil floor 2210 sin cos tan cotan asin acos atan 2211 cosh sinh tanh abs opp 2212 pi deg2rad rad2deg 2213 } { 2214 namespace export $function 2215 } 2216} 2217 2218# (AM) No "namespace import" - this should be left to the user! 2219#namespace import ::math::bigfloat::* 2220 2221package provide math::bigfloat 2.0.1 2222