mathr / blog / #

Adam Majewski set me a challenge, to find the period of the recurring decimal digits of a rational number with large numerator and denominator. In base 10, some fractions are periodic: 1/3 = 0.(3) has a period of 1 and 1/7 = 0.(142857) has a period of 6. Other fractions are preperiodic: 1/4 = 0.25(0) has a preperiod of 2 and a period of 1, and 1/14 = 0.0(714285) has a preperiod of 1 and a period of 6.

The same behaviours occur in other bases, though the preperiod and period will be different in general. For example in binary (base 2), 1/3 = 0.(01), 1/7 = 0.(001), 1/4 = 0.01(0), and 1/14 = 0.0(001).

Extracting the digits of a fraction in [0,1) is quite straightforward: keep multiplying by the base and chopping off the integer part. In the Haskell programming language:

fractionDigits :: Int -> Rational -> [Int]
fractionDigits base fraction
  = map fst
  . drop 1
  . iterate (\(_, x) -> properFraction (fromIntegral base * x))
  $ (0, fraction)

This gives an infinite stream of digits, the next challenge is to split it into preperiodic and periodic parts. It turns out that a fraction has a non-empty preperiodic part only when the denominator is divisible by a factor of the base. A bit of maths shows that we only need to check each distinct prime factor of the base, so we can check for preperiodicity like this:

divides :: Integer -> Integer -> Bool
factor `divides` number = number `mod` factor == 0

isPreperiodic :: Int -> Rational -> Bool
isPreperiodic base fraction = or
  [ fromIntegral factor `divides` denominator fraction
  | factor <- nub (primeFactors base)
  ]

Now we can find the preperiod by repeatedly multiplying by base and keeping the fractional part, until the fractional part is not preperiodic:

shift :: Int -> Rational -> Rational
shift base fraction = (fromIntegral base * fraction) `mod'` 1

preperiod :: Int -> Rational -> Int
preperiod base fraction
  = length
  . takeWhile (isPreperiodic base)
  . iterate (shift base)
  $ fraction

To find the period of a known periodic fraction, we need to find the number of shifts it takes to get back to the starting fraction. To ensure we have a periodic fraction, we drop the preperiodic part:

period :: Int -> Rational -> Int
period base fraction
  = (1 +) . length . takeWhile (/= fraction') $ fractions
  where
    fraction':fractions
      = drop (preperiod base fraction)
      . iterate (shift base)
      $ fraction

So far we've been working with fractions in [0,1). Fractions outside that range have a non-zero integer part, so we need a way to extract its digits. We do this by repeatedly dividing by the base and collecting the remainders, which gives the digits from smallest to largest.

integerDigits :: Int -> Integer -> [Int]
integerDigits base integer
  = reverse
  . map (fromInteger . snd)
  . takeWhile1 ((0 /=) . fst)
  . drop 1
  . iterate (\(i, _) -> i `divMod` fromIntegral base)
  $ (integer, 0)

takeWhile1 :: (a -> Bool) -> [a] -> [a]
takeWhile1 predicate list
  = result ++ [next]
  where
    (result, next:_) = span predicate list

So far we've been getting each digit as an Int, it remains to convert it to a character - we'll limit to maximum base 62 for now:

showDigit :: Int -> Char
showDigit digit = (['0'..'9'] ++ ['A'..'Z'] ++ ['a'..'z']) !! digit

showDigits :: [Int] -> String
showDigits = map showDigit

Finally we can put all this together to convert arbitrary fractions to strings:

showAtBase :: Int -> Rational -> String
showAtBase base rational
  =  sign
  ++ integerPart
  ++ "."
  ++ preperiodicPart
  ++ "("
  ++ periodicPart
  ++ ")"
  where
    (preperiodicPart, periodicPart')
      = splitAt
          (preperiod base fraction)
          (showDigits (fractionDigits base fraction))
    periodicPart = take (period base fraction) periodicPart'
    integerPart = showDigits (integerDigits base integer)
    (integer, fraction) = properFraction (abs rational)
    sign = if rational < 0 then "-" else ""

You can download the full source code, which uses the primes package.

Exercise: write readAtBase :: Int -> String -> Rational that converts the string representation back to the corresponding fraction.

EDIT: you can download the source code of an optimized version that inlines various loops into a single traversal.

EDIT2: I stopped the optimized version after 90+ hours of runtime, after finding a better way in IRC:

freenode/##math

14:01 < i_c-Y> ClaudiusMaximus: one can show that the period is the multiplicative order of 10 mod n where the fraction is m/n.

14:01 < i_c-Y> see http://mathworld.wolfram.com/DecimalExpansion.html for details

The large fraction in question

33877456965431938318210482471113262183356704085033125021829876006886584214655562 / 237142198758023568227473377297792835283496928595231875152809132048206089502588927

turns out to have a rather huge period in base 10:

794564201485273000257607338237654476912493997529945960250807965815440

calculated by Wolfram Alpha.

The naive program I wrote would probably never have finished in my lifetime, and even then it would have given wrong output due to Int overflow:

Prelude> 794564201485273000257607338237654476912493997529945960250807965815440 :: Int

<interactive>:7:1: Warning:
    Literal 794564201485273000257607338237654476912493997529945960250807965815440 is out of the Int range -9223372036854775808..9223372036854775807
-7502784436604661104

At least the negative result would have made the overflow error obvious!