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