Programming Praxis

import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.List (sort)

ancientSieve :: Int -> UArray Int Bool
ancientSieve n = runSTUArray $ do
    bits <- newArray (2, n) True
    forM_ [2 .. n] $ \p -> do
        isPrime <- readArray bits p
        when isPrime $ do
            forM_ [2*p, 3*p .. n] $ \i -> do
                writeArray bits i False
    return bits

ancientPrimes :: Int -> [Int]
ancientPrimes n = [p | (p, True) <-
                  assocs $ ancientSieve n]

sieve :: Int -> UArray Int Bool
sieve n = runSTUArray $ do
    let m = (n-1) `div` 2
        r = floor . sqrt $ fromIntegral n
    bits <- newArray (0, m-1) True
    forM_ [0 .. r `div` 2 - 1] $ \i -> do
        isPrime <- readArray bits i
        when isPrime $ do
            let a = 2*i*i + 6*i + 3
                b = 2*i*i + 8*i + 6
            forM_ [a, b .. (m-1)] $ \j -> do
                writeArray bits j False
    return bits

primes :: Int -> [Int]
primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]

tdPrime :: Int -> Bool
tdPrime n = prime (2:[3,5..])
    where prime (d:ds)
            | n < d * d = True
            | n `mod` d == 0 = False
            | otherwise = prime ds

tdFactors :: Int -> [Int]
tdFactors n = facts n (2:[3,5..])
    where facts n (f:fs)
            | n < f * f = [n]
            | n `mod` f == 0 =
                  f : facts (n `div` f) (f:fs)
            | otherwise = facts n fs

powmod :: Integer -> Integer -> Integer -> Integer
powmod b e m =
    let times p q = (p*q) `mod` m
        pow b e x
            | e == 0 = x
            | even e = pow (times b b)
                              (e `div` 2) x
            | otherwise = pow (times b b)
                              (e `div` 2) (times b x)
    in pow b e 1

isSpsp :: Integer -> Integer -> Bool
isSpsp n a =
    let getDandS d s =
            if even d then getDandS (d `div` 2) (s+1)
                      else (d, s)
        spsp (d, s) =
            let t = powmod a d n
            in if t == 1 then True else doSpsp t s
        doSpsp t s
            | s == 0 = False
            | t == (n-1) = True
            | otherwise = doSpsp ((t*t) `mod` n) (s-1)
    in spsp $ getDandS (n-1) 0

isPrime :: Integer -> Bool
isPrime n =
    let ps = [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]
    in n `elem` ps || all (isSpsp n) ps

rhoFactor :: Integer -> Integer -> Integer
rhoFactor n c =
    let f x = (x*x+c) `mod` n
        fact t h
            | d == 1 = fact t' h'
            | d == n = rhoFactor n (c+1)
            | isPrime d = d
            | otherwise = rhoFactor d (c+1)
            where
                t' = f t
                h' = f (f h)
                d = gcd (t' - h') n
    in fact 2 2

rhoFactors :: Integer -> [Integer]
rhoFactors n =
    let facts n
            | n == 2 = [2]
            | even n = 2 : facts (n `div` 2)
            | isPrime n = [n]
            | otherwise = let f = rhoFactor n 1
                          in f : facts (n `div` f)
    in sort $ facts n