codepad

Language:

-- programming with prime numbers
 
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed    
import Data.List (sort)
 
sieve :: Integer -> UArray Integer Bool
sieve n = runSTUArray $ do
    let m = (n-1) `div` 2
        r = ceiling . sqrt $ fromInteger n
    bits <- newArray (0, m-1) True
    forM_ [0 .. r `div` 2 - 1] $ \i -> do
        isPrime <- readArray bits i
        when isPrime $ do
            forM_ [2*i*i+6*i+3, 2*i*i+8*i+6 .. (m-1)] $ \j -> do
                writeArray bits j False
    return bits
 
primes :: Integer -> [Integer]
primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]
 
tdPrime :: Integer -> Bool
tdPrime n = prime (2:[3,5..])
    where prime (d:ds)
            | n < d * d      = True
            | n `mod` d == 0 = False
            | otherwise      = prime ds
 
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 
 
tdFactors :: Integer -> [Integer]
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
            
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
 
main = do
    print $ primes 100
    print $ length $ primes 1000000
    print $ tdPrime 600851475143
    print $ powmod 437 13 1741
    print $ isSpsp 2047 2
    print $ isPrime 600851475143
    print $ isPrime 2305843009213693951
    print $ tdFactors 600851475143
    print $ rhoFactors 600851475143

Private [?]

Run code