Problem 187

http://projecteuler.net/index.php?section=problems&id=187
まぁ，関数型ではやりにくい．
とりあえず，javaで．
非常にシンプルな，ふるいを使ったアルゴリズム．

import java.util.Arrays;
public class P187{
    public static void main(String[] args){
        int n = 100000000,c=0;
        boolean[] prime = new boolean[n];
        Arrays.fill(prime,true);
        for(int i=2,m=(int)Math.sqrt(n);i<m;i++)
            if(prime[i]) for(int j=i*i;j<n;j+=i) prime[j]=false; // sieve prime
        
        for(int i=2;i<n;i++){
            if(prime[i]) continue;
            for(int j=2*i;j<n;j+=i)prime[j]=true; // sieve semiprime
            c++;
        }
        System.out.println(c);
    }
}

あまり，速くない．
2回目のふるいの効率が良くない．
実は，公式が存在した．

import java.util.Arrays;
public class P187_1{
    public static void main(String[] args){
        int n = 100000000, m=(int)Math.sqrt(n);
        boolean[] prime = new boolean[n/2+1];
        Arrays.fill(prime,true);
        int[] pi = new int [n/2+1];
        for(int i=2 ; i<n/2+1 ; i++)
            if(prime[i]){
                if(i<m) for(int j = i*i ; j < n/2+1 ; j+=i) prime[j] = false; // sieve prime
                pi[i] = pi[i-1]+1;
            }else pi[i] = pi[i-1];
        int count = 0;
        for(int i=2 ;i <=m ;i++)
            if(prime[i]) count+=pi[n/i]-pi[i]+1;
        System.out.println(count);
    }
}

同じアルゴリズムをhaskellで実装すると…

{-# LANGUAGE BangPatterns #-}
import Control.Monad.ST (ST,runST)
import Control.Monad (when)
import Data.Array.ST (STUArray,newArray,readArray,writeArray)

p187 n = do
  q <- newArray (2,div n 2) True :: ST s (STUArray s Int Bool)
  p <- newArray (1,div n 2) 0 :: ST s (STUArray s Int Int)
       
  let m = floor.sqrt.fromIntegral $ n

      countP !i
         | i > div n 2 = return ()
         | otherwise = do
           q_i <- readArray q i
           if q_i
              then do when (i<m) $ sieve i (i*i)
                      readArray p (i-1) >>= writeArray p i.(1+)
                      countP (i+1)
              else do readArray p (i-1) >>= writeArray p i
                      countP (i+1)
      
      sieve !i !j
         | j > div n 2 = return ()
         | otherwise = writeArray q j False >> sieve i (j+i)

      semiPrime !c !i
         | i > m = return c
         | otherwise = do
           q_i <- readArray q i
           if q_i
              then do p1 <- readArray p (div n i)
                      p2 <- readArray p i
                      semiPrime (c+p1-p2+1) (i+1)
              else semiPrime c (i+1)

  countP 2
  semiPrime 0 2

main = print.runST $ p187 (10^8)

どうなんでしょうね？
全然，関数型らしくない．