-- programming with prime numbers
import Control
.Monad (forM
_, when
) import Data.Array.ST
import Data.Array.Unboxed
import Data.List (sort)
sieve n = runSTUArray $ do
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
primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]
tdPrime n = prime (2:[3,5..])
where prime (d:ds)
| n < d * d = True
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
)
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
in pow b e 1
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
in spsp $ getDandS (n-1) 0
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 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
where
t' = f t
h' = f (f h)
in fact 2 2
rhoFactors n =
let facts n
| n == 2 = [2]
| isPrime n = [n]
in sort $ facts n
main = do
print $ tdPrime
716151937 print $ powmod
437 13 1741 print $ isPrime
600851475143 print $ isPrime
2305843009213693951 print $ rhoFactors
600851475143
LS0gcHJvZ3JhbW1pbmcgd2l0aCBwcmltZSBudW1iZXJzCiAKaW1wb3J0IENvbnRyb2wuTW9uYWQgKGZvck1fLCB3aGVuKQppbXBvcnQgQ29udHJvbC5Nb25hZC5TVAppbXBvcnQgRGF0YS5BcnJheS5TVAppbXBvcnQgRGF0YS5BcnJheS5VbmJveGVkICAgIAppbXBvcnQgRGF0YS5MaXN0IChzb3J0KQogCnNpZXZlIDo6IEludCAtPiBVQXJyYXkgSW50IEJvb2wKc2lldmUgbiA9IHJ1blNUVUFycmF5ICQgZG8KICAgIGxldCBtID0gKG4tMSkgYGRpdmAgMgogICAgICAgIHIgPSBmbG9vciAuIHNxcnQgJCBmcm9tSW50ZWdyYWwgbgogICAgYml0cyA8LSBuZXdBcnJheSAoMCwgbS0xKSBUcnVlCiAgICBmb3JNXyBbMCAuLiByIGBkaXZgIDIgLSAxXSAkIFxpIC0+IGRvCiAgICAgICAgaXNQcmltZSA8LSByZWFkQXJyYXkgYml0cyBpCiAgICAgICAgd2hlbiBpc1ByaW1lICQgZG8KICAgICAgICAgICAgZm9yTV8gWzIqaSppKzYqaSszLCAyKmkqaSs4KmkrNiAuLiAobS0xKV0gJCBcaiAtPiBkbwogICAgICAgICAgICAgICAgd3JpdGVBcnJheSBiaXRzIGogRmFsc2UKICAgIHJldHVybiBiaXRzCiAKcHJpbWVzIDo6IEludCAtPiBbSW50XQpwcmltZXMgbiA9IDIgOiBbMippKzMgfCAoaSwgVHJ1ZSkgPC0gYXNzb2NzICQgc2lldmUgbl0KIAp0ZFByaW1lIDo6IEludCAtPiBCb29sCnRkUHJpbWUgbiA9IHByaW1lICgyOlszLDUuLl0pCiAgICB3aGVyZSBwcmltZSAoZDpkcykKICAgICAgICAgICAgfCBuIDwgZCAqIGQgICAgICA9IFRydWUKICAgICAgICAgICAgfCBuIGBtb2RgIGQgPT0gMCA9IEZhbHNlCiAgICAgICAgICAgIHwgb3RoZXJ3aXNlICAgICAgPSBwcmltZSBkcwogCnRkRmFjdG9ycyA6OiBJbnQgLT4gW0ludF0KdGRGYWN0b3JzIG4gPSBmYWN0cyBuICgyOlszLDUuLl0pCiAgICB3aGVyZSBmYWN0cyBuIChmOmZzKQogICAgICAgICAgICB8IG4gPCBmICogZiAgICAgID0gW25dCiAgICAgICAgICAgIHwgbiBgbW9kYCBmID09IDAgPSBmIDogZmFjdHMgKG4gYGRpdmAgZikgKGY6ZnMpCiAgICAgICAgICAgIHwgb3RoZXJ3aXNlICAgICAgPSBmYWN0cyBuIGZzCiAKcG93bW9kIDo6IEludGVnZXIgLT4gSW50ZWdlciAtPiBJbnRlZ2VyIC0+IEludGVnZXIKcG93bW9kIGIgZSBtID0KICAgIGxldCB0aW1lcyBwIHEgPSAocCpxKSBgbW9kYCBtCiAgICAgICAgcG93IGIgZSB4CiAgICAgICAgICAgIHwgZSA9PSAwICAgID0geAogICAgICAgICAgICB8IGV2ZW4gZSAgICA9IHBvdyAodGltZXMgYiBiKSAoZSBgZGl2YCAyKSB4CiAgICAgICAgICAgIHwgb3RoZXJ3aXNlID0gcG93ICh0aW1lcyBiIGIpIChlIGBkaXZgIDIpICh0aW1lcyBiIHgpCiAgICBpbiAgcG93IGIgZSAxCiAKaXNTcHNwIDo6IEludGVnZXIgLT4gSW50ZWdlciAtPiBCb29sCmlzU3BzcCBuIGEgPQogICAgbGV0IGdldERhbmRTIGQgcyA9CiAgICAgICAgICAgIGlmIGV2ZW4gZCB0aGVuIGdldERhbmRTIChkIGBkaXZgIDIpIChzKzEpIGVsc2UgKGQsIHMpCiAgICAgICAgc3BzcCAoZCwgcykgPQogICAgICAgICAgICBsZXQgdCA9IHBvd21vZCBhIGQgbgogICAgICAgICAgICBpbiAgaWYgdCA9PSAxIHRoZW4gVHJ1ZSBlbHNlIGRvU3BzcCB0IHMKICAgICAgICBkb1Nwc3AgdCBzCiAgICAgICAgICAgIHwgcyA9PSAwICAgICA9IEZhbHNlCiAgICAgICAgICAgIHwgdCA9PSAobi0xKSA9IFRydWUKICAgICAgICAgICAgfCBvdGhlcndpc2UgID0gZG9TcHNwICgodCp0KSBgbW9kYCBuKSAocy0xKQogICAgaW4gIHNwc3AgJCBnZXREYW5kUyAobi0xKSAwCiAKaXNQcmltZSA6OiBJbnRlZ2VyIC0+IEJvb2wKaXNQcmltZSBuID0KICAgIGxldCBwcyA9IFsyLDMsNSw3LDExLDEzLDE3LDE5LDIzLDI5LDMxLDM3LDQxLDQzLDQ3LDUzLDU5LDYxLDY3LDcxLDczLDc5LDgzLDg5LDk3XQogICAgaW4gIG4gYGVsZW1gIHBzIHx8IGFsbCAoaXNTcHNwIG4pIHBzIAogICAgICAgICAgICAKcmhvRmFjdG9yIDo6IEludGVnZXIgLT4gSW50ZWdlciAtPiBJbnRlZ2VyCnJob0ZhY3RvciBuIGMgPQogICAgbGV0IGYgeCA9ICh4KngrYykgYG1vZGAgbgogICAgICAgIGZhY3QgdCBoCiAgICAgICAgICAgIHwgZCA9PSAxICAgID0gZmFjdCB0JyBoJwogICAgICAgICAgICB8IGQgPT0gbiAgICA9IHJob0ZhY3RvciBuIChjKzEpCiAgICAgICAgICAgIHwgaXNQcmltZSBkID0gZAogICAgICAgICAgICB8IG90aGVyd2lzZSA9IHJob0ZhY3RvciBkIChjKzEpCiAgICAgICAgICAgIHdoZXJlCiAgICAgICAgICAgICAgICB0JyA9IGYgdAogICAgICAgICAgICAgICAgaCcgPSBmIChmIGgpCiAgICAgICAgICAgICAgICBkICA9IGdjZCAodCcgLSBoJykgbgogICAgaW4gIGZhY3QgMiAyCiAKcmhvRmFjdG9ycyA6OiBJbnRlZ2VyIC0+IFtJbnRlZ2VyXQpyaG9GYWN0b3JzIG4gPQogICAgbGV0IGZhY3RzIG4KICAgICAgICAgICAgfCBuID09IDIgICAgPSBbMl0KICAgICAgICAgICAgfCBldmVuIG4gICAgPSAyIDogZmFjdHMgKG4gYGRpdmAgMikKICAgICAgICAgICAgfCBpc1ByaW1lIG4gPSBbbl0KICAgICAgICAgICAgfCBvdGhlcndpc2UgPSBsZXQgZiA9IHJob0ZhY3RvciBuIDEKICAgICAgICAgICAgICAgICAgICAgICAgICBpbiAgZiA6IGZhY3RzIChuIGBkaXZgIGYpCiAgICBpbiAgc29ydCAkIGZhY3RzIG4KIAptYWluID0gZG8KICAgIHByaW50ICQgcHJpbWVzIDEwMAogICAgcHJpbnQgJCBsZW5ndGggJCBwcmltZXMgMTAwMDAwMAogICAgcHJpbnQgJCB0ZFByaW1lIDcxNjE1MTkzNwogICAgcHJpbnQgJCB0ZEZhY3RvcnMgODA1MQogICAgcHJpbnQgJCBwb3dtb2QgNDM3IDEzIDE3NDEKICAgIHByaW50ICQgaXNTcHNwIDIwNDcgMgogICAgcHJpbnQgJCBpc1ByaW1lIDYwMDg1MTQ3NTE0MwogICAgcHJpbnQgJCBpc1ByaW1lIDIzMDU4NDMwMDkyMTM2OTM5NTEKICAgIHByaW50ICQgcmhvRmFjdG9ycyA2MDA4NTE0NzUxNDM=
[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]