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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
module Analytics.Theory.Number.Continued.Continuations where

{-- We need to improve the straight-up prime number computation from
Analytics.Theory.Number and we do this in a couple of ways here. Firstly
we introduce a continuation to inform us where we are in our prime number
verification search (by emitting ticks for each n-steps in the computation)
(Cheesy, I know, but very reassuring to the watcher of these processes),
the second way is to introduce a stopping point to say: so far, and no
further. --}

import Control.Monad
import Control.Monad.State.Lazy
import Data.Time.Clock
import Analytics.Theory.Number            -- http://lpaste.net/107480
import Data.Stream -- for primesish       -- http://lpaste.net/107665

-- First improvement: monitoring progress for http://lpaste.net/107700

isPrimeC' :: Monad r => 
            (Integer -> Integer -> r Bool) -> Integer -> r Bool
isPrimeC' c p = liftM2 (&&) (c p 2) (primalRangeC' 3 (_sqrtInt p) p c)
   where _sqrtInt = floor . sqrt . fromIntegral
-- (p `div` 3) p c) -- increase efficiency here exponentially

-- we need to record the previously-found primes so that we test for
-- primality using only primes, speeding up this computation ... or
-- so we hope.

primalRangeC' :: Monad r => Integer -> Integer -> Integer -> 
                (Integer -> Integer -> r Bool) -> r Bool
primalRangeC' start finish p c 
   | start > finish = return True
   | otherwise      = liftM2 (&&) (c p start)
                         (primalRangeC' (start + 2) (p `div` (start + 2)) p c)

dotStepper :: Int -> Integer -> Integer -> StateT Int IO Bool
dotStepper ticks p base = get >>= \step -> 
   when (step `mod` ticks == 0) (lift (putStr ".")) >> stepper p base

stepper :: Monad m => Integer -> Integer -> StateT Int m Bool
stepper p base = get >>= \step -> put (succ step) >> return (p `mod` base > 0)

go :: Integer -> (Integer -> StateT s IO Bool) -> s -> IO (Bool, s)
go mersexp fn seed = 
   let pth = pthMersenne mersexp
   in  putStrLn ("Solving for (2 ** " ++ show mersexp ++ " - 1) mersenne "
                 ++ "prime: " ++ show pth) >>
       getCurrentTime >>= \startTime -> runStateT (fn pth) seed >>=
       \ans -> getCurrentTime >>= \endTime ->
       putStrLn ("Duration was: " ++ show (diffUTCTime endTime startTime)) >>
       return ans

gostep :: Integer -> IO (Bool, Int)
gostep mexp = go mexp (isPrimeC' (dotStepper 1000)) 1
{--

*Analytics.Theory.Number> gostep 19
Solving for (2 ** 19 - 1) mersenne prime: 524287
Duration was: 0.004946s
(True,363)

*Analytics.Theory.Number> gostep 31
Solving for (2 ** 31 - 1) mersenne prime: 2147483647
Duration was: 0.212032s
(True,23171)

*Analytics.Theory.Number> gostep 29
Solving for (2 ** 29 - 1) mersenne prime: 536870911
Duration was: 0.112159s
(False,11586)

*Analytics.Theory.Number> gostep 33
Solving for (2 ** 33 - 1) mersenne prime: 8589934591
Duration was: 0.478234s
(False,46342)

*Analytics.Theory.Number> gostep 37
Solving for (2 ** 37 - 1) mersenne prime: 137438953471
.Duration was: 2.514618s
(False,185365)

*Analytics.Theory.Number> gostep 39
Solving for (2 ** 39 - 1) mersenne prime: 549755813887
...Duration was: 6.996841s
(False,370729)

but ...

*Analytics.Theory.Number> gostep 61
Solving for (2 ** 61 - 1) mersenne prime: 2305843009213693951
......................................................<CTRL-C>

and ...

*Analytics.Theory.Number> gostep 127
Loading package transformers-0.3.0.0 ... linking ... done.
Loading package mtl-2.1.2 ... linking ... done.
.....................................................<CTRL-C>

... yeah, and ... 

... yeah. That.

 --}

-- Second improvement: a stopper! for http://lpaste.net/107701

{--

The continuation is now of the form: contC escapeC
and is in the domain StateT ((step, testedUpTo), stop) IO Bool

When step >= stop, we escape the computation.

 --}

type Step = ((Int, Stream Integer), Int)  -- (step, testedUpto), maxsteps
type StepState = StateT Step IO Bool
type IndexedTarget = Integer -- nextToTest > prime ... actually, nah! just int
type TestM = IndexedTarget -> StepState
type EscapeC = Bool -> StepState -- used to be: Monad r => r Bool

isPrimeC :: Monad r =>
            (Integer -> IndexedTarget -> TestM -> EscapeC -> r Bool) -> 
            TestM -> EscapeC -> Integer -> r Bool
isPrimeC c contC escapeC p = c (ceiling (sqrt (fromInteger p))) p contC escapeC

escape :: EscapeC
escape maxhit = get >>= \((ticks, upto :< _), maxticks) ->
   lift (putStrLn ("Tested up to " ++ show upto) >> 
   when (ticks >= maxticks) 
        (putStrLn "(still more possible factors to test)")) >> 
   return maxhit

checkNext :: TestM
checkNext p = get >>= \((ticks, start :< ps), maxTicks) ->
   put ((ticks + 1, ps), maxTicks) >> return (p `mod` start > 0)

stopper :: Int -> Integer -> IndexedTarget -> TestM -> EscapeC -> StepState
stopper ticks maxm p testM escapeC = 
   get >>= \((step, curr :< _), maxTicks) ->
   when (step `mod` ticks == 0) (lift (putStr ".")) >>
   -- ("Testing " ++ show curr))) >>
   testM p >>= \res ->
   if curr < maxm && res && step < maxTicks
   then stopper ticks maxm p testM escapeC
   else escapeC res

gostop :: Integer -> IO (Bool, Step)
gostop mexp =
   go mexp (isPrimeC (stopper 10000) checkNext escape)
           ((1, primesish), 1000000)

{--

Cool! Here's some numbers for ya!

When the max-test limit was set to 1,000:

*Analytics.Theory.Number.Continued.Continuations> gostop 31
Solving for (2 ** 31 - 1) mersenne prime: 2147483647
..........Tested up to 2999
(still more possible factors to test)
Duration was: 0.020192s
(True,((1001,Stream{ 2999 :< 3001 :< 3005 :< 3007 :< 3011 :< ... }),1000))

When the max-test limit was reset to 1,000,000:

*Analytics.Theory.Number.Continued.Continuations> gostop 31
Solving for (2 ** 31 - 1) mersenne prime: 2147483647
.Tested up to 46345
Duration was: 0.146618s
(True,((15450,Stream{ 46345 :< 46349 :< 46351 :< 46355 :< 46357 :< ... }),1000000))

*Analytics.Theory.Number.Continued.Continuations> gostop 39
Solving for (2 ** 39 - 1) mersenne prime: 549755813887
Tested up to 11
Duration was: 0.000608s
(False,((5,Stream{ 11 :< 13 :< 17 :< 19 :< 23 :< ... }),1000000))

*Analytics.Theory.Number.Continued.Continuations> gostop 61
Solving for (2 ** 61 - 1) mersenne prime: 2305843009213693951
....................................................................................................Tested up to 2999999
(still more possible factors to test)
Duration was: 13.305372s
(True,((1000001,Stream{ 2999999 :< 3000001 :< 3000005 :< 3000007 :< 3000011 :< ... }),1000000))

WOOT!

 --}