To high level refresh, the complexity of some $f(n)$ that describes an algorithm can be characterized using...

$O$ - asymptotic upper bound

$\Omega$ - asymptotic lower bound

$\Theta$ - asymptotic upper and lower bound

There are also...

$o$ - exclusive asymptotic upper bound

$\omega$ - exclusive asymptotic lower bound

The specific problem has to do with triangle numbers. All code is haskell.

"What is the value of the first triangle number to have over five hundred divisors?"

The obvious algorithm is to compute each triangle number by summing, and iterate over the integers from one to the number counting each divisor. The code for this is straightforward.

main = putStrLn $ show $ (fst . head) (filter ((>500) . snd) (zip t (map divs t))) where t = map tri [1..] tri n = sum [1..n] divs 0 = 0 divs n = length $ filter (\n' -> n `mod` n' == 0) [1..n]

It's simple, and it is correct. However, it's impractical.

Consider the complexity. The run-time to sum $1..n$ is $n$. The run-time to iterate over all possible divisors of triangle number $n$ is $\frac{n^2 + n}{2}$.

This gives a complexity of $\Theta(n^2)$. Since by definition the value being searched for has more than 500 divisors, the complexity to test each triangle number is impractical.

However, there is a practical solution. Furthermore, it only requires some optimization. Before eliminating the asymptotic bottleneck though, there is a simple strength reduction.

The $\Theta(n)$ algorithm to compute $\sum\limits_{i=1}^n i$ can easily be replaced with an $O(1)$algorithm.

tri n = (n^2 + n) `quot` 2

This equation sums $1..n$. It's often attributed to Carl Friedrich Gauss as an anecdote It reduces run-time to $\frac{n^2 + n}{2} + 1$ which is still $\Theta(n^2)$. So, again it doesn't eliminate the asymptotic bottleneck. Optimizing the algorithm to count divisors takes a bit more insight.

import Data.List (find, group) import Data.Maybe (fromMaybe) main = putStrLn $ show $ (fst . head) (filter ((>500) . snd) (zip t (map divs t))) where t = map tri [1..] facts :: Int -> [Int] facts x | (x < head pms) = [] | otherwise = h : facts (x `div` h) where h = fromMaybe x (find (\y -> ((x `mod` y) == 0)) (takeWhile (<= isqrt x) pms)) pms :: [Int] pms = sieve (2 : [3,5..]) sieve :: [Int] -> [Int] sieve (p:xs) = p : sieve [ x | x <- xs, x `mod` p > 0 ] tri :: Int -> Int tri n = (n^2 + n) `quot` 2 isqrt :: Int -> Int isqrt = truncate . sqrt . fromIntegral divs :: Int -> Int divs 1 = 1 divs x = product $ map ((+1) . length) ps where ps = group $ facts x

This algorithm is based on a simple consequence of the fundamental theorem of arithmetic. Instead of iterating over all integers below $\frac{n^2 + n}{2}$, it counts the combinations of prime factors. This replaces the $\frac{n^2 + n}{2}$ algorithm to find divisors with a $2\sqrt{\frac{n^2 + n}{2}}$ or $O(n)$ algorithm, since it only iterates over the primes below $\sqrt{\frac{n^2 + n}{2}}$

The realized run-time actually has a tighter bound, since prime numbers tend to grow more sparse further on the number line. However, knowing a bound of $O(n)$ is sufficient given the problem.

The execution time for the final version is...

real 0m0.037s

user 0m0.034s

sys 0m0.002s

There is still room for micro-optimization (starting at 14410396, better cpu utilization, etc...). However, these would likely only give constant factor speedups. More importantly though, optimizing the complexity turned the impractical algorithm, into an algorithm that finds a solution in less than a second.

## No comments:

## Post a Comment