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=; chain[1]=1;
for(int i=;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=;chain[1]=1;
for(i=;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("%lldn",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 (::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

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