Posts Tagged ‘primes’

Programming Praxis – Prime Partitions

October 19, 2012

In today’s Programming Praxis exercise, our goal is to calculate the number of prime partitions for a given number. Let’s get started, shall we?

Some imports:

import Data.List
import qualified Data.Map as M
import Data.Numbers.Primes

First we need to calculate the sum of the unique prime factors of a given number. The trivial way is to to take the list of prime factors, filter out the unique numbers and sum those. In the provided solution this is marked as the wrong way since there is a more efficient solution. On the other hand, we only need to this once each for the numbers 1 through 1000, so it’s not like it’s a massive performance bottleneck. Since performance doesn’t matter much, I personally prefer the straightforward, easy to understand implementation.

sopf :: Integer -> Integer
sopf = sum . nub . primeFactors

At first glance, the algorithm for the amount of prime partitions bears some resemblance to the Fibonacci function. The difference here is that each term refers to all of its predecessors instead of a fixed number, which means the typical zipWith solution won’t work. Instead, we build op a dictionary where for each number we store sopf(n) and k(n) so we don’t have to recalculate them.

k :: Integer -> Integer
k n = snd $ foldl' calc M.empty [1..n] M.! n where
    calc dict i = M.insert i (s, div (s + sum (map (\j -> (fst $ dict M.! j) *
        (snd $ dict M.! (i - j))) [1..i-1])) i) dict where s = sopf i

A test to see if everything is working properly:

:: IO ()
main = print $ k 1000 == 48278613741845757

Programming Praxis – Fractran

July 6, 2012

In today’s Programming Praxis exercise, our goal is to write an interpreter for the esotreric programming language Fractran and use it to execute a program that generates prime numbers. Let’s get started, shall we?

A quick import:

import Data.List

The fractran interpreter itself is pretty simple. We use tuples to represent the fractions rather than Ratios since this made for slightly more elegant code. Unfortunately the Prelude doesn’t have a function to swap the values of a tuple, or I could’ve use divMod rather than having to repeat the multiplication.

fractran :: [(Integer, Integer)] -> Integer -> [Integer]
fractran xs = unfoldr (\m -> fmap ((,) m) . lookup 0 $
    map (\(n,d) -> (mod (m*n) d, div (m*n) d)) xs)

Since we’re using tuples, we represent the primegame program as a zip.

primegame :: [(Integer, Integer)]
primegame = zip [17,78,19,23,29,77,95,77, 1,11,13,15,15,55]
                [91,85,51,38,33,29,23,19,17,13,11,14, 2, 1]

Finding the primes is a matter of finding the terms in the output of fractran primegame 2 that are powers of two.

primes :: [(Integer, Integer)]
primes = [(i,e) | (i,n) <- zip [1..] $ fractran primegame 2
                , let e = round $ log (fromIntegral n) / log 2, 2^e == n]

A test to see if everything is working properly:

main :: IO ()
main = mapM_ print $ take 26 primes

Programming Praxis – Billboard Challenge, Part 1

June 22, 2012

In today’s Programming Praxis exercise, our goal is to solve the problem posed to potential employees by a tech firm a couple of years ago: find the first 10-digit prime in the digits of e. Let’s get started, shall we?

A quick import:

import Data.List

First we reuse the unbounded spigot algorithm for calculating e from the last exercise;

stream :: Integral a => (a, a) -> (a, a, a) -> [(a, a, a)] -> [a]
stream (lo, hi) z ~(x:xs) = if lbound == approx z hi
    then lbound : stream (lo, hi) (mul (10, -10*lbound, 1) z) (x:xs)
    else stream (lo, hi) (mul z x) xs
    where lbound = approx z lo
          approx (a,b,c) n = div (a*n + b) c
          mul (a,b,c) (d,e,f) = (a*d, a*e + b*f, c*f)

streamDigits :: (Num a, Integral a, Enum a) => (a -> (a, a, a)) -> (a, a) -> [a]
streamDigits f range = stream range (1,0,1) [(n, a*d, d) | (n,d,a) <- map f [1..]]

stream_e :: [Integer]
stream_e     = streamDigits (\k -> (1, k, 1)) (1, 2)

Next we need a function to test if a number is prime. In the interest of simplicity we just use trial division. To speed things up, we use the three simplest optimizations: only testing odd numbers, only testing up to the square root of the tested number, and only testing with primes. The list of all primes is defined at root level so Haskell memoizes it.

primes :: [Integer]
primes = 2 : filter prime [3,5..]

prime :: Integer -> Bool
prime n = all ((> 0) . mod n) $ takeWhile (\p -> p * p <= n) primes

With those two building blocks out of the way, the rest is trivial: take all groups of 10 consecutive digits, convert them to numbers and find the first prime.

main :: IO ()
main = print $ find prime . map (foldl1 ((+) . (10 *)) . take 10) $ tails stream_e

Programming Praxis – Mersenne Primes

June 3, 2011

In today’s Programming Praxis exercise, our goal is to calculate the Mersenne prime exponents up to 256. Let’s get started, shall we?

A quick import:

import Data.Numbers.Primes

The Mersenne primes can be determine with a simple list comprehension. Although the Lucas-Lehmer test officially doesn’t work for M2, in practice it works just fine so there is no need to make it a special case.

mersennes :: [Int]
mersennes = [p | p <- primes,
    iterate (\n -> mod (n^2 - 2) (2^p - 1)) 4 !! p-2 == 0]

A test shows that the algorithm is working correctly. Piece of cake.

main :: IO ()
main = print $ takeWhile (<= 256) mersennes
               == [2,3,5,7,13,17,19,31,61,89,107,127]

Programming Praxis – Sieve Of Euler

February 25, 2011

In today’s Programming Praxis exercise, our goal is to implement the Sieve of Euler to generate primes. Let’s get started, shall we?

Some imports:

import Data.List.Ordered
import Data.Time.Clock

Euler’s sieve is pretty simple: just filter out all remaining multiples of each successive prime number. We add the minor optimization of only considering odd numbers.

primes_euler :: [Integer]
primes_euler = 2 : euler [3,5..] where
    euler ~(x:xs) = x : euler (minus xs $ map (* x) (x:xs))

To compare, we use the Sieve of Eratosthenes implementation found here:

primes_eratosthenes :: [Integer]
primes_eratosthenes = 2 : 3 : sieve (tail primes_eratosthenes) 3 [] where
    sieve ~(p:ps) x fs = let q=p*p in
        foldr (flip minus) [x+2,x+4..q-2] [[y+s,y+2*s..q] | (s,y) <- fs]
        ++ sieve ps q ((2*p,q):[(s,q-rem (q-y) s) | (s,y) <- fs])

A quick and dirty benchmarking function:

benchPrimes :: Num a => [a] -> IO ()
benchPrimes f = do start <- getCurrentTime
                   print . sum $ take 10000 f
                   print . flip diffUTCTime start =<< getCurrentTime

Time to run the benchmarks.

main :: IO ()
main = do benchPrimes primes_eratosthenes
          benchPrimes primes_euler

The sieve of Eratosthenes runs in about 0.02 seconds, while Euler’s sieve takes around 6.6 seconds, which means this implementation of Euler’s Sieve is about 300-400 times slower. Clearly you don’t want to use this version in practice.

Programming Praxis – Rowland’s Prime-Generating Function

November 12, 2010

In today’s Programming Praxis exercise our goal is to implement a few functions that create integer sequences related to prime generation. Let’s get started, shall we?

As in the Scheme solution, the function names refer to their number in the On-Line Encyclopedia of Integers.

The first sequence can be written in a manner similar to the typical Fibonacci sequence solution.

a106108 :: [Integer]
a106108 = 7 : zipWith (\n r -> r + gcd n r) [2..] a106108

The next two just build on the previous one and are fairly trivial.

a132199 :: [Integer]
a132199 = zipWith (-) (tail a106108) a106108

a137613 :: [Integer]
a137613 = filter (/= 1) a132199

The shortcut function is slightly longer. We store the sum of all the previous numbers so that the function runs in O(n2) rather than O(n3).

shortcut :: [Integer]
shortcut = f 0 1 where
    f s n = r : f (s + r) (n + 1) where r = lpd (6 - n + s)
    lpd n = head $ filter (\d -> mod n d == 0) [3,5..]

Some tests to see if we’re getting the correct output:

main :: IO ()
main = do print $ take 65 a106108
          print $ take 104 a132199
          print $ take 72 a137613
          print $ take 72 shortcut

Yup. Piece of cake.

Programming Praxis – Emirps

November 2, 2010

In today’s Programming Praxis, our task is to enumerate all the non-palindrome prime numbers that are still prime when their digits are reversed. Let’s get started, shall we?

Since we are to focus on maximum speed, the obvious choice is to use a library.

import Data.Numbers.Primes

I tried a couple different versions, and settled on a simple filter. This version can calculate all the emirps under one million in about 650 ms. One other version where the isPrime test was replaced by first putting all the primes in a Set and doing membership testing was a fraction faster (about 10 ms), but it is about twice as long and requires specifying the maximum up front, which is less convenient then the current infinite stream. I also tried to make it faster via parallelism, but I could only get it to go slower, so it looks like the benefits of doing everything in parallel don’t weigh up against the additional overhead. Or maybe I’m just not using the library right.

emirps :: [Integer]
emirps = [p | p <- primes, let rev = reverse (show p)
            , isPrime (read rev), show p /= rev]

All that’s left is to print all the emirps under a million.

main :: IO ()
main = print $ takeWhile (< 1000000) emirps

As mentioned, just evaluating the list takes about 650 ms, which is about the same as the Scheme version.

Programming Praxis – Goldbach’s Conjecture

March 2, 2010

In today’s Programming Praxis problem we have to test Golbach’s conjecture that every even number greater than 2 can be written as the sum of two primes. Let’s get started.

A quick import (note that both the Numbers and primes packages define this module, but only Numbers contains the isPrime function we’re using):

import Data.Numbers.Primes

The function for testing Goldbach’s conjecture is pretty simple: just look at all the primes p less than n and check if n – p is prime.

goldbach :: Integer -> (Integer, Integer)
goldbach n = head [(p, n - p) | p <- takeWhile (< n) primes, isPrime (n - p)]

A quick test reveals everything’s working correctly.

main :: IO ()
main = do print $ goldbach 28
          print $ goldbach 986332

The Scheme solution implements two versions, citing speed as a reason. I also implemented his second version, but the speed was nearly identical to the algorithm above (about 11 seconds to determine that 532 is indeed the highest first prime in the first one million numbers), and given the choice I’d much rather have one line than thirteen.