Richard Suchenwirth 2002-04-27 - Periodic decimal fractions are numbers where a sequence of digits behind the decimal point (the period) is endlessly repeated, for example:
1/7 = 0.142857142857.. 1/3 = 0.3333..
The following routine tries to detect a period in a given number and returns the period (might be 0 for integers or non-strict periods like
1/2 = 0.5000..
or (implicitly) an empty string if no period could be detected - then the input number might be irrational (not representable by a nominator/divisor expression, e.g. sqrt(2)), or it has a period longer than 7 digits, which can not be confirmed at the maximum tcl_precision of 17, a limit imposed by the underlying double representation in C. For instance, 1/17 is certainly periodic, but the period is out of sight for Tcl's expr math..
proc period x {
set frac [expr {abs(double($x)-int($x))}]
if {!$frac || [string length $frac]<10} {return 0}
set digits [string range $frac 2 end]
foreach n {1 2 3 4 5 6 7} {
foreach offset {0 1 2 3 4 5 6} {
set try [string range $digits $offset [expr $offset+$n-1]]
if {[regexp .{0,$n}$try\($try)+.{0,$n}$ $digits]} {
return $try
}
}
}
}Arjen Markus As the fraction is contained in a string, there ought to be no problem with very long periods:
set long_fraction "0.123456789012345678901234567890"
As long, of course, as you avoid interpreting the string as a number! RS admits that he does, by letting expr compute the fractional part...
KBK In the spirit of Fraction math, let's do the period of a rational number whose numerator and denominator are integers. The following has no trouble finding out that 1/17 repeats with a period of 16 places. It also handles Arjen's problem nicely, if we remember that .12345678901234567890... = 137174210/1111111111
proc period2 { n d } {
set x 0
while { 1 } {
set r [expr { $n % $d }]
if { $r == 0 } {
return 0
} elseif { [info exists seen($n)] } {
return [expr { $x - $seen($n) }]
} else {
set seen($n) $x
incr x
set n [expr { 10 * $r }]
}
}
}
for { set i 1 } { $i <= 20 } { incr i } {
puts [list $i [period2 1 $i]]
}There is some interesting number theory here, but I haven't the time to get into it; I've typed this in in about 15 minutes while waiting for a test to run.
Arjen Markus Independently of Kevin, I came up with the following proc that will return the decimal string/number of such a division with any precision you like (also known as long division, if I have translated this correctly):
# decimalFraction --
# Determine the decimal fraction of num/denom by a given precision
#
# Arguments:
# num Numerator of the fraction
# denom Denominator of the fraction
# nodec Number of decimals required
# Result:
# The decimal fraction (as a string!)
# Notes:
# The method relies entirely on integer division and can achieve
# any precision required.
#
proc decimalFraction { num denom nodec } {
#
# The integer part first
#
set fract [expr {int($num/$denom)}]
set rem [expr {$num%$denom}]
set result "$fract"
append result "."
#
# Loop over the decimals
#
for { set i 0 } { $i < $nodec } { incr i } {
set rem [expr {$rem*10}]
if { $rem < $denom } {
append result "0"
} else {
set fract [expr {int($rem/$denom)}]
set rem [expr {$rem%$denom}]
append result "$fract"
}
}
return $result
}
#
# A few test cases
#
puts [decimalFraction 1 3 10]
puts [decimalFraction 1 5 10]
puts [decimalFraction 1 13 40]
puts [decimalFraction 1 97 200]
puts [decimalFraction 1 15 10]
puts [decimalFraction 1 17 40]
puts [decimalFraction 1 7 40]The number-theoretical aspects of this problem relate (probably among other things, but I am only an amateur :-) to Fermat's little theorem:
a^(p-1)-1 is divisible by p,
if p is a prime and a is not a multiple of p.
This allows us to put an upper limit to the period of any prime, say 127, which is one of my favourites: the period for 127 can not be more than 126, and otherwise it will be a factor of 126). The latter follows from the observation that:
a^gh-1 = m^h-1 = (m-1)*(1+m+m^2+...m^h) (where m=a^g)
Yet another take by RS 2003-07-03, which needs no explicit bound:
In math books, the periodic repeating part is marked by a bar over those digits. For ASCII convenience, I insert a "p" before the period instead:
1/7 = 0.p142857
The end needs not to be marked, as it's just the end of the digit sequence - nothing else can come after that till infinity...
The following codes mimicks the way humans do division: append integer quotient to result, continue with rest*10. To be able to detect a period, it keeps track of the rests seen so far, and terminates when one repeats (which is guaranteed to happen, as for division by n there are only n-1 possible rests, plus 0 which marks the trivial non-periodic case - if num is a multiple of den, or the prime factors of den contain nothing else than 2 and 5, in our decimal system).
In periodic cases, it returns the "p" result and the list of rests. As only one digit is processed each time, this algorithm can work out periods far longer than tcl_precision allows, if you give it enough time :)
/p 1 61 0.p016393442622950819672131147540983606557377049180327868852459
proc /p {num0 den} {
set digits {}
set rests {}
set num $num0
while 1 {
set digit [expr {$num/$den}]
set rest [expr {$num%$den}]
if {$rest == 0} {
#-- trivial case, non-periodic
return [expr double($num0)/$den]
}
lappend digits $digit
set pos [lsearch $rests $rest]
if {$pos>=0} {
set res [expr $num0/$den].
append res [join [lrange $digits 1 $pos] ""] p
append res [join [lrange $digits [expr $pos+1] end] ""]
return [list $res $rests]
} else {
lappend rests $rest
set num [expr {$rest*10}]
}
}
}Cyclic numbers are the periods of primes that have only one (except for the trivial 0) of length n-1. The period repeats, but from different start points:
foreach i {1 2 3 4 5 6 7} {
puts "$i/7 = [/p $i 7]"
}
1/7 = 0.p142857 {1 3 2 6 4 5}
2/7 = 0.p285714 {2 6 4 5 1 3}
3/7 = 0.p428571 {3 2 6 4 5 1}
4/7 = 0.p571428 {4 5 1 3 2 6}
5/7 = 0.p714285 {5 1 3 2 6 4}
6/7 = 0.p857142 {6 4 5 1 3 2}
7/7 = 1.0The following code helps to detect cyclic numbers - it tortured my iPaq for a while until it reported that between 1900 and 2000, only 1913, 1949, 1979 make cyclic numbers:
proc try {from max} {
if {$from%2==0} {incr from}
for {set i $from} {$i<=$max} {incr i 2} {
set t [lindex [/p 1 $i] 1]
if {[llength $t]==$i-1} {puts $i; update}
}
}See also Remainder orbits
AM (10 july 2003) I have two remarks:
set rest [expr {$rest*10}]TV From merely a quick glance at the subject, rational numbers come to mind, not as mathematical prove mechanism but as practically dealing with simple enough arithmetical operations, where a number is then represented as {n d} where n is the numerator and d the denominator. Like with complex numbers, you'd have for instance addition:
expr ($n1*$d2 + $n2*$d1 ) / ($d1 * d2)
Multiplication:
expr $n1 * $n2 / ($d1 * $d2)
Of course you'd want a 'simplify' or 'remove common multiplier' routine which through integer largest common denominator division leaves you with a perfect but simpler answer.
For non-rational numbers of any accuracy (using integers and if needed infinite accuracy (intermedeate) results), this lets you compute complex things without continuously worrying about exonents and mantissa cutting. I graphics the technique can be used on vectors, to define a perspective projection thoughout the graphics computation pipeline. - RS: Rational math is long done - see Generic operations. TV Good! Combinations, too? I don't see that on that page, but I wouldn't know at all anymore how to do gcd in a decent and effective way... And all that with matrix computations, like a reasonable reordering gaussian sweep? I did that sort of stuff as a student to make a electronics circuit simulator on a bbc electron, that was terrific, 32 kb, and a complex matrix inversion of 30x30 you can wait for. In in decent tcl you could make proof systems with it for instance.
KBK thinks that some of TV's questions are answered in Fraction math. For combinations, see math::combinatorics in tcllib. The procedure does C(n,k) exactly up to C(34,17), and approximately (in floating point, using Lanczos's approximation to log-Gamma thereafter. (It should now probably be reworked to take advantage of 8.4's wide integers - any volunteers?) Factorial and beta distribution are there also. Reasonably careful attention has been paid to the numerical properties.
gcd is easy. My favorite implementation (fast and terse):
proc gcd { p q } {
while { $q != 0 } {
foreach { p q } [list $q [expr { $p % $q }]] break
}
return [expr { abs($p) }] ;# can be just return $p if
;# you don't need to do
;# negative numbers
}(See Greatest common denominator for more details. KPV)
There are also a few notes on continued fractions, including code to compute the nearest rational to a quadratic surd having denominator less than a specified value.
TV Cool gcd, I guess I want to spend some time browsing over tcllib man pages..