Problem 214

http://projecteuler.net/index.php?section=problems&id=214
前にHaskellで挑戦して、どうも遅くて解けなかったが、
Problem 179を解いて、ふるいを使えば効率よくΦ関数が計算できることに気がついた。
手続き型のほうが書きやすそうなのでjavaで。

public class P214{
    static int lim = 40000000;
    static int len = 25;
    public static void main(String[] args){
        int[] phi = new int[lim];
        int[] chain = new int[lim];
        long sum=0; chain[1]=1;
        for(int i=0;i<lim;phi[i]=i++);
        for(int i=2;i<lim;i++){
            if(phi[i]!=i){//if i is not prime
                chain[i]=chain[phi[i]]+1;
                continue;
            }
            phi[i]--; chain[i]=chain[i-1]+1;
            if(chain[i]==len) sum+=i;
            for(int j=2*i;j<lim;j+=i)phi[j]-=phi[j]/i;//sieve
        }
        System.out.println(sum);
    }
}

たくさんメモリを使います。chainはなくても良いが計算スピードのためにメモリを犠牲に。
Cでも書いてみた。

#include <stdio.h>
#define LIM 40000000
#define LEN 25
int main(){
  int *phi,*chain,i,j;
  long long sum;
  phi=(int*) malloc (LIM*sizeof(int));
  chain=(int*) malloc (LIM*sizeof(int));
  sum=0;chain[1]=1;
  for(i=0;i<LIM;i++)phi[i]=i;
  for(i=2;i<LIM;i++){
    if(phi[i]!=i){chain[i]=chain[phi[i]]+1;continue;}
    phi[i]--;chain[i]=chain[i-1]+1;
    if(chain[i]==LEN) sum+=i;
    for(j=2*i;j<LIM;j+=i)phi[j]-=phi[j]/i;
  }
  printf("%lld\n",sum);
}

実行時間は大して速くならなかった。
とりあえず、Haskellで書いてみる。

import Data.Array.IO (IOUArray,newArray_,readArray,writeArray)
import Data.IORef (newIORef,readIORef,modifyIORef)
import Control.Monad (forM_,when)
lim = 4*10^7;len = 25
main = do
  phi <- newArray_ (1,lim) :: IO (IOUArray Int Int)
  chain <- newArray_ (1,lim) :: IO (IOUArray Int Int)
  sum <- newIORef (0::Integer)
  writeArray chain 1 1
  forM_ [1..lim] $ \i -> writeArray phi i i
  forM_ [2..lim] $ \i -> do
         phi_i <- readArray phi i
         if phi_i /= i
            then readArray chain phi_i >>= writeArray chain i.succ
            else do writeArray phi i $ i-1
                    chain_i <- readArray chain $ i-1
                    writeArray chain i $ chain_i+1
                    when (chain_i+1==len).modifyIORef sum $ (+ toInteger i)
                    forM_ [2*i,3*i..lim] $ \j -> 
                        readArray phi j >>= writeArray phi j.((i-1)*).(`div` i)
  print =<< readIORef sum

これはもはや手続き型言語。