Problem 2

Project Euler

Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

Find the sum of all the even-valued terms in the sequence which do not exceed four million.

訳

フィボナッチ数列のそれぞれの新しい項は，それより直前の二項の和である．1，2からスタートすると，最初の10項は

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

となる．400,000を超えない偶数であるすべての項の和を求めよ．

これは数学的にこちょこちょして答えが出てくるようなものではなさそうです．
いろいろコメントいただいたり，書いたりで長くなりました．まだ増殖予定．

まずはHaskell

{-
e2.hs Project Euler Problem 2
-}
import System

fib::[Int]
fib = 0:1:zipWith (+) fib (tail fib)

{- いらないよこれ
islesseq::Int->Int->Bool
islesseq n m
    | n <= m     = True
    | otherwise = False
-}

upbfib::Int->[Int]
{- 書き直し
upbfib p = [ n | n <- takeWhile (\x -> islesseq x p) fib, n `mod`2 == 0]
-}
upbfib p = [ n | n <- takeWhile (<=p) fib, n `mod`2 == 0]

main = do getArgs >>= putStr.show.sum.list.read.head

フィボナッチ数列のfibは遅延評価を活用して一気に漸化式で作っておきます．これを単にfilterで計算させると，400,000以下という条件が成立するかしないかでどこまでも無限リスト（フィボナッチ数列）を追いかけていくので，計算が終わりません．単調増加だから最初に条件を満たさないものがあれば，あとは自動的に満たさないのですが，そこまではさすがのHaskell（というかghc?）も分かってくれません．

そこで，takeWhileで「最初の要素から条件を満たす最後の要素」までを取り出すことで，単調増加性を考えてすべての400,000以下の要素を取り出し，その和をとります．

c:\usr\local\euler>e2 4000000
e2 4000000
4613732

ちょっと数学

フィボナッチ数列は0,1,1,2,3,...ですが，フィボナッチ数列で偶数が連続することはありません．なぜなら，もし，$k$，$k+1$項がともに偶数だとすると，$k-1$項も偶数です．これを繰り返せば，最初から偶数でなければいけませんが，これは矛盾です．

したがって，フイボナッチ数列で偶数が現われるのは，その直前の二項がともに奇数の場合のみです．最初の数項をみることで，偶数項が現われるのは3の倍数の項のときのみであることが分かります．

つまり$\{F_n\}$をフィボナッチ数列とすると，偶数のフィボナッチ数列は$\{F_{3n}\}$と表せます．これより，
\begin{eqnarray}F_{3n}&=&F_{3n-1}+F_{3n-2}\\&=&F_{3n-3}+2F_{3n-2}\\&=&F_{3n-3}+2(F_{3n-3}-F_{3n-4})\\&=&3F_{3n-3}+2F_{3n-4}\\&=&3F_{3n-3}+F_{3n-4}+(F_{3n-5}+F_{3n-6})\\&=&3F_{3n-3}+(F_{3n-4}+F_{3n-5})+F_{3n-6}\\&=&4F_{3n-3}+F_{3n-6}\end{eqnarray}

となります．$F_{3n}$を$A_{n}$と置き換えて番号をずらすことにより，漸化式
$A_{n+2}=4A_{n+1}+A_{n}\quad(n=0,1,2,\ldots), A_0=0, A_1=2$
が得られます．この漸化式を解けば
A_n = \frac{1}{\sqrt{3}}\left\{(2+\sqrt{5})^n - (2-\sqrt{5})^n \right\}
です．もっともこの一般項を使って計算するのは勘弁です．

再びHaskell

上で考えたことをもとに書き直すとこんな感じです．

{-
e2-2.hs Project Euler 
-}
import System

evenfib::[Int]
evenfib = 0:2:zipWith (+) evenfib (map (\x -> 4*x) (tail evenfib))

{- いらないよこれ
isless::Int->Int->Bool
isless n m
    | n <= m     = True
    | otherwise = False
-}

upbevenfib::Int->[Int]
{- 書き直し
upbevenfib p = [ n | n <- takeWhile (\x -> isless x p) evenfib ]
-}
upbevenfib p = [ n | n <- takeWhile (<=p) evenfib ]

main = do getArgs >>= putStr.show.sum.upbevenfib.read.head

コード的には大差はないですが，漸化式部分が少し複雑になりました．計算結果は変わらないので省略します．

追記

islesseqはもともとありますよね(^^;; __さんご指摘感謝．

Perlで考える

「Perlでフィボナッチ」は前にいっぱい考えたからあっさりと処理．

#
#e2.pl Project Euler Problem 2
#
use strict;
use warnings;

{my %cache;
sub fib{
    my $n = shift;
    return $cache{$n} if defined $cache{$n};
    return 0 if $n==0;
    return 1 if $n==1;
    return $cache{$n} = fib($n-1)+fib($n-2);
    }
}

my $sum;
my $n=0;

while(1){
    my $i = fib($n);
    last if $i > $ARGV[0];
    $sum += $i if $i % 2 ==0 ;
    $n++;
}

print $sum;

計算結果は省略．下から順番にフィボナッチを計算させて，偶数だったら足すのと，$ARGV[0]を超えたらループから抜け出すというだけです．キャッシュさせてるので，計算は一瞬で終わります．もっとも自前でキャッシュしなくても，

use strict;
use warnings;

use Memoize;

memoize 'fib';

sub fib{
    my $n=shift;
    return 0 if $n==0;
    return 1 if $n==1;
    return fib($n-1)+fib($n-2);
}

my $sum;
my $n=0;

while(1){
    my $i = fib($n);
    last if $i >= $ARGV[0];
    $sum += $i if $i % 2 ==0 ;
    $n++;
}

print $sum;

とまあ，MemoizeモジュールでOKなんですが，HOP読者としては，自前で実装したいわけです．なんとなく自前キャッシュの方が速かった気がします（贔屓目？）．

「3の倍数」項を使って偶数だけを計算させるなら以下の通りです．Memoizeモジュールを使うことにします．

use strict;
use warnings;

use Memoize;

memoize 'evenfib';

sub evenfib{
    my $n=shift;
    return 0 if $n==0;
    return 2 if $n==1;
    return 4*evenfib($n-1)+evenfib($n-2);
}

my $sum;
my $n=0;

while(1){
    my $i = fib($n);
    last if $i >= $ARGV[0];
    $sum += $i;
    $n++;
}

print $sum;

さらに数学．

コメントで``別の_''さんからのご指摘による別の方法．「3の倍数」項の漸化式$A_{n+2}=4A_{n+1}+A_n (n=0,1,2,\ldots)$を$A_{n+2}-A_n=4A_{n+1}$とします．そしてこれを$n=0,1,2,\ldots,n$としてすべての和を計算してます．左辺は一個飛ばしの差分の和に，右辺が各項の和の4倍になるのです．
$\begin{array}{crcl}&A_2-A_0 &=& 4 A_1 \\&A_3-A_1 &=& 4 A_2 \\&A_4-A_2 &=& 4 A_3 \\&&\vdots&\\&A_{n-1}-A_{n-3} &=& 4 A_{n-2} \\&A_n-A_{n-2} &=& 4 A_{n-1} \\+)&A_{n+1}-A_{n-1} &=& 4 A_n \\\hline &A_{n+1}-A_n - A_0 -A_1 &=& 4\displaystyle\sum_{i=1}^n A_i\end{eqnarray}$
つまり，
$\displaystyle\sum_{i=1}^n A_i = \frac{1}{4}(A_{n+1}-A_n-2)$
が得られるわけです．

またHaskellに戻る

さて，ご指摘そのものをコードにしておきます．

{-
e2-3.hs Project Euler 
-}
import System

evenfib::[Int]
evenfib = 0:2:zipWith (+) evenfib (map (\x -> 4*x) (tail evenfib))

sumupb::Int->Int
sumupb p = (last xs + y - 2) `div` 4 where (xs,y:_) = span (<=p) evenfib

main = getArgs >>= putStr.show.sumupb.read.head

span (<=p) は p以下のevenfibの項のリストと，pより大きいevenfibの項のリストによるタプルを作って，yが大きい最初の項，last xsが「以下の最後の項」を作ります．足して2を引いて4で割れば求める値です．

またPerl

イテレータ型のPerlのコードも書いておきます．

#
#e2-2.pl Project Euler Problem 2
#
use strict;
use warnings;

my $sum;

sub earn_evenfib{
    my $upper = shift;
    my ($state,$prev,$preprev)=(0,2,0);
    return sub {
        $state++;
        return 0 if $state-1==0;
        return 2 if $state-1==1;
        my $term = 4*$prev + $preprev;
        ($prev,$preprev)=($term,$prev);
        $term < $upper ? return $term : return undef;
    };
}

my $it = earn_evenfib($ARGV[0]);

while (defined (my $i = $it->()) ){$sum += $i}
print $sum;

漸化式を使うために，項数と直前の二項だけを保持するようにして，下から順番に計算するので戻ることはできませんが，スピードは速いです．