Compare commits

...

4 Commits

Author SHA1 Message Date
Nathan McCarty 31515dc82d
Setup test and bench harness 2023-06-23 10:36:11 -04:00
Nathan McCarty 0257c9fde7
Add collie dependency 2023-06-23 04:51:17 -04:00
Nathan McCarty 8d48abed0e
Add naive prime generation code 2023-06-23 03:08:05 -04:00
Nathan McCarty a3d54c59df
Add integer square root 2023-06-23 02:24:58 -04:00
7 changed files with 328 additions and 13 deletions

View File

@ -12,10 +12,14 @@ license = "Parity Public License 7.0.0"
-- langversion
-- packages to add to search path
-- depends =
depends = contrib
, collie
-- modules to install
modules = PrimeSieve
, PrimeSieve.Util
, PrimeSieve.Trivial
, PrimeSieve.Simple
-- main file (i.e. file to load at REPL)
main = PrimeSieve

View File

@ -61,11 +61,11 @@
"rust": "rust"
},
"locked": {
"lastModified": 1686465942,
"narHash": "sha256-+6ake+PQ9zVKrKaz3g8haA/t0woq1XLmbxYtoQwLgd8=",
"lastModified": 1687504884,
"narHash": "sha256-M7eCdHpEV6WkVL2vPahM0X86eTS+zw72xeCC1LN2/I4=",
"ref": "trunk",
"rev": "b306a4aff0f183268bd11d96ab3043638097f6f0",
"revCount": 72,
"rev": "b4d806c37bdd0dc91980e4e4eb0dee5383964cf5",
"revCount": 73,
"type": "git",
"url": "https://git.sr.ht/~thatonelutenist/idr2nix"
},
@ -99,11 +99,11 @@
},
"nixpkgs": {
"locked": {
"lastModified": 1686398752,
"narHash": "sha256-nGWNQVhSw4VSL+S0D0cbrNR9vs9Bq7rlYR+1K5f5j6w=",
"lastModified": 1687488839,
"narHash": "sha256-7JDjuyHwUvGJJge9jxfRJkuYyL5G5yipspc4J3HwjGA=",
"owner": "nixos",
"repo": "nixpkgs",
"rev": "a30520bf8eabf8a5c37889d661e67a2dbcaa59e6",
"rev": "f9e94676ce6c7531c44d38da61d2669ebec0f603",
"type": "github"
},
"original": {

View File

@ -1,9 +1,21 @@
{
"sources": {},
"sorted": [],
"sources": {
"collie": {
"source": {
"url": "https://github.com/ohad/collie",
"rev": "46bff04a8d9a1598fec9b19f515541df16dc64ef",
"sha256": "sha256-0GziqzEUDs1tFd5yiu5NECV0nP1aMPe2QxS9oqhmW6I="
},
"ipkg": "collie.ipkg",
"name": "collie"
}
},
"sorted": [
"collie"
],
"idris2": {
"url": "https://github.com/idris-lang/Idris2",
"rev": "5dcf62499df5cb861d153372ef3b4386dba25c98",
"sha256": "sha256-P9fVZNtgu08cvJmanL88W4F2l0PWKbSk8h8SS5JEN/M="
"rev": "31c17ebec2e64c095076e8425637176db7658492",
"sha256": "sha256-DYwyjTdJjPDd/H5keb0tfmzR4AGwg/EaZMMQnLfksfY="
}
}

View File

@ -1,4 +1,230 @@
module PrimeSieve
import System.Clock
import Data.String
import Collie
import PrimeSieve.Trivial as Trivial
import PrimeSieve.Simple as Simple
primeSieve : Command "primesieve"
primeSieve = MkCommand
{ description = """
"""
, subcommands =
[ "--help" ::= basic "Print this help text." none
, "bench-trivial" ::= basic "Benchmark the trivial generator" none
, "bench-simple" ::= basic "Benchmark the simple sieve" none
, "test-simple" ::= basic "test the simple sieve against trivial generation" none
, "test-trivial" ::= basic "test the trival generation against trivial generation" none
]
, modifiers = []
, arguments = none }
{nm : String} -> {cmd : Command nm} -> Show (ParseTreeT f g cmd) where
show (Here x) = "\{nm} <<args>>"
show (There pos parsedSub) = "\{nm} \{show parsedSub}"
interface PrimeGen where
constructor MkPrimeGen
primesUntil : (Integral t, Ord t, Num t, Range t) => (limit : t) -> List t
primes : (Integral t, Ord t, Num t, Range t) => Maybe (Stream t)
%noinline
call : (a -> b) -> a -> IO b
call f x = fromPrim $ \w => MkIORes (f x) w
displayTime : Clock t -> String
displayTime x =
let (secs, ns) = (seconds x, nanoseconds x)
mins = secs `div` 60
secs = secs `mod` 60
nanos = ns `mod` 1_000
ns = ns `div` 1_000
micros = ns `mod` 1_000
millis = ns `div` 1_000
in """
\{padLeft 3 ' ' (show mins)}m \
\{padLeft 3 ' ' (show secs)}s \
\{padLeft 3 ' ' (show millis)}ms \
\{padLeft 3 ' ' (show micros)}μs \
\{padLeft 3 ' '(show nanos)}ns
"""
timeToDouble : Clock t -> Double
timeToDouble x =
let (secs, ns) = (seconds x, nanoseconds x)
ns_double : Double = fromInteger ns
secs_double = (fromInteger secs)
in secs_double + (ns_double / 1000000000.0)
record Timed a where
constructor MkTimed
result : a
desc : String
runCount : Nat
runsNonZero : NonZero runCount
runs : Vect runCount (Clock Duration)
%name Timed timed, timed2, timed3
(.min) : Timed a -> Clock Duration
(.min) (MkTimed result desc runCount runsNonZero runs) = minVect' runs
where
minVect : (Ord b) => (xs : Vect n b) -> b -> b
minVect [] x = x
minVect (y :: xs) x =
if y < x
then minVect xs y
else minVect xs x
minVect' : (Ord b) => {n : Nat} -> (xs : Vect n b) -> {auto prf : NonZero n} -> b
minVect' {prf} [] impossible
minVect' (x :: xs) = minVect xs x
(.max) : Timed a -> Clock Duration
(.max) (MkTimed result desc runCount runsNonZero runs) = maxVect' runs
where
maxVect : (Ord b) => (xs : Vect n b) -> b -> b
maxVect [] x = x
maxVect (y :: xs) x =
if y > x
then maxVect xs y
else maxVect xs x
maxVect' : (Ord b) => {n : Nat} -> (xs : Vect n b) -> {auto prf : NonZero n} -> b
maxVect' {prf} [] impossible
maxVect' (x :: xs) = maxVect xs x
(.avg) : Timed a -> Clock Duration
(.avg) (MkTimed result desc runCount runsNonZero runs) = avgVect runCount runs
where
totalVect : (xs : Vect n (Clock Duration)) -> (acc : (Integer, Integer)) -> (Integer, Integer)
totalVect [] acc = acc
totalVect (x :: xs) (sec_total, ns_total) =
let (sec, ns) = (seconds x, nanoseconds x)
in totalVect xs (sec_total + sec, ns_total + ns)
avgVect : (count : Nat) -> (xs : Vect count (Clock Duration)) -> Clock Duration
avgVect 0 xs = makeDuration 0 0
avgVect count@(S count') xs =
let (sec, ns) = totalVect xs (0,0)
t = ((sec * 1_000_000_000) + ns) `div` (natToInteger count)
in makeDuration 0 t
Show (Timed a) where
show t = """
\{t.desc} \
min: \{displayTime t.min} \
avg: \{displayTime t.avg} \
max: \{displayTime t.max}
"""
mergeTimed : (x : Timed a) -> (y : Timed a) -> Timed a
mergeTimed (MkTimed result desc 0 runsNonZero runs) (MkTimed x str k y xs) impossible
mergeTimed (MkTimed result desc (S j) runsNonZero runs) (MkTimed x str k y xs) =
MkTimed x str ((S j) + k) SIsNonZero (runs ++ xs)
timeIO : (prev : Maybe (Timed a))
-> (action : IO a)
-> (case prev of
Nothing => (desc : String) -> IO (Timed a)
(Just _) => IO (Timed a))
timeIO Nothing action =
(\desc =>
do start <- clockTime Monotonic
result <- action
end <- clockTime Monotonic
let time = timeDifference end start
pure (MkTimed result desc 1 SIsNonZero [time]))
timeIO (Just timed) action = do
single <- timeIO Nothing action timed.desc
pure (mergeTimed timed single)
runTimes : (times : Nat) -> {auto prf : NonZero times} -> (desc : String) -> (action : IO a) -> IO (Timed a)
runTimes 0 desc action impossible
runTimes (S times') desc action = do
timed <- timeIO Nothing action desc
runTimes' times' timed action
where
runTimes' : Nat -> (acc : Timed a) -> IO a -> IO (Timed a)
runTimes' 0 acc x = pure acc
runTimes' (S k) acc x = do
timed <- timeIO (Just acc) x
runTimes' k timed x
||| Benchmark a PrimeGen
benchmark : PrimeGen -> IO ()
benchmark x = do
putStrLn "Benching primesUntil method\n"
benchUntil 100
benchUntil 1_000
benchUntil 10_000
benchUntil 100_000
where
benchUntil : (limit : Bits64) -> IO ()
benchUntil limit =
do putStrLn " Benchmarking primes up to \{show limit} (100 times):"
time <- runTimes
100
"Primes <= \{padLeft 6 ' ' (show limit)}"
((call PrimeSieve.primesUntil) limit)
putStrLn " \{show time}"
let primesPerSec = 1.0 / (timeToDouble time.avg)
putStrLn
" Avg \{show primesPerSec} primes/s (\{show (length time.result)} primes)\n"
||| Test a PrimeGen against the default
test : PrimeGen -> IO ()
test item = do
putStrLn "Testing primesUntil method\n"
tryUntil 10
tryUntil 100
tryUntil 1_000
tryUntil 10_000
tryUntil 100_000
putStrLn ""
let stream : Maybe (Stream Bits64) = PrimeSieve.primes
case stream of
Nothing => exitSuccess
(Just x) => do
tryStream 10 x
tryStream 100 x
tryStream 1_000 x
tryStream 10_000 x
where
fail : Show t => List t -> List t -> IO ()
fail expected got = do
putStrLn "fail\n"
putStrLn "Expected: \{show expected}"
putStrLn "Got: \{show got}"
exitFailure
tryCompare : (Show t, Eq t) => List t -> List t -> IO ()
tryCompare xs ys =
if xs == ys
then putStrLn "pass"
else fail xs ys
tryUntil : Bits64 -> IO ()
tryUntil k = do
putStr " Testing up to \{show k}: "
let trivial : List Bits64 = Trivial.primesUntil k
let sample : List Bits64 = PrimeSieve.primesUntil k
tryCompare trivial sample
tryStream : Nat -> Stream Bits64 -> IO ()
tryStream k stream = do
putStr " Testing up to \{show k}: "
let trivial = take k Trivial.primes
let sample = take k stream
tryCompare trivial sample
main : IO ()
main = putStrLn "Hello!"
main = do
Right cmdParse <- primeSieve.parseArgs
| Left err => putStrLn "Error: \{err}"
case fst (lookup cmdParse) of
"bench-trivial" => benchmark (MkPrimeGen Trivial.primesUntil (Just Trivial.primes))
"bench-simple" => benchmark (MkPrimeGen Simple.primesUntil Nothing)
"test-trivial" => test (MkPrimeGen Trivial.primesUntil (Just Trivial.primes))
"test-simple" => test (MkPrimeGen Simple.primesUntil Nothing)
"--help" => putStrLn "Usage: \n\{primeSieve.usage}"
_ => putStrLn "Parse as \{show cmdParse}"

View File

@ -0,0 +1,9 @@
module PrimeSieve.Simple
||| Use a basic, unoptimized sieve to generate the primes up to a specific number
export
primesUntil : (Integral t, Ord t, Num t, Range t) => (limit : t) -> List t
||| A stream of primes, Generated by a lazily extending sieve
export
primes : (Integral t, Ord t, Num t, Range t) => Stream t

View File

@ -0,0 +1,46 @@
module PrimeSieve.Trivial
import Data.Stream
import PrimeSieve.Util
%default total
||| Test if a number is prime via trial division
export
isPrime : (Integral t, Ord t, Num t, Range t) => t -> Bool
isPrime x =
let trial_divisors = [2..((isqrt x) + 1)]
in if x <= 2
then x == 2
else isPrime' trial_divisors x
where
isPrime' : List t -> t -> Bool
isPrime' [] x = True
isPrime' (y :: xs) x =
if x `mod` y == 0
then False
else isPrime' xs x
||| A stream of primes, generated by testing via trial division
export
primes : (Integral t, Ord t, Num t, Range t) => Stream t
primes =
let naturals : Stream t = iterate (+1) 1
in unfoldr next_prime naturals
where
next_prime : Stream t -> (t, Stream t)
next_prime orig@(y :: ys) =
if isPrime y
then (y, ys)
-- We assert_smaller here as there are an infinite number of primes, so we can never
-- run out of primes, and taking one off the head of the stream will always get us
-- closer to the next prime
else next_prime (assert_smaller orig ys)
||| All the primes up until the given limit
export
primesUntil : (Integral t, Ord t, Num t, Range t) => (limit : t) -> List t
-- We assert_total here as the list of primes is infinite and strictly increasing, so this
-- Will always terminate in finite time
primesUntil limit = assert_total $ takeBefore (> limit) primes

18
src/PrimeSieve/Util.idr Normal file
View File

@ -0,0 +1,18 @@
module PrimeSieve.Util
%default total
||| Integer square root
export
isqrt : (Integral t, Ord t) => t -> t
isqrt n = sqrt' n
where
sqrt' : t -> t
sqrt' a =
let b = (a + (n `div` a)) `div` 2
in if a > b then sqrt' (assert_smaller a b) else a
||| Nat square root
export
sqrtNat : Nat -> Nat
sqrtNat = integerToNat . isqrt . natToInteger