Personal tools

Euler problems/151 to 160

From HaskellWiki

< Euler problems(Difference between revisions)
Jump to: navigation, search
m (problem_152: only need to consider primes less than 80.)
m
 
(28 intermediate revisions by 10 users not shown)
Line 1: Line 1:
[[Category:Programming exercise spoilers]]
+
== [http://projecteuler.net/index.php?section=problems&id=151 Problem 151] ==
== [http://projecteuler.net/index.php?section=view&id=151 Problem 151] ==
 
 
Paper sheets of standard sizes: an expected-value problem.
 
Paper sheets of standard sizes: an expected-value problem.
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_151 = undefined
+
problem_151 = fun (1,1,1,1)
  +
  +
fun (0,0,0,1) = 0
  +
fun (0,0,1,0) = fun (0,0,0,1) + 1
  +
fun (0,1,0,0) = fun (0,0,1,1) + 1
  +
fun (1,0,0,0) = fun (0,1,1,1) + 1
  +
fun (a,b,c,d) =
  +
(pickA + pickB + pickC + pickD) / (a + b + c + d)
  +
where
  +
pickA | a > 0 = a * fun (a-1,b+1,c+1,d+1)
  +
| otherwise = 0
  +
pickB | b > 0 = b * fun (a,b-1,c+1,d+1)
  +
| otherwise = 0
  +
pickC | c > 0 = c * fun (a,b,c-1,d+1)
  +
| otherwise = 0
  +
pickD | d > 0 = d * fun (a,b,c,d-1)
  +
| otherwise = 0
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=152 Problem 152] ==
+
== [http://projecteuler.net/index.php?section=problems&id=152 Problem 152] ==
 
Writing 1/2 as a sum of inverse squares
 
Writing 1/2 as a sum of inverse squares
   
Line 18: Line 18:
 
import Data.Ratio
 
import Data.Ratio
 
import Data.List
 
import Data.List
+
import Data.Ord (comparing)
  +
import Data.Function (on)
  +
 
invSq n = 1 % (n * n)
 
invSq n = 1 % (n * n)
 
sumInvSq = sum . map invSq
 
sumInvSq = sum . map invSq
+
 
subsets (x:xs) = let s = subsets xs in s ++ map (x :) s
 
subsets (x:xs) = let s = subsets xs in s ++ map (x :) s
 
subsets _ = [[]]
 
subsets _ = [[]]
+
 
primes = 2 : 3 : 7 : [p | p <- [11, 13..79],
 
primes = 2 : 3 : 7 : [p | p <- [11, 13..79],
 
all (\q -> p `mod` q /= 0) [3, 5, 7]]
 
all (\q -> p `mod` q /= 0) [3, 5, 7]]
+
 
-- All subsets whose sum of inverse squares,
 
-- All subsets whose sum of inverse squares,
 
-- when added to x, does not contain a factor of p
 
-- when added to x, does not contain a factor of p
 
pfree s x p = [(y, t) | t <- subsets s, let y = x + sumInvSq t,
 
pfree s x p = [(y, t) | t <- subsets s, let y = x + sumInvSq t,
 
denominator y `mod` p /= 0]
 
denominator y `mod` p /= 0]
+
-- Verify that we need not consider terms divisible by 11, or by any
+
-- prime greater than 13. Nor need we consider any term divisible
 
-- by 25, 27, 32, or 49.
 
verify = all (\p -> null $ tail $ pfree [p, 2*p..85] 0 p) $
 
11 : dropWhile (< 17) primes ++ [25, 27, 32, 49]
 
 
 
-- All pairs (x, s) where x is a rational number whose reduced
 
-- All pairs (x, s) where x is a rational number whose reduced
 
-- denominator is not divisible by any prime greater than 3;
 
-- denominator is not divisible by any prime greater than 3;
 
-- and s is all sets of numbers up to 80 divisible
 
-- and s is all sets of numbers up to 80 divisible
 
-- by a prime greater than 3, whose sum of inverse squares is x.
 
-- by a prime greater than 3, whose sum of inverse squares is x.
only23 = foldl f [(0, [[]])] [13, 7, 5]
+
only23 = foldl fun [(0, [[]])] [13, 7, 5]
where
+
where
f a p = collect $ [(y, u ++ v) | (x, s) <- a,
+
fun a p =
(y, v) <- pfree (terms p) x p,
+
collect $ [(y, u ++ v) |
u <- s]
+
(x, s) <- a,
terms p = [n * p | n <- [1..80`div`p],
+
(y, v) <- pfree (terms p) x p,
all (\q -> n `mod` q /= 0) $
+
u <- s]
11 : takeWhile (>= p) [13, 7, 5]
+
terms p =
]
+
[n * p |
collect = map (\z -> (fst $ head z, map snd z)) .
+
n <- [1..80`div`p],
groupBy fstEq . sortBy cmpFst
+
all (\q -> n `mod` q /= 0) $
fstEq (x, _) (y, _) = x == y
+
11 : takeWhile (>= p) [13, 7, 5]
cmpFst (x, _) (y, _) = compare x y
+
]
+
collect =
  +
map (\z -> (fst $ head z, map snd z)) .
  +
groupBy fstEq . sortBy cmpFst
  +
fstEq = (==) `on` fst
  +
cmpFst = comparing fst
  +
 
-- All subsets (of an ordered set) whose sum of inverse squares is x
 
-- All subsets (of an ordered set) whose sum of inverse squares is x
findInvSq x y = f x $ zip3 y (map invSq y) (map sumInvSq $ init $ tails y)
+
findInvSq x y =
where
+
fun x $ zip3 y (map invSq y) (map sumInvSq $ init $ tails y)
f 0 _ = [[]]
+
where
f x ((n, r, s):ns)
+
fun 0 _ = [[]]
| r > x = f x ns
+
fun x ((n, r, s):ns)
| s < x = []
+
| r > x = fun x ns
| otherwise = map (n :) (f (x - r) ns) ++ f x ns
+
| s < x = []
f _ _ = []
+
| otherwise = map (n :) (fun (x - r) ns) ++ fun x ns
+
fun _ _ = []
  +
 
-- All numbers up to 80 that are divisible only by the primes
 
-- All numbers up to 80 that are divisible only by the primes
 
-- 2 and 3 and are not divisible by 32 or 27.
 
-- 2 and 3 and are not divisible by 32 or 27.
 
all23 = [n | a <- [0..4], b <- [0..2], let n = 2^a * 3^b, n <= 80]
 
all23 = [n | a <- [0..4], b <- [0..2], let n = 2^a * 3^b, n <= 80]
+
solutions = if verify
+
solutions =
then [sort $ u ++ v | (x, s) <- only23,
+
[sort $ u ++ v |
u <- findInvSq (1%2 - x) all23,
+
(x, s) <- only23,
v <- s]
+
u <- findInvSq (1%2 - x) all23,
else undefined
+
v <- s
+
]
  +
 
problem_152 = length solutions
 
problem_152 = length solutions
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=153 Problem 153] ==
+
== [http://projecteuler.net/index.php?section=problems&id=153 Problem 153] ==
 
Investigating Gaussian Integers
 
Investigating Gaussian Integers
   
Solution:
+
== [http://projecteuler.net/index.php?section=problems&id=154 Problem 154] ==
<haskell>
 
problem_153 = undefined
 
</haskell>
 
 
== [http://projecteuler.net/index.php?section=view&id=154 Problem 154] ==
 
 
Exploring Pascal's pyramid.
 
Exploring Pascal's pyramid.
   
Solution:
+
{{sect-stub}}
<haskell>
 
problem_154 = undefined
 
</haskell>
 
   
== [http://projecteuler.net/index.php?section=view&id=155 Problem 155] ==
+
== [http://projecteuler.net/index.php?section=problems&id=155 Problem 155] ==
 
Counting Capacitor Circuits.
 
Counting Capacitor Circuits.
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_155 = undefined
+
--http://www.research.att.com/~njas/sequences/A051389
  +
a051389=
  +
[1, 2, 4, 8, 20, 42,
  +
102, 250, 610, 1486,
  +
3710, 9228, 23050, 57718,
  +
145288, 365820, 922194, 2327914
  +
]
  +
problem_155 = sum a051389
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=156 Problem 156] ==
+
== [http://projecteuler.net/index.php?section=problems&id=156 Problem 156] ==
 
Counting Digits
 
Counting Digits
  +
  +
== [http://projecteuler.net/index.php?section=problems&id=157 Problem 157] ==
  +
Solving the diophantine equation 1/a+1/b= p/10n
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_156 = undefined
+
-- Call (a,b,p) a primitive tuple of equation 1/a+1/b=p/10^n
</haskell>
+
-- a and b are divisors of 10^n, gcd a b == 1, a <= b and a*b <= 10^n
  +
-- I noticed that the number of variants with a primitive tuple
  +
-- is equal to the number of divisors of p.
  +
-- So I produced all possible primitive tuples per 10^n and
  +
-- summed all the number of divisors of every p
   
== [http://projecteuler.net/index.php?section=view&id=157 Problem 157] ==
+
import Data.List
Solving the diophantine equation 1/a+1/b= p/10n
+
k `divides` n = n `mod` k == 0
   
Solution:
+
divisors n
<haskell>
+
| n == 10 = [1,2,5,10]
problem_157 = undefined
+
| otherwise =
  +
[ d |
  +
d <- [1..n `div` 5],
  +
d `divides` n ]
  +
++ [n `div` 4, n `div` 2,n]
  +
fp n =
  +
[ n*(a+b) `div` ab |
  +
a <- ds,
  +
b <- dropWhile (<a) ds,
  +
gcd a b == 1,
  +
let ab = a*b,
  +
ab <= n
  +
]
  +
where
  +
ds = divisors n
  +
numDivisors :: Integer -> Integer
  +
numDivisors n = product [ toInteger (a+1) | (p,a) <- primePowerFactors n]
  +
numVgln = sum . map numDivisors . fp
  +
  +
main = do
  +
print . sum . map numVgln . takeWhile (<=10^9) . iterate (10*) $ 10
  +
primePowerFactors x = [(head a ,length a)|a<-group$primeFactors x]
  +
merge xs@(x:xt) ys@(y:yt) = case compare x y of
  +
LT -> x : (merge xt ys)
  +
EQ -> x : (merge xt yt)
  +
GT -> y : (merge xs yt)
  +
  +
diff xs@(x:xt) ys@(y:yt) = case compare x y of
  +
LT -> x : (diff xt ys)
  +
EQ -> diff xt yt
  +
GT -> diff xs yt
  +
  +
primes, nonprimes :: [Integer]
  +
primes = [2,3,5] ++ (diff [7,9..] nonprimes)
  +
nonprimes = foldr1 f . map g $ tail primes
  +
where f (x:xt) ys = x : (merge xt ys)
  +
g p = [ n*p | n <- [p,p+2..]]
  +
primeFactors n =
  +
factor n primes
  +
where
  +
factor n (p:ps)
  +
| p*p > n = [n]
  +
| n `mod` p == 0 = p : factor (n `div` p) (p:ps)
  +
| otherwise = factor n ps
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=158 Problem 158] ==
+
== [http://projecteuler.net/index.php?section=problems&id=158 Problem 158] ==
 
Exploring strings for which only one character comes lexicographically after its neighbour to the left.
 
Exploring strings for which only one character comes lexicographically after its neighbour to the left.
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_158 = undefined
+
factorial n = product [1..toInteger n]
  +
fallingFactorial x n = product [x - i | i <- [0..fromIntegral n - 1] ]
  +
choose n k = fallingFactorial n k `div` factorial k
  +
fun n=(2 ^ n - n - 1) * choose 26 n
  +
problem_158=maximum$map fun [1..26]
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=159 Problem 159] ==
+
== [http://projecteuler.net/index.php?section=problems&id=159 Problem 159] ==
 
Digital root sums of factorisations.
 
Digital root sums of factorisations.
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_159 = undefined
+
import Control.Monad
  +
import Data.Array.ST
  +
import qualified Data.Array.Unboxed as U
  +
spfArray :: U.UArray Int Int
  +
spfArray = runSTUArray (do
  +
arr <- newArray (2,m-1) 0
  +
forM_ [2 .. m-1] $ \n ->
  +
writeArray arr n (n-9*((n-1) `div` 9))
  +
forM_ [2 .. m-1] $ \x ->
  +
forM_ [2 .. m`div`n-1] $ \n ->
  +
incArray arr x n
  +
return arr
  +
)
  +
where
  +
m=10^6
  +
incArray arr x n = do
  +
a <- readArray arr x
  +
b <- readArray arr n
  +
ab <- readArray arr (x*n)
  +
when(ab<a+b) (writeArray arr (x*n) (a + b))
  +
writ x=appendFile "p159.log"$ show x ++ "\n"
  +
main=mapM_ writ $U.elems spfArray
  +
problem_159 = main
  +
  +
--at first ,make main to get file "p159.log"
  +
--then ,add all num in the file
 
</haskell>
 
</haskell>
   
== [http://projecteuler.net/index.php?section=view&id=160 Problem 160] ==
+
== [http://projecteuler.net/index.php?section=problems&id=160 Problem 160] ==
 
Factorial trailing digits
 
Factorial trailing digits
  +
  +
We use the following two facts:
  +
  +
Fact 1: <hask>(2^(d + 4*5^(d-1)) - 2^d) `mod` 10^d == 0</hask>
  +
  +
Fact 2: <hask>product [n | n <- [0..10^d], gcd n 10 == 1] `mod` 10^d == 1</hask>
  +
  +
We really only need these two facts for the special case of
  +
<hask>d == 5</hask>, and we can verify that directly by
  +
evaluating the above two Haskell expressions.
  +
  +
More generally:
  +
  +
Fact 1 follows from the fact that the group of invertible elements
  +
of the ring of integers modulo <hask>5^d</hask> has
  +
<hask>4*5^(d-1)</hask> elements.
  +
  +
Fact 2 follows from the fact that the group of invertible elements
  +
of the ring of integers modulo <hask>10^d</hask> is isomorphic to the product
  +
of a cyclic group of order 2 and another cyclic group.
   
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_160 = undefined
+
problem_160 = trailingFactorialDigits 5 (10^12)
  +
  +
trailingFactorialDigits d n = twos `times` odds
  +
where
  +
base = 10 ^ d
  +
x `times` y = (x * y) `mod` base
  +
multiply = foldl' times 1
  +
x `toPower` k = multiply $ genericReplicate n x
  +
e = facFactors 2 n - facFactors 5 n
  +
twos
  +
| e <= d = 2 `toPower` e
  +
| otherwise = 2 `toPower` (d + (e - d) `mod` (4 * 5 ^ (d - 1)))
  +
odds = multiply [odd | a <- takeWhile (<= n) $ iterate (* 2) 1,
  +
b <- takeWhile (<= n) $ iterate (* 5) a,
  +
odd <- [3, 5 .. n `div` b `mod` base],
  +
odd `mod` 5 /= 0]
  +
  +
-- The number of factors of the prime p in n!
  +
facFactors p = sum . zipWith (*) (iterate (\x -> p * x + 1) 1) .
  +
tail . radix p
  +
  +
-- The digits of n in base b representation
  +
radix p = map snd . takeWhile (/= (0, 0)) .
  +
iterate ((`divMod` p) . fst) . (`divMod` p)
 
</haskell>
 
</haskell>
  +
it have another fast way to do this .
   
[[Category:Tutorials]]
+
Solution:
[[Category:Code]]
+
<haskell>
  +
import Data.List
  +
mulMod :: Integral a => a -> a -> a -> a
  +
mulMod a b c= (b * c) `rem` a
  +
squareMod :: Integral a => a -> a -> a
  +
squareMod a b = (b * b) `rem` a
  +
pow' :: (Num a, Integral b) => (a -> a -> a) -> (a -> a) -> a -> b -> a
  +
pow' _ _ _ 0 = 1
  +
pow' mul sq x' n' = f x' n' 1
  +
where
  +
f x n y
  +
| n == 1 = x `mul` y
  +
| r == 0 = f x2 q y
  +
| otherwise = f x2 q (x `mul` y)
  +
where
  +
(q,r) = quotRem n 2
  +
x2 = sq x
  +
powMod :: Integral a => a -> a -> a -> a
  +
powMod m = pow' (mulMod m) (squareMod m)
  +
  +
productMod =foldl (mulMod (10^5)) 1
  +
hFacial 0=1
  +
hFacial a
  +
|gcd a 5==1=(a*hFacial(a-1)) `mod` (5^5)
  +
|otherwise=hFacial(a-1)
  +
fastFacial a= hFacial $a `mod` 6250
  +
numPrime x p=takeWhile(>0) [x `div` (p^a)|a<-[1..]]
  +
p160 x=mulMod t5 a b
  +
where
  +
t5=10^5
  +
lst=numPrime x 5
  +
a=powMod t5 1563 $c `mod` 2500
  +
b=productMod c6
  +
c=sum lst
  +
c6=map fastFacial $x:lst
  +
problem_160 = p160 (10^12)
  +
  +
</haskell>

Latest revision as of 08:22, 23 February 2010

Contents

[edit] 1 Problem 151

Paper sheets of standard sizes: an expected-value problem.

Solution:

problem_151 = fun (1,1,1,1)
 
fun (0,0,0,1) = 0
fun (0,0,1,0) = fun (0,0,0,1) + 1
fun (0,1,0,0) = fun (0,0,1,1) + 1
fun (1,0,0,0) = fun (0,1,1,1) + 1
fun (a,b,c,d) = 
    (pickA + pickB + pickC + pickD) / (a + b + c + d)
    where
    pickA | a > 0 = a * fun (a-1,b+1,c+1,d+1)
          | otherwise = 0
    pickB | b > 0 = b * fun (a,b-1,c+1,d+1)
          | otherwise = 0
    pickC | c > 0 = c * fun (a,b,c-1,d+1)
          | otherwise = 0
    pickD | d > 0 = d * fun (a,b,c,d-1)
          | otherwise = 0

[edit] 2 Problem 152

Writing 1/2 as a sum of inverse squares

Note that if p is an odd prime, the sum of inverse squares of all terms divisible by p must have reduced denominator not divisible by p.

Solution:

import Data.Ratio
import Data.List
import Data.Ord (comparing)
import Data.Function (on)
 
invSq n = 1 % (n * n)
sumInvSq = sum . map invSq
 
subsets (x:xs) = let s = subsets xs in s ++ map (x :) s
subsets _      = [[]]
 
primes = 2 : 3 : 7 : [p | p <- [11, 13..79],
                          all (\q -> p `mod` q /= 0) [3, 5, 7]]
 
-- All subsets whose sum of inverse squares,
-- when added to x, does not contain a factor of p
pfree s x p = [(y, t) | t <- subsets s, let y =  x + sumInvSq t,
                        denominator y `mod` p /= 0]
 
 
-- All pairs (x, s) where x is a rational number whose reduced
-- denominator is not divisible by any prime greater than 3;
-- and s is all sets of numbers up to 80 divisible
-- by a prime greater than 3, whose sum of inverse squares is x.
only23 = foldl fun [(0, [[]])] [13, 7, 5]
    where
    fun a p = 
        collect $ [(y, u ++ v) |
        (x, s) <- a,
        (y, v) <- pfree (terms p) x p,
        u <- s]
    terms p = 
        [n * p | 
        n <- [1..80`div`p],
        all (\q -> n `mod` q /= 0) $
        11 : takeWhile (>= p) [13, 7, 5]
        ]
    collect = 
        map (\z -> (fst $ head z, map snd z)) .
        groupBy fstEq . sortBy cmpFst
    fstEq = (==) `on` fst
    cmpFst = comparing fst
 
-- All subsets (of an ordered set) whose sum of inverse squares is x
findInvSq x y = 
    fun x $ zip3 y (map invSq y) (map sumInvSq $ init $ tails y)
    where
    fun 0 _        = [[]]
    fun x ((n, r, s):ns)
        | r > x     = fun x ns
        | s < x     = []
        | otherwise = map (n :) (fun (x - r) ns) ++ fun x ns
    fun _ _        = []
 
-- All numbers up to 80 that are divisible only by the primes
-- 2 and 3 and are not divisible by 32 or 27.
all23 = [n | a <- [0..4], b <- [0..2], let n = 2^a * 3^b, n <= 80]
 
solutions = 
    [sort $ u ++ v |
    (x, s) <- only23,
    u <- findInvSq (1%2 - x) all23,
    v <- s
    ]
 
problem_152 = length solutions

[edit] 3 Problem 153

Investigating Gaussian Integers

[edit] 4 Problem 154

Exploring Pascal's pyramid.

[edit] 5 Problem 155

Counting Capacitor Circuits.

Solution:

--http://www.research.att.com/~njas/sequences/A051389
a051389=
   [1, 2, 4, 8, 20, 42,
    102, 250, 610, 1486, 
    3710, 9228, 23050, 57718,
    145288, 365820, 922194, 2327914
   ]
problem_155 = sum a051389

[edit] 6 Problem 156

Counting Digits

[edit] 7 Problem 157

Solving the diophantine equation 1/a+1/b= p/10n

Solution:

-- Call (a,b,p) a primitive tuple of equation 1/a+1/b=p/10^n
-- a and b are divisors of 10^n, gcd a b == 1, a <= b and a*b <= 10^n
-- I noticed that the number of variants with a primitive tuple
-- is equal to the number of divisors of p.
-- So I produced all possible primitive tuples per 10^n and
-- summed all the number of divisors of every p
 
import Data.List
k `divides` n = n `mod` k == 0
 
divisors n 
    | n == 10 = [1,2,5,10]
    | otherwise = 
        [ d |
        d <- [1..n `div` 5],
        d `divides` n ] 
        ++ [n `div` 4, n `div` 2,n]
fp n = 
    [ n*(a+b) `div` ab |
    a <- ds,
    b <- dropWhile (<a) ds,
    gcd a b == 1,
    let ab = a*b,
    ab <= n 
    ]
    where
    ds = divisors n
numDivisors :: Integer -> Integer
numDivisors n = product [ toInteger (a+1) | (p,a) <- primePowerFactors n]
numVgln = sum . map numDivisors . fp
 
main = do
    print . sum . map numVgln . takeWhile (<=10^9) . iterate (10*) $ 10 
primePowerFactors x = [(head a ,length a)|a<-group$primeFactors x] 
merge xs@(x:xt) ys@(y:yt) = case compare x y of
    LT -> x : (merge xt ys)
    EQ -> x : (merge xt yt)
    GT -> y : (merge xs yt)
 
diff  xs@(x:xt) ys@(y:yt) = case compare x y of
    LT -> x : (diff xt ys)
    EQ -> diff xt yt
    GT -> diff xs yt
 
primes, nonprimes :: [Integer]
primes    = [2,3,5] ++ (diff [7,9..] nonprimes) 
nonprimes = foldr1 f . map g $ tail primes
    where f (x:xt) ys = x : (merge xt ys)
          g p = [ n*p | n <- [p,p+2..]]
primeFactors n =
    factor n primes
    where
    factor n (p:ps) 
        | p*p > n        = [n]
        | n `mod` p == 0 = p : factor (n `div` p) (p:ps)
        | otherwise      = factor n ps

[edit] 8 Problem 158

Exploring strings for which only one character comes lexicographically after its neighbour to the left.

Solution:

factorial n = product [1..toInteger n]
fallingFactorial x n = product [x - i | i <- [0..fromIntegral n - 1] ]
choose n k = fallingFactorial n k `div` factorial k
fun n=(2 ^ n - n - 1) * choose 26 n 
problem_158=maximum$map fun [1..26]

[edit] 9 Problem 159

Digital root sums of factorisations.

Solution:

import Control.Monad
import Data.Array.ST
import qualified Data.Array.Unboxed as U
spfArray :: U.UArray Int Int
spfArray  = runSTUArray (do
  arr <- newArray (2,m-1) 0 
  forM_ [2 .. m-1] $ \n ->
    writeArray arr n (n-9*((n-1) `div` 9))
  forM_ [2 .. m-1] $ \x ->
    forM_ [2 .. m`div`n-1] $ \n ->
      incArray arr x n
  return arr 
  )
  where
  m=10^6
  incArray arr x n = do
      a <- readArray arr x
      b <- readArray arr n
      ab <- readArray arr (x*n)
      when(ab<a+b) (writeArray arr (x*n) (a + b))
writ x=appendFile "p159.log"$ show x ++ "\n"
main=mapM_ writ $U.elems spfArray
problem_159 = main
 
--at first ,make main to get file "p159.log"
--then ,add all num in the file

[edit] 10 Problem 160

Factorial trailing digits

We use the following two facts:

Fact 1:
(2^(d + 4*5^(d-1)) - 2^d) `mod` 10^d == 0
Fact 2:
product [n | n <- [0..10^d], gcd n 10 == 1] `mod` 10^d == 1

We really only need these two facts for the special case of

d == 5
, and we can verify that directly by

evaluating the above two Haskell expressions.

More generally:

Fact 1 follows from the fact that the group of invertible elements

of the ring of integers modulo
5^d
has
4*5^(d-1)
elements.

Fact 2 follows from the fact that the group of invertible elements

of the ring of integers modulo
10^d
is isomorphic to the product

of a cyclic group of order 2 and another cyclic group.

Solution:

problem_160 = trailingFactorialDigits 5 (10^12)
 
trailingFactorialDigits d n = twos `times` odds
  where
    base = 10 ^ d
    x `times` y = (x * y) `mod` base
    multiply = foldl' times 1
    x `toPower` k = multiply $ genericReplicate n x
    e = facFactors 2 n - facFactors 5 n
    twos
     | e <= d    = 2 `toPower` e
     | otherwise = 2 `toPower` (d + (e - d) `mod` (4 * 5 ^ (d - 1)))
    odds = multiply [odd | a <- takeWhile (<= n) $ iterate (* 2) 1,
                           b <- takeWhile (<= n) $ iterate (* 5) a,
                           odd <- [3, 5 .. n `div` b `mod` base],
                           odd `mod` 5 /= 0]
 
-- The number of factors of the prime p in n!
facFactors p = sum . zipWith (*) (iterate (\x -> p * x + 1) 1) .
               tail . radix p
 
-- The digits of n in base b representation
radix p = map snd . takeWhile (/= (0, 0)) .
          iterate ((`divMod` p) . fst) . (`divMod` p)

it have another fast way to do this .

Solution:

import Data.List
mulMod :: Integral a => a -> a -> a -> a
mulMod a b c= (b * c) `rem` a
squareMod :: Integral a => a -> a -> a
squareMod a b = (b * b) `rem` a
pow' :: (Num a, Integral b) => (a -> a -> a) -> (a -> a) -> a -> b -> a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
    where
    f x n y
        | n == 1 = x `mul` y
        | r == 0 = f x2 q y
        | otherwise = f x2 q (x `mul` y)
        where
            (q,r) = quotRem n 2
            x2 = sq x
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)
 
productMod =foldl (mulMod (10^5)) 1
hFacial 0=1
hFacial a
    |gcd a 5==1=(a*hFacial(a-1)) `mod` (5^5)
    |otherwise=hFacial(a-1)
fastFacial a= hFacial $a `mod` 6250
numPrime x p=takeWhile(>0) [x `div` (p^a)|a<-[1..]]
p160 x=mulMod t5 a b
    where
    t5=10^5
    lst=numPrime x 5
    a=powMod t5 1563 $c `mod` 2500
    b=productMod  c6 
    c=sum lst
    c6=map fastFacial $x:lst
problem_160 = p160 (10^12)