######################################################################## # BigFloat for Tcl # Copyright (C) 2003-2005 ARNOLD Stephane # It is published with the terms of tcllib's BSD-style license. # See the file named license.terms. ######################################################################## package require Tcl 8.5 # this line helps when I want to source this file again and again catch {namespace delete ::math::bigfloat} # private namespace # this software works only with Tcl v8.4 and higher # it is using the package math::bignum namespace eval ::math::bigfloat { # cached constants # ln(2) with arbitrary precision variable Log2 # Pi with arb. precision variable Pi variable _pi0 } ################################################################################ # procedures that handle floating-point numbers # these procedures are sorted by name (after eventually removing the underscores) # # BigFloats are internally represented as a list : # {"F" Mantissa Exponent Delta} where "F" is a character which determins # the datatype, Mantissa and Delta are two big integers and Exponent another integer. # # The BigFloat value equals to (Mantissa +/- Delta)*2^Exponent # So the internal representation is binary, but trying to get as close as possible to # the decimal one when converted to a string. # When calling [fromstr], the Delta parameter is set to the value of 1 at the position # of the last decimal digit. # Example : 1.50 belongs to [1.49,1.51], but internally Delta may not equal to 1. # Because of the binary representation, it is between 1 and 1+(2^-15). # # So Mantissa and Delta are not limited in size, but in practice Delta is kept under # 2^32 by the 'normalize' procedure, to avoid a never-ended growth of memory used. # Indeed, when you perform some computations, the Delta parameter (which represent # the uncertainty on the value of the Mantissa) may increase. # Exponent, as an integer, is limited to 32 bits, and this limit seems fair. # The exponent is indeed involved in logarithmic computations, so it may be # a mistake to give it a too large value. # Retrieving the parameters of a BigFloat is often done with that command : # foreach {dummy int exp delta} $bigfloat {break} # (dummy is not used, it is just used to get the "F" marker). # The isInt, isFloat, checkNumber and checkFloat procedures are used # to check data types # # Taylor development are often used to compute the analysis functions (like exp(),log()...) # To learn how it is done in practice, take a look at ::math::bigfloat::_asin # While doing computation on Mantissas, we do not care about the last digit, # because if we compute correctly Deltas, the digits that remain will be exact. ################################################################################ ################################################################################ # returns the absolute value ################################################################################ proc ::math::bigfloat::abs {number} { checkNumber $number if {[isInt $number]} { # set sign to positive for a BigInt return [expr {abs($number)}] } # set sign to positive for a BigFloat into the Mantissa (index 1) lset number 1 [expr {abs([lindex $number 1])}] return $number } ################################################################################ # arccosinus of a BigFloat ################################################################################ proc ::math::bigfloat::acos {x} { # handy proc for checking datatype checkFloat $x foreach {dummy entier exp delta} $x {break} set precision [expr {($exp<0)?(-$exp):1}] # acos(0.0)=Pi/2 # 26/07/2005 : changed precision from decimal to binary # with the second parameter of pi command set piOverTwo [floatRShift [pi $precision 1]] if {[iszero $x]} { # $x is too close to zero -> acos(0)=PI/2 return $piOverTwo } # acos(-x)= Pi/2 + asin(x) if {$entier<0} { return [add $piOverTwo [asin [abs $x]]] } # we always use _asin to compute the result # but as it is a Taylor development, the value given to [_asin] # has to be a bit smaller than 1 ; by using that trick : acos(x)=asin(sqrt(1-x^2)) # we can limit the entry of the Taylor development below 1/sqrt(2) if {[compare $x [fromstr 0.7071]]>0} { # x > sqrt(2)/2 : trying to make _asin converge quickly # creating 0 and 1 with the same precision as the entry set fzero [list F 0 -$precision 1] # 1.000 with $precision zeros set fone [list F [expr {1<<$precision}] -$precision 1] # when $x is close to 1 (acos(1.0)=0.0) if {[equal $fone $x]} { return $fzero } if {[compare $fone $x]<0} { # the behavior assumed because acos(x) is not defined # when |x|>1 error "acos on a number greater than 1" } # acos(x) = asin(sqrt(1 - x^2)) # since 1 - cos(x)^2 = sin(x)^2 # x> sqrt(2)/2 so x^2 > 1/2 so 1-x^2<1/2 set x [sqrt [sub $fone [mul $x $x]]] # the parameter named x is smaller than sqrt(2)/2 return [_asin $x] } # acos(x) = Pi/2 - asin(x) # x$expB} { set diff [expr {$expA-$expB}] set integerA [expr {$integerA<<$diff}] set deltaA [expr {$deltaA<<$diff}] incr integerA $integerB incr deltaA $deltaB return [normalize [list F $integerA $expB $deltaA]] } elseif {$expA==$expB} { # nothing to shift left return [normalize [list F [incr integerA $integerB] $expA [incr deltaA $deltaB]]] } else { error "internal error" } } ################################################################################ # returns the sum A(BigFloat) + B(BigInt) # the greatest advantage of this method is that the uncertainty # of the result remains unchanged, in respect to the entry's uncertainty (deltaA) ################################################################################ proc ::math::bigfloat::addInt2Float {a b} { # type checking checkFloat $a if {![isInt $b]} { error "second argument is not an integer" } # retrieving data from $a foreach {dummy integerA expA deltaA} $a {break} # to add an int to a BigFloat,... if {$expA>0} { # we have to put the integer integerA # to the level of zero exponent : 1e8 --> 100000000e0 set shift $expA set integerA [expr {($integerA<<$shift)+$b}] set deltaA [expr {$deltaA<<$shift}] # we have to normalize, because we have shifted the mantissa # and the uncertainty left return [normalize [list F $integerA 0 $deltaA]] } elseif {$expA==0} { # integerA is already at integer level : float=(integerA)e0 return [normalize [list F [incr integerA $b] \ 0 $deltaA]] } else { # here we have something like 234e-2 + 3 # we have to shift the integer left by the exponent |$expA| incr integerA [expr {$b<<(-$expA)}] return [normalize [list F $integerA $expA $deltaA]] } } ################################################################################ # arcsinus of a BigFloat ################################################################################ proc ::math::bigfloat::asin {x} { # type checking checkFloat $x foreach {dummy entier exp delta} $x {break} if {$exp>-1} { error "not enough precision on input (asin)" } set precision [expr {-$exp}] # when x=0, return 0 at the same precision as the input was if {[iszero $x]} { return [list F 0 -$precision 1] } # asin(-x)=-asin(x) if {$entier<0} { return [opp [asin [abs $x]]] } # 26/07/2005 : changed precision from decimal to binary set piOverTwo [floatRShift [pi $precision 1]] # now a little trick : asin(x)=Pi/2-asin(sqrt(1-x^2)) # so we can limit the entry of the Taylor development # to 1/sqrt(2)~0.7071 # the comparison is : if x>0.7071 then ... if {[compare $x [fromstr 0.7071]]>0} { set fone [list F [expr {1<<$precision}] -$precision 1] # asin(1)=Pi/2 (with the same precision as the entry has) if {[equal $fone $x]} { return $piOverTwo } if {[compare $x $fone]>0} { error "asin on a number greater than 1" } # asin(x)=Pi/2-asin(sqrt(1-x^2)) set x [sqrt [sub $fone [mul $x $x]]] return [sub $piOverTwo [_asin $x]] } return [normalize [_asin $x]] } ################################################################################ # _asin : arcsinus of numbers between 0 and +1 ################################################################################ proc ::math::bigfloat::_asin {x} { # Taylor development # asin(x)=x + 1/2 x^3/3 + 3/2.4 x^5/5 + 3.5/2.4.6 x^7/7 + ... # into this iterative form : # asin(x)=x * (1 + 1/2 * x^2 * (1/3 + 3/4 *x^2 * (... # ...* (1/(2n-1) + (2n-1)/2n * x^2 / (2n+1))...))) # we show how is really computed the development : # we don't need to set a var with x^n or a product of integers # all we need is : x^2, 2n-1, 2n, 2n+1 and a few variables foreach {dummy mantissa exp delta} $x {break} set precision [expr {-$exp}] if {$precision+1<[bits $mantissa]} { error "sinus greater than 1" } # precision is the number of after-dot digits set result $mantissa set delta_final $delta # resultat is the final result, and delta_final # will contain the uncertainty of the result # square is the square of the mantissa set square [expr {$mantissa*$mantissa>>$precision}] # dt is the uncertainty of Mantissa set dt [expr {$mantissa*$delta>>($precision-1)}] incr dt set num 1 # two will be used into the loop set i 3 set denom 2 # the nth factor equals : $num/$denom* $mantissa/$i set delta [expr {$delta*$square + $dt*($delta+$mantissa)}] set delta [expr {($delta*$num)/ $denom >>$precision}] incr delta # we do not multiply the Mantissa by $num right now because it is 1 ! # but we have Mantissa=$x # and we want Mantissa*$x^2 * $num / $denom / $i set mantissa [expr {($mantissa*$square>>$precision)/$denom}] # do not forget the modified Taylor development : # asin(x)=x * (1 + 1/2*x^2*(1/3 + 3/4*x^2*(...*(1/(2n-1) + (2n-1)/2n*x^2/(2n+1))...))) # all we need is : x^2, 2n-1, 2n, 2n+1 and a few variables # $num=2n-1 $denom=2n $square=x^2 and $i=2n+1 set mantissa_temp [expr {$mantissa/$i}] set delta_temp [expr {1+$delta/$i}] # when the Mantissa increment is smaller than the Delta increment, # we would not get much precision by continuing the development while {$mantissa_temp!=0} { # Mantissa = Mantissa * $num/$denom * $square # Add Mantissa/$i, which is stored in $mantissa_temp, to the result incr result $mantissa_temp incr delta_final $delta_temp # here we have $two instead of [fromstr 2] (optimization) # num=num+2,i=i+2,denom=denom+2 # because num=2n-1 denom=2n and i=2n+1 incr num 2 incr i 2 incr denom 2 # computes precisly the future Delta parameter set delta [expr {$delta*$square+$dt*($delta+$mantissa)}] set delta [expr {($delta*$num)/$denom>>$precision}] incr delta set mantissa [expr {$mantissa*$square>>$precision}] set mantissa [expr {($mantissa*$num)/$denom}] set mantissa_temp [expr {$mantissa/$i}] set delta_temp [expr {1+$delta/$i}] } return [normalize [list F $result $exp $delta_final]] } ################################################################################ # arctangent : returns atan(x) ################################################################################ proc ::math::bigfloat::atan {x} { checkFloat $x foreach {dummy mantissa exp delta} $x {break} if {$exp>=0} { error "not enough precision to compute atan" } set precision [expr {-$exp}] # atan(0)=0 if {[iszero $x]} { return [list F 0 -$precision $delta] } # atan(-x)=-atan(x) if {$mantissa<0} { return [opp [atan [abs $x]]] } # now x is strictly positive # at this moment, we are trying to limit |x| to a fair acceptable number # to ensure that Taylor development will converge quickly set float1 [list F [expr {1<<$precision}] -$precision 1] if {[compare $float1 $x]<0} { # compare x to 2.4142 if {[compare $x [fromstr 2.4142]]<0} { # atan(x)=Pi/4 + atan((x-1)/(x+1)) # as 10} { # atan(x)=Pi/4 + atan((x-1)/(x+1)) # x>0.420 so (x-1)/(x+1)=1 - 2/(x+1) > 1-2/1.414 # > -0.414 # x<1 so (x-1)/(x+1)<0 set pi_sur_quatre [floatRShift [pi $precision 1] 2] return [add $pi_sur_quatre [atan \ [div [sub $x $float1] [add $x $float1]]]] } # precision increment : to have less uncertainty # we add a little more precision so that the result would be more accurate # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) # when we have n steps in Taylor development : the nth term is : # x^(2n-1)/(2n-1) # and the loss of precision is of 2n (n sums and n divisions) # this command is called with x(precision-3/2)*log(2)-log(2n-1) # hence log(2n-1)<2n-1 # n*sqrt(2)>(precision-1.5)*log(2)+1-2n # n*(sqrt(2)+2)>(precision-1.5)*log(2)+1 set n [expr {int((log(2)*($precision-1.5)+1)/(sqrt(2)+2)+1)}] incr precision $n set mantissa [expr {$mantissa<<$n}] set delta [expr {$delta<<$n}] # end of adding precision increment # now computing Taylor development : # atan(x)=x - x^3/3 + x^5/5 - x^7/7 ... + (-1)^n*x^(2n+1)/(2n+1) # atan(x)=x * (1 - x^2 * (1/3 - x^2 * (1/5 - x^2 * (...*(1/(2n-1) - x^2 / (2n+1))...)))) # what do we need to compute this ? # x^2 ($square), 2n+1 ($divider), $result, the nth term of the development ($t) # and the nth term multiplied by 2n+1 ($temp) # then we do this (with care keeping as much precision as possible): # while ($t <>0) : # $result=$result+$t # $temp=$temp * $square # $divider = $divider+2 # $t=$temp/$divider # end-while set result $mantissa set delta_end $delta # we store the square of the integer (mantissa) # Delta of Mantissa^2 = Delta * 2 = Delta << 1 set delta_square [expr {$delta<<1}] set square [expr {$mantissa*$mantissa>>$precision}] # the (2n+1) divider set divider 3 # computing precisely the uncertainty set delta [expr {1+($delta_square*$mantissa+$delta*$square>>$precision)}] # temp contains (-1)^n*x^(2n+1) set temp [expr {-$mantissa*$square>>$precision}] set t [expr {$temp/$divider}] set dt [expr {1+$delta/$divider}] while {$t!=0} { incr result $t incr delta_end $dt incr divider 2 set delta [expr {1+($delta_square*abs($temp)+$delta*($delta_square+$square)>>$precision)}] set temp [expr {-$temp*$square>>$precision}] set t [expr {$temp/$divider}] set dt [expr {1+$delta/$divider}] } # we have to normalize because the uncertainty might be greater than 2**16 # moreover it is the most often case return [normalize [list F $result [expr {$exp-$n}] $delta_end]] } ################################################################################ # compute atan(1/integer) at a given precision # this proc is only used to compute Pi # it is using the same Taylor development as [atan] ################################################################################ proc ::math::bigfloat::_atanfract {integer precision} { # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) # when we have n steps in Taylor development : the nth term is : # 1/denom^(2n+1)/(2n+1) # and the loss of precision is of 2n (n sums and n divisions) # this command is called with integer>=5 # # We do not want to compute the Delta parameter, so we just # can increment precision (with lshift) in order for the result to be precise. # Remember : we compute atan2(1,$integer) with $precision bits # $integer has no Delta parameter as it is a BigInt, of course, so # theorically we could compute *any* number of digits. # # if we add an increment to the precision, say n: # (1/5)^(2n-1)/(2n-1) has to be lower than (1/2)^(precision+n-1) # Calculus : # log(left term) < log(right term) # log(1/left term) > log(1/right term) # (2n-1)*log(5)+log(2n-1)>(precision+n-1)*log(2) # n(2log(5)-log(2))>(precision-1)*log(2)-log(2n-1)+log(5) # -log(2n-1)>-(2n-1) # n(2log(5)-log(2)+2)>(precision-1)*log(2)+1+log(5) set n [expr {int((($precision-1)*log(2)+1+log(5))/(2*log(5)-log(2)+2)+1)}] incr precision $n # first term of the development : 1/integer set a [expr {(1<<$precision)/$integer}] # 's' will contain the result set s $a # Taylor development : x - x^3/3 + x^5/5 - ... + (-1)^(n+1)*x^(2n-1)/(2n-1) # equals x (1 - x^2 * (1/3 + x^2 * (... * (1/(2n-3) + (-1)^(n+1) * x^2 / (2n-1))...))) # all we need to store is : 2n-1 ($denom), x^(2n+1) and x^2 ($square) and two results : # - the nth term => $u # - the nth term * (2n-1) => $t # + of course, the result $s set square [expr {$integer*$integer}] set denom 3 # $t is (-1)^n*x^(2n+1) set t [expr {-$a/$square}] set u [expr {$t/$denom}] # we break the loop when the current term of the development is null while {$u!=0} { incr s $u # denominator= (2n+1) incr denom 2 # div $t by x^2 set t [expr {-$t/$square}] set u [expr {$t/$denom}] } # go back to the initial precision return [expr {$s>>$n}] } # # bits : computes the number of bits of an integer, approx. # proc ::math::bigfloat::bits {int} { set l [string length [set int [expr {abs($int)}]]] # int<10**l -> log_2(int)=l*log_2(10) set l [expr {int($l*log(10)/log(2))+1}] if {$int>>$l!=0} { error "bad result: $l bits" } while {($int>>($l-1))==0} { incr l -1 } return $l } ################################################################################ # returns the integer part of a BigFloat, as a BigInt # the result is the same one you would have # if you had called [expr {ceil($x)}] ################################################################################ proc ::math::bigfloat::ceil {number} { checkFloat $number set number [normalize $number] if {[iszero $number]} { return 0 } foreach {dummy integer exp delta} $number {break} if {$exp>=0} { error "not enough precision to perform rounding (ceil)" } # saving the sign ... set sign [expr {$integer<0}] set integer [expr {abs($integer)}] # integer part set try [expr {$integer>>(-$exp)}] if {$sign} { return [opp $try] } # fractional part if {($try<<(-$exp))!=$integer} { return [incr try] } return $try } ################################################################################ # checks each variable to be a BigFloat # arguments : each argument is the name of a variable to be checked ################################################################################ proc ::math::bigfloat::checkFloat {number} { if {![isFloat $number]} { error "BigFloat expected" } } ################################################################################ # checks if each number is either a BigFloat or a BigInt # arguments : each argument is the name of a variable to be checked ################################################################################ proc ::math::bigfloat::checkNumber {x} { if {![isFloat $x] && ![isInt $x]} { error "input is not an integer, nor a BigFloat" } } ################################################################################ # returns 0 if A and B are equal, else returns 1 or -1 # accordingly to the sign of (A - B) ################################################################################ proc ::math::bigfloat::compare {a b} { if {[isInt $a] && [isInt $b]} { set diff [expr {$a-$b}] if {$diff>0} {return 1} elseif {$diff<0} {return -1} return 0 } checkFloat $a checkFloat $b if {[equal $a $b]} {return 0} if {[lindex [sub $a $b] 1]<0} {return -1} return 1 } ################################################################################ # gets cos(x) # throws an error if there is not enough precision on the input ################################################################################ proc ::math::bigfloat::cos {x} { checkFloat $x foreach {dummy integer exp delta} $x {break} if {$exp>-2} { error "not enough precision on floating-point number" } set precision [expr {-$exp}] # cos(2kPi+x)=cos(x) foreach {n integer} [divPiQuarter $integer $precision] {break} # now integer>=0 and exp (multiplied by -1) #idebug break lset l 1 [expr {-([lindex $l 1])}] # set the sign if {$signe} { lset l 0 [expr {-[lindex $l 0]}] } #idebug break return [normalize [linsert $l 0 F]] } ################################################################################ # compute cos(x) where 0<=x>1}] set pis4 [expr {$pis2>>1}] if {$x>=$pis4} { # cos(Pi/2-x)=sin(x) set x [expr {$pis2-$x}] incr delta return [_sin $x $precision $delta] } #idebug break return [_cos $x $precision $delta] } ################################################################################ # compute cos(x) where 0<=x>$precision)}] # x=x^2 (because in this Taylor development, there are only even powers of x) set x [expr {$x*$x>>$precision}] set denom1 1 set denom2 2 set t [expr {-($x>>1)}] set dt $d while {$t!=0} { incr s $t incr delta $dt incr denom1 2 incr denom2 2 set dt [expr {$x*$dt+($t+$dt)*$d>>$precision}] incr dt set t [expr {$x*$t>>$precision}] set t [expr {-$t/($denom1*$denom2)}] } return [list $s $precision $delta] } ################################################################################ # cotangent : the trivial algorithm is used ################################################################################ proc ::math::bigfloat::cotan {x} { return [::math::bigfloat::div [::math::bigfloat::cos $x] [::math::bigfloat::sin $x]] } ################################################################################ # converts angles from degrees to radians # deg/180=rad/Pi ################################################################################ proc ::math::bigfloat::deg2rad {x} { checkFloat $x set xLen [expr {-[lindex $x 2]}] if {$xLen<3} { error "number too loose to convert to radians" } set pi [pi $xLen 1] return [div [mul $x $pi] 180] } ################################################################################ # private proc to get : x modulo Pi/2 # and the quotient (x divided by Pi/2) # used by cos , sin & others ################################################################################ proc ::math::bigfloat::divPiQuarter {integer precision} { incr precision 2 set integer [expr {$integer<<1}] #idebug break set P [_pi $precision] # modulo 2Pi set integer [expr {$integer%$P}] # end modulo 2Pi # 2Pi>>1 = Pi of course! set P [expr {$P>>1}] set n [expr {$integer/$P}] set integer [expr {$integer%$P}] # now divide by Pi/2 # multiply n by 2 set n [expr {$n<<1}] # pi/2=Pi>>1 set P [expr {$P>>1}] return [list [incr n [expr {$integer/$P}]] [expr {($integer%$P)>>1}]] } ################################################################################ # divide A by B and returns the result # throw error : divide by zero ################################################################################ proc ::math::bigfloat::div {a b} { checkNumber $a checkNumber $b # dispatch to an appropriate procedure if {[isInt $a]} { if {[isInt $b]} { return [expr {$a/$b}] } error "trying to divide an integer by a BigFloat" } if {[isInt $b]} {return [divFloatByInt $a $b]} foreach {dummy integerA expA deltaA} $a {break} foreach {dummy integerB expB deltaB} $b {break} # computes the limits of the doubt (or uncertainty) interval set BMin [expr {$integerB-$deltaB}] set BMax [expr {$integerB+$deltaB}] if {$BMin>$BMax} { # swap BMin and BMax set temp $BMin set BMin $BMax set BMax $temp } # multiply by zero gives zero if {$integerA==0} { # why not return any number or the integer 0 ? # because there is an exponent that might be different between two BigFloats # 0.00 --> exp = -2, 0.000000 -> exp = -6 return $a } # test of the division by zero if {$BMin*$BMax<0 || $BMin==0 || $BMax==0} { error "divide by zero" } # shift A because we need accuracy set l [bits $integerB] set integerA [expr {$integerA<<$l}] set deltaA [expr {$deltaA<<$l}] set exp [expr {$expA-$l-$expB}] # relative uncertainties (dX/X) are added # to give the relative uncertainty of the result # i.e. 3% on A + 2% on B --> 5% on the quotient # d(A/B)/(A/B)=dA/A + dB/B # Q=A/B # dQ=dA/B + dB*A/B*B # dQ is "delta" set delta [expr {($deltaB*abs($integerA))/abs($integerB)}] set delta [expr {([incr delta]+$deltaA)/abs($integerB)}] set quotient [expr {$integerA/$integerB}] if {$integerB*$integerA<0} { incr quotient -1 } return [normalize [list F $quotient $exp [incr delta]]] } ################################################################################ # divide a BigFloat A by a BigInt B # throw error : divide by zero ################################################################################ proc ::math::bigfloat::divFloatByInt {a b} { # type check checkFloat $a if {![isInt $b]} { error "second argument is not an integer" } foreach {dummy integer exp delta} $a {break} # zero divider test if {$b==0} { error "divide by zero" } # shift left for accuracy ; see other comments in [div] procedure set l [bits $b] set integer [expr {$integer<<$l}] set delta [expr {$delta<<$l}] incr exp -$l set integer [expr {$integer/$b}] # the uncertainty is always evaluated to the ceil value # and as an absolute value set delta [expr {$delta/abs($b)+1}] return [normalize [list F $integer $exp $delta]] } ################################################################################ # returns 1 if A and B are equal, 0 otherwise # IN : a, b (BigFloats) ################################################################################ proc ::math::bigfloat::equal {a b} { if {[isInt $a] && [isInt $b]} { return [expr {$a==$b}] } # now a & b should only be BigFloats checkFloat $a checkFloat $b foreach {dummy aint aexp adelta} $a {break} foreach {dummy bint bexp bdelta} $b {break} # set all Mantissas and Deltas to the same level (exponent) # with lshift set diff [expr {$aexp-$bexp}] if {$diff<0} { set diff [expr {-$diff}] set bint [expr {$bint<<$diff}] set bdelta [expr {$bdelta<<$diff}] } elseif {$diff>0} { set aint [expr {$aint<<$diff}] set adelta [expr {$adelta<<$diff}] } # compute limits of the number's doubt range set asupInt [expr {$aint+$adelta}] set ainfInt [expr {$aint-$adelta}] set bsupInt [expr {$bint+$bdelta}] set binfInt [expr {$bint-$bdelta}] # A & B are equal # if their doubt ranges overlap themselves if {$bint==$aint} { return 1 } if {$bint>$aint} { set r [expr {$asupInt>=$binfInt}] } else { set r [expr {$bsupInt>=$ainfInt}] } return $r } ################################################################################ # returns exp(X) where X is a BigFloat ################################################################################ proc ::math::bigfloat::exp {x} { checkFloat $x foreach {dummy integer exp delta} $x {break} if {$exp>=0} { # shift till exp<0 with respect to the internal representation # of the number incr exp set integer [expr {$integer<<$exp}] set delta [expr {$delta<<$exp}] set exp -1 } # add 8 bits of precision for safety set precision [expr {8-$exp}] set integer [expr {$integer<<8}] set delta [expr {$delta<<8}] set Log2 [_log2 $precision] set new_exp [expr {$integer/$Log2}] set integer [expr {$integer%$Log2}] # $new_exp = integer part of x/log(2) # $integer = remainder # exp(K.log(2)+r)=2^K.exp(r) # so we just have to compute exp(r), r is small so # the Taylor development will converge quickly incr delta $new_exp foreach {integer delta} [_exp $integer $precision $delta] {break} set delta [expr {$delta>>8}] incr precision -8 # multiply by 2^K , and take care of the sign # example : X=-6.log(2)+0.01 # exp(X)=exp(0.01)*2^-6 # if {abs($new_exp)>>30!=0} { # error "floating-point overflow due to exp" # } set exp [expr {$new_exp-$precision}] incr delta return [normalize [list F [expr {$integer>>8}] $exp $delta]] } ################################################################################ # private procedure to compute exponentials # using Taylor development of exp(x) : # exp(x)=1+ x + x^2/2 + x^3/3! +...+x^n/n! # input : integer (the mantissa) # precision (the number of decimals) # delta (the doubt limit, or uncertainty) # returns a list : 1. the mantissa of the result # 2. the doubt limit, or uncertainty ################################################################################ proc ::math::bigfloat::_exp {integer precision delta} { if {$integer==0} { # exp(0)=1 return [list [expr {1<<$precision}] $delta] } set s [expr {(1<<$precision)+$integer}] set d [expr {1+$delta/2}] incr delta $delta # dt = uncertainty on x^2 set dt [expr {1+($d*$integer>>$precision)}] # t= x^2/2 = x^2>>1 set t [expr {$integer*$integer>>$precision+1}] set denom 2 while {$t!=0} { # the sum is called 's' incr s $t incr delta $dt # we do not have to keep trace of the factorial, we just iterate divisions incr denom # add delta set d [expr {1+$d/$denom}] incr dt $d # get x^n from x^(n-1) set t [expr {($integer*$t>>$precision)/$denom}] } return [list $s $delta] } ################################################################################ # divide a BigFloat by 2 power 'n' ################################################################################ proc ::math::bigfloat::floatRShift {float {n 1}} { return [lset float 2 [expr {[lindex $float 2]-$n}]] } ################################################################################ # procedure floor : identical to [expr floor($x)] in functionality # arguments : number IN (a BigFloat) # returns : the floor value as a BigInt ################################################################################ proc ::math::bigfloat::floor {number} { checkFloat $number if {[iszero $number]} { # returns the BigInt 0 return 0 } foreach {dummy integer exp delta} $number {break} if {$exp>=0} { error "not enough precision to perform rounding (floor)" } # floor(n.xxxx)=n when n is positive if {$integer>0} {return [expr {$integer>>(-$exp)}]} set integer [expr {abs($integer)}] # integer part set try [expr {$integer>>(-$exp)}] # floor(-n.xxxx)=-(n+1) when xxxx!=0 if {$try<<(-$exp)!=$integer} { incr try } return [expr {-$try}] } ################################################################################ # returns a list formed by an integer and an exponent # x = (A +/- C) * 10 power B # return [list "F" A B C] (where F is the BigFloat tag) # A and C are BigInts, B is a raw integer # return also a BigInt when there is neither a dot, nor a 'e' exponent # # arguments : -base base integer # or integer # or float # or float trailingZeros ################################################################################ proc ::math::bigfloat::fromstr {number {addzeros 0}} { if {$addzeros<0} { error "second argument has to be a positive integer" } # eliminate the sign problem # added on 05/08/2005 # setting '$signe' to the sign of the number set number [string trimleft $number +] if {[string index $number 0]=="-"} { set signe 1 set string [string range $number 1 end] } else { set signe 0 set string $number } # integer case (not a floating-point number) if {[string is digit $string]} { if {$addzeros!=0} { error "second argument not allowed with an integer" } # we have completed converting an integer to a BigInt # please note that most math::bigfloat procs accept BigInts as arguments return $number } # floating-point number : check for an exponent # scientific notation set tab [split $string e] if {[llength $tab]>2} { # there are more than one 'e' letter in the number error "syntax error in number : $string" } if {[llength $tab]==2} { set exp [lindex $tab 1] # now exp can look like +099 so you need to handle octal numbers # too bad... # find the sign (if any?) regexp {^[\+\-]?} $exp expsign # trim the number with left-side 0's set found [string length $expsign] set exp $expsign[string trimleft [string range $exp $found end] 0] set mantissa [lindex $tab 0] } else { set exp 0 set mantissa [lindex $tab 0] } # a floating-point number may have a dot set tab [split [string trimleft $mantissa 0] .] if {[llength $tab]>2} {error "syntax error in number : $string"} if {[llength $tab]==2} { set mantissa [join $tab ""] # increment by the number of decimals (after the dot) incr exp -[string length [lindex $tab 1]] } # this is necessary to ensure we can call fromstr (recursively) with # the mantissa ($number) if {![string is digit $mantissa]} { error "$number is not a number" } # take account of trailing zeros incr exp -$addzeros # multiply $number by 10^$trailingZeros append mantissa [string repeat 0 $addzeros] # add the sign # here we avoid octal numbers by trimming the leading zeros! # 2005-10-28 S.ARNOLD if {$signe} {set mantissa [expr {-[string trimleft $mantissa 0]}]} # the F tags a BigFloat # a BigInt is like any other integer since Tcl 8.5, # because expr now supports arbitrary length integers return [_fromstr $mantissa $exp] } ################################################################################ # private procedure to transform decimal floats into binary ones # IN : # - number : a BigInt representing the Mantissa # - exp : the decimal exponent (a simple integer) # OUT : # $number * 10^$exp, as the internal binary representation of a BigFloat ################################################################################ proc ::math::bigfloat::_fromstr {number exp} { set number [string trimleft $number 0] if {$number==""} { return [list F 0 [expr {int($exp*log(10)/log(2))-15}] [expr {1<<15}]] } if {$exp==0} { return [list F $number 0 1] } if {$exp>0} { # mul by 10^exp, then normalize set power [expr {10**$exp}] set number [expr {$number*$power}] return [normalize [list F $number 0 $power]] } # now exp is negative or null # the closest power of 2 to the 'exp'th power of ten, but greater than it # 10**$exp<2**$binaryExp # $binaryExp>$exp*log(10)/log(2) set binaryExp [expr {int(-$exp*log(10)/log(2))+1+16}] # then compute n * 2^binaryExp / 10^(-exp) # (exp is negative) # equals n * 2^(binaryExp+exp) / 5^(-exp) set diff [expr {$binaryExp+$exp}] if {$diff<0} { error "internal error" } set power [expr {5**(-$exp)}] set number [expr {($number<<$diff)/$power}] set delta [expr {(1<<$diff)/$power}] return [normalize [list F $number [expr {-$binaryExp}] [incr delta]]] } ################################################################################ # fromdouble : # like fromstr, but for a double scalar value # arguments : # double - the number to convert to a BigFloat # exp (optional) - the total number of digits ################################################################################ proc ::math::bigfloat::fromdouble {double {exp {}}} { set mantissa [lindex [split $double e] 0] # line added by SArnold on 05/08/2005 set mantissa [string trimleft [string map {+ "" - ""} $mantissa] 0] set precision [string length [string map {. ""} $mantissa]] if { $exp != {} && [incr exp]>$precision } { return [fromstr $double [expr {$exp-$precision}]] } else { # tests have failed : not enough precision or no exp argument return [fromstr $double] } } ################################################################################ # converts a BigInt into a BigFloat with a given decimal precision ################################################################################ proc ::math::bigfloat::int2float {int {decimals 1}} { # it seems like we need some kind of type handling # very odd in this Tcl world :-( if {![isInt $int]} { error "first argument is not an integer" } if {$decimals<1} { error "non-positive decimals number" } # the lowest number of decimals is 1, because # [tostr [fromstr 10.0]] returns 10. # (we lose 1 digit when converting back to string) set int [expr {$int*10**$decimals}] return [_fromstr $int [expr {-$decimals}]] } ################################################################################ # multiplies 'leftop' by 'rightop' and rshift the result by 'shift' ################################################################################ proc ::math::bigfloat::intMulShift {leftop rightop shift} { return [::math::bignum::rshift [::math::bignum::mul $leftop $rightop] $shift] } ################################################################################ # returns 1 if x is a BigFloat, 0 elsewhere ################################################################################ proc ::math::bigfloat::isFloat {x} { # a BigFloat is a list of : "F" mantissa exponent delta if {[llength $x]!=4} { return 0 } # the marker is the letter "F" if {[string equal [lindex $x 0] F]} { return 1 } return 0 } ################################################################################ # checks that n is a BigInt (a number create by math::bignum::fromstr) ################################################################################ proc ::math::bigfloat::isInt {n} { if {[llength $n]>1} { return 0 } # if {[string is digit $n]} { # return 1 # } return 1 } ################################################################################ # returns 1 if x is null, 0 otherwise ################################################################################ proc ::math::bigfloat::iszero {x} { if {[isInt $x]} { return [expr {$x==0}] } checkFloat $x # now we do some interval rounding : if a number's interval englobs 0, # it is considered to be equal to zero foreach {dummy integer exp delta} $x {break} if {$delta>=abs($integer)} {return 1} return 0 } ################################################################################ # compute log(X) ################################################################################ proc ::math::bigfloat::log {x} { checkFloat $x foreach {dummy integer exp delta} $x {break} if {$integer<=0} { error "zero logarithm error" } if {[iszero $x]} { error "number equals zero" } set precision [bits $integer] # uncertainty of the logarithm set delta [_logOnePlusEpsilon $delta $integer $precision] incr delta # we got : x = 1xxxxxx (binary number with 'precision' bits) * 2^exp # we need : x = 0.1xxxxxx(binary) *2^(exp+precision) incr exp $precision foreach {integer deltaIncr} [_log $integer] {break} incr delta $deltaIncr # log(a * 2^exp)= log(a) + exp*log(2) # result = log(x) + exp*log(2) # as x<1 log(x)<0 but 'integer' (result of '_log') is the absolute value # that is why we substract $integer to log(2)*$exp set integer [expr {[_log2 $precision]*$exp-$integer}] incr delta [expr {abs($exp)}] return [normalize [list F $integer -$precision $delta]] } ################################################################################ # compute log(1-epsNum/epsDenom)=log(1-'epsilon') # Taylor development gives -x -x^2/2 -x^3/3 -x^4/4 ... # used by 'log' command because log(x+/-epsilon)=log(x)+log(1+/-(epsilon/x)) # so the uncertainty equals abs(log(1-epsilon/x)) # ================================================ # arguments : # epsNum IN (the numerator of epsilon) # epsDenom IN (the denominator of epsilon) # precision IN (the number of bits after the dot) # # 'epsilon' = epsNum*2^-precision/epsDenom ################################################################################ proc ::math::bigfloat::_logOnePlusEpsilon {epsNum epsDenom precision} { if {$epsNum>=$epsDenom} { error "number is null" } set s [expr {($epsNum<<$precision)/$epsDenom}] set divider 2 set t [expr {$s*$epsNum/$epsDenom}] set u [expr {$t/$divider}] # when u (the current term of the development) is zero, we have reached our goal # it has converged while {$u!=0} { incr s $u # divider = order of the term = 'n' incr divider # t = (epsilon)^n set t [expr {$t*$epsNum/$epsDenom}] # u = t/n = (epsilon)^n/n and is the nth term of the Taylor development set u [expr {$t/$divider}] } return $s } ################################################################################ # compute log(0.xxxxxxxx) : log(1-epsilon)=-eps-eps^2/2-eps^3/3...-eps^n/n ################################################################################ proc ::math::bigfloat::_log {integer} { # the uncertainty is nbSteps with nbSteps<=nbBits # take nbSteps=nbBits (the worse case) and log(nbBits+increment)=increment set precision [bits $integer] set n [expr {int(log($precision+2*log($precision)))}] set integer [expr {$integer<<$n}] incr precision $n set delta 3 # 1-epsilon=integer set integer [expr {(1<<$precision)-$integer}] set s $integer # t=x^2 set t [expr {$integer*$integer>>$precision}] set denom 2 # u=x^2/2 (second term) set u [expr {$t/$denom}] while {$u!=0} { # while the current term is not zero, it has not converged incr s $u incr delta # t=x^n set t [expr {$t*$integer>>$precision}] # denom = n (the order of the current development term) # u = x^n/n (the nth term of Taylor development) set u [expr {$t/[incr denom]}] } # shift right to restore the precision set delta return [list [expr {$s>>$n}] [expr {($delta>>$n)+1}]] } ################################################################################ # computes log(num/denom) with 'precision' bits # used to compute some analysis constants with a given accuracy # you might not call this procedure directly : it assumes 'num/denom'>4/5 # and 'num/denom'<1 ################################################################################ proc ::math::bigfloat::__log {num denom precision} { # Please Note : we here need a precision increment, in order to # keep accuracy at $precision digits. If we just hold $precision digits, # each number being precise at the last digit +/- 1, # we would lose accuracy because small uncertainties add to themselves. # Example : 0.0001 + 0.0010 = 0.0011 +/- 0.0002 # This is quite the same reason that made tcl_precision defaults to 12 : # internally, doubles are computed with 17 digits, but to keep precision # we need to limit our results to 12. # The solution : given a precision target, increment precision with a # computed value so that all digits of he result are exacts. # # p is the precision # pk is the precision increment # 2 power pk is also the maximum number of iterations # for a number close to 1 but lower than 1, # (denom-num)/denum is (in our case) lower than 1/5 # so the maximum nb of iterations is for: # 1/5*(1+1/5*(1/2+1/5*(1/3+1/5*(...)))) # the last term is 1/n*(1/5)^n # for the last term to be lower than 2^(-p-pk) # the number of iterations has to be # 2^(-pk).(1/5)^(2^pk) < 2^(-p-pk) # log(1/5).2^pk < -p # 2^pk > p/log(5) # pk > log(2)*log(p/log(5)) # now set the variable n to the precision increment i.e. pk set n [expr {int(log(2)*log($precision/log(5)))+1}] incr precision $n # log(num/denom)=log(1-(denom-num)/denom) # log(1+x) = x + x^2/2 + x^3/3 + ... + x^n/n # = x(1 + x(1/2 + x(1/3 + x(...+ x(1/(n-1) + x/n)...)))) set num [expr {$denom-$num}] # $s holds the result set s [expr {($num<<$precision)/$denom}] # $t holds x^n set t [expr {$s*$num/$denom}] set d 2 # $u holds x^n/n set u [expr {$t/$d}] while {$u!=0} { incr s $u # get x^n * x set t [expr {$t*$num/$denom}] # get n+1 incr d # then : $u = x^(n+1)/(n+1) set u [expr {$t/$d}] } # see head of the proc : we return the value with its target precision return [expr {$s>>$n}] } ################################################################################ # computes log(2) with 'precision' bits and caches it into a namespace variable ################################################################################ proc ::math::bigfloat::__logbis {precision} { set increment [expr {int(log($precision)/log(2)+1)}] incr precision $increment # ln(2)=3*ln(1-4/5)+ln(1-125/128) set a [__log 125 128 $precision] set b [__log 4 5 $precision] set r [expr {$b*3+$a}] set ::math::bigfloat::Log2 [expr {$r>>$increment}] # formerly (when BigFloats were stored in ten radix) we had to compute log(10) # ln(10)=10.ln(1-4/5)+3*ln(1-125/128) } ################################################################################ # retrieves log(2) with 'precision' bits ; the result is cached ################################################################################ proc ::math::bigfloat::_log2 {precision} { variable Log2 if {![info exists Log2]} { __logbis $precision } else { # the constant is cached and computed again when more precision is needed set l [bits $Log2] if {$precision>$l} { __logbis $precision } } # return log(2) with 'precision' bits even when the cached value has more bits return [_round $Log2 $precision] } ################################################################################ # returns A modulo B (like with fmod() math function) ################################################################################ proc ::math::bigfloat::mod {a b} { checkNumber $a checkNumber $b if {[isInt $a] && [isInt $b]} {return [expr {$a%$b}]} if {[isInt $a]} {error "trying to divide an integer by a BigFloat"} set quotient [div $a $b] # examples : fmod(3,2)=1 quotient=1.5 # fmod(1,2)=1 quotient=0.5 # quotient>0 and b>0 : get floor(quotient) # fmod(-3,-2)=-1 quotient=1.5 # fmod(-1,-2)=-1 quotient=0.5 # quotient>0 and b<0 : get floor(quotient) # fmod(-3,2)=-1 quotient=-1.5 # fmod(-1,2)=-1 quotient=-0.5 # quotient<0 and b>0 : get ceil(quotient) # fmod(3,-2)=1 quotient=-1.5 # fmod(1,-2)=1 quotient=-0.5 # quotient<0 and b<0 : get ceil(quotient) if {[sign $quotient]} { set quotient [ceil $quotient] } else { set quotient [floor $quotient] } return [sub $a [mul $quotient $b]] } ################################################################################ # returns A times B ################################################################################ proc ::math::bigfloat::mul {a b} { checkNumber $a checkNumber $b # dispatch the command to appropriate commands regarding types (BigInt & BigFloat) if {[isInt $a]} { if {[isInt $b]} { return [expr {$a*$b}] } return [mulFloatByInt $b $a] } if {[isInt $b]} {return [mulFloatByInt $a $b]} # now we are sure that 'a' and 'b' are BigFloats foreach {dummy integerA expA deltaA} $a {break} foreach {dummy integerB expB deltaB} $b {break} # 2^expA * 2^expB = 2^(expA+expB) set exp [expr {$expA+$expB}] # mantissas are multiplied set integer [expr {$integerA*$integerB}] # compute precisely the uncertainty set delta [expr {$deltaA*(abs($integerB)+$deltaB)+abs($integerA)*$deltaB+1}] # we have to normalize because 'delta' may be too big return [normalize [list F $integer $exp $delta]] } ################################################################################ # returns A times B, where B is a positive integer ################################################################################ proc ::math::bigfloat::mulFloatByInt {a b} { checkFloat $a foreach {dummy integer exp delta} $a {break} if {$b==0} { return [list F 0 $exp $delta] } # Mantissa and Delta are simply multplied by $b set integer [expr {$integer*$b}] set delta [expr {$delta*$b}] # We normalize because Delta could have seriously increased return [normalize [list F $integer $exp $delta]] } ################################################################################ # normalizes a number : Delta (accuracy of the BigFloat) # has to be limited, because the memory use increase # quickly when we do some computations, as the Mantissa and Delta # increase together # The solution : limit the size of Delta to 16 bits ################################################################################ proc ::math::bigfloat::normalize {number} { checkFloat $number foreach {dummy integer exp delta} $number {break} set l [bits $delta] if {$l>16} { incr l -16 # $l holds the supplementary size (in bits) # now we can shift right by $l bits # always round upper the Delta set delta [expr {$delta>>$l}] incr delta set integer [expr {$integer>>$l}] incr exp $l } return [list F $integer $exp $delta] } ################################################################################ # returns -A (the opposite) ################################################################################ proc ::math::bigfloat::opp {a} { checkNumber $a if {[iszero $a]} { return $a } if {[isInt $a]} { return [expr {-$a}] } # recursive call lset a 1 [expr {-[lindex $a 1]}] return $a } ################################################################################ # gets Pi with precision bits # after the dot (after you call [tostr] on the result) ################################################################################ proc ::math::bigfloat::pi {precision {binary 0}} { if {![isInt $precision]} { error "'$precision' expected to be an integer" } if {!$binary} { # convert decimal digit length into bit length set precision [expr {int(ceil($precision*log(10)/log(2)))}] } return [list F [_pi $precision] -$precision 1] } # # Procedure that resets the stored cached Pi constant # proc ::math::bigfloat::reset {} { variable _pi0 if {[info exists _pi0]} {unset _pi0} } proc ::math::bigfloat::_pi {precision} { # the constant Pi begins with 3.xxx # so we need 2 digits to store the digit '3' # and then we will have precision+2 bits in the mantissa variable _pi0 if {![info exists _pi0]} { set _pi0 [__pi $precision] } set lenPiGlobal [bits $_pi0] if {$lenPiGlobal<$precision} { set _pi0 [__pi $precision] } return [expr {$_pi0 >> [bits $_pi0]-2-$precision}] } ################################################################################ # computes an integer representing Pi in binary radix, with precision bits ################################################################################ proc ::math::bigfloat::__pi {precision} { set safetyLimit 8 # for safety and for the better precision, we do so ... incr precision $safetyLimit # formula found in the Math litterature (on Wikipedia # Pi/4 = 44.atan(1/57) + 7.atan(1/239) - 12.atan(1/682) + 24.atan(1/12943) set a [expr {[_atanfract 57 $precision]*44}] incr a [expr {[_atanfract 239 $precision]*7}] set a [expr {$a - [_atanfract 682 $precision]*12}] incr a [expr {[_atanfract 12943 $precision]*24}] return [expr {$a>>$safetyLimit-2}] } ################################################################################ # shift right an integer until it haves $precision bits # round at the same time ################################################################################ proc ::math::bigfloat::_round {integer precision} { set shift [expr {[bits $integer]-$precision}] if {$shift==0} { return $integer } # $result holds the shifted integer set result [expr {$integer>>$shift}] # $shift-1 is the bit just rights the last bit of the result # Example : integer=1000010 shift=2 # => result=10000 and the tested bit is '1' if {$integer & (1<<($shift-1))} { # we round to the upper limit return [incr result] } return $result } ################################################################################ # returns A power B, where B is a positive integer ################################################################################ proc ::math::bigfloat::pow {a b} { checkNumber $a if {$b<0} { error "pow : exponent is not a positive integer" } # case where it is obvious that we should use the appropriate command # from math::bignum (added 5th March 2005) if {[isInt $a]} { return [expr {$a**$b}] } # algorithm : exponent=$b = Sum(i=0..n) b(i)2^i # $a^$b = $a^( b(0) + 2b(1) + 4b(2) + ... + 2^n*b(n) ) # we have $a^(x+y)=$a^x * $a^y # then $a^$b = Product(i=0...n) $a^(2^i*b(i)) # b(i) is boolean so $a^(2^i*b(i))= 1 when b(i)=0 and = $a^(2^i) when b(i)=1 # then $a^$b = Product(i=0...n and b(i)=1) $a^(2^i) and 1 when $b=0 if {$b==0} {return 1} # $res holds the result set res 1 while {1} { # at the beginning i=0 # $remainder is b(i) set remainder [expr {$b&1}] # $b 'rshift'ed by 1 bit : i=i+1 # so next time we will test bit b(i+1) set b [expr {$b>>1}] # if b(i)=1 if {$remainder} { # mul the result by $a^(2^i) # if i=0 we multiply by $a^(2^0)=$a^1=$a set res [mul $res $a] } # no more bits at '1' in $b : $res is the result if {$b==0} { return [normalize $res] } # i=i+1 : $a^(2^(i+1)) = square of $a^(2^i) set a [mul $a $a] } } ################################################################################ # converts angles for radians to degrees ################################################################################ proc ::math::bigfloat::rad2deg {x} { checkFloat $x set xLen [expr {-[lindex $x 2]}] if {$xLen<3} { error "number too loose to convert to degrees" } # $rad/Pi=$deg/180 # so result in deg = $radians*180/Pi return [div [mul $x 180] [pi $xLen 1]] } ################################################################################ # retourne la partie entière (ou 0) du nombre "number" ################################################################################ proc ::math::bigfloat::round {number} { checkFloat $number #set number [normalize $number] # fetching integers (or BigInts) from the internal representation foreach {dummy integer exp delta} $number {break} if {$integer==0} { return 0 } if {$exp>=0} { error "not enough precision to round (in round)" } set exp [expr {-$exp}] # saving the sign, ... set sign [expr {$integer<0}] set integer [expr {abs($integer)}] # integer part of the number set try [expr {$integer>>$exp}] # first bit after the dot set way [expr {$integer>>($exp-1)&1}] # delta is shifted so it gives the integer part of 2*delta set delta [expr {$delta>>($exp-1)}] # when delta is too big to compute rounded value ( if {$delta!=0} { error "not enough precision to round (in round)" } if {$way} { incr try } # ... restore the sign now if {$sign} {return [expr {-$try}]} return $try } ################################################################################ # round and divide by 10^n ################################################################################ proc ::math::bigfloat::roundshift {integer n} { # $exp= 10^$n incr n -1 set exp [expr {10**$n}] set toround [expr {$integer/$exp}] if {$toround%10>=5} { return [expr {$toround/10+1}] } return [expr {$toround/10}] } ################################################################################ # gets the sign of either a bignum, or a BitFloat # we keep the bignum convention : 0 for positive, 1 for negative ################################################################################ proc ::math::bigfloat::sign {n} { if {[isInt $n]} { return [expr {$n<0}] } checkFloat $n # sign of 0=0 if {[iszero $n]} {return 0} # the sign of the Mantissa, which is a BigInt return [sign [lindex $n 1]] } ################################################################################ # gets sin(x) ################################################################################ proc ::math::bigfloat::sin {x} { checkFloat $x foreach {dummy integer exp delta} $x {break} if {$exp>-2} { error "sin : not enough precision" } set precision [expr {-$exp}] # sin(2kPi+x)=sin(x) # $integer is now the modulo of the division of the mantissa by Pi/4 # and $n is the quotient foreach {n integer} [divPiQuarter $integer $precision] {break} incr delta $n set d [expr {$n%4}] # now integer>=0 # x = $n*Pi/4 + $integer and $n belongs to [0,3] # sin(2Pi-x)=-sin(x) # sin(Pi-x)=sin(x) # sin(Pi/2+x)=cos(x) set sign 0 switch -- $d { 0 {set l [_sin2 $integer $precision $delta]} 1 {set l [_cos2 $integer $precision $delta]} 2 {set sign 1;set l [_sin2 $integer $precision $delta]} 3 {set sign 1;set l [_cos2 $integer $precision $delta]} default {error "internal error"} } # $l is a list : {Mantissa Precision Delta} # precision --> the opposite of the exponent # 1.000 = 1000*10^-3 so exponent=-3 and precision=3 digits lset l 1 [expr {-([lindex $l 1])}] # the sign depends on the switch statement below #::math::bignum::setsign integer $sign if {$sign} { lset l 0 [expr {-[lindex $l 0]}] } # we insert the Bigfloat tag (F) and normalize the final result return [normalize [linsert $l 0 F]] } proc ::math::bigfloat::_sin2 {x precision delta} { set pi [_pi $precision] # shift right by 1 = divide by 2 # shift right by 2 = divide by 4 set pis2 [expr {$pi>>1}] set pis4 [expr {$pis2>>1}] if {$x>=$pis4} { # sin(Pi/2-x)=cos(x) incr delta set x [expr {$pis2-$x}] return [_cos $x $precision $delta] } return [_sin $x $precision $delta] } ################################################################################ # sin(x) with 'x' lower than Pi/4 and positive # 'x' is the Mantissa - 'delta' is Delta # 'precision' is the opposite of the exponent ################################################################################ proc ::math::bigfloat::_sin {x precision delta} { # $s holds the result set s $x # sin(x) = x - x^3/3! + x^5/5! - ... + (-1)^n*x^(2n+1)/(2n+1)! # = x * (1 - x^2/(2*3) * (1 - x^2/(4*5) * (...* (1 - x^2/(2n*(2n+1)) )...))) # The second expression allows us to compute the less we can # $double holds the uncertainty (Delta) of x^2 : 2*(Mantissa*Delta) + Delta^2 # (Mantissa+Delta)^2=Mantissa^2 + 2*Mantissa*Delta + Delta^2 set double [expr {$x*$delta>>$precision-1}] incr double [expr {1+$delta*$delta>>$precision}] # $x holds the Mantissa of x^2 set x [expr {$x*$x>>$precision}] set dt [expr {$x*$delta+$double*($s+$delta)>>$precision}] incr dt # $t holds $s * -(x^2) / (2n*(2n+1)) # mul by x^2 set t [expr {$s*$x>>$precision}] set denom2 2 set denom3 3 # mul by -1 (opp) and divide by 2*3 set t [expr {-$t/($denom2*$denom3)}] while {$t!=0} { incr s $t incr delta $dt # incr n => 2n --> 2n+2 and 2n+1 --> 2n+3 incr denom2 2 incr denom3 2 # $dt is the Delta corresponding to $t # $double "" "" "" "" $x (x^2) # ($t+$dt) * ($x+$double) = $t*$x + ($dt*$x + $t*$double) + $dt*$double # Mantissa^ ^--------Delta-------------------^ set dt [expr {$x*$dt+($t+$dt)*$double>>$precision}] set t [expr {$t*$x>>$precision}] # removed 2005/08/31 by sarnold75 #set dt [::math::bignum::add $dt $double] set denom [expr {$denom2*$denom3}] # now computing : div by -2n(2n+1) set dt [expr {1+$dt/$denom}] set t [expr {-$t/$denom}] } return [list $s $precision $delta] } ################################################################################ # procedure for extracting the square root of a BigFloat ################################################################################ proc ::math::bigfloat::sqrt {x} { checkFloat $x foreach {dummy integer exp delta} $x {break} # if x=0, return 0 if {[iszero $x]} { # return zero, taking care of its precision ($exp) return [list F 0 $exp $delta] } # we cannot get sqrt(x) if x<0 if {[lindex $integer 0]<0} { error "negative sqrt input" } # (1+epsilon)^p = 1 + epsilon*(p-1) + epsilon^2*(p-1)*(p-2)/2! + ... # + epsilon^n*(p-1)*...*(p-n)/n! # sqrt(1 + epsilon) = (1 + epsilon)^(1/2) # = 1 - epsilon/2 - epsilon^2*3/(4*2!) - ... # - epsilon^n*(3*5*..*(2n-1))/(2^n*n!) # sqrt(1 - epsilon) = 1 + Sum(i=1..infinity) epsilon^i*(3*5*...*(2i-1))/(i!*2^i) # sqrt(n +/- delta)=sqrt(n) * sqrt(1 +/- delta/n) # so the uncertainty on sqrt(n +/- delta) equals sqrt(n) * (sqrt(1 - delta/n) - 1) # sqrt(1+eps) < sqrt(1-eps) because their logarithm compare as : # -ln(2)(1+eps) < -ln(2)(1-eps) # finally : # Delta = sqrt(n) * Sum(i=1..infinity) (delta/n)^i*(3*5*...*(2i-1))/(i!*2^i) # here we compute the second term of the product by _sqrtOnePlusEpsilon set delta [_sqrtOnePlusEpsilon $delta $integer] set intLen [bits $integer] # removed 2005/08/31 by sarnold75, readded 2005/08/31 set precision $intLen # intLen + exp = number of bits before the dot #set precision [expr {-$exp}] # square root extraction set integer [expr {$integer<<$intLen}] incr exp -$intLen incr intLen $intLen # there is an exponent 2^$exp : when $exp is odd, we would need to compute sqrt(2) # so we decrement $exp, in order to get it even, and we do not need sqrt(2) anymore ! if {$exp&1} { incr exp -1 set integer [expr {$integer<<1}] incr intLen incr precision } # using a low-level (taken from math::bignum) root extraction procedure # using binary operators set integer [_sqrt $integer] # delta has to be multiplied by the square root set delta [expr {$delta*$integer>>$precision}] # round to the ceiling the uncertainty (worst precision, the fastest to compute) incr delta # we are sure that $exp is even, see above return [normalize [list F $integer [expr {$exp/2}] $delta]] } ################################################################################ # compute abs(sqrt(1-delta/integer)-1) # the returned value is a relative uncertainty ################################################################################ proc ::math::bigfloat::_sqrtOnePlusEpsilon {delta integer} { # sqrt(1-x) - 1 = x/2 + x^2*3/(2^2*2!) + x^3*3*5/(2^3*3!) + ... # = x/2 * (1 + x*3/(2*2) * ( 1 + x*5/(2*3) * # (...* (1 + x*(2n-1)/(2n) ) )...))) set l [bits $integer] # to compute delta/integer we have to shift left to keep the same precision level # we have a better accuracy computing (delta << lg(integer))/integer # than computing (delta/integer) << lg(integer) set x [expr {($delta<<$l)/$integer}] # denom holds 2n set denom 4 # x/2 set result [expr {$x>>1}] # x^2*3/(2!*2^2) # numerator holds 2n-1 set numerator 3 set temp [expr {($result*$delta*$numerator)/($integer*$denom)}] incr temp while {$temp!=0} { incr result $temp incr numerator 2 incr denom 2 # n = n+1 ==> num=num+2 denom=denom+2 # num=2n+1 denom=2n+2 set temp [expr {($temp*$delta*$numerator)/($integer*$denom)}] } return $result } # # Computes the square root of an integer # Returns an integer # proc ::math::bigfloat::_sqrt {n} { set i [expr {(([bits $n]-1)/2)+1}] set b [expr {$i*2}] ; # Bit to set to get 2^i*2^i set r 0 ; # guess set x 0 ; # guess^2 set s 0 ; # guess^2 backup set t 0 ; # intermediate result for {} {$i >= 0} {incr i -1; incr b -2} { set x [expr {$s+($t|(1<<$b))}] if {abs($x)<= abs($n)} { set s $x set r [expr {$r|(1<<$i)}] set t [expr {$t|(1<<$b+1)}] } set t [expr {$t>>1}] } return $r } ################################################################################ # substracts B to A ################################################################################ proc ::math::bigfloat::sub {a b} { checkNumber $a checkNumber $b if {[isInt $a] && [isInt $b]} { # the math::bignum::sub proc is designed to work with BigInts return [expr {$a-$b}] } return [add $a [opp $b]] } ################################################################################ # tangent (trivial algorithm) ################################################################################ proc ::math::bigfloat::tan {x} { return [::math::bigfloat::div [::math::bigfloat::sin $x] [::math::bigfloat::cos $x]] } ################################################################################ # returns a power of ten ################################################################################ proc ::math::bigfloat::tenPow {n} { return [expr {10**$n}] } ################################################################################ # converts a BigInt to a double (basic floating-point type) # with respect to the global variable 'tcl_precision' ################################################################################ proc ::math::bigfloat::todouble {x} { global tcl_precision set precision $tcl_precision if {$precision==0} { # this is a cheat, I must admit, for Tcl 8.5 set precision 16 } checkFloat $x # get the string repr of x without the '+' sign # please note: here we call math::bigfloat::tostr set result [string trimleft [tostr $x] +] set minus "" if {[string index $result 0]=="-"} { set minus - set result [string range $result 1 end] } set l [split $result e] set exp 0 if {[llength $l]==2} { # exp : x=Mantissa*2^Exp set exp [lindex $l 1] } # caution with octal numbers : we have to remove heading zeros # but count them as digits regexp {^0*} $result zeros incr exp -[string length $zeros] # Mantissa = integerPart.fractionalPart set l [split [lindex $l 0] .] set integerPart [lindex $l 0] set integerLen [string length $integerPart] set fractionalPart [lindex $l 1] # The number of digits in Mantissa, excluding the dot and the leading zeros, of course set integer [string trimleft $integerPart$fractionalPart 0] if {$integer eq ""} { set integer 0 } set len [string length $integer] # Now Mantissa is stored in $integer if {$len>$precision} { set lenDiff [expr {$len-$precision}] # true when the number begins with a zero set zeroHead 0 if {[string index $integer 0]==0} { incr lenDiff -1 set zeroHead 1 } set integer [roundshift $integer $lenDiff] if {$zeroHead} { set integer 0$integer } set len [string length $integer] if {$len<$integerLen} { set exp [expr {$integerLen-$len}] # restore the true length set integerLen $len } } # number = 'sign'*'integer'*10^'exp' if {$exp==0} { # no scientific notation set exp "" } else { # scientific notation set exp e$exp } # place the dot just before the index $integerLen in the Mantissa set result [string range $integer 0 [expr {$integerLen-1}]] append result .[string range $integer $integerLen end] # join the Mantissa with the sign before and the exponent after return $minus$result$exp } ################################################################################ # converts a number stored as a list to a string in which all digits are true ################################################################################ proc ::math::bigfloat::tostr {args} { if {[llength $args]==2} { if {![string equal [lindex $args 0] -nosci]} {error "unknown option: should be -nosci"} set nosci yes set number [lindex $args 1] } else { if {[llength $args]!=1} {error "syntax error: should be tostr ?-nosci? number"} set nosci no set number [lindex $args 0] } if {[isInt $number]} { return $number } checkFloat $number foreach {dummy integer exp delta} $number {break} if {[iszero $number]} { # we do matter how much precision $number has : # it can be 0.0000000 or 0.0, the result is not the same zero #return 0 } if {$exp>0} { # the power of ten the closest but greater than 2^$exp # if it was lower than the power of 2, we would have more precision # than existing in the number set newExp [expr {int(ceil($exp*log(2)/log(10)))}] # 'integer' <- 'integer' * 2^exp / 10^newExp # equals 'integer' * 2^(exp-newExp) / 5^newExp set binExp [expr {$exp-$newExp}] if {$binExp<0} { # it cannot happen error "internal error" } # 5^newExp set fivePower [expr {5**$newExp}] # 'lshift'ing $integer by $binExp bits is like multiplying it by 2^$binExp # but much, much faster set integer [expr {($integer<<$binExp)/$fivePower}] # $integer is the Mantissa - Delta should follow the same operations set delta [expr {($delta<<$binExp)/$fivePower}] set exp $newExp } elseif {$exp<0} { # the power of ten the closest but lower than 2^$exp # same remark about the precision set newExp [expr {int(floor(-$exp*log(2)/log(10)))}] # 'integer' <- 'integer' * 10^newExp / 2^(-exp) # equals 'integer' * 5^(newExp) / 2^(-exp-newExp) set binShift [expr {-$exp-$newExp}] set fivePower [expr {5**$newExp}] # rshifting is like dividing by 2^$binShift, but faster as we said above about lshift set integer [expr {$integer*$fivePower>>$binShift}] set delta [expr {$delta*$fivePower>>$binShift}] set exp -$newExp } # saving the sign, to restore it into the result set result [expr {abs($integer)}] set sign [expr {$integer<0}] # rounded 'integer' +/- 'delta' set up [expr {$result+$delta}] set down [expr {$result-$delta}] if {($up<0 && $down>0)||($up>0 && $down<0)} { # $up>0 and $down<0 or vice-versa : then the number is considered equal to zero set isZero yes # delta <= 2**n (n = bits(delta)) # 2**n <= 10**exp , then # exp >= n.log(2)/log(10) # delta <= 10**(n.log(2)/log(10)) incr exp [expr {int(ceil([bits $delta]*log(2)/log(10)))}] set result 0 } else { # iterate until the convergence of the rounding # we incr $shift until $up and $down are rounded to the same number # at each pass we lose one digit of precision, so necessarly it will success for {set shift 1} { [roundshift $up $shift]!=[roundshift $down $shift] } { incr shift } {} incr exp $shift set result [roundshift $up $shift] set isZero no } set l [string length $result] # now formatting the number the most nicely for having a clear reading # would'nt we allow a number being constantly displayed # as : 0.2947497845e+012 , would we ? if {$nosci} { if {$exp >= 0} { append result [string repeat 0 $exp]. } elseif {$l + $exp > 0} { set result [string range $result 0 end-[expr {-$exp}]].[string range $result end-[expr {-1-$exp}] end] } else { set result 0.[string repeat 0 [expr {-$exp-$l}]]$result } } else { if {$exp>0} { # we display 423*10^6 as : 4.23e+8 # Length of mantissa : $l # Increment exp by $l-1 because the first digit is placed before the dot, # the other ($l-1) digits following the dot. incr exp [incr l -1] set result [string index $result 0].[string range $result 1 end] append result "e+$exp" } elseif {$exp==0} { # it must have a dot to be a floating-point number (syntaxically speaking) append result . } else { set exp [expr {-$exp}] if {$exp < $l} { # we can display the number nicely as xxxx.yyyy* # the problem of the sign is solved finally at the bottom of the proc set n [string range $result 0 end-$exp] incr exp -1 append n .[string range $result end-$exp end] set result $n } elseif {$l==$exp} { # we avoid to use the scientific notation # because it is harder to read set result "0.$result" } else { # ... but here there is no choice, we should not represent a number # with more than one leading zero set result [string index $result 0].[string range $result 1 end]e-[expr {$exp-$l+1}] } } } # restore the sign : we only put a minus on numbers that are different from zero if {$sign==1 && !$isZero} {set result "-$result"} return $result } ################################################################################ # PART IV # HYPERBOLIC FUNCTIONS ################################################################################ ################################################################################ # hyperbolic cosinus ################################################################################ proc ::math::bigfloat::cosh {x} { # cosh(x) = (exp(x)+exp(-x))/2 # dividing by 2 is done faster by 'rshift'ing return [floatRShift [add [exp $x] [exp [opp $x]]] 1] } ################################################################################ # hyperbolic sinus ################################################################################ proc ::math::bigfloat::sinh {x} { # sinh(x) = (exp(x)-exp(-x))/2 # dividing by 2 is done faster by 'rshift'ing return [floatRShift [sub [exp $x] [exp [opp $x]]] 1] } ################################################################################ # hyperbolic tangent ################################################################################ proc ::math::bigfloat::tanh {x} { set up [exp $x] set down [exp [opp $x]] # tanh(x)=sinh(x)/cosh(x)= (exp(x)-exp(-x))/2/ [(exp(x)+exp(-x))/2] # =(exp(x)-exp(-x))/(exp(x)+exp(-x)) # =($up-$down)/($up+$down) return [div [sub $up $down] [add $up $down]] } # exporting public interface namespace eval ::math::bigfloat { foreach function { add mul sub div mod pow iszero compare equal fromstr tostr fromdouble todouble int2float isInt isFloat exp log sqrt round ceil floor sin cos tan cotan asin acos atan cosh sinh tanh abs opp pi deg2rad rad2deg } { namespace export $function } } # (AM) No "namespace import" - this should be left to the user! #namespace import ::math::bigfloat::* package provide math::bigfloat 2.0.1