Primarily Primes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
module Analytics.Theory.Number.Prime where

-- Okay, scrap all this. Complete rewrite

-- isPrimeC: checks if x is a prime by accumulating primes

-- How do we know if it's in there?

-- A set (partitioned?) contains the primes for fast lookup. Solved

-- How do we know what is max prime? separate index into the set, so we
-- don't have to do a full set scan.

import Control.Monad
import Control.Monad.State
import Data.Set (Set)
import qualified Data.Set as Set

import Analytics.Theory.Number            -- http://lpaste.net/107480
import Control.Logic.Frege (mand)         -- http://lpaste.net/111101
import Data.Stream                        -- http://lpaste.net/107665

-- This was originally motivated by http://lpaste.net/107700

type Primes = (Set Prime, Prime)

type PrimesS = StateT Primes

-- checks if p is prime, and, along the way, adds all primes up to the
-- checking value for p's primality.

isPrimeC :: Monad r => Integer -> PrimesS r Bool
isPrimeC p = get >>= \(primes, max) ->
   if p < max then return (Set.member p primes) else verify p

-- gets you the prime after x, x may possibly not be prime, nor in the
-- current state for that matter ... deal

primeAfterC :: Monad r => Integer -> PrimesS r Prime
primeAfterC x = get >>= \(set, max) ->
   if x < max then return (head (dropWhile (<= x) (Set.toList set)))
   else mapM_ verify (primeCandidates max x) >>
        primeAfterC (x + 2) >> -- first we ensure we DO have a prime after
        primeAfterC x          -- then we actually return that prime (after)

-- prime candidates gives us a list of primes(ish) numbers outside our
-- known primes

primeCandidates :: Prime -> Integer -> [Integer]
primeCandidates max p = takeWhilePlus1S (< p) (dropWhileS (<= max) primesish)

-- takes from a stream while pred is true, then takes just one more

takeWhilePlus1S :: (Integer -> Bool) -> Stream Integer -> [Integer]
takeWhilePlus1S pred (s :< tream) =
   if pred s then s : takeWhilePlus1S pred tream else [s]

verify :: Monad r => Integer -> PrimesS r Bool
verify p = get >>= \(primes, max) ->
   verifyC p max (Set.toList primes ++ primeCandidates max p)

verifyC :: Monad r => Integer -> Prime -> [Integer] -> PrimesS r Bool
verifyC p _ [] = get >>= \(primes, _) -> addC p >> return True
verifyC p max (h:t) = if p `mod` h == 0 then return False
                      else when (h > max) (checkC h) >> verifyC p max t

addC :: Monad r => Integer -> PrimesS r ()
addC p = get >>= \(set, _) -> put (Set.insert p set, p)

checkC :: Monad r => Integer -> PrimesS r ()
checkC p = get >>= \(set, _) -> 
   when (hardCheck p (floor (sqrt (fromIntegral p))) (Set.toList set)) (addC p)

-- See, we need to hard-check the values going in greater than (already-found)
-- max-prime, because primesish is as it is named: it's prime-..ish, but not
-- guaranteed prime (it includes 25, 55, 85, ...) so we need a hard-check here.

hardCheck :: Integer -> Integer -> [Prime] -> Bool
hardCheck p max [] = True
hardCheck p max (h:t) = 
   if max < h then True else p `mod` h > 0 && hardCheck p max t

-- hardcheck is linear on the square root of p
-- that is to say: it's NOT linear on p but more like logarhythmic

{--
*Main> runState (isPrimeC 1000) (Set.empty , 0) ~> (False,(fromList [],0))

fast-fail: 1000 not prime based off of 2.

runState (isPrimeC 10001) (Set.empty , 0) ~>
(False,(fromList [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71],71))

We get that 10001 is not prime AND a set of primes up to its composite, 71.

Sweet!

Fun fact: there are 78498 primes < 1,000,000

How do I know this? "Ask me about my grandprimeren"

And this newly rewritten Analytics module
--}
81:4: Error: Redundant if
Found:
if max < h then True else p `mod` h > 0 && hardCheck p max t
Why not:
(max < h) || (p `mod` h > 0 && hardCheck p max t)