Calculating multiples in Haskell (conversion from C)? [closed]

I'm not really into understanding your C, so I implemented a solution afresh using the algorithm discussed here. The N in the linked algorithm is the product of the primes up to a in your problem description.

So first we'll need a list of primes. There's a standardish trick for getting a list of primes that is at once very idiomatic and relatively efficient:

primes :: [Integer]
primes = 2:filter isPrime [3..]

-- Doesn't work right for n<2, but we never call it there, so who cares?
isPrime :: Integer -> Bool
isPrime n = go primes n where
    go (p:ps) n | p*p>n = True
                | otherwise = n `rem` p /= 0 && go ps n

Next up: we want a way to iterate over the positive square-free divisors of N. This can be achieved by iterating over the subsets of the primes less than a. There's a standard idiomatic way to get a powerset, namely:

-- import Control.Monad
-- powerSet :: [a] -> [[a]]
-- powerSet = filterM (const [False, True])

That would be a fine component to use, but since at the end of the day we only care about the product of each powerset element and the value of the Mobius function of that product, we would end up duplicating a lot of multiplications and counting problems. It's cheaper to compute those two things directly while producing the powerset. So:

-- Given the prime factorization of a square-free number, produce a list of
-- its divisors d together with mu(d).
divisorsWithMu :: Num a => [a] -> [(a, a)]
divisorsWithMu [] = [(1, 1)]
divisorsWithMu (p:ps) = rec ++ [(p*d, -mu) | (d, mu) <- rec] where
    rec = divisorsWithMu ps

With that in hand, we can just iterate and do a little arithmetic.

f :: Integer -> Integer -> Integer
f a b = b - sum
    [ mu * (b `div` d)
    | (d, mu) <- divisorsWithMu (takeWhile (<=a) primes)
    ]

And that's all the code. Crunched 137 lines of C down to 15 lines of Haskell -- not bad! Try it out in ghci:

> f 3 30
20

As an additional optimization, one could consider modifying divisorsWithMu to short-circuit when its divisor is bigger than b, as we know such terms will not contribute to the final sum. This makes a noticeable difference for large a, as without it there are exponentially many elements in the powerset. Here's how that modification looks:

-- Given an upper bound and the prime factorization of a square-free number,
-- produce a list of its divisors d that are no larger than the upper bound
-- together with mu(d).
divisorsWithMuUnder :: (Ord a, Num a) => a -> [a] -> [(a, a)]
divisorsWithMuUnder n [] = [(1, 1)]
divisorsWithMuUnder n (p:ps) = rec ++ [(p*d, -mu) | (d, mu) <- rec, p*d<=n]
    where rec = divisorsWithMuUnder n ps

f' :: Integer -> Integer -> Integer
f' a b = b - sum
    [ mu * (b `div` d)
    | (d, mu) <- divisorsWithMuUnder b (takeWhile (<=a) primes)
    ]

Not much more complicated; the only really interesting difference is that there's now a condition in the list comprehension. Here's an example of f' finishing quickly for inputs that would take infeasibly long with f:

> f' 100 100000
88169

With data-ordlist package mentioned by Daniel Wagner in the comments, it is just

f a b = length $ unionAll [ [p,p+p..b] | p <- takeWhile (<= a) primes]

That is all. Some timings, for non-compiled code run inside GHCi:

~> f 100 (10^5)
88169
(0.05 secs, 48855072 bytes)

~> f 131 (3*10^6)
2659571
(0.55 secs, 1493586480 bytes)

~> f 131 132
131
(0.00 secs, 0 bytes)

~> f 500 300000
274055
(0.11 secs, 192704760 bytes)

Compiling will surely make the memory consumption a non-issue, by converting the length to a counting loop.

You'll have to use recursion in place of loops.

In (most) procedural or object-orientated languages, you should hardly ever (never?) be using recursion. It is horribly inefficient, as a new stack frame must be created each time the recursive function is called.

However, in a functional language, like Haskell, the compiler is often able to optimize the recursion away into a loop, which makes it much faster then its procedural counterparts.

I've converted your sieve function into a set of recursive functions in C. I'll leave it to you to convert it into Haskell:

int main(void) {
    //...
    int root = sqrt(PRIME_RANGE);
    primes[primeCnt++] = 2;

    sieve(3, PRIME_RANGE, root);
    //...
}

void sieve(int i, int end, int root) {
    if(i > end) {
        return;
    }

    if(!primeFlag[i]) {
        primes[primeCnt++] = i;

        if(root >= i) {
            markMultiples(i * i, PRIME_RANGE, i);
        }
    }

    i += 2;
    sieve(i, end, root);
}

void markMultiples(int j, int end, int prime) {
    if(j > end) {
        return;
    }

    primeFlag[j] = 1;
    j += i << 1;

    markMultiples(j, end, prime);
}

The point of recursion is that the same function is called repeatedly, until a condition is met. The results of one recursive call are passed onto the next call, until the condition is met.

Also, why are you bit-fiddling instead of just multiplying or dividing by 2? Any half-decent compiler these days can convert most multiplications and divisions by 2 into a bit-shift.