*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

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:

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

So now we can combine 'groupsOf' and 'matches' to write a function that gets us closer to what we want:

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.

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:

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

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 }

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.

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

### Homework:

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 }

#### Feedback

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:Hmm, where to start on that one?

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.

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.

Here's an example:

> ghciOk, no problem spotting the pattern there. It's 6 digits long.

Prelude> 10^32 `div` 7

14285714285714285714285714285714

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]]Now we can convert the number to a list of characters and group it like this:

groupsOf n [] = []

groupsOf n xs =

let n' = fromIntegral n

in (take n' xs) : groupsOf n (drop n' xs)

groupsOf 5 . show . div (10^32) $ 7No way. 5 is not the right grouping -- let's try going up to 6...

["14285","71428","57142","85714","28571","42857","14"]

groupsOf 6 . show . div (10^32) $ 7Cool, 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.

["142857","142857","142857","142857","142857","14"]

matches :: Eq(a) => [a] -> [a] -> BoolThis 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.

matches xs ys = xs == ys || lys < lxs && take lys xs == ys

where (lxs, lys) = (length xs, length ys)

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] -> BoolLet's try it out ...

isPeriodicWith n xs =

let (g:gs) = groupsOf n xs

in all (matches g) gs

isPeriodicWith 5 (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.

False

isPeriodicWith 6 (show . div (10^32) $ 7)

True

period :: Eq(a) => [a] -> Maybe IntegerNice! If it's not periodic, it returns Nothing, otherwise it returns Just (the period).

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]

Nothing

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]Let's try it out.

intInvertedDecimal p n =

let p' = fromIntegral p

in show . div (10^p') $ n

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

"14285714285714285714285714285714"

period $ intInvertedDecimal 32 7

Just 6

intInvertedDecimal 32 16Hmm... 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.

"6250000000000000000000000000000"

period $ intInvertedDecimal 32 16

Nothing

nonTrivials :: Integer -> [Integer]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.

nonTrivials u = [x | x <- [3..(u-1)], x `rem` 2 /= 0, x `rem` 5 /= 0]

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 23Not enough digits, ok so let's try 10^2 digits.

Nothing

period $ intInvertedDecimal 100 23Great, if 10^2 hadn't worked, we could have tried 10^3 and so on. So we can try something like this:

Just 22

periodForNumber :: Integer -> IntegerLet's try it

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

periodForNumber 11Ok, 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.

2

periodForNumber 23

22

periodForNumber 7

6

map (\x -> (x, periodForNumber x)) $ nonTrivials 100Seems reasonable -- Let's get write a function that gives us the max pair for an upper bound

[(3,1),(7,6),(9,1),(11,2),(13,6),(17,16),(19,18),(21,6),(23,22),(27,3),(29,28),(31,15),(33,2),(37,3),(39,6),(41,5),(43,21),(47,46),(49,42),(51,16),(53,13),(57,18),(59,58),(61,60),(63,6),(67,33),(69,22),(71,35),(73,8),(77,6),(79,13),(81,9),(83,41),(87,28),(89,44),(91,6),(93,15),(97,96),(99,2)]

maxPeriod :: Integer -> (Integer, Integer)So there's our solution to the Project Euler problem that started us off.

maxPeriod u =

let pairs = map (\x -> (x, periodForNumber x)) $ nonTrivials u

in maximumBy (\(_,p1) (_,p2) -> compare p1 p2) $ pairs

maxPeriod 1000

(983,982)

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 $ 1000Hmm, 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 $ 1000Interesting, they all look like they could be primes. Let's check using the Sieve of Eratosthenes. We can implement it quickly like this:

maxExps

[7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499,503,509,541,571,577,593,619,647,659,701,709,727,743,811,821,823,857,863,887,937,941,953,971,977,983]

sieve :: [Integer] -> [Integer]If you get tired of checking each element of maxExps against the above list of primes, you can do this

sieve [] = []

sieve (x:xs) = x:(sieve . filter (\y -> y `rem` x /= 0) $ xs)

let primes = sieve [2..1000]

primes

[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997]

:m + Data.SetOkay, so they're all primes. So here's the next conjecture...

maxExps \\ primes

[]

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.

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]It's just a more general version of our decimal-specific function

intInvertedForBase base p n =

let result = (base^p) `div` n

in showIntAtBase base intToDigit result ""

intInvertedDecimal :: Integer -> Integer -> [Char].. 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.

intInvertedDecimal p n =

let p' = fromIntegral p

in show . div (10^p') $ n

periodForNumber :: (Integer -> Integer -> [Char]) -> Integer -> IntegerWe'll need to make some corresponding changes to the consumers of periodForNumber

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

maxPeriod :: (Integer -> Integer -> [Char]) -> Integer -> (Integer, Integer)To generalize our list of nonTrivials to handle an arbitrary base, we can do this:

maxPeriod f u =

let pairs = map (\x -> (x, periodForNumber f x)) $ nonTrivials u

in maximumBy (\(_,p1) (_,p2) -> compare p1 p2) $ pairs

nonTrivialsForBase :: Integer -> Integer -> [Integer]Ok, let's test out our generalized versions for base 10 to make sure we didn't break anything...

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

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

maxExps

[7,17,19,23,29,47,59,61,97,109,113,131,149,167,179,181,193,223,229,233,257,263,269,313,337,367,379,383,389,419,433,461,487,491,499,503,509,541,571,577,593,619,647,659,701,709,727,743,811,821,823,857,863,887,937,941,953,971,977,983]

let maxExps = map fst . filter (\(n, p) -> p == n-1) . map (\x -> (x, periodForNumber (intInvertedForBase 8) x)) . nonTrivialsForBase 8 $ 1000Awesome! We can see that we're still getting a subset of primes, but it's different!

maxExps

[3,5,11,29,53,59,83,101,107,131,149,173,179,197,227,269,293,317,347,389,419,443,461,467,491,509,557,563,587,653,659,677,701,773,797,821,827,941,947]

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]Nice! The (un-polished) code is at Github if you feel like playing around with it.

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

[11,13,17,41,59,61,79,83,89,103,107,109,113,127,131,137,151,157,179,199,223,227,229,233,251,257,271,277,347,367,373,397,401,419,443,449,467,487,491,521,563,569,587,593,613,641,659,661,683,709,733,757,761,809,823,827,829,853,857,881,929,947,953,967,971,977,991]

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 PostsAged Man Consulting

dowdrake@msn.com

3754 SE Madison St.

Portland, OR 97214

Phone: (503) 236-8881

my Bitbucket page

my GitHub page

- Ruby-Doc.org
- Programming Ruby.org
- Learn Rails by Example
- Objects on Rails
- Ruby on Rails API
- Rails Guides
- Ruby Rogues
- Rubies in the Rough
- Ryan Bates Railscasts

- the polymath project
- Galaxy Zoo
- Fold it
- GenBank
- Innocentive
- Mathworks competition
- Sloan Digital Sky Survey
- Open Dinosaur Project
- ArXiv
- PLOS
- Terence Tao Blog