A Mathematical Exploration with Haskell

August 16, 2014

When people ask me what I like to do for fun and I say, "Just studying math and computer science", they usually give me a funny look and wander off shaking their heads.

So here's an attempt to explain why this such an enjoyable pastime for me:

I was working today on this problem from Project Euler 
A unit fraction contains 1 in the numerator. The decimal representation of the unit fractions with denominators 2 to 10 are given:

    1/2 =   0.5
    1/3 =   0.(3)
    1/4 =   0.25
    1/5 =   0.2
    1/6 =   0.1(6)
    1/7 =   0.(142857)
    1/8 =   0.125
    1/9 =   0.(1)
    1/10  =   0.1

Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be seen that 1/7 has a 6-digit recurring cycle.

Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.
Hmm, where to start on that one? 

The first thing to consider is that a floating point representation (like on a calculator) is not going to give us enough accuracy.  Instead, let's use Haskell Integers which can be as big as we want.  If we take a big power of 10 and divide it by our integer value, we'll get back an integer with lots of digits that will help us spot the patterns.

Aside: If you're new to Haskell, I'd recommend starting with the excellent book "Learn You a Haskell for Great Good!"

Here's an example:
> ghci

Prelude> 10^32 `div` 7
Ok, no problem spotting the pattern there.  It's 6 digits long.

But what if the pattern was 100 digits long.  Counting the lengths of the patterns by hand could get to be a pain.

So to start off with we need some way to get Haskell to spot recurring patterns in a list of elements.  Here's an idea -- how about if we try constructing different groupings of the digits and look for the smallest sized one where each of the groups (except possibly the last one) match the first group.

So here's a simple function to help us with that
groupsOf :: Eq(a) => Integer -> [a] -> [[a]]
groupsOf n [] = []
groupsOf n xs =
  let n' = fromIntegral n
  in (take n' xs) : groupsOf n (drop n' xs)
Now we can convert the number to a list of characters and group it like this:
groupsOf 5 . show . div (10^32) $ 7
No way. 5 is not the right grouping -- let's try going up to 6...
groupsOf 6 . show . div (10^32) $ 7
Cool, that works.  So now, assuming that our repeating pattern is always going to start at the beginning of the list (can we assume that?  more on this later...) we can write a function that tells us if our groups match up.
matches :: Eq(a) => [a] -> [a] -> Bool
matches xs ys = xs == ys || lys < lxs && take lys xs == ys
  where (lxs, lys) = (length xs, length ys)
This should do just what we want.  It basically says that two lists match if either they're equal or the second list is shorter and matches the first part of the first list.

So now we can combine 'groupsOf' and 'matches' to write a function that gets us closer to what we want:
isPeriodicWith :: Eq(a) => Integer -> [a] -> Bool
isPeriodicWith n xs =
  let (g:gs) = groupsOf n xs
  in all (matches g) gs
Let's try it out ...
isPeriodicWith 5 (show . div (10^32) $ 7)
isPeriodicWith 6 (show . div (10^32) $ 7)
Ok, but we don't want to have to check isPeriodicWith for each possibility, so let's write a function that does that for us.
period :: Eq(a) => [a] -> Maybe Integer
period xs =
  let lmax = fromIntegral $ length xs `div` 2
      period' n mx xs
        | n > mx = Nothing
        | otherwise = if isPeriodicWith n xs then Just n else period' (n+1) mx xs
  in period' 1 lmax xs

period (show . div (10^32) $ 7)
Just 6
period [1,3,4,3,5,4,6,4,32,2,5,6,7,4]
Nice!  If it's not periodic, it returns Nothing, otherwise it returns Just (the period).

Next let's make a helper function that takes a power to raise 10 (which will control how many digits we'll get in the expansion) and the number we want to test.  It returns a list of the digits obtained by dividing that big power of 10 by our number.
intInvertedDecimal :: Integer -> Integer -> [Char]
intInvertedDecimal p n =
  let p' = fromIntegral p
  in show . div (10^p') $ n 
Let's try it out.
intInvertedDecimal 32 7

period $ intInvertedDecimal 32 7
Just 6
Ok, great, just like before.  But what if we try it with the inverse of 16?
intInvertedDecimal 32 16

period $ intInvertedDecimal 32 16
Hmm... So we can't really assume that the repeating part of the expansion starts with the first number in our quotient.  But maybe we can -- if we exclude any multiples of 2 and 5.  We know those will have a repeating expansion of zeros (or nines) with period 1, so let's remove them from consideration.  Will all the other numbers work the way we want?  Let's find out...  Here's a function that lets us specify an upper bound and returns a list of all the 'non-trivial' integers up to that bound.
nonTrivials :: Integer -> [Integer]
nonTrivials u = [x | x <- [3..(u-1)], x `rem` 2 /= 0, x `rem` 5 /= 0]
Let's try putting this all together.  We've got some awesome functions that will tell us if a number has a reciprocal cycle and what its period is.  The confusing part is that we don't know up front how many digits we're going to need to look at before we can spot the pattern.  We could use some huge number, but that probably wouldn't be too efficient. 

If we were going to work it out by hand for, say 23, we'd probably start with a small number of digits like 10.  So we do this:
period $ intInvertedDecimal 10 23
Not enough digits, ok so let's try 10^2 digits.
period $ intInvertedDecimal 100 23
Just 22
Great, if 10^2 hadn't worked, we could have tried 10^3 and so on.  So we can try something like this:
periodForNumber :: Integer -> Integer
periodForNumber n =
    let periodForNumber' p n =
          case period (intInvertedDecimal (10^p) n) of
              Nothing -> periodForNumber' (p+1) n
              (Just per) -> per
    in periodForNumber' 1 n
Let's try it
periodForNumber 11
periodForNumber 23
periodForNumber 7
Ok, seems to work.  Let's try mapping it over some non-trivial integers.  Let's use pairs like (number, period) so we can keep track of which numbers we're looking at.
map (\x -> (x, periodForNumber x)) $ nonTrivials 100
Seems reasonable -- Let's get write a function that gives us the max pair for an upper bound
maxPeriod :: Integer -> (Integer, Integer)
maxPeriod u =
  let pairs = map (\x -> (x, periodForNumber x)) $ nonTrivials u
  in maximumBy (\(_,p1) (_,p2) -> compare p1 p2) $ pairs

maxPeriod 1000
So there's our solution to the Project Euler problem that started us off. 
Yup, it's correct (Hey -- I was the 42,000th person to solve it there -- cool!)

But I notice something funny... do you see it?  983 is the number (less than 1000) whose reciprocal cycle has the longest period.  It's period is 982.  One less.  That seems vaguely familiar.  Oh yeah, 7 had a period of 6 and 23 had a period of 22.  Interesting, I wonder if a number could ever have a period greater than or equal to itself.  Let's see...
filter (\(n, p) -> p >= n) . map (\x -> (x, periodForNumber x)) . nonTrivials $ 1000
Hmm, absence of a counter-example is not a proof, but here's a conjecture:

Conjecture 1: For any given number n, the length of the repeating decimal expansion of 1/n cannot exceed n-1

So I'm wondering whether there could be something special about these numbers with cycles of length one less than themselves.  Let's tak a look at them.  First, though, let's give them a name:  Define 'maxExps' to be { n <- N | the repeating decimal expansion of 1/n has length equal to n-1 }
let maxExps = map fst . filter (\(n, p) -> p == n-1) . map (\x -> (x, periodForNumber x)) . nonTrivials $ 1000
Interesting, they all look like they could be primes.  Let's check using the Sieve of Eratosthenes.  We can implement it quickly like this:
sieve :: [Integer] -> [Integer]
sieve [] = []
sieve (x:xs) = x:(sieve . filter (\y -> y `rem` x /= 0) $ xs)

let primes = sieve [2..1000]
If you get tired of checking each element of maxExps against the above list of primes, you can do this
:m + Data.Set
maxExps \\ primes
Okay, so they're all primes.  So here's the next conjecture...
Conjecture 2: maxExps is a subset of the primes

This is getting interesting.  If we can prove this conjecture, we have discovered a way of finding primes without having to do any factoring.  All we have to do is divide two numbers and count the digits in the repeating expansion. 

I'm also wondering how much of this behavior is special to the fact that we're looking at decimal expansions versus some other base, say octal.

Aside: At this point, you may be asking yourself -- "Hey is this for real?  I mean, are we really figuring this stuff out or was this all discovered a thousand years ago by some ancient Greek?"  Here's my honest answer: "Yes, and probably yes".  On one hand, we are becoming aware of interesting facets of a landscape that we've never visited before.  It would be cool to think that we might be the first explorers to stumble upon this place, but that would be pretty improbable given the zillions of smart mathematician types who have come before us.  So we probably aren't going to get our names on any theorems or anything, but one nice thing about math is that no matter how many feet have trodden on our path before us, the landscape remains as pristine and beautiful as it was for the first explorer...

Ok, back to work.  Let's try doing the same sort of thing in octal (base 8).  So how would that work?  Well, I'm thinking that instead of dividing a power of 10 by our number, we'll want to divide a power of 8 by it.  And then when we go looking for patterns, we'll want to look at an octal representation of the quotient so it will only use the symbols 0 - 7 to represent the number.

By importing a couple more modules: Numeric and Data.Char, we can define a more general function that takes a base, a power to raise it to (to control how many digits we get in the result) and our number -- and returns a list of digits in the base that we specified.
intInvertedForBase :: Integer -> Integer -> Integer -> [Char]
intInvertedForBase base p n =
  let result = (base^p) `div` n
  in showIntAtBase base intToDigit result ""
It's just a more general version of our decimal-specific function
intInvertedDecimal :: Integer -> Integer -> [Char]
intInvertedDecimal p n =
  let p' = fromIntegral p
  in show . div (10^p') $ n 
.. but just in a more general way.  So how can we swap it in to our existing system?  Well, we could re-write our code to use this new function, but instead let's make this function an argument.   Let's change our periodForNumber function so it takes an additional argument -- the function to use to convert a number into a list of symbols.  In Haskell and several other languages with a functional streak, functions are first class objects so we can pass them around to swap out algorithms.
periodForNumber :: (Integer -> Integer -> [Char]) -> Integer -> Integer
periodForNumber f n =
    let periodForNumber' p n =
          case period (f (10^p) n) of
              Nothing -> periodForNumber' (p+1) n
              (Just per) -> per
    in periodForNumber' 1 n
We'll need to make some corresponding changes to the consumers of periodForNumber
maxPeriod :: (Integer -> Integer -> [Char]) -> Integer -> (Integer, Integer)
maxPeriod f u =
  let pairs = map (\x -> (x, periodForNumber f x)) $ nonTrivials u
  in maximumBy (\(_,p1) (_,p2) -> compare p1 p2) $ pairs
To generalize our list of nonTrivials to handle an arbitrary base, we can do this:
nonTrivialsForBase :: Integer -> Integer -> [Integer]
nonTrivialsForBase base u =
  let divs = nub $ divisors base
      prd x = all (\z -> rem x z /= 0)
  in [x | x <- [2..(u-1)], prd x divs]

divisors :: Integer -> [Integer]
divisors n =
  let divisors' _ 1 = []
      divisors' start n =
        let (q, r) = quotRem n start
        in if r == 0 then start:divisors' start q
           else divisors' (start+1) n
  in divisors' 2 n
Ok, let's test out our generalized versions for base 10 to make sure we didn't break anything...
let maxExps = map fst . filter (\(n, p) -> p == n-1) . map (\x -> (x, periodForNumber (intInvertedForBase 10) x)) . nonTrivialsForBase 10 $ 1000

Great! Now let's see what we get for base 8
let maxExps = map fst . filter (\(n, p) -> p == n-1) . map (\x -> (x, periodForNumber (intInvertedForBase 8) x)) . nonTrivialsForBase 8 $ 1000

Awesome!  We can see that we're still getting a subset of primes, but it's different!
It appears that 3, 5 and 11 are now MaxExps, while 7, 17, 19 and 23 aren't anymore.
13 isn't a MaxExp for either base and 59 is for both bases.  Wow, now those are what I call 'odd' numbers!

Alright, one final step before we get started on the proofs.  Let's make a nicer function to get the maxExps for a given base so we can play around with different bases
maxExpsForBase :: Integer -> Integer -> [Integer]
maxExpsForBase base u =
  let pairs = map (\x -> (x, periodForNumber (intInvertedForBase base) x)) . nonTrivialsForBase base $ u
  in map fst . filter (\(n, p) -> p == n-1) $ pairs

maxExpsForBase 6 1000
Nice!  The (un-polished) code is at Github  if you feel like playing around with it.


Prove or provide a counterexample for each of the following statements:

Let n be a natural number greater than 1 and let b be a base for representing numbers.  Then if n and b are relatively prime (no common divisors except 1) the reciprocal cycle of 1/n begins with the first non-zero digit of the quotient.

Conjecture 1: For any given number n, the length of the repeating decimal expansion of 1/n cannot exceed n-1

Conjecture 2: maxExps is a subset of the primes
where 'maxExps' is defined to be { n <- N | the repeating decimal expansion of 1/n has length equal to n-1 }


Your feedback is welcome! If you find any errors in this post or have any additional pointers or insights, please take a moment to register and share your thoughts.

No Comments Yet

You must be logged in to comment

All Posts