Project Euler

Problem 187

posted Sep 29, 2009, 7:31 PM by Wei Hu   [ updated Sep 29, 2009, 7:32 PM ]

http://projecteuler.net/index.php?section=problems&id=187

import NumberTheory.Sieve.ONeill

process [] p2 acc = acc
process (p:p1) p2 acc = process p1 (tail p2) $ acc + calc p p2

calc p p2 = let maxx = 10^8
                limit = maxx `div` p
                go [] acc = acc
                go (p:ps) acc = if p>limit then acc
                                else go ps (acc+1)
            in go p2 0

main = let p2 = primesToLimit (5*10^7)
           p1 = primesToLimit (10^4)
           ans = process p1 p2 0
       in print ans

Problem 207

posted Feb 4, 2009, 6:41 PM by Wei Hu   [ updated Feb 4, 2009, 7:50 PM ]


import Data.Int
import Data.Ratio
import Data.List

{-                                                                                                              
2^t = [2..]                                                                                                     
k = (2^t)^2 - (2^t)                                                                                             
-}
k :: Int -> Int64
k i = x*x - x
    where x = toEnum (i+2)

{-
Perfect?  Y     Y           Y
2^t  ==  [2, 3, 4, 5, 6, 7, 8, ..]
-}
p :: [Rational]
p = snd $ unzip $ tail $ scanl f ((2,0,0),0) [2..]
    where f ((power,perfect,total),_) x = ((power',perfect',total'), perfect'%total')
              where power' = if power == x then power*2 else power
                    perfect' = if power == x then perfect + 1 else perfect
                    total' = total + 1

main = print ans
    where ans = k index
          Just index = findIndex (< 1%12345) p

A much faster solution by daniel.is.fischer

{-# LANGUAGE BangPatterns #-}
 
-- highest bit or integer log2 of a positive integer
-- since GHC implements bit-shifting for Integer via multiplication and division,
-- we do it directly
-- the overhead of memoising powers of 2 and shifts seems smaller than that of doing divisions on the way down
highestBit :: Integer -> Integer
highestBit n
    | n < 1     = error "highestBit needs positive n"
    | n == 1    = 0
    | n < 4     = 1
    | n < 8     = 2
    | n < 16    = 3
    | n < 32    = 4
    | otherwise = go (0 :: Int) [(2,1)] n
      where
        go shft _ 1 = fromIntegral shft
        go !shft lst@((p1,s1):tl) m
            | m < p1    = go shft tl m
            | m < p0    = go (shft+s1) tl (m `quot` p1)
            | otherwise = go (shft+s0) ((p0,s0):lst) (m `quot` p0)
              where
                p0 = p1*p1
                s0 = s1+s1
 
-- find the smallest n for which the proportion of perfect partitions is below
-- tnum/tden
firstBelow :: Integer -> Integer -> [(Integer,Integer)]
firstBelow tnum tden
    | tnum < 1 || tden < 1  = error "expecting positive arguments"
    | tden < tnum           = [(2,2)]
    | tden == tnum          = [(3,6)]
    | otherwise             = go 2 1
      where
        g = gcd tnum tden
        tn = tnum `quot` g
        td = tden `quot` g
        go u r
            | nu == u   = [(u,u*(u-1))]
            | u < nu    = (r,u):go nu nr
            | otherwise = (r,u):(nr,nu):error "Heck"
          where
            nu = (r*td) `quot` tn + 2
            nr = highestBit nu

main = print $ snd $ last $ firstBelow 1 12345

Problem 204

posted Feb 1, 2009, 11:35 PM by Wei Hu   [ updated Feb 2, 2009, 12:01 AM ]

http://projecteuler.net/index.php?section=problems&id=204

A simple recursive solution.

primes = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]

-- number of ways to combine the primes not exceeding x, excluding 1
count :: [Int] -> Int -> Integer
count (p:ps) x = go x
    where go x
              | x >= p  = 1 + go (div x p) + count ps x
              | otherwise = 0
count [] _ = 0

main = print $ (1+) $ count primes (10^9)

{- For debugging

count' :: [Int] -> Int -> [Int]
count' (p:ps) x = go p x
    where go prod x
              | x >= p  = [prod] ++ go (prod*p) (div x p) ++ (map ((`div` p) . (*prod)) $ count' ps x)
              | otherwise = []
count' [] _ = []
-}

Some crazy code that computes the list directly, by suspect42

merge xs [] = xs 
merge [] ys = ys 
merge xss@(x:xs) yss@(y:ys) = if x == y then x : merge xs ys
                                else if x < y
                                     then x : merge xs yss
                                     else y : merge xss ys   
s = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]   
h = 1 : foldr1 merge [map (*n) h | n <- s]   
main = print $ length $ takeWhile (<=10^9) h

Problem 230

posted Jan 31, 2009, 10:55 PM by Wei Hu

http://projecteuler.net/index.php?section=problems&id=230

import Prelude hiding ((!!))
import Data.List (genericIndex)
(!!) = genericIndex

fib :: Integer -> Integer
fib = (100 *) . (map fib' [0..] !!)
    where fib' = fst . fib2
          fib2 0 = (1, 1)
          fib2 1 = (1, 2)
          fib2 n
              | even n    = (a*a + b*b, c*c - a*a)
              | otherwise = (c*c - a*a, b*b + c*c)
              where (a,b) = fib2 (n `div` 2 - 1)
                    c     = a + b

-- Given a number x, locate which fib will be at least x.
locate x = let go i
                   | fib i >= x = i
                   | otherwise = go (i+1)
           in
             go 0

a = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
b = "8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"

get :: Integer -> Integer
get x = read [go x $ locate x]
    where go x 0 = a !! (x-1)
          go x 1 = b !! (x-1)
          -- Find the xth number in the yth fib
          go x y
              | u >= x = go x (y-2)
              | otherwise = go (x-u) (y-1)
              where u = fib (y-2)

main = print $ sum $ zipWith (*) xs ys
    where g n = (127 + 19*n) * 7^n
          xs = map get $ map g [0..17]
          ys = 1:(map (*10) ys)

Problem 203

posted Jan 7, 2009, 6:42 PM by Wei Hu   [ updated Jan 7, 2009, 7:54 PM ]


import Data.List

-- Dropped the common 1's                                                                                       
rows :: [[Integer]]
rows = [[1], [], [], [2]] ++ map computeRow [4..]

-- A better way to compute the next row is something like zipWith (+) (0:list) (list ++ [0])
computeRow n
    | even n = this
    | otherwise = this ++ [2 * last above]
    where this = tail $ fst $ unzip intermediate
          intermediate = scanl f (1,1) above
          f (_, y) z = (y+z, z)
          above = rows !! (n - 1)

-- Instead of foldr union, we can also nub . concat
-- Or, convert to a Set and then back
fiftyOne = sort $ foldr union [] $ take 52 rows

-- This was my guess that turned out to be correct
primess = map (\x -> x * x) [2,3,5,7,11,13,17,19]

-- Conservatively, we can sieve all the primes up to the 4th root of the maximum value in the distinct list.
-- Better yet, any prime dividing C(n,k) is ≤ n, so you need only the primes up to 47.
-- Even sweeter, we only need to sieve [2,3,5,7] since 112 > 51.
answer = foldr f fiftyOne primess
    where f pp = filter (\x -> x `mod` pp /= 0)

main = print answer

Problem 220

posted Dec 14, 2008, 5:26 PM by Wei Hu   [ updated Jan 31, 2009, 8:41 PM ]

http://projecteuler.net/index.php?section=problems&id=220

UPDATE: See this post for code drawing the curve.

Computing the string

D0 = Fa                                                                                                         
D1 = FaRbFR                                                                                                     
                                                                                                                
a -> aRbFR                                                                                                      
b -> LFaLb                                                                                                      
                                                                                                                
z = FaRb -> (FaRb)(FRRLFaLb) = (FaRb)(FR)(FaLb) = z(FR)u                                                        
u = FaLb -> (FaRb)(FRLLFaLb) = (FaRb)(FL)(FaLb) = z(FL)u                                                        
                                                                                                                
X(1) = FR                                                                                                        
Y(1) = FL                                                                                                        
X(n) = X(n-1) ++ FR ++ Y(n-1)                                                                                   
Y(n) = X(n-1) ++ FL ++ Y(n-1)                                                                                   

Steps(1) = 1
Steps(n) = 2 * Steps(n-1) + 1

Expanding Steps, we get Steps(n) = 2^n - 1
                                                                                                               
Dn = X(n) ++ FR            

We're asked to calculate the position of the cursor after 1012 steps in D50, so we don't have to compute the whole string. The problem is equivalent to calculate the position of the cursor after 1012 steps in D40. Still, the problem is not simplified much.

After looking up wikipedia, I found that I overlooked a very important fact: Y(n) is the reverse of X(n) with L and R swapped.

import Data.Bits

data Turn = FL | FR
            deriving Show

data Direction = U | D | L | R
                 deriving Show

class DeepSeq a where
    deepSeq :: a -> b -> b
instance  DeepSeq Integer   where  deepSeq = seq
instance DeepSeq Direction where deepSeq = seq
instance  (DeepSeq a,DeepSeq b,DeepSeq c) => DeepSeq (a,b,c)  where
    deepSeq (a,b,c)         y = deepSeq a $ deepSeq b $ deepSeq c y

move :: (Integer,Integer,Direction) -> Turn -> (Integer,Integer,Direction)
move (x,y,U) FL = (x,y+1,L)
move (x,y,D) FL = (x,y-1,R)
move (x,y,L) FL = (x-1,y,D)
move (x,y,R) FL = (x+1,y,U)
move (x,y,U) FR = (x,y+1,R)
move (x,y,D) FR = (x,y-1,L)
move (x,y,L) FR = (x-1,y,U)
move (x,y,R) FR = (x+1,y,D)

decide :: Integer -> Turn
decide n = if turn then FR else FL
    where turn = (((n .&. (-n)) `shiftL` 1) .&. n) == 0

{- A hand-written tail-recursive function
count :: Integer -> (Integer,Integer,Direction)
count n = go 1 initial
    where initial = (0,0,U)
          go i pos = if i == n then pos'
                     else go (i+1) pos'
                         where pos' = move pos turn
                               turn = decide i
-}

foldl' f a []     = a
foldl' f a (x:xs) = let a' = f a x in a' `deepSeq` foldl' f a' xs

count n = foldl' move initial turns
    where turns = map decide [1..n]
          initial = (0,0,U)

main = print $ count (10^12)

At first I simply used foldl' from Data.List that uses seq to force evaluation. Turns out that it's still too lazy: it builds up too many thunks that eats my memory away. Why? Because (a,b,c) is already in WHNF, so seq is not going to evaluate a, b, or c. It's a shame that ghc -O2 didn't do the strictness analysis deeply enough. Although deepSeq fixes our memory problem, the program is still very inefficient: it computes the positions sequentially.

To be able to interrupt the program and pick up the progress later, I wrote the following ugly imperative program:
import Data.Bits
import Data.IORef
import System.Posix.Signals
import System.Environment
import System.IO.Unsafe

data Turn = FL | FR
            deriving Show

data Direction = U | D | L | R
                 deriving (Read, Show)

deepSeq (d,(a,b,c)) y = seq a $ seq b $ seq c $ seq d y

move :: (Integer,Integer,Direction) -> Turn -> (Integer,Integer,Direction)
move (x,y,U) FL = (x,y+1,L)
move (x,y,D) FL = (x,y-1,R)
move (x,y,L) FL = (x-1,y,D)
move (x,y,R) FL = (x+1,y,U)
move (x,y,U) FR = (x,y+1,R)
move (x,y,D) FR = (x,y-1,L)
move (x,y,L) FR = (x-1,y,U)
move (x,y,R) FR = (x+1,y,D)

decide :: Integer -> Turn
decide n = if turn then FR else FL
    where turn = (((n .&. (-n)) `shiftL` 1) .&. n) == 0

count n = do
  (i, p) <- readIORef progress
  if i == n then print p
   else do
     let x = (i+1, move p $ decide (i+1))
     x `deepSeq` writeIORef progress x
     count n

progress :: IORef (Integer, (Integer, Integer, Direction))
progress = unsafePerformIO $ newIORef (0, (0,0,U))

handler = do
  p <- readIORef progress
  print p

main = do
  args <- getArgs
  if length args > 0
   then writeIORef progress (read $ head args)
   else return ()
  installHandler sigINT (CatchOnce handler) Nothing
  count (10^12)

Finally, after about 10 days, the computer worked out the answer. With the correct answer, I can now find out more efficient algorithms at the discussion board. Turns out that I overlooked yet another important fact: at nth iteration, a duplicate of the graph at (n-1)th stage is rotated anti-clockwise by 90 degrees with the center of rotation at the end point of the (n-1)th graph. This observation reduces O(n) to O(log n). One user even posted the solution at the Wikipedia page. I find that Ruby program pretty easy to follow, once I understood the theory. Here's a solution by BrianHuffman, using slightly different properties:
-- dragon' n = (dragon n, dragon (n+1))
dragon' 0 = ((0, 0), (0, 1))
dragon' n
  | even n = ((x+y, y-x), (x+y, v-u))
  | otherwise = ((x+y, v-u), (u+v, v-u))
  where
    ((x, y), (u, v)) = dragon' (n `div` 2)
 
dragon :: Integer -> (Integer, Integer)
dragon n = fst (dragon' n)

main = print $ dragon (10^12)

Problem 205

posted Sep 23, 2008, 9:33 PM by Wei Hu   [ updated Sep 23, 2008, 9:52 PM ]


-- number of ways to sum to n using 9x4 dices                                                                                
way94' n
    | n < 9 || n > 36 = 0
    | otherwise = sum [1 | a<-x, b<-x, c<-x, d<-x, e<-x, f<-x, g<-x, h<-x, i<-[n-a-b-c-d-e-f-g-h], i>0, i<5]
    where x = [1..4]

way94 n = (map way94' [0..]) !! n

-- number of ways to sum to n using 6x6 dices                                                                                
way66' n
    | n < 6 || n > 36 = 0
    | otherwise = sum [1 | a<-x, b<-x, c<-x, d<-x, e<-x, f<-[n-a-b-c-d-e], f>0, f<7]
    where x = [1..6]

way66 n = (map way66' [0..]) !! n

total = 4^9 * 6^6

-- number of ways to throw 9x4 dices to be higher than x                                                                     
higher x = sum $ map way94 [x .. 36]

ans = sum $ zipWith (*) (map way66 [1..35]) (map higher [2..36])

main = print (ans/total)

A smarter solution can use some math.

Problem 209

posted Sep 21, 2008, 11:47 PM by Wei Hu   [ updated Sep 21, 2008, 11:55 PM ]

http://projecteuler.net/index.php?section=problems&id=209

Although I solved this problem, I didn't formally prove its correctness. A forum member pointed out an important fact, that the map sending (a,b,c,d,e,f) -> (b,c,d,e,f,a xor (b and c)) 
is invertible, so it gives a permutation of the 64 possible inputs. Thus, we can partition the 64 inputs into groups. Here gives many ways of calculating Fib.

import Data.Bits
import Data.List

set :: [(Int, Int, Int, Int, Int, Int)]
set = [(a,b,c,d,e,f) | a<-x, b<-x, c<-x, d<-x, e<-x, f<-x]
    where x = [0, 1]

-- we divide the whole set into loops                                                                                        
loops = go set
    where go [] = []
          go (x:xs) = y : (go $ xs \\ y)
              where y = count x

-- for a given 6-tuple, find what it transforms to until we hit a loop                                                       
count x = go (transform x) [x]
    where transform (a,b,c,d,e,f) = (b,c,d,e,f, a `xor` (b .&. c))
          go start xs = if start == head xs then xs
                        else go (transform start) (xs ++ [start])

-- for a loop, the constraint is that AND of every adjacent pair is 0                                                        
loop 1 = 1
loop 2 = 3
loop n = path' (n-3) + -- pick one and assign 1                                                                              
         path' (n-1)   -- pick one and assign 0      

-- the path has n inner nodes, and its ends are both 0                                                                       
path 0 = 1
path 1 = 2
path n = path (n-1) + -- pick one and assign 0                                                                               
         path (n-2)   -- pick one and assign 1                                                                               

-- a more efficient impl                                                                                                     
path' n = fibs !! n
    where fibs = scanl (+) 1 (1:fibs)

main = print $ product $ map (loop . length) loops

Problem 206

posted Sep 21, 2008, 10:51 PM by Wei Hu   [ updated Sep 21, 2008, 10:59 PM ]


The first answer took me more time to work out. But it's very fast (0.5s). Most of other people solved it using brute force. So I quickly wrote the second answer, but the running time is 25s.

1.hs

import Data.List

square x = x * x
nth x n = (x `div` (10^n)) `mod` 10

-- n digits can safely only determine (n+1) digits squared                                                                   
-- the lowest two digits must be 30 or 70, so we start from there.                                                           

-- before this function, we got n digits from pool, to determine the (n+1)th bit.                                                               
-- to determine the (n+3)th digit squared, we need to prefix the pool with two more digits                                   
-- after this function, we have (n+2) digits that determine the (n+3)th digit
compute :: Integer -> [Integer] -> [Integer]
compute n pool = get (map f pool) pool
    where get [] [] = []
          get ([]:highs) (low:lows) = get highs lows
          get (hh:highs) (low:lows) = xs ++ (get highs lows)
              where xs = map my hh
                    my h = h*(10^n) + low
          f low = filter g [0..99]
              where g x = nth (square (x*(10^n) + low)) (n+2) == 9 - (div n 2)

-- to determine the 19th digit, we need to compute 16 $ compute 14 $ ...                                                     
go 0 = [30, 70]
go n = compute n $ go (n-2)

-- but, since the final answer must be 10 digits, we don't need to go that far                                               
-- go 8 only gives us assurance about 11 digits: 5 _ 6 _ 7 _ 8 _ 9 _ 0                                                       
main = print $ filter (\x -> let y = square x in
                             nth y 18 == 1 &&
                             nth y 16 == 2 &&
                             nth y 14 == 3 &&
                             nth y 12 == 4 &&
                             y < 2 * 10^18
                      )
       $ go 8

2.hs

f ['1',_, '2',_, '3',_, '4',_, '5',_, '6',_, '7',_, '8',_, '9',_, '0'] = True
f _ = False

main = print $ filter (f . show . square) $ filter g [low, low+10..high]
    where low = 1000000030
          high = 1500000070
          square x = x * x
          g x = y == 30 || y == 70
              where y = x `mod` 100

1-9 of 9