Personal tools

Euler problems/181 to 190

From HaskellWiki

< Euler problems(Difference between revisions)
Jump to: navigation, search
(another solution to problem 185)
(Problem 186: efficient solution using STUArray)
 
(7 intermediate revisions by 3 users not shown)
Line 44: Line 44:
 
Triangles containing the origin.
 
Triangles containing the origin.
   
Solution:
+
{{sect-stub}}
<haskell>
 
problem_184 = undefined
 
</haskell>
 
   
 
== [http://projecteuler.net/index.php?section=problems&id=185 Problem 185] ==
 
== [http://projecteuler.net/index.php?section=problems&id=185 Problem 185] ==
Line 133: Line 133:
   
 
bestGuess :: Mind -> [[Char]] -> [Int]
 
bestGuess :: Mind -> [[Char]] -> [Int]
bestGuess mind guesses =
+
bestGuess mind guesses = bestGuesses
let scores = map (scoreMind mind) guesses
+
where scores = map (scoreMind mind) guesses
bestScore = minimum scores
+
bestScore = minimum scores
bestGuesses = findIndices (==bestScore) scores
+
bestGuesses = elemIndices bestScore scores
in bestGuesses
 
   
 
iterateGuesses :: Mind -> [Char] -> [Char]
 
iterateGuesses :: Mind -> [Char] -> [Char]
iterateGuesses mind value =
+
iterateGuesses mind value = nextguess value mins
let allguesses = map (guesses value) [0..(length value)-1]
+
where allguesses = map (guesses value) [0..(length value)-1]
mins = map (bestGuess mind) allguesses
+
mins = map (bestGuess mind) allguesses
in nextguess value mins
 
   
 
nextguess :: [Char] -> [[Int]] -> [Char]
 
nextguess :: [Char] -> [[Int]] -> [Char]
nextguess prev mins =
+
nextguess = zipWith choose
let choose x = if length (snd x) == 1 then intToDigit ((snd x)!!0) else fst x
+
where choose x [y] = intToDigit y
both = zip prev mins
+
choose x _ = x
in foldr (\x y -> (choose x) : y) "" both
 
   
   
 
iterateMind :: Mind -> [Char] -> [([Char], Int)]
 
iterateMind :: Mind -> [Char] -> [([Char], Int)]
iterateMind mind n =
+
iterateMind mind n = zip b c
let a = drop 2 $ inits $ iterate (iterateGuesses mind) n
+
where a = drop 2 $ inits $ iterate (iterateGuesses mind) n
b = last $ takeWhile (\x -> (last x) `notElem` (init x)) a
+
b = last $ takeWhile (\x -> (last x) `notElem` (init x)) a
c = map (scoreMind mind) b
+
c = map (scoreMind mind) b
in zip b c
 
 
 
 
 
randomStart :: (Num a, Enum a) => a -> IO [Char]
+
randomStart :: Int -> IO [Char]
randomStart n = mapM (\_ -> getStdRandom (randomR ('0','9'))) [1..n]
+
randomStart n = replicateM n (getStdRandom (randomR ('0','9')))
   
 
main :: IO ()
 
main :: IO ()
Line 164: Line 164:
 
x <- randomStart (length (fst (head ex)))
 
x <- randomStart (length (fst (head ex)))
 
let y = iterateMind ex x
 
let y = iterateMind ex x
let done = (snd (last y) == 0)
+
done = snd (last y) == 0
when done (putStrLn $ (fst.last) y)
+
if done then (putStrLn $ (fst.last) y)
unless done main
+
else main
   
 
</haskell>
 
</haskell>
Line 206: Line 206:
   
 
solve e@(Environment _ _ 0) = [[]]
 
solve e@(Environment _ _ 0) = [[]]
solve e@(Environment [] r _) = sequence $ map (S.toList . (domain S.\\)) r
+
solve e@(Environment [] r _) = mapM (S.toList . (domain S.\\)) r
 
solve e' = do
 
solve e' = do
 
is <- m
 
is <- m
Line 242: Line 242:
 
problem_185 = head . solve $ initial
 
problem_185 = head . solve $ initial
 
</haskell>
 
</haskell>
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=186 Problem 186] ==
  +
Connectedness of a network.
  +
  +
This implements the known tree structure for imperative [http://en.wikipedia.org/wiki/Disjoint-set_data_structure union/find].
  +
On my machine it runs in about 3 seconds.
  +
  +
<haskell>
  +
module ProjectEuler186 where
  +
  +
import qualified Control.Monad.ST.Strict as StrictST
  +
import qualified Control.Monad.ST.Lazy as LazyST
  +
import Data.Array.ST (STUArray, newArray, newListArray, readArray, writeArray, )
  +
import Control.Monad (liftM2, when, )
  +
  +
  +
data Network s =
  +
Network
  +
(STUArray s Telephone Telephone)
  +
(STUArray s Telephone Size)
  +
  +
  +
type Telephone = Int
  +
type Size = Int
  +
  +
  +
initialize :: StrictST.ST s (Network s)
  +
initialize =
  +
let bnds = (0,999999)
  +
in liftM2 Network
  +
(newListArray bnds [0..])
  +
(newArray bnds 1)
  +
  +
union :: Network s -> Telephone -> Telephone -> StrictST.ST s ()
  +
union net@(Network repr size) m n = do
  +
mroot <- find net m
  +
nroot <- find net n
  +
when (mroot/=nroot) $ do
  +
msize <- readArray size mroot
  +
nsize <- readArray size nroot
  +
let totalSize = msize+nsize
  +
if msize < nsize
  +
then do
  +
writeArray repr mroot nroot
  +
writeArray size nroot totalSize
  +
else do
  +
writeArray repr nroot mroot
  +
writeArray size mroot totalSize
  +
  +
find :: Network s -> Telephone -> StrictST.ST s Telephone
  +
find (Network repr _size) =
  +
let go n = do
  +
m <- readArray repr n
  +
if n==m
  +
then return m
  +
else go m
  +
in go
  +
  +
  +
callersStart :: [Telephone]
  +
callersStart =
  +
map (\k -> fromInteger $ mod (100003 - 200003*k + 300007*k^(3::Int)) 1000000) [1..]
  +
  +
callers :: [Telephone]
  +
callers =
  +
take 55 callersStart ++
  +
zipWith (\s24 s55 -> mod (s24+s55) 1000000) (drop 31 callers) callers
  +
  +
pairs :: [a] -> [(a,a)]
  +
pairs (x0:x1:xs) = (x0,x1) : pairs xs
  +
pairs _ = error "list of callers must be infinite"
  +
  +
answer :: Int
  +
answer =
  +
(1+) . length . takeWhile (<990000) $
  +
LazyST.runST (do
  +
net@(Network _repr size) <- LazyST.strictToLazyST initialize
  +
mapM
  +
(\(from,to) ->
  +
LazyST.strictToLazyST $
  +
union net from to >>
  +
find net 524287 >>=
  +
readArray size) $
  +
filter (uncurry (/=)) $
  +
pairs $ callers)
  +
  +
main :: IO ()
  +
main = print answer
  +
</haskell>
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=187 Problem 187] ==
  +
Semiprimes
  +
  +
Solution:
  +
  +
The solution to this problem isn't terribly difficult, once you
  +
know that the numbers the problem is referring to are called
  +
semiprimes. In fact there is an excellent write-up at:
  +
http://mathworld.wolfram.com/Semiprime.html which provides an
  +
explicit formula for the number of semiprimes less than a given
  +
number.
  +
  +
The problem with this formula is the use of pi(n) where pi is
  +
the number of primes less than n. For Mathematica users this is
  +
no problem since the function is built in, but for us poor
  +
Haskeller's we have to build it ourselves. This is what took
  +
the most time for me, to come up with a pi(n) function that can
  +
compute pi(50000000) in less than the lifetime of the universe.
  +
Thus I embarked on a long and circuitous journey that eventually
  +
led me to the PI module below, which does little more than read
  +
a 26MB file of prime numbers into memory and computes a map from
  +
which we can calculate pi(n). I am sure there must be a better
  +
way of doing this, and I look forward to this entry being
  +
amended (or replaced) with a more reasonable solution.
  +
HenryLaxen Mar 24, 2008
  +
  +
  +
<haskell>
  +
  +
module PI (prime,primes,pI) where
  +
import Control.Monad.ST
  +
import System.IO.Unsafe
  +
import qualified Data.Map as M
  +
import Data.ByteString.Char8 (readFile,lines,readInt)
  +
import Data.Maybe
  +
import Data.Array.ST
  +
import Data.Array.IArray
  +
  +
r = runSTUArray
  +
(do a <- newArray (0,3001134) 0 :: ST s (STUArray s Int Int)
  +
writeArray a 0 1
  +
zipWithM_ (writeArray a) [0..] primes
  +
return a)
  +
  +
{-# NOINLINE s #-}
  +
s = M.fromDistinctAscList $ map (\(x,y) -> (y,x)) (assocs r)
  +
  +
prime n = r!n
  +
  +
pI :: Int -> Int
  +
pI n =
  +
case M.splitLookup n s of
  +
(_,Just x,_) -> (x+1)
  +
(_,_,b) -> snd (head (M.assocs b))
  +
  +
{-# NOINLINE primes #-}
  +
primes :: [Int]
  +
primes = unsafePerformIO $ do
  +
l <- Data.ByteString.Char8.readFile "primes_to_5M.txt"
  +
let p = Data.ByteString.Char8.lines l
  +
r = map (fst . fromJust . readInt) p::[Int]
  +
return r
  +
  +
  +
import PI
  +
import Data.List
  +
  +
s n k = ( n `div` (prime (k-1))) - k + 1
  +
  +
semiPrimeCount n =
  +
let last = pI $ floor $ sqrt (fromIntegral n)
  +
s k = pI ( n `div` (prime (k-1))) - k + 1
  +
pI2 = foldl' (\x y -> s y + x) 0 [1..last]
  +
in pI2
  +
  +
main = print (semiPrimeCount 100000000)
  +
  +
</haskell>
  +
  +
The file primes_to_5M.txt is:
  +
2
  +
3
  +
5 ..
  +
49999991
  +
50000017
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=188 Problem 188] ==
  +
hyperexponentiation
  +
  +
The idea here is actually very simple. Euler's theorem tells us that a^phi(n) mod n = 1, so the exponent never need grow above 40000000 in this case, which is phi(10^8). Henry Laxen April 6, 2008
  +
  +
<haskell>
  +
fastPower x 1 modulus = x `mod` modulus
  +
fastPower x n modulus
  +
| even n = fastPower ((x*x) `mod` modulus) (n `div` 2) modulus
  +
| otherwise = (x * (fastPower x (n-1) modulus)) `mod` modulus
  +
  +
modulus = 10^8
  +
phi = 40000000 -- eulerTotient of 10^8
  +
  +
a :: Int -> Int -> Int
  +
a n 1 = n
  +
-- a n k = n^(a n (k-1))
  +
a n k = fastPower n (a n (k-1) `mod` phi) modulus
  +
  +
problem_188 = a 1777 1855
  +
  +
  +
</haskell>
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=190 Problem 190] ==
  +
Maximising a weighted product
  +
  +
For those of us who can reach all the back into our sophomore
  +
calculus class memories, we will find the dreaded Lagrange
  +
Multipliers, which say that if you want to max(min)imize f
  +
subject to a constraint g, then the gradients of f and g have to
  +
be parallel, or grad f = L * grad g.
  +
  +
Happy happy, for in our case f = x * y^2 * .. z^n and even better
  +
g = x + y + .. z - n, so grad g = (1,1,..1) and
  +
grad f = (f/x, 2f/y, .. nf/z) That leads us to:
  +
  +
x = f/L, y = 2f/L .. z = nf/L and
  +
sum (1..n) x == n so x = 2/(n+1)
  +
  +
Thus the Haskell program is nothing more than:
  +
  +
<haskell>
  +
f l = product $ zipWith (^) l [1..]
  +
  +
fmin n = floor $ f [ i*2/(n+1) | i<-[1..n]]
  +
  +
problem_190 = sum $ map fmin [2..15]
  +
</haskell>
  +
  +
Aren't you glad you didn't skip that math class now? HenryLaxen Apr 18,2008

Latest revision as of 13:39, 25 November 2011

Contents

[edit] 1 Problem 181

Investigating in how many ways objects of two different colours can be grouped.

[edit] 2 Problem 182

RSA encryption.

Solution:

fun a1 b1 = sum [ e | e <- [2..a*b-1],
                      gcd e (a*b) == 1,
                      gcd (e-1) a == 2,
                      gcd (e-1) b == 2 ]
  where a = a1-1
        b = b1-1
 
problem_182 = fun 1009 3643

[edit] 3 Problem 183

Maximum product of parts.

Solution:

-- Does the decimal expansion of p/q terminate?
terminating p q = 1 == reduce [2,5] (q `div` gcd p q)
        where reduce   []   n = n
              reduce (x:xs) n | n `mod` x == 0 = reduce (x:xs) (n `div` x)
                              | otherwise      = reduce xs n
 
-- The expression (round $ fromIntegral n / e) computes the integer k
-- for which (n/k)^k is at a maximum. Also note that, given a rational number
-- r and a natural number k, the decimal expansion of r^k terminates if
-- and only if the decimal expansion of r does.
answer = sum [if terminating n (round $ fromIntegral n / e) then -n else n
                | n <- [5 .. 10^4]]
        where e = exp 1
 
main = print answer

[edit] 4 Problem 184

Triangles containing the origin.

[edit] 5 Problem 185

Number Mind

Solution:

This approach does NOT solve the problem in under a minute, unless of course you are extremely lucky. The best time I've seen so far has been about 76 seconds. Before I came up with this code, I tried to search for the solution by generating a list of all possible solutions based on the information given in the guesses. This was feasible with the 5 digit problem, but completely intractable with the 16 digit problem. The approach here, simple yet effective, is to make a random guess, and then vary each digit in the guess from [0..9], generating a score of how well the guess matched the given numbermind clues. You then improve the guess by selecting those digits that had a unique low score. It turns out this approach converges rather quickly, but can often be stuck in cycles, so we test for this and try a differenct random first guess if a cycle is detected. Once you run the program, you might have time for a cup of coffee, or maybe even a dinner. HenryLaxen 2008-03-12

import Data.List
import Control.Monad
import Data.Char
import System.Random
 
type Mind = [([Char],Int)]
 
values :: [Char]
values = "0123456789"
 
 
score :: [Char] -> [Char] -> Int
score guess answer = foldr (\(a,b) y -> if a == b then y+1 else y) 0 
      (zip guess answer)
 
scores :: Mind -> [Char] -> [Int]
scores m g = map (\x -> abs ((snd x) - score (fst x) g)) m
 
scoreMind :: Mind -> [Char] -> Int
scoreMind m g = sum $ scores m g
 
 
ex1 :: Mind
ex1 = 
  [("90342",2),
   ("39458",2),
   ("51545",2),
   ("34109",1),
   ("12531",1),
   ("70794",0)]
 
ex2 :: Mind   
ex2 = 
  [
   ("5616185650518293",2),
   ("3847439647293047",1),
   ("5855462940810587",3),
   ("9742855507068353",3),
   ("4296849643607543",3),
   ("3174248439465858",1),
   ("4513559094146117",2),
   ("7890971548908067",3),
   ("8157356344118483",1),
   ("2615250744386899",2),
   ("8690095851526254",3),
   ("6375711915077050",1),
   ("6913859173121360",1),
   ("6442889055042768",2),
   ("2321386104303845",0),
   ("2326509471271448",2),
   ("5251583379644322",2),
   ("1748270476758276",3),
   ("4895722652190306",1),
   ("3041631117224635",3),
   ("1841236454324589",3),
   ("2659862637316867",2)]
 
 
guesses :: [Char] -> Int -> [[Char]]
guesses str pos = [ left ++ n:(tail right) | n<-values]
  where (left,right) = splitAt pos str
 
bestGuess :: Mind -> [[Char]] -> [Int]
bestGuess mind guesses = bestGuesses
  where scores = map (scoreMind mind) guesses
        bestScore = minimum scores
        bestGuesses = elemIndices bestScore scores
 
iterateGuesses :: Mind -> [Char] -> [Char]
iterateGuesses mind value = nextguess value mins
   where allguesses = map (guesses value) [0..(length value)-1]
         mins = map (bestGuess mind) allguesses
 
nextguess :: [Char] -> [[Int]] -> [Char]
nextguess = zipWith choose
  where choose x [y] = intToDigit y
        choose x _   = x
 
 
iterateMind :: Mind -> [Char] -> [([Char], Int)]
iterateMind mind n = zip b c
  where a = drop 2 $ inits $ iterate (iterateGuesses mind) n
        b = last $ takeWhile (\x -> (last x) `notElem` (init x)) a
        c = map (scoreMind mind) b
 
 
randomStart :: Int -> IO [Char]
randomStart n = replicateM n (getStdRandom (randomR ('0','9')))
 
main :: IO ()
main = do
  let ex = ex1
  x <- randomStart (length (fst (head ex)))
  let y = iterateMind ex x
      done = snd (last y) == 0
  if done then (putStrLn $ (fst.last) y)
          else main

Here's another solution, and this one squeaks by in just under a minute on my machine. The basic idea is to just do a back-tracking search, but with some semi-smart pruning and guess ordering. The code is in pretty much the order I wrote it, so most prefixes of this code should also compile. This also means you should be able to figure out what each function does one at a time.

import Control.Monad
import Data.List
import qualified Data.Set as S
 
ensure p x = guard (p x) >> return x
selectDistinct 0 _  = [[]]
selectDistinct n [] = []
selectDistinct n (x:xs) = map (x:) (selectDistinct (n - 1) xs) ++ selectDistinct n xs
 
data Environment a = Environment { guesses      :: [(Int, [a])]
                                 , restrictions :: [S.Set a]
                                 , assignmentsLeft :: Int
                                 } deriving (Eq, Ord, Show)
 
reorder e = e { guesses = sort . guesses $ e }
domain  = S.fromList "0123456789"
initial = Environment gs (replicate a S.empty) a where
    a   = length . snd . head $ gs
    gs  = [(2, "5616185650518293"), (1, "3847439647293047"), (3, "5855462940810587"), (3, "9742855507068353"), (3, "4296849643607543"), (1, "3174248439465858"), (2, "4513559094146117"), (3, "7890971548908067"), (1, "8157356344118483"), (2, "2615250744386899"), (3, "8690095851526254"), (1, "6375711915077050"), (1, "6913859173121360"), (2, "6442889055042768"), (0, "2321386104303845"), (2, "2326509471271448"), (2, "5251583379644322"), (3, "1748270476758276"), (1, "4895722652190306"), (3, "3041631117224635"), (3, "1841236454324589"), (2, "2659862637316867")]
 
acceptableCounts e = small >= 0 && big <= assignmentsLeft e where
    ns    = (0:) . map fst . guesses $ e
    small = minimum ns
    big   = maximum ns
 
positions s = map fst . filter (not . snd) . zip [0..] . zipWith S.member s
acceptableRestriction  r (n, s) = length (positions s r) >= n
acceptableRestrictions e = all (acceptableRestriction (restrictions e)) (guesses e)
 
firstGuess = head . guesses
sanityCheck e = acceptableRestrictions e && acceptableCounts e
 
solve e@(Environment _  _ 0) = [[]]
solve e@(Environment [] r _) = mapM (S.toList . (domain S.\\)) r
solve e' = do
    is   <- m
    newE <- f is
    rest <- solve newE
    return $ interleaveAscIndices is (l is) rest
        where
    f = ensure sanityCheck . update e
    m = selectDistinct n (positions g (restrictions e)) 
    e = reorder e'
    l = fst . flip splitAscIndices g
    (n, g) = firstGuess e
 
splitAscIndices = indices 0 where
    indices _ [] xs = ([], xs)
    indices n (i:is) (x:xs)
        | i == n = let (b, e) = indices (n + 1) is     xs in (x:b, e)
        | True   = let (b, e) = indices (n + 1) (i:is) xs in (b, x:e)
 
interleaveAscIndices = indices 0 where
    indices _ [] [] ys = ys
    indices n (i:is) (x:xs) ys
        | i == n = x : indices (n + 1) is xs ys
        | True   = head ys : indices (n + 1) (i:is) (x:xs) (tail ys)
 
update (Environment ((_, a):gs) r l) is = Environment newGs restriction (l - length is) where
    (assignment, newRestriction)    = splitAscIndices is a
    (_, oldRestriction)             = splitAscIndices is r
    restriction                     = zipWith S.insert newRestriction oldRestriction
    newGs                           = map updateEntry gs
    updateEntry (n', a')            = (newN, newA) where
        (dropped, newA) = splitAscIndices is a'
        newN            = n' - length (filter id $ zipWith (==) assignment dropped)
 
problem_185 = head . solve $ initial

[edit] 6 Problem 186

Connectedness of a network.

This implements the known tree structure for imperative union/find. On my machine it runs in about 3 seconds.

module ProjectEuler186 where
 
import qualified Control.Monad.ST.Strict as StrictST
import qualified Control.Monad.ST.Lazy as LazyST
import Data.Array.ST (STUArray, newArray, newListArray, readArray, writeArray, )
import Control.Monad (liftM2, when, )
 
 
data Network s =
   Network
      (STUArray s Telephone Telephone)
      (STUArray s Telephone Size)
 
 
type Telephone = Int
type Size = Int
 
 
initialize :: StrictST.ST s (Network s)
initialize =
   let bnds = (0,999999)
   in  liftM2 Network
          (newListArray bnds [0..])
          (newArray bnds 1)
 
union :: Network s -> Telephone -> Telephone -> StrictST.ST s ()
union net@(Network repr size) m n = do
   mroot <- find net m
   nroot <- find net n
   when (mroot/=nroot) $ do
      msize <- readArray size mroot
      nsize <- readArray size nroot
      let totalSize = msize+nsize
      if msize < nsize
        then do
           writeArray repr mroot nroot
           writeArray size nroot totalSize
        else do
           writeArray repr nroot mroot
           writeArray size mroot totalSize
 
find :: Network s -> Telephone -> StrictST.ST s Telephone
find (Network repr _size) =
   let go n = do
          m <- readArray repr n
          if n==m
            then return m
            else go m
   in  go
 
 
callersStart :: [Telephone]
callersStart =
   map (\k -> fromInteger $ mod (100003 - 200003*k + 300007*k^(3::Int)) 1000000) [1..]
 
callers :: [Telephone]
callers =
   take 55 callersStart ++
   zipWith (\s24 s55 -> mod (s24+s55) 1000000) (drop 31 callers) callers
 
pairs :: [a] -> [(a,a)]
pairs (x0:x1:xs) = (x0,x1) : pairs xs
pairs _ = error "list of callers must be infinite"
 
answer :: Int
answer =
   (1+) . length . takeWhile (<990000) $
   LazyST.runST (do
      net@(Network _repr size) <- LazyST.strictToLazyST initialize
      mapM
         (\(from,to) ->
            LazyST.strictToLazyST $
               union net from to >>
               find net 524287 >>=
               readArray size) $
         filter (uncurry (/=)) $
         pairs $ callers)
 
main :: IO ()
main = print answer

[edit] 7 Problem 187

Semiprimes

Solution:

The solution to this problem isn't terribly difficult, once you know that the numbers the problem is referring to are called semiprimes. In fact there is an excellent write-up at: http://mathworld.wolfram.com/Semiprime.html which provides an explicit formula for the number of semiprimes less than a given number.

The problem with this formula is the use of pi(n) where pi is the number of primes less than n. For Mathematica users this is no problem since the function is built in, but for us poor Haskeller's we have to build it ourselves. This is what took the most time for me, to come up with a pi(n) function that can compute pi(50000000) in less than the lifetime of the universe. Thus I embarked on a long and circuitous journey that eventually led me to the PI module below, which does little more than read a 26MB file of prime numbers into memory and computes a map from which we can calculate pi(n). I am sure there must be a better way of doing this, and I look forward to this entry being amended (or replaced) with a more reasonable solution. HenryLaxen Mar 24, 2008


module PI (prime,primes,pI) where 
import Control.Monad.ST
import System.IO.Unsafe
import qualified Data.Map as M
import Data.ByteString.Char8 (readFile,lines,readInt) 
import Data.Maybe
import Data.Array.ST
import Data.Array.IArray
 
r = runSTUArray 
            (do a <- newArray (0,3001134) 0 :: ST s (STUArray s Int Int)
                writeArray a 0 1
                zipWithM_ (writeArray a) [0..] primes
                return a) 
 
{-# NOINLINE s #-}
s = M.fromDistinctAscList $ map (\(x,y) -> (y,x)) (assocs r)
 
prime n = r!n
 
pI :: Int -> Int
pI n = 
    case M.splitLookup n s of
       (_,Just x,_) -> (x+1)
       (_,_,b) -> snd (head (M.assocs b))
 
{-# NOINLINE primes #-}
primes :: [Int]
primes = unsafePerformIO $ do
  l <- Data.ByteString.Char8.readFile "primes_to_5M.txt"
  let p = Data.ByteString.Char8.lines l
      r = map (fst . fromJust . readInt) p::[Int]
  return r  
 
 
import PI
import Data.List
 
s n k =       ( n `div` (prime (k-1))) - k + 1
 
semiPrimeCount n =
  let last = pI $ floor $ sqrt (fromIntegral n)
      s k = pI ( n `div` (prime (k-1))) - k + 1
      pI2 = foldl' (\x y -> s y + x) 0 [1..last]
  in pI2
 
main = print (semiPrimeCount 100000000)

The file primes_to_5M.txt is: 2 3 5 .. 49999991 50000017

[edit] 8 Problem 188

hyperexponentiation

The idea here is actually very simple. Euler's theorem tells us that a^phi(n) mod n = 1, so the exponent never need grow above 40000000 in this case, which is phi(10^8). Henry Laxen April 6, 2008

fastPower x 1 modulus = x `mod` modulus
fastPower x n modulus
  | even n = fastPower ((x*x) `mod` modulus) (n `div` 2) modulus
  | otherwise = (x * (fastPower x (n-1) modulus)) `mod` modulus
 
modulus = 10^8
phi = 40000000 -- eulerTotient of 10^8
 
a :: Int ->  Int  -> Int
a n 1 = n
-- a n k = n^(a n (k-1))
a n k = fastPower n (a n (k-1) `mod` phi) modulus
 
problem_188 = a 1777 1855

[edit] 9 Problem 190

Maximising a weighted product

For those of us who can reach all the back into our sophomore calculus class memories, we will find the dreaded Lagrange Multipliers, which say that if you want to max(min)imize f subject to a constraint g, then the gradients of f and g have to be parallel, or grad f = L * grad g.

Happy happy, for in our case f = x * y^2 * .. z^n and even better g = x + y + .. z - n, so grad g = (1,1,..1) and grad f = (f/x, 2f/y, .. nf/z) That leads us to:

x = f/L, y = 2f/L .. z = nf/L and sum (1..n) x == n so x = 2/(n+1)

Thus the Haskell program is nothing more than:

f l = product $ zipWith (^) l [1..]
 
fmin n = floor $ f [ i*2/(n+1) | i<-[1..n]]
 
problem_190 = sum $ map fmin [2..15]

Aren't you glad you didn't skip that math class now? HenryLaxen Apr 18,2008