castle-bravo February 2016

Totient-function generating sieve consuming too much memory

I have written a sieve-based generator for the list of totients, and want to take the sum up to 1000000.

applyEvery :: (a -> a) -> Int -> [a] -> [a]
applyEvery f n xs = xf ++ (\(y:ys) -> f y : applyEvery f n ys) xb
  where
    (xf, xb) = splitAt (n - 1) xs

totients :: [Int]
totients = 1 : sieve [2..] [2..]
  where 
    sieve (x:xs) (y:ys) 
        | x == y    = (y - 1) : sieve xs (propagatePrime x ys)  
        | otherwise = y : sieve xs ys
    propagatePrime j = applyEvery (\x -> (quot x j)*(j - 1)) j

totientSum :: Int -> Int
totientSum n = sum $ take n totients

When computing totientSum n for n above 40000, ghci takes ages to evaluate and starts consuming tremendous amounts of memory. Compiling to an executable doesn't help. I assume that this has something to do with the way Haskell handles lazy evaluation.

I would like to know if it's possible to selectively apply strictness to improve memory consumption of the above functions so that I can compute the totient sum up to 1000000. I would also like to know if there's a better way to do this using lists, or if I should use a random-access data structure. If you know of a faster algorithm for computing the totient sum, please share a reference.

I thought that the definition of applyEvery might make a difference, so I tried the following other implementations, but they both seemed to consume more memory than the definition used above.

-- <https://www.reddit.com/2tpqip/>
applyEvery' :: (a -> a) -> Int -> [a] -> [a]
applyEvery' f n = zipWith ($) (cycle (replicate (n - 1) id ++ [f]))

applyEvery'' :: (a -> a) -> Int -> [a] -> [a]
applyEvery'' f n xs = xf ++ (\(y:ys) -> f y : applyEvery'' f n ys) xb
  where
    xf = take (n - 1) xs
    xb = drop (n - 1) xs

Answers


behzad.nouri February 2016

In implementing Euler product formula:

Euler product formula

you can take advantage of the fact that you are calculating Euler Phi numbers for all the numbers in the range [1..n]

Doing so, you may first find all the primes in the range [1..n] and then instead of finding prime divisors of each number, find all the multiples of each prime number. Obviously, the latter can be done much more efficiently.

A possible implementation would be:

import Data.Int (Int64)
import Control.Applicative ((<$>))
import Data.Array.Unboxed (UArray, elems, accum, listArray)

primes :: Integral a => [a]
primes = 2: 3: filter pred (chain [5,11..] [7,13..])
    where
    chain (x:xs) (y:ys) = x: y: chain xs ys
    pred a = all (\i -> a `mod` i /= 0) $ takeWhile (\i -> i * i <= a) primes

euler_phi :: Int64 -> [Int64]
euler_phi n = elems $ accum (\a p -> a - a `div` p) arr inc
    where
    val = takeWhile (<= n) primes
    idx = map (\i -> takeWhile (<= n) [i,2 * i..]) val

    inc = concat $ zipWith (\i j -> ($j) <$> (,) <$> i) idx val
    arr = listArray (1, n) [1..n] :: UArray Int64 Int64

main = getLine >>= print . sum . euler_phi . read

then:

\> euler_phi 20
[1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,8]

would be Euler totient function for the first 20 numbers; and if you compile with -O2 flag you may calculate the cumulative sums with a pretty decent performance:

$ ghc --make -O2 e 

Post Status

Asked in February 2016
Viewed 3,936 times
Voted 9
Answered 1 times

Search




Leave an answer