labunix's blog

labunixのラボUnix

フィボナッチ数列、リュカ数列、黄金比、おまけにletコマンドの最大値

■フィボナッチ関連のネタは過去にもある。

 「LLまつりでシェル芸人をジェネレートしてきた」を参加していないのに解いてみた。
 http://labunix.hateblo.jp/entry/20130915/1379249444

 bashとclispとawkで小学三年生レベルの等差数列の総和
 http://d.hatena.ne.jp/labunix/20130114

■フィボナッチ数列を30個。黄金比付き。
 初期値F{0}=0,F{1}=1,F{n+2}=F{n}+F{n+1},n>=0
 F{0}=a,F{1}=b,F{n+2}=cと置く。

$ a=0;b=1; \
  for n in `seq 0 30`;do \
    let c=$a+$b;echo "$a $b" | \
      awk '{if($1!=0)printf "%8f %8d\n",$2/$1,$1;else printf "%17d\n",$1}'; \
    a=$b;b=$c; \
  done
                0
1.000000        1
2.000000        1
1.500000        2
1.666667        3
1.600000        5
1.625000        8
1.615385       13
1.619048       21
1.617647       34
1.618182       55
1.617978       89
1.618056      144
1.618026      233
1.618037      377
1.618033      610
1.618034      987
1.618034     1597
1.618034     2584
1.618034     4181
1.618034     6765
1.618034    10946
1.618034    17711
1.618034    28657
1.618034    46368
1.618034    75025
1.618034   121393
1.618034   196418
1.618034   317811
1.618034   514229
1.618034   832040

■リュカ数列30個。黄金比付き。
 初期値L{0}=2,L{1}=1,L{n+2}=L{n}+L{n+1},n>=0
 L{0}=a,L{1}=b,L{n+2}=cと置く。

$ a=2;b=1; \
  for n in `seq 0 30`;do \
    let c=$a+$b;echo "$a $b" | awk '{printf "%8f %8d\n",$2/$1,$1}'; \
    a=$b;b=$c; \
  done
0.500000        2
3.000000        1
1.333333        3
1.750000        4
1.571429        7
1.636364       11
1.611111       18
1.620690       29
1.617021       47
1.618421       76
1.617886      123
1.618090      199
1.618012      322
1.618042      521
1.618031      843
1.618035     1364
1.618034     2207
1.618034     3571
1.618034     5778
1.618034     9349
1.618034    15127
1.618034    24476
1.618034    39603
1.618034    64079
1.618034   103682
1.618034   167761
1.618034   271443
1.618034   439204
1.618034   710647
1.618034  1149851
1.618034  1860498

■黄金比は無理数なので、フィボナッチ数からclispで桁を増やした黄金比を求めてみる。

$ a=1;b=1; \
  for n in `seq 0 30`;do \
    let c=$a+$b;echo "(setf (EXT:LONG-FLOAT-DIGITS) 240)(/ $b.0L0 $a)" | \
      clisp -q | tail -1 | sed s/"L0"//g; \
     a=$b;b=$c; \
  done
1.0
2.0
1.5
1.66666666666666666666666666666666666666666666666666666666666666666666666666666
1.6
1.625
1.61538461538461538461538461538461538461538461538461538461538461538461538461539
1.61904761904761904761904761904761904761904761904761904761904761904761904761905
1.61764705882352941176470588235294117647058823529411764705882352941176470588235
1.61818181818181818181818181818181818181818181818181818181818181818181818181819
1.6179775280898876404494382022471910112359550561797752808988764044943820224719
1.61805555555555555555555555555555555555555555555555555555555555555555555555555
1.61802575107296137339055793991416309012875536480686695278969957081545064377683
1.61803713527851458885941644562334217506631299734748010610079575596816976127321
1.6180327868852459016393442622950819672131147540983606557377049180327868852459
1.61803444782168186423505572441742654508611955420466058763931104356636271529889
1.61803381340012523481527864746399499060738885410144020037570444583594239198496
1.61803405572755417956656346749226006191950464396284829721362229102167182662539
1.61803396316670652953838794546759148529060033484812245874192776847644104281273
1.6180339985218033998521803399852180339985218033998521803399852180339985218034
1.61803398501735793897314087337840306961447103964918691759546866435227480358122
1.6180339901755970865563773925808819377787815481903901530122522725989498052058
1.61803398820532505147084481976480441078968489374323899919740377569180304986566
1.61803398895790200138026224982746721877156659765355417529330572808833678398895
1.61803398867044318560479840053315561479506831056314561812729090303232255914696
1.61803398878024268285650737686687041262675772079114940729696110978392493801125
1.61803398873830300685273243796393405899662963679499842173324237086214094431264
1.6180339887543225376088304054925726296446630229916522713184880321952355330684
1.61803398874820362134379819107829391185639082976650480622446419785737482716845
1.61803398875054083938272198451997500120186529493774337772222489303398875054084
1.618033988749648101530971893432887483853524072826455931169773648505610691474

■(5+1)/2で求めた黄金比

$ echo "(setf (EXT:LONG-FLOAT-DIGITS) 240)(/ (+ (sqrt 5.0L0) 1) 2)" | \
  clisp -q | tail -1 | sed s/"L0"//g
1.61803398874989484820458683436563811772030917980576286213544862270526046281891

■四捨五入でn番目のリュカ数を予測。
 まずは10番目。

$ n=10;echo "(- (expt 1.618033 $n) (expt -1.618033 -$n))" | clisp -q | tail -1 | \
  awk '{printf "%d\n",$1+0.5}'
123

$ a=2;b=1; \
  for n in `seq 0 10`;do \
    let c=$a+$b;echo "$a"; \
    a=$b;b=$c; \
  done | tail -1
123100番目。
 letがオーバーフローしてる。

$ n=100;echo "(- (expt 1.618033 $n) (expt -1.618033 -$n))" | clisp -q | tail -1 | awk '{printf "%d\n",$1+0.5}'
792025799999999967232

$ a=2;b=1; \
  for n in `seq 0 100`;do \
    let c=$a+$b;echo "$a"; \
    a=$b;b=$c; \
  done | tail -1
-1139155321138466361

■ちなみにletがオーバーフローするタイミングは2^63-1を超えたところ。
 いわゆる64bitの「INTMAX_MAX」の値です。

$ let n="$((0x7FFFFFFFFFFFFFFF))";while true;do echo $n;test "$n" -le 0 && break;let n=$n+1;done
9223372036854775807
-9223372036854775808

$ echo "(- (expt 2 63) 1)" | clisp -q | tail -1
9223372036854775807

■簡単に浮動小数点を扱うようclispに変更したが、
 今度は黄金比の有効桁数が不足している様子。

$ a=2;b=1; \
  for n in `seq 0 100`;do \
    c=`echo "(+ $a $b)" | clisp -q |tail -1`;echo "$a"; \
    a=$b;b=$c; \
  done | tail -1 
792070839848372253127

■桁数を増やしてみるというより、黄金比の計算式をそのまま入れてみる。
 うーむ。ちょっと残念な結果だけど、分からなくもない。。。

$ n=100;echo "(- (expt (/ (+ (sqrt 5.0L0) 1) 2) $n) \
              (expt (/ (+ (sqrt 5.0L0) 1) 2) -$n))" | \
          clisp -q | tail -1 | awk '{printf "%d\n",$1+0.5}'
792072340000000049152

$ n=100;echo "(setf (EXT:LONG-FLOAT-DIGITS) 240) \
              (- (expt (/ (+ (sqrt 5.0L0) 1) 2) $n) \
              (expt (/ (+ (sqrt 5.0L0) 1) 2) -$n))" | \
          clisp -q | tail -1 | sed s/"L20"//g 
7.9207083984837225312699999999999999999999747497332387231411562195458792078607

■n=1000にしてみる。
 なるほど。惜しいw。まあ、予測の範疇なので十分かと。。。
 「INTMAX_MAX」を超える整数値には注意が必要なようです。。。

$ a=2;b=1; \
  for n in `seq 0 1000`;do \
    c=`echo "(+ $a $b)" | clisp -q |tail -1`;echo "$a"; \
    a=$b;b=$c; \
  done | tail -1 | sed s/"1..80"/"&\n"/g
97194177735908175207981982079326473737797879155345685082728081084772518818444815269
08061914904596829767957830540320934740116303690766057397174086246375180164120149028
4097309096322681531675707666695323797578127

$ n=1000;  echo "(setf (EXT:LONG-FLOAT-DIGITS) 1200) \
        (- (expt (/ (+ (sqrt 5.0L0) 1) 2) $n) \
        (expt (/ (+ (sqrt 5.0L0) 1) 2) -$n))" |     clisp -q | tail -1 | sed s/"L20"//g 
9.7194177735908175207981982079326473737797879155345685082728081084772518818444815269
08061914904596829767957830540320934740116303690766057397174086246375180164120149028
40973090963226815316757076666953237975781270000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000005488