Some prime factors for ya!

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
224
225
226
import Control.Monad.State
import Control.Arrow

import Data.Time.Clock

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

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

import Analytics.Theory.Number   -- http://lpaste.net/107480

{--

From the P-99 Prolog problem set, exercises 35, 36: prime factors of n.

A solution-set to the Haskell exercise stated at http://lpaste.net/107878

first we solve the bonus-bonus problem as it is the basis of the other
solutions, simplifying them significantly: --}

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))

{--

-- the above takes an extraordinarily long time if we recompute primes
-- from scrate each time. [ed: previous statement made before the isPrime
-- sanity check was in place] Let's leverage state here:

baggedPrimeFactorsS' :: Integer -> PrimeS (Bag Prime)
baggedPrimeFactorsS' n = primeAfter 1 >>= b' n emptyBag
--    primesUpto (n `div` 2) >> get >>= \(primes, _) -> -- BLECH! TOO EAGER!
   where b' 1 bag _ = return bag
         b' n bag prime =
                 let (d, m) = n `divMod` prime
                 in  if m == 0 then b' d (add prime bag) prime
                               else primeAfter prime >>= b' n bag

This was how far I tolerated before I implemented scanPrimes

*Main> process' "big.rnds" 
(554621,Bag.fromList [31,17891])
(229558,Bag.fromList [2,7,19,863])
(280884,Bag.fromList [2,2,3,89,263])
(513462,Bag.fromList [2,3,85577])
(551418,Bag.fromList [2,3,7,19,691])
(96467,Bag.fromList [7,13781])
(319263,Bag.fromList [3,7,23,661])
(480355,Bag.fromList [5,23,4177])
(882602,Bag.fromList [2,7,23,2741])
(278410,Bag.fromList [2,5,11,2531])
(964438,Bag.fromList [2,127,3797])
(460282,Bag.fromList [2,373,617])
(469102,Bag.fromList [2,79,2969])
(806036,Bag.fromList [2,2,7,11,2617])
(416922,Bag.fromList [2,3,11,6317])
(630363,Bag.fromList [3,19,11059])
(788380,Bag.fromList ^CInterrupted.

 --}

-- adding scanPrimes (then isPrime) for a sanity check:

baggedPrimeFactorsS :: Integer -> PrimeS (Bag Prime)
baggedPrimeFactorsS n = primeAfter 1 >>= b' n emptyBag
   where b' 1 bag _ = return bag
         b' n bag prime = get >>= \(primes, _) ->
            if isPrime n -- third enhancement instead of: scanPrimes n primes 
            then return (add n bag)
            else cont bag (n `divMod` prime) n prime
         cont bag (d, m) n prime =
            if m == 0 then b' d (add prime bag) prime
                      else primeAfter prime >>= \p ->
                           cont bag (n `divMod` p) n p

{-- the below was trial and error code. Did not work; moved on.

-- a quick scan to see if x is a prime already computed

scanPrimes :: Integer -> Deque Integer -> Bool
scanPrimes x Empty = False
scanPrimes x (FirstAndLast y) = x == y
scanPrimes x (Deck a d c) = x == a || x == c || (x < c && scanPrimes x d)

-- the above presupposes an ordered deque
-- nope, not get primes because this does not advance state. :/
getPrimes :: PrimeS (Prime, Deque Prime)
getPrimes = get >>= \(primes, _) -> g' primes
   where g' Empty = nextPrime >> getPrimes
         g' primes = return (nab primes)
 --}

runit :: Integer -> IO (Bag Prime, (Vector Prime, Stream Prime))
runit n = getCurrentTime >>= \start ->
          return (runState (baggedPrimeFactorsS n) (newVector, primesish)) >>=
          \ans -> getCurrentTime >>= \stop ->
          putStrLn (show (diffUTCTime stop start)) >> return ans

process' :: FilePath -> IO (NominalDiffTime, Int)
-- PrimeS [(Integer, Bag Prime)]) 
-- StateT (Deque Prime, Stream Prime) (,) IO
process' file = readFile file >>= return . map read . lines >>= \nums ->
   getCurrentTime >>= \start ->
   return (runState (mapBagM nums emptyDL) (newVector, primesish)) >>=
   \(dlist, (primes, _)) -> mapM_ (putStrLn . show) (dlToList dlist) >>
   getCurrentTime >>= \end -> return (diffUTCTime end start, lengthV primes)

{--

Final result with vectors and isPrime check:

*Main> process' "big.rnds" ~>

(554621,Bag.fromList [31,17891])
(229558,Bag.fromList [2,7,19,863])
(280884,Bag.fromList [2,2,3,89,263])
(513462,Bag.fromList [2,3,85577])
(551418,Bag.fromList [2,3,7,19,691])
(96467,Bag.fromList [7,13781])
(319263,Bag.fromList [3,7,23,661])
(480355,Bag.fromList [5,23,4177])
(882602,Bag.fromList [2,7,23,2741])
(278410,Bag.fromList [2,5,11,2531])
(964438,Bag.fromList [2,127,3797])
(460282,Bag.fromList [2,373,617])
(469102,Bag.fromList [2,79,2969])
(806036,Bag.fromList [2,2,7,11,2617])
(416922,Bag.fromList [2,3,11,6317])
(630363,Bag.fromList [3,19,11059])
(788380,Bag.fromList [2,2,5,39419])
(461331,Bag.fromList [3,3,13,3943])
(419818,Bag.fromList [2,7,157,191])
(997638,Bag.fromList [2,3,166273])
(588634,Bag.fromList [2,294317])
(51549,Bag.fromList [3,17183])
(848759,Bag.fromList [17,49927])
(686765,Bag.fromList [5,137353])
(739023,Bag.fromList [3,181,1361])
(582112,Bag.fromList [2,2,2,2,2,18191])
(645056,Bag.fromList [2,2,2,2,2,2,10079])
(861956,Bag.fromList [2,2,229,941])
(692063,Bag.fromList [692063])
(751578,Bag.fromList [2,3,229,547])
(831612,Bag.fromList [2,2,3,37,1873])
(648665,Bag.fromList [5,129733])
(760310,Bag.fromList [2,5,76031])
(564409,Bag.fromList [564409])
(377543,Bag.fromList [377543])
(576664,Bag.fromList [2,2,2,11,6553])
(208616,Bag.fromList [2,2,2,89,293])
(785950,Bag.fromList [2,5,5,11,1429])
(813692,Bag.fromList [2,2,11,18493])
(988871,Bag.fromList [13,29,43,61])
(445098,Bag.fromList [2,3,31,2393])
(298833,Bag.fromList [3,99611])
(811843,Bag.fromList [389,2087])
(363767,Bag.fromList [363767])
(365218,Bag.fromList [2,7,19,1373])
(914264,Bag.fromList [2,2,2,13,59,149])
(798554,Bag.fromList [2,399277])
(976592,Bag.fromList [2,2,2,2,67,911])
(55115,Bag.fromList [5,73,151])
(240255,Bag.fromList [3,3,5,19,281])
(3.846796s,77)

Three seconds to give the prime factors of a bunch of numbers. Not bad!

So, at the end of the day ... (and into tomorrow :/), we have
not only our prime factors but also, at each application, a growing list
(vector, actually) of known primes against which we can do ... whatever
with. YAY!

 --}

type IndexedPrimeFactors = DList (Integer, Bag Prime)

mapBagM :: [Integer] -> IndexedPrimeFactors -> PrimeS IndexedPrimeFactors
mapBagM [] ans = return ans
mapBagM (h : t) accum = baggedPrimeFactorsS h >>= \bag ->
                        mapBagM t (accum <| (h, bag))

{--

So, against the rnd set in big.rnds:

readFile "big.rnds" >>= mapM read . lines >>= \nums -> 

With the above, P36 is simply the toList of the above bagged solution.

P36: Now, group these prime factors. One way do to this is as was demonstrated
in the P-99 problem set: --}

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

{--

So, for example:

*Main> groupedPrimeFactors 315 ~> [(3, 2), (5, 1), (7, 1)]
 becomes an unwind of the above solution. (So, we try that).

After that, P35 is simply the unwind of P36. Or, in otherwords, P35 is
toList baggedPrimeFactors

P35: Find the prime factors of a given positive integer greater than 1:

 --}

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

-- That was easy!
-- ... eventually ;)
27: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))
28:39: Warning: Redundant bracket
Found:
n `divMod` (headS primus)
Why not:
n `divMod` headS primus
34:44: Warning: Redundant bracket
Found:
n `divMod` (headS ps)
Why not:
n `divMod` headS ps
107:11: Error: Monad law, left identity
Found:
return (runState (baggedPrimeFactorsS n) (newVector, primesish))
>>=
\ ans ->
getCurrentTime >>=
\ stop -> putStrLn (show (diffUTCTime stop start)) >> return ans
Why not:
((\ ans ->
getCurrentTime >>=
\ stop -> putStrLn (show (diffUTCTime stop start)) >> return ans)
(runState (baggedPrimeFactorsS n) (newVector, primesish)))
109:11: Error: Use print
Found:
putStrLn (show (diffUTCTime stop start))
Why not:
print (diffUTCTime stop start)
114:17: Warning: Use liftM
Found:
readFile file >>= return . map read . lines
Why not:
Control.Monad.liftM (map read . lines) (readFile file)
116:4: Error: Monad law, left identity
Found:
return (runState (mapBagM nums emptyDL) (newVector, primesish)) >>=
\ (dlist, (primes, _)) ->
mapM_ (putStrLn . show) (dlToList dlist) >> getCurrentTime >>=
\ end -> return (diffUTCTime end start, lengthV primes)
Why not:
((\ (dlist, (primes, _)) ->
mapM_ (putStrLn . show) (dlToList dlist) >> getCurrentTime >>=
\ end -> return (diffUTCTime end start, lengthV primes))
(runState (mapBagM nums emptyDL) (newVector, primesish)))
117:36: Error: Use print
Found:
putStrLn . show
Why not:
print