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.0

The 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:

- If the fraction can be expressed in a finite number of digits, then the above code will give the
*exact*answer only when that number is lower than 15 or 16. That can easily be helped by extending the code for this case a bit :) - There is only one place where the decimal character of the result is determined - so, the script can be adapted to generate fractions in other number systems without even frowning your eye brows. The line in question:

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..