Problem 179

http://projecteuler.net/index.php?section=problems&id=179
たいていの場合、素因数分解は重い。
エラトステネスのふるいとDPを組み合わせたようなアルゴリズム。
単純な繰り返しでかける。Cのコード。

#include <stdio.h>
#define LIM 10000001
#define SQR 3163

int main(){
  int *d,p[SQR],i,j,c,t,s;
  d=(int*) malloc (LIM*sizeof(int));
  for(i=1;i<LIM;d[i++]=0);
  for(i=2;i<SQR;p[i++]=1);
  for(i=2;i<SQR;i++){
    if(p[i]){//if  i is prime
      for(j=i*i;j<SQR;j+=i) p[j]=0;//j is not prime
      for(j=i*i;j<LIM;j+=i) d[j]=i;//j's prime factor (less than sqrt(j)) is i
    }
  }
  d[1]=1;
  for(i=2;i<LIM;i++){
    if(d[i]==0) d[i]=2;//i is prime
    else {
      for(j=d[i],t=i/j,c=2;t%j==0;c++,t/=j);//delete prime factor j from i
      d[i]=d[t]*c;
    }
  }
  s=0;
  for(i=2;i<LIM-1;i++) if(d[i]==d[i+1]) s++;
  printf("%d\n",s);
  free(d);
}

javaで

import java.util.Arrays;
public class P179_2{
    public static void main(String[] args){
        int lim = 10000001,sqr=3163;
        int[] d = new int[lim];
        boolean[] p = new boolean[lim];
        Arrays.fill(d,0);
        Arrays.fill(p,true);
        for(int i=2;i<sqr;i++){
            if(p[i]){
                for(int j=i*i;j<sqr;j+=i) p[j]=false;
                for(int j=i*i;j<lim;j+=i) d[j]=i;
            }
        }
        d[1]=1;
        for(int i=2;i<lim;i++){
            if(d[i]==0) d[i]=2;
            else {
                int t=i/d[i],c=2;
                for(int j=d[i];t%j==0;c++,t/=j);
                d[i]=d[t]*c;
            }
        }
        int s=0;
        for(int i=2;i<lim-1;i++) if(d[i]==d[i+1]) s++;
        System.out.println(s);
    }
}

同じアルゴリズムをHaskellで。

import Data.IORef (IORef,newIORef,readIORef,modifyIORef)
import Data.Array.IO (IOUArray,newArray,readArray,writeArray)
import Control.Monad (forM_,when)

lim = 10^7
sqr = floor.sqrt.fromIntegral $ lim

delPrime p (n,q) | mod q p == 0 = delPrime p (n+1,div q p)
                 | otherwise = (n,q)

main = do d <- newArray (1,lim) 0 :: IO (IOUArray Int Int)
          p <- newArray (1,sqr) True :: IO (IOUArray Int Bool)

          forM_ [2..sqr] $ \i -> do
            pi <- readArray p i
            when pi $ do
              forM_ [i*i,i*(i+1)..sqr] $ \j -> writeArray p j False
              forM_ [i*i,i*(i+1)..lim] $ \j -> writeArray d j i

          writeArray d 1 1
          forM_ [2..lim] $ \i -> do
            di <- readArray d i
            if di==0 then writeArray d i 2
               else do case delPrime di (1,i) of
                         (c,t) -> readArray d t >>= writeArray d i.(c*)

          s <- newIORef (0::Int)
          forM_ [2..lim-1] $ \i -> do
            a <- readArray d i
            b <- readArray d $ i+1
            when (a==b) $ modifyIORef s (+1)

          readIORef s >>= print

関数型らしくない。