## Sunday, March 29, 2009

### Generating Circular Primes

I've been having way too much fun with Project Euler (and Haskell) lately. Below is a listing of my solution to problem 35: "How many circular primes are there below a million?"

It solves in about 2 seconds on my laptop (using -O2). This isn't great, but not bad given that the only optimizations in the code are: a memoized prime number generator, and the use of Data.Set for quick primality tests.

`-- Project Euler: Problem 35---- The number, 197, is called a circular prime because all rotations of the-- digits: 197, 971, and 719, are themselves prime.---- There are thirteen such primes below 100: 2, 3, 5, 7, 11, 13, 17, 31, 37,-- 71, 73, 79, and 97.---- How many circular primes are there below one million?import Data.Set (Set)import qualified Data.Set as Setimport System (getArgs)-- How many circular primes do you want?limit = 1000000-- Memoized prime number generator.primes :: [Integer]primes = 2 : 3 : 5 : filter (not . hasFactor) [7..]  where hasFactor n = any (divides n) \$ takeWhile (<= lowestFactor n) primes        divides n m = n `mod` m == 0        lowestFactor = ceiling . sqrt . fromIntegral-- Given an integer, return the number of digits in it.digits :: Integer -> Intdigits n = ceiling \$ logBase 10 \$ fromIntegral (n + 1)-- Rotate number right:--   rotate 456 = 645--   rotate 3498 = 8349---- Note that this does not work for numbers that have zeroes in them. But-- that's okay because any number with a zero in it is not a circular prime.-- Also, right rotation catches the zero before it disappears into the -- significant digit.rotate :: Integer -> Integerrotate n = strip_last_digit + (last_digit * tens_place)  where last_digit = n `mod` 10        strip_last_digit = (n - last_digit) `div` 10        tens_place = (10 ^ (digits n - 1))-- Set of million primes. This enables fast lookup in isCircularPrime.primeSet :: Set IntegerprimeSet = Set.fromList \$ takeWhile (< limit) primes-- Returns True if number is prime.isPrime :: Integer -> BoolisPrime n = Set.member n primeSet-- Returns True if the number is a circular prime.isCircularPrime :: Integer -> BoolisCircularPrime n = all isPrime \$ take (digits n) \$ iterate rotate n-- Go.main :: IO ()main = print \$ length \$ filter isCircularPrime \$ takeWhile (< limit) primes`

So... how can I optimize this some more?

1. Well, you could adjust your prime number generator to generate two or more digit prime numbers only containing the digits 1, 3, 7 and 9. Otherwise the rotated numbers will be divisible by 2 or 5.

2. Divisibility by 3 and 11 should be done separately.

http://home.comcast.net/~babdulbaki/Circular_Primes.html

3. @ts That's a pretty good observation. Thanks!

4. @bassam Thanks for the pointer.

5. I've got this sitting around my \$HOME, possibly useful. Also, nice mach-o article.

(looks like the indentation is gonna be butchered..)

{-# OPTIONS_GHC -O2 #-}

import Data.List(foldl')
import Data.Monoid(mempty)
import Data.IntMap(IntMap) --Data.Map(Map)
import qualified Data.IntMap as M --Data.Map as M
import System.Environment(getArgs)

type Z = Int --Integer
type ZMap = IntMap [Z] --Map Z [Z]

main = print
. (sieve !!)

sieve :: [Z]
sieve = go [2..] mempty
where go [] _ = []
go (x:xs) m = maybe
(x : go xs (insert x m))
(go xs . fixup x m)
(M.lookup x m)

insert :: Z -> ZMap -> ZMap
insert x m = M.insert (x*x) [x] m

fixup :: Z -> ZMap -> [Z] -> ZMap
fixup x m = flip foldl' (M.delete x m)
(\m p -> M.insertWith (++) (x+p) [p] m)