Personal tools

Euler problems/141 to 150

From HaskellWiki

< Euler problems(Difference between revisions)
Jump to: navigation, search
m
 
(17 intermediate revisions by 9 users not shown)
Line 1: Line 1:
== [http://projecteuler.net/index.php?section=view&id=141 Problem 141] ==
+
== [http://projecteuler.net/index.php?section=problems&id=141 Problem 141] ==
 
Investigating progressive numbers, n, which are also square.
 
Investigating progressive numbers, n, which are also square.
  
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_141 = undefined
+
import Data.List
 +
intSqrt :: Integral a => a -> a
 +
intSqrt n
 +
    | n < 0 = error "intSqrt: negative n"
 +
    | otherwise = f n
 +
    where
 +
        f x = if y < x then f y else x
 +
            where y = (x + (n `quot` x)) `quot` 2
 +
isSqrt n = n==((^2).intSqrt) n
 +
takec a b =
 +
    two++takeWhile (<=e12)
 +
    [sq| c1<-[1..], let c=c1*c1,let sq=(c^2*a^3*b+b^2*c) ]
 +
    where
 +
    e12=10^12
 +
    two=[sq|c<-[b,2*b],let sq=(c^2*a^3*b+b^2*c) ]
 +
problem_141=
 +
    sum$nub[c|
 +
    (a,b)<-takeWhile (\(a,b)->a^3*b+b^2<e12)
 +
        [(a,b)|
 +
        a<-[2..e4],
 +
        b<-[1..(a-1)]
 +
        ],
 +
    gcd a b==1,
 +
    c<-takec a b,
 +
    isSqrt c
 +
    ]
 +
    where
 +
    e4=120
 +
    e12=10^12
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=142 Problem 142] ==
+
== [http://projecteuler.net/index.php?section=problems&id=142 Problem 142] ==
 
Perfect Square Collection
 
Perfect Square Collection
  
Line 16: Line 44:
 
aToX (a,b,c)=[x,y,z]
 
aToX (a,b,c)=[x,y,z]
 
     where
 
     where
     x=div (a+b) 2
+
     x=(a+b)`div`2
     y=div (a-b) 2
+
     y=(a-b)`div`2
 
     z=c-x
 
     z=c-x
 
{-
 
{-
Line 57: Line 85:
 
     let n=(a2+b2)*(a2*b2+1),
 
     let n=(a2+b2)*(a2*b2+1),
 
     isSquare n,
 
     isSquare n,
     let t=div n 4,
+
     let t=n`div`4,
 
     let t2=a2*b2,
 
     let t2=a2*b2,
     let t3=div (a2*(b2+1)^2) 4
+
     let t3=(a2*(b2+1)^2)`div`4
 
     ]
 
     ]
  
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=143 Problem 143] ==
+
== [http://projecteuler.net/index.php?section=problems&id=143 Problem 143] ==
 
Investigating the Torricelli point of a triangle
 
Investigating the Torricelli point of a triangle
 +
 +
== [http://projecteuler.net/index.php?section=problems&id=144 Problem 144] ==
 +
Investigating multiple reflections of a laser beam.
  
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
import Data.List
+
type Point = (Double, Double)
import Data.Map((!),fromList,member)
+
type Vector = (Double, Double)
merge xs@(x:xt) ys@(y:yt) = case compare x y of
+
type Normal = (Double, Double)
    LT -> x : (merge xt ys)
+
    EQ -> x : (merge xt yt)
+
sub :: Vector -> Vector -> Vector
    GT -> y : (merge xs yt)
+
sub (x,y) (a,b) = (x-a, y-b)
   
+
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]
+
mull :: Double -> Vector -> Vector
primes    = [2,3,5] ++ (diff [7,9..] nonprimes)  
+
mull s (x,y) = (s*x, s*y)
nonprimes = foldr1 f . map g $ tail primes
+
    where f (x:xt) ys = x : (merge xt ys)
+
mulr :: Vector -> Double -> Vector
          g p = [ n*p | n <- [p,p+2..]]
+
mulr v s = mull s v
primeFactors :: Integer -> [Integer]
+
primeFactors n = factor n primes
+
dot :: Vector -> Vector -> Double
 +
dot (x,y) (a,b) = x*a + y*b
 +
 +
normSq :: Vector -> Double
 +
normSq v = dot v v
 +
 +
normalize :: Vector -> Vector
 +
normalize v
 +
    |len /= 0 =mulr v (1.0/len)
 +
    |otherwise=error "Vettore nullo.\n
 
     where
 
     where
        factor _ [] = []
+
    len = (sqrt . normSq) v
        factor m (p:ps) | p*p > m        = [m]
+
                        | m `mod` p == 0 = p : factor (m `div` p) (p:ps)
+
                        | otherwise      = factor m ps
+
 
   
 
   
fstfac x = [(head a ,length a)|a<-group$primeFactors x]
+
proj :: Vector -> Vector -> Vector
fac [(x,y)]=[x^a|a<-[0..y]]
+
proj a b = mull ((dot a b)/normSq b) b
fac (x:xs)=[a*b|a<-fac [x],b<-fac xs]
+
factors x=fac$fstfac x
+
intSqrt :: Integral a => a -> a
+
intSqrt n
+
    | n < 0 = error "intSqrt: negative n"
+
    | otherwise = f n
+
    where
+
        f x = if y < x then f y else x
+
            where y = (x + (n `quot` x)) `quot` 2
+
prim40=tail$take 40 primes
+
primeSqr=
+
    fromList[(a,fromList$zip b [1..])|
+
    a<-prim40,
+
    let b=nub[t|c<-[0..a-1],
+
    let t=mod (c*c) a]
+
    ]
+
  
isSqrt n
+
reflect :: Vector -> Normal -> Vector
    |k= n==((^2).intSqrt) n
+
reflect i n = sub i $ mulr (proj i n) 2.0
    |otherwise=False
+
    where
+
type Ray = (Point, Vector)
    k=foldl (&&) True [member k ma |
+
        a<-prim40,
+
makeRay :: Point -> Vector -> Ray
        let ma=(primeSqr !a),
+
makeRay p v = (p, v)
        let k=mod n a
+
   
        ]
+
getPoint :: Ray -> Double -> Point
getOne a = [c|
+
getPoint ((px,py),(vx,vy)) t = (px + t*vx, py + t*vy)
    x<-factors t,
+
    a>x,
+
type Ellipse = (Double, Double)
    let y=(a-x)*(3*a+x),
+
    let k=4*x,
+
getNormal :: Ellipse -> Point -> Normal
    let (c,m)=divMod y k,
+
getNormal (a,b) (x,y) = ((-b/a)*x, (-a/b)*y)
    m==0
+
    ]
+
    where
+
    t=(3*a^2)
+
  
getThree a = [[a,m,n]|
+
rayFromPoint :: Ellipse -> Vector -> Point -> Ray
    m<-t,
+
rayFromPoint e v p = makeRay p (reflect v (getNormal e p))
    n<-[k|k<-t,mod k 5/=0],
+
    let z=(2*m+n)^2+3*n*n,
+
test :: Point -> Bool
     isSqrt z
+
test (x,y) = y > 0 && x >= -0.01 && x <= 0.01
    ]
+
 +
intersect :: Ellipse -> Ray -> Point
 +
intersect (e@(a,b)) (r@((px,py),(vx,vy))) =
 +
     getPoint r t1
 
     where
 
     where
     t=getOne a
+
     c0 = normSq (vx/a, vy/b)
gcdlst [x,y]=gcd x y
+
     c1 = 2.0 * dot (vx/a, vy/b) (px/a, py/b)
gcdlst (x:xs)=gcd x$gcdlst xs
+
     c2 = (normSq (px/a, py/b)) - 1.0
      
+
     (t0, t1) = quadratic c0 c1 c2
p143 k=[c|
+
    a<-[1+k*groups..groups*(k+1)],
+
quadratic :: Double -> Double -> Double -> (Double, Double)
     c<-getThree (a*5),
+
quadratic a b c
     gcdlst c==1
+
     |d < 0= error "Discriminante minore di zero"
    ]
+
     |otherwise= if (t0 < t1) then (t0, t1) else (t1, t0)
-- run test find test==[],so one of a b c is 5*x
+
test=[(a,b,c)|a<-t,b<-t,c<-t,
+
     f a b,
+
     f b c,
+
    f c a]
+
 
     where
 
     where
     t=[1..4]
+
     d = b * b - 4.0 * a * c
     f a b=elem (mod (a^2+b^2+a*b) 5) [0,1,4]
+
     sqrtD = sqrt d
 
+
    q = if b < 0 then -0.5*(b - sqrtD) else 0.5*(b + sqrtD)
groups=200
+
    t0 = q / a
 +
    t1 = c / q
 
   
 
   
google num
+
calculate :: Ellipse -> Ray -> Int -> IO ()
-- write file to change bignum to small num
+
calculate e (r@(o,d)) n
  =if (num>33)
+
    |test p=print n
      then return()
+
    |otherwise=do
      else do let k=p143 num
+
        putStrLn $ "\rHit " ++ show n
              appendFile "file.log" $(show$k)  ++" "++(show num) ++"\n"
+
        calculate e (rayFromPoint e d p) (n+1)
              appendFile "files.log" $(show$map sum k) ++"  "++(show num) ++"\n"
+
    where
              google (num+1)
+
    p = intersect e r
-- first use main to make file.log
+
-- then run problem_143
+
main=google 0
+
split :: Char -> String -> [String]
+
split = unfoldr . split'
+
 
   
 
   
split' :: Char -> String -> Maybe (String, String)
+
origin = (0.0,10.1)
split' c l
+
direction = sub (1.4,-9.6) origin
  | null l = Nothing
+
ellipse = (5.0,10.0)
  | otherwise = Just (h, drop 1 t)
+
  where (h, t) = span (/=c) l
+
 
   
 
   
sToInt x=((++[-1]).read) $head$split ' ' x
+
problem_144 = do
+
     calculate ellipse (makeRay origin direction) 0
filer x
+
    |x<0=False
+
    |x>100000=False
+
    |otherwise=True
+
problem_143=do
+
     x<-readFile "files.log"
+
    let y=concat$map sToInt $lines x
+
    let z= filter filer y
+
    let t=[b|a<-z,b<-takeWhile (<=100000) [a*b|b<-[1..]]]
+
    print$ sum$nub t
+
 
+
</haskell>
+
 
+
== [http://projecteuler.net/index.php?section=view&id=144 Problem 144] ==
+
Investigating multiple reflections of a laser beam.
+
 
+
Solution:
+
<haskell>
+
problem_144 = undefined
+
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=145 Problem 145] ==
+
== [http://projecteuler.net/index.php?section=problems&id=145 Problem 145] ==
 
How many reversible numbers are there below one-billion?
 
How many reversible numbers are there below one-billion?
  
Line 228: Line 213:
 
     y=floor$logBase 10 $fromInteger x
 
     y=floor$logBase 10 $fromInteger x
 
     ten=10^y
 
     ten=10^y
     s=mod x 10
+
     s=x`mod`10
     h=div x ten
+
     h=x`div`ten
  
 
a2=[i|i<-[10..99],isOdig i]
 
a2=[i|i<-[10..99],isOdig i]
Line 254: Line 239:
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=146 Problem 146] ==
+
== [http://projecteuler.net/index.php?section=problems&id=146 Problem 146] ==
 
Investigating a Prime Pattern
 
Investigating a Prime Pattern
  
Line 260: Line 245:
 
<haskell>
 
<haskell>
 
import List
 
import List
find2km :: Integral a => a -> (a,a)
 
find2km n = f 0 n
 
    where
 
        f k m
 
            | r == 1 = (k,m)
 
            | otherwise = f (k+1) q
 
            where (q,r) = quotRem m 2       
 
 
millerRabinPrimality :: Integer -> Integer -> Bool
 
millerRabinPrimality n a
 
    | a <= 1 || a >= n-1 =
 
        error $ "millerRabinPrimality: a out of range ("
 
              ++ show a ++ " for "++ show n ++ ")"
 
    | n < 2 = False
 
    | even n = False
 
    | b0 == 1 || b0 == n' = True
 
    | otherwise = iter (tail b)
 
    where
 
        n' = n-1
 
        (k,m) = find2km n'
 
        b0 = powMod n a m
 
        b = take (fromIntegral k) $ iterate (squareMod n) b0
 
        iter [] = False
 
        iter (x:xs)
 
            | x == 1 = False
 
            | x == n' = True
 
            | otherwise = iter xs
 
 
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
 
 
mulMod :: Integral a => a -> a -> a -> a
 
mulMod a b c = (b * c) `mod` a
 
squareMod :: Integral a => a -> a -> a
 
squareMod a b = (b * b) `rem` a
 
powMod :: Integral a => a -> a -> a -> a
 
powMod m = pow' (mulMod m) (squareMod m)
 
 
isPrime x=millerRabinPrimality x 2
 
isPrime x=millerRabinPrimality x 2
--isPrime x=foldl  (&& )True [millerRabinPrimality x y|y<-[2,3,7,61,24251]]
+
--isPrime x=all (millerRabinPrimality x) [2,3,7,61,24251]
 
six=[1,3,7,9,13,27]
 
six=[1,3,7,9,13,27]
allPrime x=foldl (&&) True [isPrime k|a<-six,let k=x^2+a]
+
allPrime x=all (\a -> isPrime (x^2+a)) six
 
linkPrime [x]=filterPrime x
 
linkPrime [x]=filterPrime x
 
linkPrime (x:xs)=[y|
 
linkPrime (x:xs)=[y|
Line 315: Line 254:
 
     b<-[0..(x-1)],
 
     b<-[0..(x-1)],
 
     let y=b*prxs+a,
 
     let y=b*prxs+a,
     let c=mod y x,
+
     let c=y`mod`x,
 
     elem c d]
 
     elem c d]
 
     where
 
     where
Line 324: Line 263:
 
     [a|
 
     [a|
 
     a<-[0..(p-1)],
 
     a<-[0..(p-1)],
     length[b|b<-six,mod (a^2+b) p/=0]==6
+
     length[b|b<-six,(a^2+b)`mod`p/=0]==6
 
     ]
 
     ]
 
testPrimes=[2,3,5,7,11,13,17,23]
 
testPrimes=[2,3,5,7,11,13,17,23]
Line 334: Line 273:
 
     allPrime (y)
 
     allPrime (y)
 
     ]==1242490
 
     ]==1242490
p146 =[y|y<-linkPrime primes,y<150000000,allPrime (y)]
+
p146 =[y|y<-linkPrime primes,y<150000000,allPrime y]
 
problem_146=[a|a<-p146, allNext a]
 
problem_146=[a|a<-p146, allNext a]
 
allNext x=
 
allNext x=
Line 340: Line 279:
 
     where
 
     where
 
     a=[x^2+b|b<-six]
 
     a=[x^2+b|b<-six]
     b=head a:(map nextPrime a)
+
     b=head a:map nextPrime a
 
nextPrime x=head [a|a<-[(x+1)..],isPrime a]
 
nextPrime x=head [a|a<-[(x+1)..],isPrime a]
 
main=writeFile "p146.log" $show $sum problem_146
 
main=writeFile "p146.log" $show $sum problem_146
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=147 Problem 147] ==
+
== [http://projecteuler.net/index.php?section=problems&id=147 Problem 147] ==
 
Rectangles in cross-hatched grids
 
Rectangles in cross-hatched grids
 +
 +
== [http://projecteuler.net/index.php?section=problems&id=148 Problem 148] ==
 +
Exploring Pascal's triangle.
  
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
problem_147 = undefined
+
triangel 0 = 0
 +
triangel n
 +
    |n <7 =n+triangel (n-1) 
 +
    |n==k7 =28^k
 +
    |otherwise=(triangel i) + j*(triangel (n-i))
 +
    where
 +
    i=k7*((n-1)`div`k7)
 +
    j= -(n`div`(-k7))
 +
    k7=7^k
 +
    k=floor . logBase 7 . fromIntegral $ n
 +
problem_148=triangel (10^9)
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=148 Problem 148] ==
+
== [http://projecteuler.net/index.php?section=problems&id=149 Problem 149] ==
Exploring Pascal's triangle.
+
Searching for a maximum-sum subsequence.
  
 
Solution:
 
Solution:
 
<haskell>
 
<haskell>
import List
+
import Data.Array
digits n
+
import Data.List (foldl')
{- 123->[3,2,1]
+
   
-  -}
+
n = 2000
    |n<7=[n]
+
    |otherwise= y:digits x
+
res = maximum' $ concat [rows, cols, diags, diags']
    where
+
    (x,y)=divMod n 7
+
 
+
notDivX x=product$map (+1) $digits x
+
 
+
array::[Integer]
+
array=
+
    [a*b*c*d*e*f|
+
    let t=[1..7],
+
    a<-t,
+
    b<-t,
+
    c<-t,
+
    d<-t,
+
    e<-t,
+
    f<-t
+
    ]
+
 
+
fastNotDivX::Integer->Integer
+
fastNotDivX x=sum[k*a|a<-array]
+
 
     where
 
     where
    k=product$map (+1) $digits x
+
        rows  = map (maxSumInRow . getRow  laggedFibArray) [0 .. n-1]
 
+
        cols  = map (maxSumInRow . getCol  laggedFibArray) [0 .. n-1]
sumNotDivX x=sum[notDivX a|a<-[0..x]]
+
        diags  = map (maxSumInRow . getDiag  laggedFibArray) [-(n-2) .. (n-2)]
 
+
        diags' = map (maxSumInRow . getDiag' laggedFibArray) [-(n-2) .. (n-2)]
-- sum[fastNotDivX x|x<-[0..b]]=sumNotDivX ((b+1)*7^6-1)
+
 
   
 
   
moreNotDivX =sum[notDivX a|a<-[1000000000.. 1000016499 ]]
 
 
google num
 
-- write file to change bignum to small num
 
  =if (num>8499)
 
      then return()
 
      else do appendFile "file.log" $(show$fastNotDivX num)  ++"  "++(show num) ++"\n"
 
              google (num+1)
 
-- first use main to make file.log
 
-- then run problem_148
 
main=google 0
 
 
split :: Char -> String -> [String]
 
split = unfoldr . split'
 
 
   
 
   
split' :: Char -> String -> Maybe (String, String)
+
laggedFibArray :: Array Integer Integer
split' c l
+
laggedFibArray = listArray (0, n^2-1) $ map f [1..n^2]
   | null l = Nothing
+
    where
   | otherwise = Just (h, drop 1 t)
+
        f k = norm $ if k < 56
  where (h, t) = span (/=c) l
+
              then 100003 - (200003*k) + (300007*(k^3))
 +
              else (laggedFibArray ! (k-25)) + (laggedFibArray ! (k-56)) + (10^6)
 +
 +
        norm x = mod x (10^6) - 500000
 +
 +
 +
getRow   a i = map (a!) [i*n .. (i+1)*n-1]
 +
getCol   a i = map (a!) [i,n+i .. n*(n-1)+i]
 +
getDiag  a i = map (a!) $
 +
    if i >= 0
 +
    then [(i*n) + (k*(n+1)) | k <- [0..n-i-1]]
 +
    else [k + n*(k+i) | k <- [-i .. n-1]]
 +
getDiag' a i = map (a!) $
 +
    if i >= 0
 +
    then [(n*k) + n-k-i-1 | k <- [0..n-i-1]]
 +
    else [n*(k-i) + n-k-1 | k <- [0..n+i-1]]
 +
 +
 +
maxSumInRow = snd . foldl' f (0,0)
 +
    where
 +
        f (line_sum, line_max) x = (line_sum', max line_max line_sum')
 +
            where line_sum' = max (line_sum+x) 0
  
sToInt x=((+0).read) $head$split ' ' x
+
-- strict version of maximum
 +
maximum' (x:xs) = foldl' max x xs
  
problem_148=do
+
main = print res
    x<-readFile "file.log"
+
    let y=sum$map sToInt $lines x
+
    print (  y-(fromInteger moreNotDivX))
+
 
</haskell>
 
</haskell>
  
== [http://projecteuler.net/index.php?section=view&id=149 Problem 149] ==
+
== [http://projecteuler.net/index.php?section=problems&id=150 Problem 150] ==
Searching for a maximum-sum subsequence.
+
 
+
Solution:
+
<haskell>
+
problem_149 = undefined
+
</haskell>
+
 
+
== [http://projecteuler.net/index.php?section=view&id=150 Problem 150] ==
+
 
Searching a triangular array for a sub-triangle having minimum-sum.
 
Searching a triangular array for a sub-triangle having minimum-sum.
  
Solution:
+
{{sect-stub}}
<haskell>
+
problem_150 = undefined
+
</haskell>
+

Latest revision as of 10:51, 12 February 2010

Contents

[edit] 1 Problem 141

Investigating progressive numbers, n, which are also square.

Solution:

import Data.List
intSqrt :: Integral a => a -> a
intSqrt n
    | n < 0 = error "intSqrt: negative n"
    | otherwise = f n
    where
        f x = if y < x then f y else x
            where y = (x + (n `quot` x)) `quot` 2
isSqrt n = n==((^2).intSqrt) n
takec a b =
    two++takeWhile (<=e12) 
    [sq| c1<-[1..], let c=c1*c1,let sq=(c^2*a^3*b+b^2*c) ]
    where
    e12=10^12
    two=[sq|c<-[b,2*b],let sq=(c^2*a^3*b+b^2*c) ]
problem_141=
    sum$nub[c|
    (a,b)<-takeWhile (\(a,b)->a^3*b+b^2<e12) 
        [(a,b)|
        a<-[2..e4],
        b<-[1..(a-1)]
        ],
    gcd a b==1,
    c<-takec a b,
    isSqrt c
    ]
    where
    e4=120
    e12=10^12

[edit] 2 Problem 142

Perfect Square Collection

Solution:

import List
isSquare n = (round . sqrt $ fromIntegral n) ^ 2 == n
aToX (a,b,c)=[x,y,z]
    where
    x=(a+b)`div`2
    y=(a-b)`div`2
    z=c-x
{-
 -                                2    2    2
 -                               a  = c  + d
 -                                2    2    2
 -                               a  = e  + f
 -                                2    2    2
 -                               c  = e  + b
 -   let b=x*y  then 
 -                                             (y + xb)
 -                                          c= ---------
 -                                                 2
 -                                             (-y + xb)
 -                                          e= ---------
 -                                                 2
 -                                             (-x + yb)
 -                                          d= ---------
 -                                                 2
 -                                             (x + yb)
 -                                          f= ---------
 -                                                 2
 -
 - and 
 -                                2    2    2
 -                               a  = c  + d
 - then 
 -                                   2    2    2  2
 -                              2  (y  + x ) (x  y  + 1)
 -                             a = ---------------------
 -                                           4
 -
 -}
problem_142 = sum$head[aToX(t,t2 ,t3)|
    a<-[3,5..50],
    b<-[(a+2),(a+4)..50],
    let a2=a^2,
    let b2=b^2,
    let n=(a2+b2)*(a2*b2+1),
    isSquare n,
    let t=n`div`4,
    let t2=a2*b2,
    let t3=(a2*(b2+1)^2)`div`4
    ]

[edit] 3 Problem 143

Investigating the Torricelli point of a triangle

[edit] 4 Problem 144

Investigating multiple reflections of a laser beam.

Solution:

type Point = (Double, Double)
type Vector = (Double, Double)
type Normal = (Double, Double)
 
sub :: Vector -> Vector -> Vector
sub (x,y) (a,b) = (x-a, y-b)
 
mull :: Double -> Vector -> Vector
mull s (x,y) = (s*x, s*y)
 
mulr :: Vector -> Double -> Vector
mulr v s = mull s v
 
dot :: Vector -> Vector -> Double
dot (x,y) (a,b) = x*a + y*b
 
normSq :: Vector -> Double
normSq v = dot v v
 
normalize :: Vector -> Vector
normalize v 
    |len /= 0 =mulr v (1.0/len)
    |otherwise=error "Vettore nullo.\n"  
    where
    len = (sqrt . normSq) v 
 
proj :: Vector -> Vector -> Vector
proj a b = mull ((dot a b)/normSq b) b
 
reflect :: Vector -> Normal -> Vector
reflect i n = sub i $ mulr (proj i n) 2.0
 
type Ray = (Point, Vector)
 
makeRay :: Point -> Vector -> Ray
makeRay p v = (p, v)
 
getPoint :: Ray -> Double -> Point
getPoint ((px,py),(vx,vy)) t = (px + t*vx, py + t*vy)
 
type Ellipse = (Double, Double)
 
getNormal :: Ellipse -> Point -> Normal
getNormal (a,b) (x,y) = ((-b/a)*x, (-a/b)*y)
 
rayFromPoint :: Ellipse -> Vector -> Point -> Ray
rayFromPoint e v p = makeRay p (reflect v (getNormal e p))
 
test :: Point -> Bool
test (x,y) = y > 0 && x >= -0.01 && x <= 0.01
 
intersect :: Ellipse -> Ray -> Point
intersect (e@(a,b)) (r@((px,py),(vx,vy))) =
    getPoint r t1
    where
    c0 = normSq (vx/a, vy/b)
    c1 = 2.0 * dot (vx/a, vy/b) (px/a, py/b)
    c2 = (normSq (px/a, py/b)) - 1.0
    (t0, t1) = quadratic c0 c1 c2 
 
quadratic :: Double -> Double -> Double -> (Double, Double)
quadratic a b c  
    |d < 0= error "Discriminante minore di zero"
    |otherwise= if (t0 < t1) then (t0, t1) else (t1, t0)
    where
    d = b * b - 4.0 * a * c
    sqrtD = sqrt d
    q = if b < 0 then -0.5*(b - sqrtD) else 0.5*(b + sqrtD)
    t0 = q / a
    t1 = c / q 
 
calculate :: Ellipse -> Ray -> Int -> IO ()
calculate e (r@(o,d)) n 
    |test p=print n 
    |otherwise=do
         putStrLn $ "\rHit " ++ show n
         calculate e (rayFromPoint e d p) (n+1)
    where
    p = intersect e r 
 
origin = (0.0,10.1)
direction = sub (1.4,-9.6) origin
ellipse = (5.0,10.0)
 
problem_144 = do
    calculate ellipse (makeRay origin direction) 0

[edit] 5 Problem 145

How many reversible numbers are there below one-billion?

Solution:

import List
 
digits n 
{-  123->[3,2,1]
 -}
    |n<10=[n]
    |otherwise= y:digits x 
    where
    (x,y)=divMod n 10
-- 123 ->321
dmm=(\x y->x*10+y)
palind n=foldl dmm 0 (digits n) 
 
isOdd x=(length$takeWhile odd x)==(length x)
isOdig x=isOdd m && s<=h
    where
    k=x+palind x
    m=digits k
    y=floor$logBase 10 $fromInteger x
    ten=10^y
    s=x`mod`10
    h=x`div`ten
 
a2=[i|i<-[10..99],isOdig i]
aa2=[i|i<-[10..99],isOdig i,mod i 10/=0]
a3=[i|i<-[100..999],isOdig i]
m5=[i|i1<-[0..99],i2<-[0..99],
      let i3=i1*1000+3*100+i2,
      let i=10^6*   8+i3*10+5,
      isOdig i
   ]
 
fun i
    |i==2  =2*le aa2
    |even i=(fun 2)*d^(m-1)
    |i==3  =2*le a3
    |i==7  =fun 3*le m5
    |otherwise=0
    where
    le=length
    m=div i 2
    d=2*le a2
 
problem_145 = sum[fun a|a<-[1..9]]

[edit] 6 Problem 146

Investigating a Prime Pattern

Solution:

import List
isPrime x=millerRabinPrimality x 2
--isPrime x=all (millerRabinPrimality x) [2,3,7,61,24251]
six=[1,3,7,9,13,27]
allPrime x=all (\a -> isPrime (x^2+a)) six
linkPrime [x]=filterPrime x
linkPrime (x:xs)=[y|
    a<-linkPrime xs,
    b<-[0..(x-1)],
    let y=b*prxs+a,
    let c=y`mod`x,
    elem c d]
    where
    prxs=product xs
    d=filterPrime x
 
filterPrime p=
    [a|
    a<-[0..(p-1)],
    length[b|b<-six,(a^2+b)`mod`p/=0]==6
    ]
testPrimes=[2,3,5,7,11,13,17,23]
primes=[2,3,5,7,11,13,17,23,29]
test =
    sum[y|
    y<-linkPrime testPrimes,
    y<1000000,
    allPrime (y)
    ]==1242490
p146 =[y|y<-linkPrime primes,y<150000000,allPrime y]
problem_146=[a|a<-p146, allNext a]
allNext x=
    sum [1|(x,y)<-zip a b,x==y]==6
    where
    a=[x^2+b|b<-six]
    b=head a:map nextPrime a
nextPrime x=head [a|a<-[(x+1)..],isPrime a]
main=writeFile "p146.log" $show $sum problem_146

[edit] 7 Problem 147

Rectangles in cross-hatched grids

[edit] 8 Problem 148

Exploring Pascal's triangle.

Solution:

triangel 0 = 0
triangel n 
    |n <7 =n+triangel (n-1)  
    |n==k7 =28^k 
    |otherwise=(triangel i) + j*(triangel (n-i))
    where
    i=k7*((n-1)`div`k7)
    j= -(n`div`(-k7))
    k7=7^k
    k=floor . logBase 7 . fromIntegral $ n
problem_148=triangel (10^9)

[edit] 9 Problem 149

Searching for a maximum-sum subsequence.

Solution:

import Data.Array
import Data.List (foldl')
 
n = 2000
 
res = maximum' $ concat [rows, cols, diags, diags']
    where
        rows   = map (maxSumInRow . getRow   laggedFibArray) [0 .. n-1]
        cols   = map (maxSumInRow . getCol   laggedFibArray) [0 .. n-1]
        diags  = map (maxSumInRow . getDiag  laggedFibArray) [-(n-2) .. (n-2)]
        diags' = map (maxSumInRow . getDiag' laggedFibArray) [-(n-2) .. (n-2)]
 
 
laggedFibArray :: Array Integer Integer
laggedFibArray = listArray (0, n^2-1) $ map f [1..n^2]
    where
        f k = norm $ if k < 56
              then 100003 - (200003*k) + (300007*(k^3))
              else (laggedFibArray ! (k-25)) + (laggedFibArray ! (k-56)) + (10^6)
 
        norm x = mod x (10^6) - 500000
 
 
getRow   a i = map (a!) [i*n .. (i+1)*n-1]
getCol   a i = map (a!) [i,n+i .. n*(n-1)+i]
getDiag  a i = map (a!) $
    if i >= 0
    then [(i*n) + (k*(n+1)) | k <- [0..n-i-1]]
    else [k + n*(k+i) | k <- [-i .. n-1]]
getDiag' a i = map (a!) $
    if i >= 0
    then [(n*k) + n-k-i-1 | k <- [0..n-i-1]]
    else [n*(k-i) + n-k-1 | k <- [0..n+i-1]]
 
 
maxSumInRow = snd . foldl' f (0,0)
    where
        f (line_sum, line_max) x = (line_sum', max line_max line_sum')
            where line_sum' = max (line_sum+x) 0
 
-- strict version of maximum
maximum' (x:xs) = foldl' max x xs
 
main = print res

[edit] 10 Problem 150

Searching a triangular array for a sub-triangle having minimum-sum.