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++]=);
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]=;//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]==) d[i]=2;//i is prime
else {
for(j=d[i],t=i/j,c=2;t%j==;c++,t/=j);//delete prime factor j from i
d[i]=d[t]*c;
}
}
s=;
for(i=2;i<LIM-1;i++) if(d[i]==d[i+1]) s++;
printf("%dn",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,);
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]==) d[i]=2;
else {
int t=i/d[i],c=2;
for(int j=d[i];t%j==;c++,t/=j);
d[i]=d[t]*c;
}
}
int s=;
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 ==  = delPrime p (n+1,div q p)
| otherwise = (n,q)
main = do d <- newArray (1,lim)  :: 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== then writeArray d i 2
else do case delPrime di (1,i) of
(c,t) -> readArray d t >>= writeArray d i.(c*)
s <- newIORef (::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

関数型らしくない。