### 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=187import NumberTheory.Sieve.ONeillprocess [] p2 acc = accprocess (p:p1) p2 acc = process p1 (tail p2) \$ acc + calc p p2calc 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 0main = 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=230import Prelude hiding ((!!))import Data.List (genericIndex)(!!) = genericIndexfib :: Integer -> Integerfib = (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 0a = "1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"b = "8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"get :: Integer -> Integerget 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 ]

 http://projecteuler.net/index.php?section=problems&id=203import 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 backfiftyOne = sort \$ foldr union [] \$ take 52 rows-- This was my guess that turned out to be correctprimess = 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

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
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
print p

main = do
args <- getArgs
if length args > 0
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 ]

 http://projecteuler.net/index.php?section=problems&id=205 -- 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..]) !! ntotal = 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=209Although 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.Bitsimport Data.Listset :: [(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 = 1loop 2 = 3loop 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 = 1path 1 = 2path 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