[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}]
       }
    }
 }
if 0 {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}
   }
 }

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

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.
----
[Category Mathematics] | [Arts and crafts of Tcl-Tk programming]