# Euler problems/191 to 200

### From HaskellWiki

Henrylaxen (Talk | contribs) (solution to problem 193) |
(→Problem 193) |
||

(3 intermediate revisions by one user not shown) | |||

Line 37: | Line 37: | ||

| n == 0 = sum [ hasM_LsAndEndsInN_As (m-1) i (k-1) | i<-[0..2]] -- +L | | n == 0 = sum [ hasM_LsAndEndsInN_As (m-1) i (k-1) | i<-[0..2]] -- +L | ||

+ sum [ hasM_LsAndEndsInN_As m i (k-1) | i<-[0..2]] -- +O | + sum [ hasM_LsAndEndsInN_As m i (k-1) | i<-[0..2]] -- +O | ||

− | | n 0 = hasM_LsAndEndsInN_As m (n-1) (k-1) -- + A | + | | n > 0 = hasM_LsAndEndsInN_As m (n-1) (k-1) -- + A |

-- Count awards of length k that have "m" L's in them, and end in "n" A's | -- Count awards of length k that have "m" L's in them, and end in "n" A's | ||

problem191 n = do | problem191 n = do | ||

let p a b c d = "hasM_LsAndEndsInN_As " ++ | let p a b c d = "hasM_LsAndEndsInN_As " ++ | ||

− | + | unwords (map show [a,b,c]) ++ | |

− | "= " ++ | + | " = " ++ show d |

− | putStrLn $ "award " ++ | + | putStrLn $ "award " ++ show n ++ " = " ++ show (award n) |

mapM_ (\(i,j) -> putStrLn $ p i j n (hasM_LsAndEndsInN_As i j n)) | mapM_ (\(i,j) -> putStrLn $ p i j n (hasM_LsAndEndsInN_As i j n)) | ||

[ (i,j) | i<-[0..1], j<-[0..2]] | [ (i,j) | i<-[0..1], j<-[0..2]] | ||

Line 57: | Line 57: | ||

Before going through the code below, it is important to have a | Before going through the code below, it is important to have a | ||

good understanding of continued fractions. Have a look at the | good understanding of continued fractions. Have a look at the | ||

− | wikipedia article below. Pay particular attention to the | + | wikipedia article below. Pay particular attention to the section |

on semiconvergents, for that is the key to the code in closest2, | on semiconvergents, for that is the key to the code in closest2, | ||

which calculates '''the other''' candidate for closest rational. | which calculates '''the other''' candidate for closest rational. | ||

Line 74: | Line 74: | ||

squareRootPeriod r n d = m : rest | squareRootPeriod r n d = m : rest | ||

where | where | ||

− | m = truncate | + | m = (truncate (sqrt (fromIntegral r)) + n) `div` d |

a = n - d * m | a = n - d * m | ||

− | rest | + | rest | d == 1 && n /= 0 = [] |

− | + | | otherwise = squareRootPeriod r (-a) ((r - a ^ 2) `div` d) | |

− | + | ||

-- Turn the period into an infinite stream | -- Turn the period into an infinite stream | ||

continuedFraction :: [Int] -> [Integer] | continuedFraction :: [Int] -> [Integer] | ||

− | continuedFraction | + | continuedFraction (x:xs) = map fromIntegral $ x : cycle xs |

-- calculate successive convergents as a ratio | -- calculate successive convergents as a ratio | ||

Line 97: | Line 96: | ||

-- calculated as p3%q3 | -- calculated as p3%q3 | ||

closest2 :: Integer -> Integer -> [Ratio Integer] | closest2 :: Integer -> Integer -> [Ratio Integer] | ||

− | closest2 bound n = | + | closest2 bound n = [p,(p3%q3)] |

− | + | where | |

− | b = takeWhile ( | + | a = convergents $ continuedFraction $ squareRootPeriod (fromIntegral n) 0 1 |

+ | b = takeWhile ((<= bound) . denominator) a | ||

c = reverse b | c = reverse b | ||

(p:q:_) = c | (p:q:_) = c | ||

Line 106: | Line 106: | ||

p3 = ((bound-q2) `div` q1) * p1 + p2 | p3 = ((bound-q2) `div` q1) * p1 + p2 | ||

q3 = ((bound-q2) `div` q1) * q1 + q2 | q3 = ((bound-q2) `div` q1) * q1 + q2 | ||

− | |||

-- pick the ratio returned from closest2 which is | -- pick the ratio returned from closest2 which is | ||

-- actually closer to the square root, and return the denominator | -- actually closer to the square root, and return the denominator | ||

denomClosest :: Integer -> Integer -> Integer | denomClosest :: Integer -> Integer -> Integer | ||

− | denomClosest bound n = | + | denomClosest bound n = denominator c1 |

− | + | where [l,r] = closest2 bound n | |

− | + | c1 | abs (l*l - (n%1)) < abs (r*r - (n%1)) = l | |

− | + | | otherwise = r | |

isSquare :: Integer -> Bool | isSquare :: Integer -> Bool | ||

Line 176: | Line 175: | ||

</haskell> | </haskell> | ||

− | which makes f evaluate its | + | which makes f evaluate its arguments. Again, I think the Earth |

would have been cold and dark before that program finished. | would have been cold and dark before that program finished. | ||

That is one complaint I have about Haskell, namely that (at | That is one complaint I have about Haskell, namely that (at | ||

Line 192: | Line 191: | ||

generated. Surely we (haskellers) can do better than the C# | generated. Surely we (haskellers) can do better than the C# | ||

guys, no? HenryLaxen July 2,2008 | guys, no? HenryLaxen July 2,2008 | ||

+ | |||

+ | * [[Prime_numbers#Tree_merging_with_Wheel]] calculates the 2 mln-th prime in just above 3 seconds on [http://ideone.com/p0e81 Ideone.com]. A straightforward [http://ideone.com/fapob C++ code] does the same in 0.16 seconds. [[User:WillNess|WillNess]] 13:14, 1 July 2011 (UTC) | ||

+ | |||

<haskell> | <haskell> | ||

Line 202: | Line 204: | ||

primesSquared :: Array Int Integer | primesSquared :: Array Int Integer | ||

primesSquared = listArray (0,length p2) p2 | primesSquared = listArray (0,length p2) p2 | ||

− | where p2 = | + | where p2 = takeWhile (<= maxN) |

− | (map ( | + | (map (^2) $ 2:[p | p<-[3,5..], isPrime p]) ++ [maxN+1] |

f (result,partialProduct,i,isOdd) | f (result,partialProduct,i,isOdd) | ||

| seq result $ seq partialProduct $ seq i $ seq isOdd $ False = undefined | | seq result $ seq partialProduct $ seq i $ seq isOdd $ False = undefined | ||

− | + | | nextProduct < maxN = | |

− | + | (nextPrimeSquared . descendLevel) (result+adjustCount) | |

+ | | otherwise = result+adjustCount | ||

+ | where | ||

+ | nextProduct = ((primesSquared ! i) * partialProduct) | ||

numberOfSquareFull = maxN `div` nextProduct | numberOfSquareFull = maxN `div` nextProduct | ||

− | adjustCount | + | adjustCount | isOdd = numberOfSquareFull |

− | + | | otherwise = -numberOfSquareFull | |

descendLevel x = f (x, nextProduct, i+1, not isOdd) | descendLevel x = f (x, nextProduct, i+1, not isOdd) | ||

− | nextPrimeSquared x = f (x, partialProduct, i+1, | + | nextPrimeSquared x = f (x, partialProduct, i+1, isOdd) |

− | + | ||

− | + | ||

− | + | ||

main = do | main = do | ||

a <- seq primesSquared $ getCPUTime | a <- seq primesSquared $ getCPUTime | ||

− | let b = maxN - | + | let b = maxN - f (0,1,0,True) |

putStrLn $ "Result = " ++ show b | putStrLn $ "Result = " ++ show b | ||

c <- getCPUTime | c <- getCPUTime | ||

putStrLn $ "Prime Calculation took " ++ show (a `div` 10^12) ++ " seconds" | putStrLn $ "Prime Calculation took " ++ show (a `div` 10^12) ++ " seconds" | ||

− | putStrLn $ "Moebius Calculation took " ++ show ((c-a) `div` 10^12) ++ " seconds" | + | putStrLn $ "Moebius Calculation took " ++ show ((c-a) `div` 10^12) |

+ | ++ " seconds" | ||

</haskell> | </haskell> |

## Latest revision as of 13:14, 1 July 2011

## [edit] 1 Problem 191

Prize Strings

A couple of notes. I was too lazy to memoize this, so I just ran it twice, once with 15 and then again with 30. I pasted the output of the 15 run into the code. The way to get a handle on this is to just case it out. Ask yourself what can I add to award (n-1) to generate award (n). You can add an O to the end of all of award (n-1). You can add an L to any award (n-1) that doesn't contain an L, and you can add an A to award (n-1) provided it doesn't end with two A's. So the function hasM_LsAndEndsInN_As is just what is needed to cover all of the cases. Henry Laxen April 29, 2008

award 1 = 3 award 15 = 107236 award k = award (k-1) -- + O + sum [ hasM_LsAndEndsInN_As 0 i (k-1) | i<-[0..2] ] -- +L + sum [ hasM_LsAndEndsInN_As i j (k-1) | i<-[0,1], j<-[0,1] ] -- +A hasM_LsAndEndsInN_As 0 0 1 = 1 -- O hasM_LsAndEndsInN_As 1 0 1 = 1 -- L hasM_LsAndEndsInN_As 0 1 1 = 1 -- A hasM_LsAndEndsInN_As _ _ 1 = 0 hasM_LsAndEndsInN_As 0 0 15 = 5768 hasM_LsAndEndsInN_As 0 1 15 = 3136 hasM_LsAndEndsInN_As 0 2 15 = 1705 hasM_LsAndEndsInN_As 1 0 15 = 54736 hasM_LsAndEndsInN_As 1 1 15 = 27820 hasM_LsAndEndsInN_As 1 2 15 = 14071 hasM_LsAndEndsInN_As m n k | m < 0 || n < 0 = 0 | n == 0 = sum [ hasM_LsAndEndsInN_As (m-1) i (k-1) | i<-[0..2]] -- +L + sum [ hasM_LsAndEndsInN_As m i (k-1) | i<-[0..2]] -- +O | n > 0 = hasM_LsAndEndsInN_As m (n-1) (k-1) -- + A -- Count awards of length k that have "m" L's in them, and end in "n" A's problem191 n = do let p a b c d = "hasM_LsAndEndsInN_As " ++ unwords (map show [a,b,c]) ++ " = " ++ show d putStrLn $ "award " ++ show n ++ " = " ++ show (award n) mapM_ (\(i,j) -> putStrLn $ p i j n (hasM_LsAndEndsInN_As i j n)) [ (i,j) | i<-[0..1], j<-[0..2]]

A brief tutorial on solving this problem is available here

## [edit] 2 Problem 192

Best Approximations

Before going through the code below, it is important to have a
good understanding of continued fractions. Have a look at the
wikipedia article below. Pay particular attention to the section
on semiconvergents, for that is the key to the code in closest2,
which calculates **the other** candidate for closest rational.
HenryLaxen May 5, 2008

http://en.wikipedia.org/wiki/Continued_fraction

http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion

import Data.List import Data.Ratio -- Continued fraction representations of square roots are periodic -- Here we calculate the periodic expansion squareRootPeriod :: Int -> Int -> Int -> [Int] squareRootPeriod r n d = m : rest where m = (truncate (sqrt (fromIntegral r)) + n) `div` d a = n - d * m rest | d == 1 && n /= 0 = [] | otherwise = squareRootPeriod r (-a) ((r - a ^ 2) `div` d) -- Turn the period into an infinite stream continuedFraction :: [Int] -> [Integer] continuedFraction (x:xs) = map fromIntegral $ x : cycle xs -- calculate successive convergents as a ratio convergents :: [Integer] -> [Ratio Integer] convergents l = zipWith (%) (drop 2 $ hn) (drop 2 $ kn) where hn = 0:1:zipWith3 (\x y z -> x*y+z) l (tail hn) hn kn = 1:0:zipWith3 (\x y z -> x*y+z) l (tail kn) kn -- here are the guts of the solution -- we calculate convergents until the size of the denominator exceeds -- the given bound. This is one candidate for the closest rational -- approximation. The other candidate is a semiconvergent, which is -- calculated as p3%q3 closest2 :: Integer -> Integer -> [Ratio Integer] closest2 bound n = [p,(p3%q3)] where a = convergents $ continuedFraction $ squareRootPeriod (fromIntegral n) 0 1 b = takeWhile ((<= bound) . denominator) a c = reverse b (p:q:_) = c (p1,q1) = (numerator p, denominator p) (p2,q2) = (numerator q, denominator q) p3 = ((bound-q2) `div` q1) * p1 + p2 q3 = ((bound-q2) `div` q1) * q1 + q2 -- pick the ratio returned from closest2 which is -- actually closer to the square root, and return the denominator denomClosest :: Integer -> Integer -> Integer denomClosest bound n = denominator c1 where [l,r] = closest2 bound n c1 | abs (l*l - (n%1)) < abs (r*r - (n%1)) = l | otherwise = r isSquare :: Integer -> Bool isSquare n = n `elem` takeWhile (<= n) [n*n | n <- [x..] ] where x = floor . sqrt . fromIntegral $ n nonSquares :: Integer -> [Integer] nonSquares k = [ n | n<-[2..k] , (not . isSquare) n] bound = (10^12) problem_192 :: Integer problem_192 = sum $ map (denomClosest bound) (nonSquares 100000)

## [edit] 3 Problem 193

Square free numbers Before trying to understand the solution below, you may want to have a look at

The Moebius function:

http://en.wikipedia.org/wiki/Moebius_function and

The Inclusion-exclusion principle:

http://en.wikipedia.org/wiki/Inclusion-exclusion_principle

It turns out they are basically the same thing. To count the number of square free numbers less than N, we instead count the number of numbers less than N that have a square factor (are squareFull) , and then subtract that from N. The number of squareFull numbers less than N that have 2^2 = 4 as a factor is just floor (N/4). To that we add (N/9) , (N/25), (N/49) ... The problem is that we have counted some squareFull numbers twice, namely (N/4*9), (N/4*25), (N/9*25) ... We need to exclude those numbers. This goes on an on, we need to include (or re-include) divisors with an odd number of prime square factors, and we need to exclude divisors with an even number of odd prime factors. Not co-incidentally, this is exactly what the Moebius function does. So, another way of writing the answer is:

N = sum floor(N/k^2) * moebius(k) where the sum runs from 1 to sqrt N

Unfortunately, if you try to compute it this way using the mobius function found in the wonderful Maths library ArithmeticFunctions, the Earth will be a cold dark cinder before it finally finishes.

Instead the calculation below computes all primes squared less than 2^50 and stuffs them into an array for O(1) access. It then goes through the inclusion-exclusion calculation described above, fortunately terminating early once the divisor exceed 2^50.

A couple of comments about the code below. The first time I wrote it I didn't include the lines:

f (result,partialProduct,i,isOdd) | seq result $ seq partialProduct $ seq i $ seq isOdd $ False = undefined

which makes f evaluate its arguments. Again, I think the Earth would have been cold and dark before that program finished. That is one complaint I have about Haskell, namely that (at least to me) it isn't at all obvious what goes on "under the hood" and the performance of programs can change very drastically with seemingly minor changes in the source code.

My second comment is that someone more skilled than me needs to write a really good prime number generator. Reading through the Euler Project Forums I see that the C/C++/C# guys can generate all primes < 2^25 in about 8 seconds. On my machine using the wonderful "Haskell for Maths" library by David Amos, it takes about 260 seconds to generate those primes. The actual Moebius calculation then only take about 47 seconds once the primes are generated. Surely we (haskellers) can do better than the C# guys, no? HenryLaxen July 2,2008

- Prime_numbers#Tree_merging_with_Wheel calculates the 2 mln-th prime in just above 3 seconds on Ideone.com. A straightforward C++ code does the same in 0.16 seconds. WillNess 13:14, 1 July 2011 (UTC)

import Primes (isPrime) import Data.Array.Unboxed (Array, listArray, (!)) import System.CPUTime (getCPUTime) maxN = 2^50 primesSquared :: Array Int Integer primesSquared = listArray (0,length p2) p2 where p2 = takeWhile (<= maxN) (map (^2) $ 2:[p | p<-[3,5..], isPrime p]) ++ [maxN+1] f (result,partialProduct,i,isOdd) | seq result $ seq partialProduct $ seq i $ seq isOdd $ False = undefined | nextProduct < maxN = (nextPrimeSquared . descendLevel) (result+adjustCount) | otherwise = result+adjustCount where nextProduct = ((primesSquared ! i) * partialProduct) numberOfSquareFull = maxN `div` nextProduct adjustCount | isOdd = numberOfSquareFull | otherwise = -numberOfSquareFull descendLevel x = f (x, nextProduct, i+1, not isOdd) nextPrimeSquared x = f (x, partialProduct, i+1, isOdd) main = do a <- seq primesSquared $ getCPUTime let b = maxN - f (0,1,0,True) putStrLn $ "Result = " ++ show b c <- getCPUTime putStrLn $ "Prime Calculation took " ++ show (a `div` 10^12) ++ " seconds" putStrLn $ "Moebius Calculation took " ++ show ((c-a) `div` 10^12) ++ " seconds"