fork download
copy 
  1. -- programming with prime numbers
  2.  
  3. import Control.Monad (forM_, when)
  4. import Control.Monad.ST
  5. import Data.Array.ST
  6. import Data.Array.Unboxed    
  7. import Data.List (sort)
  8.  
  9. sieve :: Int -> UArray Int Bool
  10. sieve n = runSTUArray $ do
  11.     let m = (n-1) `div` 2
  12.         r = floor . sqrt $ fromIntegral n
  13.     bits <- newArray (0, m-1) True
  14.     forM_ [0 .. r `div` 2 - 1] $ \i -> do
  15.         isPrime <- readArray bits i
  16.         when isPrime $ do
  17.             forM_ [2*i*i+6*i+3, 2*i*i+8*i+6 .. (m-1)] $ \j -> do
  18.                 writeArray bits j False
  19.     return bits
  20.  
  21. primes :: Int -> [Int]
  22. primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]
  23.  
  24. tdPrime :: Int -> Bool
  25. tdPrime n = prime (2:[3,5..])
  26.     where prime (d:ds)
  27.             | n < d * d      = True
  28.             | n `mod` d == 0 = False
  29.             | otherwise      = prime ds
  30.  
  31. tdFactors :: Int -> [Int]
  32. tdFactors n = facts n (2:[3,5..])
  33.     where facts n (f:fs)
  34.             | n < f * f      = [n]
  35.             | n `mod` f == 0 = f : facts (n `div` f) (f:fs)
  36.             | otherwise      = facts n fs
  37.  
  38. powmod :: Integer -> Integer -> Integer -> Integer
  39. powmod b e m =
  40.     let times p q = (p*q) `mod` m
  41.         pow b e x
  42.             | e == 0    = x
  43.             | even e    = pow (times b b) (e `div` 2) x
  44.             | otherwise = pow (times b b) (e `div` 2) (times b x)
  45.     in  pow b e 1
  46.  
  47. isSpsp :: Integer -> Integer -> Bool
  48. isSpsp n a =
  49.     let getDandS d s =
  50.             if even d then getDandS (d `div` 2) (s+1) else (d, s)
  51.         spsp (d, s) =
  52.             let t = powmod a d n
  53.             in  if t == 1 then True else doSpsp t s
  54.         doSpsp t s
  55.             | s == 0     = False
  56.             | t == (n-1) = True
  57.             | otherwise  = doSpsp ((t*t) `mod` n) (s-1)
  58.     in  spsp $ getDandS (n-1) 0
  59.  
  60. isPrime :: Integer -> Bool
  61. isPrime n =
  62.     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]
  63.     in  n `elem` ps || all (isSpsp n) ps 
  64.  
  65. rhoFactor :: Integer -> Integer -> Integer
  66. rhoFactor n c =
  67.     let f x = (x*x+c) `mod` n
  68.         fact t h
  69.             | d == 1    = fact t' h'
  70.             | d == n    = rhoFactor n (c+1)
  71.             | isPrime d = d
  72.             | otherwise = rhoFactor d (c+1)
  73.             where
  74.                 t' = f t
  75.                 h' = f (f h)
  76.                 d  = gcd (t' - h') n
  77.     in  fact 2 2
  78.  
  79. rhoFactors :: Integer -> [Integer]
  80. rhoFactors n =
  81.     let facts n
  82.             | n == 2    = [2]
  83.             | even n    = 2 : facts (n `div` 2)
  84.             | isPrime n = [n]
  85.             | otherwise = let f = rhoFactor n 1
  86.                           in  f : facts (n `div` f)
  87.     in  sort $ facts n
  88.  
  89. main = do
  90.     print $ primes 100
  91.     print $ length $ primes 1000000
  92.     print $ tdPrime 716151937
  93.     print $ tdFactors 8051
  94.     print $ powmod 437 13 1741
  95.     print $ isSpsp 2047 2
  96.     print $ isPrime 600851475143
  97.     print $ isPrime 2305843009213693951
  98.     print $ rhoFactors 600851475143
Success #stdin #stdout 0.04s 4644KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[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]
78498
False
[83,97]
819
True
False
True
[71,839,1471,6857]
https://ideone.com/uc0M8
language:
Haskell (ghc 8.4.4)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!