Euler problems/151 to 160
Jump to navigation
Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
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
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
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]
-- 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
-- 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 f [(0, [[]])] [13, 7, 5]
where
f 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 (x, _) (y, _) = x == y
cmpFst (x, _) (y, _) = compare x y
-- 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)
where
f 0 _ = [[]]
f x ((n, r, s):ns)
| r > x = f x ns
| s < x = []
| otherwise = map (n :) (f (x - r) ns) ++ f x ns
f _ _ = []
-- 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 = if verify
then [sort $ u ++ v | (x, s) <- only23,
u <- findInvSq (1%2 - x) all23,
v <- s]
else undefined
problem_152 = length solutions
Problem 153
Investigating Gaussian Integers
Solution:
#include <stdio.h>
#include <math.h>
typedef long long lolo;
static const lolo sumTo( lolo n ) { return n * ( n + 1 ) / 2; }
#define LL (1000)
lolo ssTab[ LL ];
int gcd(int a, int b) {
if (b==0) return a;
return gcd(b, a%b);
}
static const lolo sumSigma( lolo n ) {
lolo a, r, s;
if( n == 0 ) return 0;
if( n < LL ) { r = ssTab[ n ]; if( r ) return r; }
s = floor(sqrt( n ));
r = 0;
for( a = 1; a <= s; ++a ) r += a * ( n / a );
for( a = 1; a <= s; ++a ) r += ( sumTo( n / a ) - sumTo ( n / ( a + 1 ) ) ) * a;
if( n / s == s ) r -= s * s;
if( n < LL ) ssTab[ n ] = r;
return r;
}
int main() {
const lolo m = 100000000;
lolo t;
int a, b;
long ab;
t = sumSigma(m);
for( a = 1; a <=floor(sqrt(m)); ++a ) {
for( b = 1; b <= a && a * a + b * b <= m; ++b ) {
ab=(a*a+b*b);
if( ( a | b ) & 1 && gcd( a, b ) == 1 ) {
t += 2 * sumSigma( m / ab) * ( a == b ? a : a + b );
}
}
}
printf( "t = %lld\n", t );
return 1;
}
problem_153 = main
Problem 154
Exploring Pascal's pyramid.
Solution:
#include <stdio.h>
int main(){
int bound = 200000;
long long sum = 0;
int val2[bound+1], val5[bound+1]; // number of factors 2/5 in i!
int v2 = 0, v5 = 0;
int i;
int n;
for(n=0;n<=bound;n++)
{val5[n]=n/5+n/25+n/125+n/625+n/3125+n/15625+n/78125;
val2[n]=n/2+n/4+n/8+n/16+n/32+n/64+n/128+n/256+n/512+n/1024
+n/2048+n/4096+n/8192+n/16384+n/32768+n/65536+n/131072;}
v2 =val2[bound]- 11;
v5 = val5[bound]-11;
int j,k,vi2,vi5;
for(i = 2; i < 65625; i++){
if (!(i&1023)){
// look how many we got so far
printf("%d:\t%lld\n",i,sum);
}
vi5 = val5[i];
vi2 = val2[i];
int jb = ((bound - i) >> 1)+1;
// I want i <= j <= k
// by carry analysis, I know that if i < 4*5^5+2, then
// j must be at least 2*5^6+2
for(j = (i < 12502) ? 31252 : i; j < jb; j++){
k = bound - i - j;
if (vi5 + val5[j] + val5[k] < v5
&& vi2 + val2[j] + val2[k] < v2){
if (j == k || i == j){
sum += 3;
} else {
sum += 6;
}
}
}
}
printf("Total:\t%lld\n",sum);
return 0;
}
problem_154 = main
Problem 155
Counting Capacitor Circuits.
Solution:
problem_155 = undefined
Problem 156
Counting Digits
Solution:
digits =reverse.digits'
where
digits' n
|n<10=[n]
|otherwise= y:digits' x
where
(x,y)=divMod n 10
digitsToNum n=foldl dmm 0 n
where
dmm=(\x y->x*10+y)
countA :: Int -> Integer
countA 0 = 0
countA k = fromIntegral k * (10^(k-1))
countFun :: Integer -> Integer -> Integer
countFun _ 0 = 0
countFun d n = countL ds k
where
ds = digits n
k = length ds - 1
countL [a] _
| a < d = 0
| otherwise = 1
countL (a:tl) m
| a < d = a*countA m + countL tl (m-1)
| a == d = a*countA m + digitsToNum tl + 1 + countL tl (m-1)
| otherwise = a*countA m + 10^m + countL tl (m-1)
fixedPoints :: Integer -> [Integer]
fixedPoints d
= [a*10^10+b | a <- [0 .. d-1], b <- findFrom 0 (10^10-1)]
where
fun = countFun d
good r = r == fun r
findFrom lo hi
| hi < lo = []
| good lo = lgs ++ findFrom (last lgs + 2) hi
| good hi = findFrom lo (last hgs - 2) ++ reverse hgs
| h1 < l1 = []
| l1 == h1 = if good l1 then [l1] else []
| m0 == m1 = findFrom l1 (head mgs - 2) ++ mgs
++ findFrom (last mgs + 2) h1
| m0 < m1 = findFrom l1 (m0-1) ++ findFrom (goUp h1 m1) h1
| otherwise = findFrom l1 (goDown l1 m1) ++ findFrom (m0+1) h1
where
l1 = goUp hi lo
h1 = goDown l1 hi
goUp bd k
| k < k1 && k < bd = goUp bd k1
| otherwise = k
where
k1 = fun k
goDown bd k
| k1 < k && bd < k = goDown bd k1
| otherwise = k
where
k1 = fun k
m0 = (l1 + h1) `div` 2
m1 = fun m0
lgs = takeWhile good [lo .. hi]
hgs = takeWhile good [hi,hi-1 .. lo]
mgs = reverse (takeWhile good [m0,m0-1 .. l1])
++ takeWhile good [m0+1 .. h1]
problem_156=sum[sum $fixedPoints a|a<-[1..9]]
Problem 157
Solving the diophantine equation 1/a+1/b= p/10n
Solution:
n `splitWith` p = doSplitWith 0 n
where doSplitWith s t
| p `divides` t = doSplitWith (s+1) (t `div` p)
| otherwise = (s, t)
primesTo100 :: [Integer]
primesTo100 = [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]
trialDivision ps n = doTrialDivision ps
where doTrialDivision (p:ps) = let (q,r) = n `quotRem` p in if r == 0 then False else if q < p then True else doTrialDivision ps
doTrialDivision [] = True
primesTo10000 = primesTo100 ++ filter (trialDivision primesTo100) [101,103..9999]
d `divides` n = n `mod` d == 0
power (idG,multG) x n = doPower idG x n
where
doPower y _ 0 = y
doPower y x n =
let y' = if odd n then (y `multG` x) else y
x' = x `multG` x
n' = n `div` 2
in doPower y' x' n'
isStrongPseudoPrime :: Integer -> (Int,Integer) -> Integer -> Bool
isStrongPseudoPrime n (s,t) b =
let b' = power (1, \x y -> x*y `mod` n) b t
in if b' == 1 then True else doSquaring s b'
where
doSquaring 0 x = False
doSquaring s x
| x == n-1 = True
| x == 1 = False
| otherwise = doSquaring (s-1) (x*x `mod` n)
isMillerRabinPrime :: Integer -> Bool
isMillerRabinPrime n
| n < 100 = n `elem` primesTo100
| otherwise = all (isStrongPseudoPrime n (s,t)) primesTo100
where (s,t) = (n-1) `splitWith` 2
primePowerFactors :: Integer -> [(Integer,Int)]
primePowerFactors n | n > 0 = takeOutFactors n primesTo10000
where
takeOutFactors n (p:ps)
| p*p > n = finish n
| otherwise =
let (s,n') = n `splitWith` p
in if s > 0 then (p,s) : takeOutFactors n' ps else takeOutFactors n ps
takeOutFactors n [] = finish n
finish 1 = []
finish n =
if n < 100000000 || isMillerRabinPrime n
then [(n,1)]
else error ("primePowerFactors: unable to factor " ++ show n)
primeFactors :: Integer -> [Integer]
primeFactors n = concat (map (\(p,a) -> replicate a p) (primePowerFactors n))
fromPrimePowerFactors :: [(Integer,Int)] -> Integer
fromPrimePowerFactors factors = product [p^a | (p,a) <- factors]
numDivisors :: Integer -> Integer
numDivisors n = product [ toInteger (a+1) | (p,a) <- primePowerFactors n]
k `deelt` n = n `mod` k == 0
delers n | n == 10 = [1,2,5,10]
| otherwise = [ d | d <- [1..n `div` 5], d `deelt` 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 = delers n
ttLn = sum . map numDivisors . fp
problem_157 = do print . sum . map ttLn . takeWhile (<=10^9) . iterate (10*) $ 10
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 - fromInteger i | i <- [0..toInteger 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]
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 (0,m-1) 0
loop arr 2
forM_ [2 .. m - 1] $ \ x ->
loop2 arr x 2
return arr
)
where
m=10^6
loop arr n
|n>=m=return ()
|otherwise=do writeArray arr n (n-9*(div (n-1) 9))
loop arr (n+1)
loop2 arr x n
|n*x>=m=return ()
|otherwise=do incArray arr x n
loop2 arr x (n+1)
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"$foldl (++) "" [show x,"\n"]
main=do
mapM_ writ $U.elems spfArray
problem_159 = main
--at first ,make main to get file "p159.log"
--then ,add all num in the file
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=mod (a*hFacial(a-1)) (5^5)
|otherwise=hFacial(a-1)
fastFacial a= hFacial $mod a 6250
numPrime x p=takeWhile(>0) [div x (p^a)|a<-[1..]]
p160 x=mulMod t5 a b
where
t5=10^5
lst=numPrime x 5
a=powMod t5 1563 $mod c 2500
b=productMod c6
c=sum lst
c6=map fastFacial $x:lst
problem_160 = p160 (10^12)