KBK says:
I liked RS's fraction procedures in Bag of algorithms, but I didn't like the comment that "precision must be specified." I saw it as a challenge to come up with something that would do:
% num2fract 0.125 1/8 % num2fract 0.44 11/25 % num2fract 0.3333333333 1/3 % num2fract 0.5714285714 4/7
In short, I wanted to come as close as possible to converting ANY float to a fraction. This actually isn't as hard as might seem. The following code does it reasonably well (and is much improved [2002-01-11] over an earlier version that lived on the Wiki for over a year).
# Convert a floating point number to a mixed number, e.g.,
# 3.75 -> 3 3/4. The optional arg "maxint" is the largest
# integer that will appear in the denominator.
proc num2mixed { float { maxint 2147483647 } } {
set numer [expr { int( $float ) }]
set denom [expr { $float - $numer }]
return [list $numer [num2fract $denom $maxint]]
}
# Convert a floating point number to a fraction, e.g.,
# 3.375 -> 27/8. The optional arg "maxint" is the largest
# integer that will appear in the fraction, e.g.,
# [num2fract 3.14159 1000] -> 355/113
proc num2fract { float { maxint 2147483647 } } {
return [join [quotient_rep $float $maxint] "/"]
}KBK (11 January 2002) - I've replaced quotient_rep with an 'improved' version, which, alas, doesn't use the tricky [catch]. The old version is at Tricky catch.
# Convert a floating point number to a pair of integers
# that are the numerator and denominator of its fractional
# representation. The optional arg "maxint" is the largest
# integer that will appear in the fraction, e.g.,
# [quotient_rep 3.14159 1000] -> {355 113}
proc quotient_rep { num { maxint 2147483647 } } {
set i 0
set p_i 1
set q_i 0
set p_im1 0
set q_im1 1
while { 1 } {
incr i
set a_i [expr { int($num) }]
set fract [expr { $num - $a_i }]
if { 1.0 * $a_i * $p_i + $p_im1 > $maxint } {
break
}
if { 1.0 * $a_i * $q_i + $q_im1 > $maxint } {
break
}
foreach { p_i p_im1 } \
[list [expr { $a_i * $p_i + $p_im1 }] $p_i] break
foreach { q_i q_im1 } \
[list [expr { $a_i * $q_i + $q_im1 }] $q_i] break
if { abs( $fract * $maxint ) < 1 } {
break
}
set num [expr { 1.0 / $fract }]
}
return [list $p_i $q_i]
}Lars H: The two foreach commands above are unnecessary. An alternative method for doing the same thing is
set p_i [expr {$p_im1 + $a_i*[set p_im1 $p_i]}]
set q_i [expr {$q_im1 + $a_i*[set q_im1 $q_i]}]which works because expr looks up the value of e.g. p_im1 for $p_im1 before the set p_im1 $p_i is evaluated.
The reader who wants to know the maths behind the above (regular continued fractions and how this algorithm computes them) can take a look at [L1 ], where the necessary formulae are derived from elementary mathematics. It is an excerpt from the source of a PS type 1 font linker I've written. (Type 1 fonts probably provide one of the few practical applications of algorithms such as the above, since the Type 1 charstring format [L2 ] offers no other way of encoding a non-integer than as the fraction of two integers.)
KBK Exercise (for the mathematically inclined): Prove that the [while] loop terminates no later than $i == 41 and find an input to [quotient_rep] for which the performance is that bad. Fraction Math - Answer to exercise
Arjen Markus I have not yet analysed the above (not even the original one, though both are intriguing), but I do know a two-step approximation that does seem nice as well. I will post it here as soon as convenient. The only drawback: no guarantee that is "best", only that it is better than the trivial first step.
MS Note that this algorithm does produce an approx p/q that is "best" approximator to num in the following sense: no other rational r/s with s<=q is closer. But it is not "best" in the sense: ''no other rational involving integers <= $maxint is closer": consider the counterexample
quotient_rep 3.1416305 500 --> 355/113
However, the fraction 377/120 is a better approximation with both numerator and denominator <=500.
KBK If you're curious what's going on here, and where else you can go with it, check out the wonderful musings at [L3 ]! MS Some of the claims in there seem (at this point, to me) to be wrong - look at Notes on continued fractions for counterexamples and notes. [This paragraph has been slightly edited.]
RS: Good job, as far I've seen it under time pressure. But just to justify why I had the mandatory denominator: I wanted to emulate the behavior on fractions of e.g. inches or stock quotes, where denominator was always a power of two. I think 2-1/3" is not really idiomatic ;-)
AM: Richard's original code might be useful to get a nice-looking decimal representation - my problem is: how to determine the number of decimals that will accurately represent a bunch of numbers with the least number of decimals:
1.05000, 1.10000 will not do - I want: 1.05, 1.10 for instance.
JFR: Try this:
set return_val [format "%0.9f" $return_val]
set return_val [string trimright $return_val 0]
set return_val [string trimright $return_val .]DKF: Just fixed a minor documentation typo in num2fract.
Now that we've got code to convert arbitrary decimal expansions to fractions, we can put together some simple quotient arithmetic code. Note that it tries to give you the same kind of fractions out as you put in. Cannot handle negative numbers, decimal input or very large quotients. Take this and build on it!
namespace eval quotientmath {
# First, some code to handle "p q/r" and "p/q" input and output
# formats. I'll use a list format internally.
proc Qscan {fraction} {
set bits [split $fraction " "]
if {[llength $bits] < 1 || [llength $bits] > 2} {
return -code error "\"$fraction\" is not a fraction"
} elseif {[llength $bits] == 1} {
set int 0
set kind 0
} elseif {![string is integer [lindex $bits 0]]} {
return -code error "\"$fraction\" is not a fraction"
} else {
set int [lindex $bits 0]
set fraction [lindex $bits 1]
set kind 1
}
set bits [split $fraction "/"]
if {
[llength $bits] != 2 ||
![string is integer [lindex $bits 0]] ||
![string is integer [lindex $bits 1]]
} then {
return -code error "\"$fraction\" is not a fraction"
}
set numer [lindex $bits 0]
set denom [lindex $bits 1]
list $kind [list [expr {$numer+$int*$denom}] $denom]
}
proc Qformat {kind fraction} {
foreach {numer denom} $fraction {}
if {$kind && $numer>$denom} {
set int [expr {$numer/$denom}]
set numer [expr {$numer%$denom}]
return [format "%d %d/%d" $int $numer $denom]
} else {
return [format "%d/%d" $numer $denom]
}
}
# Now the basic arithmetic operators.
proc Qmult {a b} {
foreach {an ad} $a {}; foreach {bn bd} $b {}
list [expr {$an*$bn}] [expr {$ad*$bd}]
}
proc Qdiv {a b} {
foreach {an ad} $a {}; foreach {bn bd} $b {}
list [expr {$an*$bd}] [expr {$ad*$bn}]
}
proc Qadd {a b} {
foreach {an ad} $a {}; foreach {bn bd} $b {}
list [expr {$an*$bd+$bn*$ad}] [expr {$ad*$bd}]
}
proc Qsub {a b} {
foreach {an ad} $a {}; foreach {bn bd} $b {}
list [expr {$an*$bd-$bn*$ad}] [expr {$ad*$bd}]
}
# A normalisation function
proc Qnorm {x} {
foreach {p q} $x {}
# Handle canonical forms
if {$p == 0} {return [list 0 1]}
if {$q == 0} {return [list 1 0]}
# Calculate GCD
while {1} {
if {$q == 0} {
# Termination case
break
} elseif {$p>$q} {
# Swap p and q
set t $p
set p $q
set q $t
}
set q [expr {$q%$p}]
}
# Divide both numerator and denominator by the GCD
foreach {xn xd} $x {}
list [expr {$xn/$p}] [expr {$xd/$p}]
}
# The public functions
proc add {x y} {
foreach {kx x} [Qscan $x] {}
foreach {ky y} [Qscan $y] {}
Qformat [expr {$kx||$ky}] [Qnorm [Qadd $x $y]]
}
proc sub {x y} {
foreach {kx x} [Qscan $x] {}
foreach {ky y} [Qscan $y] {}
Qformat [expr {$kx||$ky}] [Qnorm [Qsub $x $y]]
}
proc mult {x y} {
foreach {kx x} [Qscan $x] {}
foreach {ky y} [Qscan $y] {}
Qformat [expr {$kx||$ky}] [Qnorm [Qmult $x $y]]
}
proc div {x y} {
foreach {kx x} [Qscan $x] {}
foreach {ky y} [Qscan $y] {}
Qformat [expr {$kx||$ky}] [Qnorm [Qdiv $x $y]]
}
namespace export add sub mult div
}RS built on it (or reinvented the wheel) when playing with rationals, and has this tiny "rationalizer":
proc dbl2frac {dbl {eps 0.001}} {
for {set den 1} {$den<1024} {incr den} {
set num [expr {round($dbl*$den)}]
if {abs(double($num)/$den - $dbl) < $eps} break
}
list $num $den
}
% dbl2frac 0.333
1 3
% dbl2frac 1.23
16 13
% dbl2frac -0.42
-13 31AM (21 june 2004) I used the quotientmath module that Donal implemented with modifications to take care of negative numbers to produce the so-called Bernoulli numbers. These numbers arise in combinatorics, the gamma function and sundry other places. There exists a simple recurrence relation:
B0 = 1
n /n+1\
Sum | | Bk = 0
k=0 \ k /Because of the large numbers involved this is difficult to get numerically stable, unless you use symbolic manipulation or indeed rational numbers. Here is the script (it implicitly uses my changes to the above module and I did not have the math module from Tcllib available, so I reprogrammed the binomial coefficients):
# bernoulli.tcl --
# Compute the first N Bernoulli numbers:
#
# n /n+1\
# sum | | B = 0, n = 1, ...
# k=0 \ k / k
# binomial_coeffs --
# Determine the binomial coefficients for
# the given argument
#
# Arguments:
# number The power of the binomial
# Result:
# List of coefficients
#
proc binomial_coeffs {number} {
if { $number == 1 } {
return {1 1}
} else {
set prev [binomial_coeffs [expr {$number-1}]]
set result {}
foreach p1 [concat 0 $prev] p2 [concat $prev 0] {
lappend result [expr {$p1+$p2}]
}
}
return $result
}
# sum_bernoulli_nm1 --
# Compute the sum of the first n-1 numbers
#
# Arguments:
# bernoulli List of the first n-1 Bernoulli numbers
# Result:
# Sum of the numbers (as a rational number)
#
proc sum_bernoulli_nm1 {bernoulli} {
set coeffs [binomial_coeffs [expr {[llength $bernoulli]+1}]]
set sum "0/1"
foreach c [lrange $coeffs 0 end-2] b $bernoulli {
set sum [add $sum [mult $b $c/1]]
}
return $sum
}
# main --
# Compute the first N Bernoulli numbers
#
source quotientmath.tcl
namespace import quotientmath::*
set b0 {1/1}
set bernoulli [list $b0]
for { set i 1 } { $i < 50 } { incr i } {
set bi [div [sum_bernoulli_nm1 $bernoulli] [expr {-($i+1)}]/1]
puts "New number $i: $bi"
update
lappend bernoulli $bi
}The result:
New number 1: -1/2 New number 2: 1/6 New number 3: 0/1 New number 4: -1/30 New number 5: 0/1 New number 6: 1/42 New number 7: 0/1 New number 8: -1/30 New number 9: 0/1 New number 10: 5/66 New number 11: 0/1 New number 12: -691/2730 New number 13: 0/1 New number 14: 7/6 New number 15: 0/1 New number 16: -3617/510 New number 17: 0/1 New number 18: 43867/798 New number 19: 0/1 New number 20: -174611/330 New number 21: 0/1 New number 22: 854513/138
Two remarks:
DKF: Here's a version that's focused on just the computation of Bernoulli numbers; it's quite happy to calculate really large values in 8.5 or 8.6 (provided you've got the memory, of course).
proc bernoulli {n} {
for {set m 0} {$m <= $n} {incr m} {
lappend A [list 1 [expr {$m + 1}]]
for {set j $m} {[set i $j] >= 1} {} {
lassign [lindex $A [incr j -1]] a1 b1
lassign [lindex $A $i] a2 b2
set x [set p [expr {$i * ($a1*$b2 - $a2*$b1)}]]
set y [set q [expr {$b1 * $b2}]]
while {$q} {set q [expr {$p % [set p $q]}]}
lset A $j [list [expr {$x/$p}] [expr {$y/$p}]]
}
}
return [lindex $A 0]
}Though this is going to be quite fast, it's not as clear what the code is doing as I might've hoped.