Posts Tagged ‘primality’

Programming Praxis – Proving Primality

February 2, 2010

In today’s Programming Praxis exercise we have to implement an algorithm to prove the primality of a number. Let’s get started, shall we?

Some imports:

import Data.Bits
import Data.List

the tdFactors and expm functions have both been featured in previous exercises.

tdFactors :: Integer -> [Integer]
tdFactors n = f n [2..floor . sqrt $ fromIntegral n] where
    f _ []     = []
    f r (x:xs) | mod r x == 0 = x : f (div r x) (x:xs)
               | otherwise    = f r xs

expm :: Integer -> Integer -> Integer -> Integer
expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
             filter (flip testBit 0 . snd) .
             zip (iterate (flip mod m . (^ 2)) b) .
             takeWhile (> 0) $ iterate (`shiftR` 1) e

Unfortunately, these math-heavy problems typically resist shortening, so this is pretty similar to the Scheme solution.

isPrime :: Integer -> Bool
isPrime n = f 2 $ tdFactors (n - 1) where
    f _ []     = False
    f b (q:qs) | expm b (n - 1)         n /= 1 = f (b + 1) (q:qs)
               | expm b (n - 1 `div` q) n /= 1 = True
               | otherwise                     = f b qs

All that’s left is a quick test:

main :: IO ()
main = print . isPrime $ 2^89 - 1

Programming Praxis – Monte Carlo factorization

June 19, 2009

In today‚Äôs Programming Praxis problem we have to implement John Pollard’s factorization algorithm. Our target is 16 lines (I’m not counting the code for the primality test, since we did that already).

First, we’re going to need to reuse the code from the Miller-Rabin primality test exercise, since we need to determine whether or not the number is prime:

import Control.Arrow
import Data.Bits
import Data.List
import System.Random

isPrime :: Integer -> StdGen -> Bool
isPrime n g =
    let (s, d) = (length *** head) . span even $ iterate (`div` 2) (n-1)
        xs = map (expm n d) . take 50 $ randomRs (2, n - 2) g
    in all (\x -> elem x [1, n - 1] ||
                  any (== n-1) (take s $ iterate (expm n 2) x)) xs

expm :: Integer -> Integer -> Integer -> Integer
expm m e b = foldl' (\r (b', _) -> mod (r * b') m) 1 .
             filter (flip testBit 0 . snd) .
             zip (iterate (flip mod m . (^ 2)) b) $
             takeWhile (> 0) $ iterate (`shiftR` 1) e

There’s not much more to the factor function than simply writing down the algorithm. It’s just math with Haskell syntax, really. The only thing you have to take care of is not to calculate gcd (x-y) n when x and y are still 2, since that will give you an incorrect result.

factor :: Integer -> Integer -> Integer
factor c n = factor' 2 2 1 where
    f x = mod (x * x + c) n
    factor' x y 1 = factor' x' y' (gcd (x' - y') n) where
                        (x', y') = (f x, f $ f y)
    factor' _ _ d = if d == n then factor (c + 1) n else d

And factors does little more than recursively calling factor, while filtering out factors of two and making sure not to call factor on a prime (since the algorithm is not designed for that).

factors :: Integer -> StdGen -> [Integer]
factors n g = sort $ fs n where
    fs x | even x      = 2 : fs (div x 2)
         | isPrime x g = [x]
         | otherwise   = f : fs (div x f) where f = factor 1 x

And to test:

main :: IO ()
main = print . factors (2^98 - 1) =<< getStdGen

The end result is 9 lines of code, which was to be expected, given that it’s just a question of writing math in your language’s syntax.