Problem50 - Project Euler

30分プログラム、その324。Problem50 - Project Euler。

素数41は6つの連続する素数の和として表せる:

41 = 2 + 3 + 5 + 7 + 11 + 13.

100未満の素数を連続する素数の和で表したときにこれが最長になる.
同様に, 連続する素数の和で1000未満の素数を表したときに最長になるのは953で21項を持つ.
100万未満の素数を連続する素数の和で表したときに最長になるのはどの素数か?

とりあえず、Erlangで作ってみた。けど、ファンがものすごい音がしてる。
アルゴリズムは、

1. 条件を満すNの上限を決める
2. 素数列のうちN個足して、素数になるやつが存在するかチェックする
3. 存在しないなら、N-1についてチェックする

そもそも、100万までの素数の計算が終了しない・・・。

使い方

1> problem50:solve(100).
6
2> problem50:solve(1000).
21
3> problem50:solve(1000000).
?

ソースコード

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

%% Prime generator
sieve([X|Xs])->
    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    -> 
	    is_prime(Y) orelse satisfy(To,N,Xs)
    end;
satisfy(To,_,[]) ->
    false.

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

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

参考

• 過去の30分プログラム
• Problem 50 - PukiWiki