Problem 193

Problem 193 - Project Euler
平方数で割り切れない2^50以下の数を数える．
メビウス関数と包除原理を利用．

import Control.Monad.ST (ST,runST)
import Control.Monad (when)
import Data.Array.ST (STUArray,newArray,readArray,writeArray)

p193 :: Integer -> ST s Integer
p193 n = do let n' = floor.sqrt.fromIntegral $ n
            m <- newArray (2,n') 2 :: ST s (STUArray s Integer Int)
                        
            let count !i !c
                    | i > n'    = return c
                    | otherwise = do
                  m_i <- readArray m i
                  when (m_i==2) $ sieve1 i i >> sieve2 (i*i) (i*i)
                  m_i <- readArray m i
                  count (i+1) $ div n (i*i)*toInteger m_i+c

                sieve1 !i !j = when (j <= n') $ do
                                 m_j <- readArray m j
                                 writeArray m j $ if m_j==2 then -1 else -m_j
                                 sieve1 i (j+i)
                                                  
                sieve2 !i !j = when (j <= n') $ do
                                 writeArray m j 0
                                 sieve2 i (j+i)
            count 2 n

main :: IO ()
main = print.runST $ p193 (2^50)

これも，関数型らしくないなぁ…
javaで実装すると，下のようになった．

import java.util.Arrays;
public class P193{
    public static void main(String[] args){
        long n = 1L<<50, c = n;
        int s = 1<<25;
        int[] m = new int[s+1];
        Arrays.fill(m,2);
        for(int i=2;i<=s;i++){
            long k = ((long)i)*i;
            if(m[i]==2){
                for(int j=i;j<=s;j+=i) m[j] = m[j]==2?-1:-m[j];
                for(long j=k;j<=s;j+=k) m[(int)j] = 0;
            }
            c += (n/k)*m[i];
        }
        System.out.println(c);
    }
}