The start of a primal inquiry

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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
module Analytics.Theory.Number where

import Control.Monad
import Control.Monad.State.Lazy

import Data.Foldable hiding (toList)
import Data.Time.Clock

import Control.DList               -- http://lpaste.net/107607

import Data.Bag hiding (addn)      -- http://lpaste.net/107815 for the bonus
import Data.Deque                  -- http://lpaste.net/106780
import Data.Peano                  -- http://lpaste.net/107204
import Data.Stream                 -- http://lpaste.net/107665
import Data.Vector hiding ((<|))   -- http://lpaste.net/106865

-- just some math functions

isPrime :: Integer -> Bool
isPrime n = null (factors n)

primes :: Stream Integer
primes = filterS (\p -> Prelude.all (\q -> p `mod` q > 0)
                            (takeWhileS (\r -> r < p) primesish))
                 primesish  -- ugh! Linear on top of exponential cost!

primeNumberTheorem :: Stream (Integer, Int)
primeNumberTheorem = p' 10
   where p' n = (n, length (takeWhileS (< succ n) primes)) :< p' (n * 10)

-- so how do we represent a stream of primes? primesish is pretty-good
-- but a 'pretty-good' prime may not be a prime, and less so, given the
-- Prime Number Theorem.

type PrimeS a = State (Deque Integer, Stream Integer) a
type Prime = Integer

nprimes :: Int -> PrimeS [Prime]
nprimes n = get >>= \(primes, primeses) ->
   let ans = takeD n primes
       ps  = fromInt (n - length ans)
   in  addn ps ans

primesUpto :: Prime -> PrimeS [Prime]
primesUpto prime = get >>= \(primes, primese) ->
   if needMorePrimes primes prime
   then nextPrime >> primesUpto prime
   else return (takeWhileD (<= prime) primes)

-- e.g.: runState (primesUpto 1000) (Empty, primesish)

needMorePrimes :: Deque Prime -> Prime -> Bool
needMorePrimes Empty _ = True
needMorePrimes deck prime = Data.Deque.last deck < prime

addn :: Peano -> [Prime] -> PrimeS [Prime]
addn Z ans = return ans
addn (S n) partial = a' (S n) (dl partial)
   where a' Z dl = return (dlToList dl)
         a' (S n) dl = nextPrime >>= \prime -> a' n (dl <| prime)

nextPrime :: PrimeS Prime
nextPrime = get >>= \(primes, p :< ps) ->
   if Data.Foldable.foldr (\n b -> mod p n == 0 || b) False primes
   then put (primes, ps) >> nextPrime
   else put (primes |< p, ps) >> return p

-- so, for example:
-- execState (nprimes 50) (Empty, primesish)

{-- 

e.g.:

*Main> getCurrentTime >>= \t1 -> return (takeS 5000 primes) >>= \a -> 
       return (last a) >>= \b -> putStrLn (show b) >> getCurrentTime >>= \te ->
       return (b, diffUTCTime te t1)
7919
(7919,0.835521s)

*Analytics.Theory.Number> getCurrentTime >>= \t1 -> 
     return (execState (nprimes 1000) (Empty, primesish)) >>= \(a, b) -> 
     getCurrentTime >>= \te -> return ((first a, b), diffUTCTime te t1) 
((2,Stream{ 7921 :< 7925 :< 7927 :< 7931 :< 7933 :< ... }),0.409623s)

*Main> getCurrentTime >>= \t1 -> return (takeS 5000 primes) >>= \a -> 
       return (last a) >>= \b -> putStrLn (show b) >> getCurrentTime >>= \te ->
       return (b, diffUTCTime te t1)
48611
(48611,26.840493s)

*Analytics.Theory.Number> getCurrentTime >>= \t1 -> 
  return (execState (nprimes 5000) (Empty, primesish)) >>= \(a, b) -> 
  getCurrentTime >>= \te -> return ((Data.Deque.last a, b), diffUTCTime te t1) 
((48611,Stream{ 48613 :< 48617 :< 48619 :< 48623 :< 48625 :< ... }),10.640899s)

*Main> getCurrentTime >>= \t1 -> return (takeS 10000 primes) >>= \a -> 
       return (last a) >>= \b -> putStrLn (show b) >> getCurrentTime >>= \te ->
       return (b, diffUTCTime te t1)
104729
(104729,117.449069s)

*Analytics.Theory.Number> getCurrentTime >>= \t1 ->
  return (execState (nprimes 10000) (Empty, primesish)) >>= \(a, b) -> 
  getCurrentTime >>= \te -> return ((Data.Deque.last a, b), diffUTCTime te t1) 
((104729,Stream{ 104731 :< 104735 :< 104737 :< 104741 :< 104743 :< ... }),43.079645s)

 --}

factors' :: Integer -> [Integer]
factors' n = filter ((== 0) . mod n) [2 .. n `div` 2]

{--

So:

*Main Data.Time.Clock> getCurrentTime >>= \t1 ->
          putStrLn (show (factors 53463515)) >> 
          getCurrentTime >>= return . diffUTCTime t1
[5,7,35,1527529,7637645,10692703]
-11.885973s

and: factors 53463517 ~> [491,108887]

in about the same time. Changing the upper bound to (n `div` 2) cuts
the time in half to return in just under 6 seconds: -5.969551s

Of course, there's a much more efficient implementation of factors:

 --}

factors :: Integer -> [Integer]
factors n = factorRange 2 (n `div` 2) n

factorRange :: Integer -> Integer -> Integer -> [Integer]
factorRange start finish n = 
   nextFactor start finish n >>= \h ->
      h : factorRange (succ h) (n `div` h) n ++ [n `div` h]

nextFactor :: Integer -> Integer -> Integer -> [Integer]
nextFactor start finish n
   | start > finish = []    -- I do the maybeToList by hand here.
   | otherwise      = if n `mod` start == 0 then [start]
                      else nextFactor (succ start) finish n

{--

This gives the result in under 4 seconds:

*Main Data.Time.Clock> getCurrentTime >>= \t1 -> 
           putStrLn (show (factors 53463515)) >> 
           getCurrentTime >>= return . diffUTCTime t1
[5,7,35,1527529,35,1527529,7637645,10692703]
-3.751684s

 --}

pow :: Integer -> Integer -> Integer
pow x toDaY | toDaY <= 0 = 1
            | otherwise  = x * pow x (pred toDaY)

pthMersenne :: Integer -> Integer
pthMersenne prime = 2 `pow` prime - 1

{--

The first 10 mersenne primes:

##	p        digits
    (exponent)	Mp      Pp     year	discoverer
1	 2	 1	 1	----	----	 
2	 3	 1	 2	----	----	 
3	 5	 2	 3	----	----	 
4	 7	 3	 4	----	----	 
5	 13	 4	 8	 1456	 anonymous	 
6	 17	 6	 10	 1588	 Cataldi	 
7	 19	 6	 12	 1588	 Cataldi	 
8	31	 10	 19	 1772	 Euler	 
9	61	 19	 37	 1883	 Pervushin	 
10	 89	 27	 54	 1911	 Powers

from http://primes.utm.edu/mersenne/index.html

 --}

-- PRIME FACTORS -------------------------------------------------------

-- three different views of prime factors:

baggedPrimeFactors :: Integer -> Bag Prime
baggedPrimeFactors n = factors' n primes emptyBag
   where factors' 1 _ bag = bag
         factors' n primus bag
            = if isPrime n then (add n bag)
              else cont bag n primus (n `divMod` (headS primus))

         -- a continued function that rechecks primality of new numbers or
         -- just iterates (not checking primality) over non-factored primes
         cont bag n primus@(p :< ps) (d, m)
              | m == 0    = factors' d primus (add p bag)
              | otherwise = cont bag n ps (n `divMod` (headS ps))

groupedPrimeFactors :: Integer -> [(Prime, Int)]
groupedPrimeFactors = toAssocList . baggedPrimeFactors

primeFactors :: Integer -> [Prime]
primeFactors = toList . baggedPrimeFactors

-- COPRIMALITY ---------------------------------------------------------

data CoPrime = CoPrime Integer Integer
   deriving Show
data CoComposite =
    CoComposite { a :: Integer, b :: Integer,  exempliGratia :: Integer }
   deriving Show

coprime :: Integer -> Integer -> Either CoPrime CoComposite
coprime a b = co' (gcd a b)
   where co' 1 = Left $ CoPrime a b
         co' n = Right $ CoComposite a b n

isCoprime :: Integer -> Integer -> Bool
isCoprime a b = either (const True) (const False) (coprime a b)
24:42: Warning: Avoid lambda
Found:
\ r -> r < p
Why not:
(< p)
194:15: Warning: Redundant bracket
Found:
if isPrime n then (add n bag) else
cont bag n primus (n `divMod` (headS primus))
Why not:
if isPrime n then add n bag else
cont bag n primus (n `divMod` (headS primus))
195:39: Warning: Redundant bracket
Found:
n `divMod` (headS primus)
Why not:
n `divMod` headS primus
201:44: Warning: Redundant bracket
Found:
n `divMod` (headS ps)
Why not:
n `divMod` headS ps