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 Set

import 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 -> Int

digits 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 -> Integer

rotate 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 Integer

primeSet = Set.fromList $ takeWhile (< limit) primes

-- Returns True if number is prime.

isPrime :: Integer -> Bool

isPrime n = Set.member n primeSet

-- Returns True if the number is a circular prime.

isCircularPrime :: Integer -> Bool

isCircularPrime 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?

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.

ReplyDeleteDivisibility by 3 and 11 should be done separately.

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

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

ReplyDelete@bassam Thanks for the pointer.

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

ReplyDelete(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 !!)

. read . head =<< getArgs

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)