Euler problems/181 to 190
From HaskellWiki
(→Problem 186: efficient solution using STUArray) 

(14 intermediate revisions by 8 users not shown)  
Line 1:  Line 1:  
== [http://projecteuler.net/index.php?section=problems&id=181 Problem 181] == 
== [http://projecteuler.net/index.php?section=problems&id=181 Problem 181] == 

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

+  
+  {{sectstub}} 

+  
+  == [http://projecteuler.net/index.php?section=problems&id=182 Problem 182] == 

+  RSA encryption. 

Solution: 
Solution: 

<haskell> 
<haskell> 

−  import Data.Map ((!),Map) 
+  fun a1 b1 = sum [ e  e < [2..a*b1], 
−  import qualified Data.Map as M 
+  gcd e (a*b) == 1, 
+  gcd (e1) a == 2, 

+  gcd (e1) b == 2 ] 

+  where a = a11 

+  b = b11 

+  
+  problem_182 = fun 1009 3643 

+  </haskell> 

+  
+  == [http://projecteuler.net/index.php?section=problems&id=183 Problem 183] == 

+  Maximum product of parts. 

+  
+  Solution: 

+  <haskell> 

+   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 

+  </haskell> 

+  
+  == [http://projecteuler.net/index.php?section=problems&id=184 Problem 184] == 

+  Triangles containing the origin. 

+  
+  {{sectstub}} 

+  
+  == [http://projecteuler.net/index.php?section=problems&id=185 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 20080312 

+  
+  <haskell> 

+  
import Data.List 
import Data.List 

import Control.Monad 
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 :: IO () 

main = do 
main = do 

−  let es = [40,60] 
+  let ex = ex1 
−  dg = sum es 
+  x < randomStart (length (fst (head ex))) 
−  mon = Mon dg es 
+  let y = iterateMind ex x 
−  Poly mp = partitionPol mon 
+  done = snd (last y) == 0 
−  print $ mp!mon 
+  if done then (putStrLn $ (fst.last) y) 
−  +  else main 

−  data Monomial 

−  = Mon 

−  { degree :: !Int 

−  , expos :: [Int] 

−  } 

−  
−  infixl 7 <*>, *> 

−  
−  (<*>) :: Monomial > Monomial > Monomial 

−  (Mon d1 e1) <*> (Mon d2 e2) 

−  = Mon (d1+d2) (zipWithZ (+) e1 e2) 

−  
−  unit :: Monomial 

−  unit = Mon 0 [] 

−  
−  (<<) :: Monomial > Monomial > Bool 

−  (Mon d1 e1) << (Mon d2 e2) 

−  = d1 <= d2 && and (zipWithZ (<=) e1 e2) 

−  
−  upTo :: Monomial > [Monomial] 

−  upTo (Mon 0 _) = [unit] 

−  upTo (Mon d es) = 

−  sort $ go 0 [] es 

−  where 

−  go dg acc [] = return (Mon dg $ reverse acc) 

−  go dg acc (n:ns) = do 

−  k < [0 .. n] 

−  go (dg+k) (k:acc) ns 

−  
−  newtype Polynomial = 

−  Poly { mapping :: (Map Monomial Integer) } 

−  deriving (Eq, Ord) 

−  
−  (*>) :: Integer > Monomial > Polynomial 

−  n *> m = Poly $ M.singleton m n 

−  
−   

−   The hard stuff  

−   

−  
−  one :: Map Monomial Integer 

−  one = M.singleton unit 1 

−  reciprocal :: Monomial > Polynomial 
+  </haskell> 
−  reciprocal m = 

−  Poly . foldl' extend one . reverse . drop 1 . upTo $ m 

−  where 

−  extend mp mon = 

−  M.filter (/= 0) $ 

−  foldl' (flip (uncurry $ M.insertWith' (+))) mp list 

−  where 

−  list = filter ((<< m) . fst) [(mon <*> mn, c)  

−  (mn,c) < M.assocs mp] 

−  
−  partitionPol :: Monomial > Polynomial 

−  partitionPol m = 

−  Poly . foldl' update one $ sliced m 

−  where 

−  Poly rec = reciprocal m 

−  sliced mon = sortBy (comparing expos) . drop 1 $ upTo mon 

−  comparing f x y = compare (f x) (f y) 

−  update mp mon@(Mon d es) 

−   es /= ses = M.insert mon (mp!(Mon d ses)) mp 

−   otherwise = M.insert mon (negate clc) mp 

−  where 

−  ses = sort es 

−  clc = sum $ do 

−  mn@(Mon dg xs) < sliced mon 

−  let cmn = Mon (ddg) (zipWithZ () es xs) 

−  case M.lookup mn rec of 

−  Nothing > [] 

−  Just c > return $ c*(mp!(Mon (ddg) 

−  (zipWithZ () es xs))) 

−   
+  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 backtracking search, but with some semismart 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. 
−   Auxiliary Functions  
+  
−   
+  <haskell> 
−  +  import Control.Monad 

−  zipWithZ :: (Int > Int > a) > [Int] > [Int] > [a] 
+  import Data.List 
−  zipWithZ _ [] [] = [] 
+  import qualified Data.Set as S 
−  zipWithZ f [] ys = map (f 0) ys 
+  
−  zipWithZ f xs [] = map (flip f 0) xs 
+  ensure p x = guard (p x) >> return x 
−  zipWithZ f (x:xs) (y:ys) = f x y:zipWithZ f xs ys 
+  selectDistinct 0 _ = [[]] 
−  +  selectDistinct n [] = [] 

−  unknowns :: [String] 
+  selectDistinct n (x:xs) = map (x:) (selectDistinct (n  1) xs) ++ selectDistinct n xs 
−  unknowns = ['X':show i  i < [1 .. ]] 
+  
−  +  data Environment a = Environment { guesses :: [(Int, [a])] 

−  instance Show Monomial where 
+  , restrictions :: [S.Set a] 
−  showsPrec _ (Mon 0 _) = showString "1" 
+  , assignmentsLeft :: Int 
−  showsPrec _ (Mon _ es) = foldr (.) id $ intersperse (showString "*") us 
+  } 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 
where 

−  ps = filter ((/= 0) . snd) $ zip unknowns es 
+  f = ensure sanityCheck . update e 
−  us = map (\(s,e) > showString s . showString "^" 
+  m = selectDistinct n (positions g (restrictions e)) 
−  . showParen (e < 0) (shows e)) ps 
+  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) 

−  instance Eq Monomial where 
+  update (Environment ((_, a):gs) r l) is = Environment newGs restriction (l  length is) where 
−  (Mon d1 e1) == (Mon d2 e2) 
+  (assignment, newRestriction) = splitAscIndices is a 
−  = d1 == d2 && (d1 == 0  e1 == e2) 
+  (_, oldRestriction) = splitAscIndices is r 
−  +  restriction = zipWith S.insert newRestriction oldRestriction 

−  instance Ord Monomial where 
+  newGs = map updateEntry gs 
−  compare (Mon d1 e1) (Mon d2 e2) 
+  updateEntry (n', a') = (newN, newA) where 
−  = case compare d1 d2 of 
+  (dropped, newA) = splitAscIndices is a' 
−  EQ  d1 == 0 > EQ 
+  newN = n'  length (filter id $ zipWith (==) assignment dropped) 
−   otherwise > compare e2 e1 

−  other > other 

−  
−  instance Show Polynomial where 

−  showsPrec p (Poly m) 

−  = showP p . filter ((/= 0) . snd) $ M.assocs m 

−  showP :: Int > [(Monomial,Integer)] > ShowS 
+  problem_185 = head . solve $ initial 
−  showP _ [] = showString "0" 

−  showP p cs = 

−  showParen (p > 6) showL 

−  where 

−  showL = foldr (.) id $ intersperse (showString " + ") ms 

−  ms = map (\(m,c) > showParen (c < 0) (shows c) 

−  . showString "*" . shows m) cs 

−  
−  instance Num Polynomial where 

−  (Poly m1) + (Poly m2) = Poly (M.filter (/= 0) $ addM m1 m2) 

−  p1  p2 = p1 + (negate p2) 

−  (Poly m1) * (Poly m2) = Poly (mulM (M.assocs m1) (M.assocs m2)) 

−  negate (Poly m) = Poly $ M.map negate m 

−  abs = id 

−  signum = id 

−  fromInteger n 

−   n == 0 = Poly (M.empty) 

−   otherwise = Poly (M.singleton unit n) 

−  
−  addM :: Map Monomial Integer > Map Monomial Integer > Map Monomial Integer 

−  addM p1 p2 = 

−  foldl' (flip (uncurry (M.insertWith' (+)))) p1 $ 

−  M.assocs p2 

−  
−  mulM :: [(Monomial,Integer)] > [(Monomial,Integer)] > Map Monomial Integer 

−  mulM p1 p2 = 

−  M.filter (/= 0) . 

−  foldl' (flip (uncurry (M.insertWith' (+)))) M.empty $ 

−  liftM2 (\(e1,c1) (e2,c2) > (e1 <*> e2,c1*c2)) p1 p2 

−  problem_181 = main 

</haskell> 
</haskell> 

−  == [http://projecteuler.net/index.php?section=problems&id=182 Problem 182] == 
+  == [http://projecteuler.net/index.php?section=problems&id=186 Problem 186] == 
−  RSA encryption. 
+  Connectedness of a network. 
+  
+  This implements the known tree structure for imperative [http://en.wikipedia.org/wiki/Disjointset_data_structure union/find]. 

+  On my machine it runs in about 3 seconds. 

−  Solution: 

<haskell> 
<haskell> 

−  fun a1 b1 = 
+  module ProjectEuler186 where 
−  sum [ e  
+  
−  e < [2..a*b1], 
+  import qualified Control.Monad.ST.Strict as StrictST 
−  gcd e (a*b) == 1, 
+  import qualified Control.Monad.ST.Lazy as LazyST 
−  gcd (e1) a == 2, 
+  import Data.Array.ST (STUArray, newArray, newListArray, readArray, writeArray, ) 
−  gcd (e1) b == 2 
+  import Control.Monad (liftM2, when, ) 
−  ] 
+  
−  where 
+  
−  a=a11 
+  data Network s = 
−  b=b11 
+  Network 
−  problem_182=fun 1009 3643 
+  (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> 
</haskell> 

−  == [http://projecteuler.net/index.php?section=problems&id=183 Problem 183] == 
+  == [http://projecteuler.net/index.php?section=problems&id=187 Problem 187] == 
−  Maximum product of parts. 
+  Semiprimes 
Solution: 
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 writeup 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> 
<haskell> 

−  pmax x a=a*(log xlog a) 
+  
−  tofloat x=encodeFloat x 0 
+  module PI (prime,primes,pI) where 
−  fun x= 
+  import Control.Monad.ST 
−  div n1 $gcd n1 x 
+  import System.IO.Unsafe 
−  where 
+  import qualified Data.Map as M 
−  e=exp 1 
+  import Data.ByteString.Char8 (readFile,lines,readInt) 
−  n=floor(fromInteger x/e) 
+  import Data.Maybe 
−  n1=snd.maximum$[(b,a)a<[n..n+1],let b=pmax (tofloat x) (tofloat a)] 
+  import Data.Array.ST 
−  n `splitWith` p = doSplitWith 0 n 
+  import Data.Array.IArray 
−  where doSplitWith s t 
+  
−   p `divides` t = doSplitWith (s+1) (t `div` p) 
+  r = runSTUArray 
−   otherwise = (s, t) 
+  (do a < newArray (0,3001134) 0 :: ST s (STUArray s Int Int) 
−  d `divides` n = n `mod` d == 0 
+  writeArray a 0 1 
−  funD x 
+  zipWithM_ (writeArray a) [0..] primes 
−  is25 k=(x) 
+  return a) 
−  otherwise =x 
+  
−  where 
+  {# NOINLINE s #} 
−  k=fun x 
+  s = M.fromDistinctAscList $ map (\(x,y) > (y,x)) (assocs r) 
−  is25 x 
+  
−  s==1=True 
+  prime n = r!n 
−  otherwise=False 
+  
−  where 
+  pI :: Int > Int 
−  s=snd(splitWith (snd (splitWith x 2)) 5) 
+  pI n = 
−  problem_183 =sum[funD aa< [5..10000]] 
+  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 (k1)))  k + 1 

+  
+  semiPrimeCount n = 

+  let last = pI $ floor $ sqrt (fromIntegral n) 

+  s k = pI ( n `div` (prime (k1)))  k + 1 

+  pI2 = foldl' (\x y > s y + x) 0 [1..last] 

+  in pI2 

+  
+  main = print (semiPrimeCount 100000000) 

+  
</haskell> 
</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 (n1) 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 (k1)) 

+  a n k = fastPower n (a n (k1) `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*b1], gcd e (a*b) == 1, gcd (e1) a == 2, gcd (e1) b == 2 ] where a = a11 b = b11 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 20080312
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 backtracking search, but with some semismart 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 writeup 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 (k1)))  k + 1 semiPrimeCount n = let last = pI $ floor $ sqrt (fromIntegral n) s k = pI ( n `div` (prime (k1)))  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 (n1) 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 (k1)) a n k = fastPower n (a n (k1) `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