Problem 193

平方数で割り切れない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 
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] = ;
}
c += (n/k)*m[i];
}
System.out.println(c);
}
}

コンテンツ