Euler Solution 234

From ProgSoc Wiki

Revision as of 12:59, 17 May 2009 by Sanguinev (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Solutions for Problem 234

For an integer n:

  • the lower prime square root lps(n) is the largest prime <= sqrt(n)
  • the upper prime square root ups(n) is the smallest prime >= sqrt(n)
  • an integer k is semidivisble for n if k <= n and k is divisible by lps(n) or ups(n) but not both.

Find the sum of all semidivisible numbers in the range 4 to 999966663333.

Haskell by SanguineV

Runtime: 5.178 seconds

{- Generate primes as usual -}
primes :: [Integer]
primes = 2:3:primes'
  where
    1:p:candidates  = [6*k+r | k <- [0..], r <- [1,5]]
    primes'        = p : filter isPrime candidates
    isPrime n       = all (not . divides n) $ takeWhile (\p -> p*p <= n) primes'
    divides n p     = mod n p == 0

{- The disjoint union of two (sorted) lists -}
dunion :: [Integer] -> [Integer] -> [Integer]
dunion xxs@(x:xs) yys@(y:ys) =
  case compare x y of
    LT -> x : dunion xs yys
    EQ -> dunion xs ys
    GT -> y : dunion xxs ys

{- Given a divisor and an inital number "n", find the next integer divisble by
 - the divisor >= n -}
nextDiv :: Integer -> Integer -> Integer
nextDiv div n | mod n div == 0 = n
nextDiv div n = n + div - (mod n div)

{- Given two numbers and an initial point, find all numbers divisible by either
 - (but not both) of the starting numbers beyond the limit. -}
divs :: Integer -> Integer -> Integer -> [Integer]
divs m n init = dunion [m',m'+m..] [n',n'+n..]
  where
    m' = nextDiv m init
    n' = nextDiv n init

{- Given a limit, find the sum of all semidivisible numbers up to that limit.
 - The trick is to work with consecutive primes and remember the current "d" or
 - position. -}
semiDivs :: Integer -> Integer
semiDivs lim = inner 0 2 3 (drop 2 primes) (divs 2 3 5)
  where
    inner :: Integer -> Integer -> Integer -> [Integer] -> [Integer] -> Integer
    inner acc p1 p2 ps (d:ds) | d > lim = acc
    inner acc p1 p2 (p3:ps) (d:ds) | d >= p2 * p2 = inner acc p2 p3 ps (divs p2 p3 (d + 1))
    inner acc p1 p2 ps (d:ds) = inner (acc + d) p1 p2 ps ds

{- Find the semidivisible numbers up to 999966663333 -}
main = print (semiDivs 999966663333)
Personal tools