Problem50 - Project Euler

Erlang 30分プログラム

30分プログラム、その325。Problem50 - Project Euler。昨日のやつ(id:mzp:20080620:euler)の改良版。ログメッセージを出すようにしたのと、ちゃんと答えを表示するようにした。

使い方

1> problem50:solve(1000000).
sive prime:2
sive prime:3
...
sive prime:999983
checking N=546
checking N=545
checking N=544
checking N=543
{543,
 997651}

ソースコード

-module(problem50).
-compile(export_all).

%% Prime generator
sieve([X|Xs])->
    io:format("sive prime:~p~n",[X]),
    Ys = lists:filter(fun (N)->N rem X =/= 0 end,Xs),
    [X|sieve(Ys)];
sieve([]) ->
    [].

make_primes(N)->
    sieve(lists:seq(2,N)).

is_prime(N,[X|Xs])->
    N rem X =/= 0 andalso is_prime(N,Xs);
is_prime(_,[]) ->
    true.
is_prime(N)->
    To = round(math:sqrt(N))+1,
    is_prime(N,lists:seq(2,To)).
    
%% estimate upper bound for N
up_bound(To,[X|Xs])->
    if 
	To-X < 0 -> 0;
	true     -> 1 + up_bound(To-X,Xs)
    end;
up_bound(To,[]) ->
    To.
    
%% check N is satified
take(0,_Xs,Ys) ->
    Ys;
take(N,[X|Xs],Ys) ->
    take(N-1,Xs,[X|Ys]);
take(_N,[],_Ys) ->
    [].

take(N,Xs)->
    take(N,Xs,[]).

satisfy(To,N,[_|Xs]=Ys)->
    Y = lists:sum(take(N,Ys)),
    if
	Y =:= 0;Y > To ->
	    false;
	true -> 
	    case is_prime(Y) of
		true -> 
		    Y;
		_ ->
		    satisfy(To,N,Xs)
	    end
    end;
satisfy(_To,_,[]) ->
    false.

%% find satisfy N
loop(To,N,Primes) when N > 0->
    io:format("checking N=~p~n",[N]),
    case satisfy(To,N,Primes) of
	false ->
	    loop(To,N-1,Primes);
	P ->
	    {N,P}
    end;
loop(_,_,_) ->
    unexpected.

%% run
solve(To)->
    Primes = make_primes(To),
    N = up_bound(To,Primes),
    loop(To,N,Primes).

参考