「課題: フィボナッチ数を出力せよ」の実装例

# これは何？

https://qiita.com/Nabetani/items/2441696a019317173f61

の実装例と軽く解説。

# 課題について

一見平易に見えるように作ったつもりだけど、どうだろう。

課題は以下のように分解できるんじゃないかな。

## フィボナッチ数を計算する

軽く考えれば定義通りに普通に計算するだけとなる。

C++ なんかだと 64bit でオーバーフローしたりするけど、 ruby なんかだとぼんやりしていても大丈夫。

出力例で f(800万以上) = 175万桁 のフィボナッチ数を計算しているので、これに追いつくのは定義通りの手順では厳しいと思う。

とはいえ、f(10万), f(20万), f(30万), ... と、10万間隔ぐらいで 2000万ぐらいまで事前に計算した値を持っておいたら、用意しておいた数値ぐらいまでは行けるかもね。それはそれでアリだと思う。試してないけど。

f(800万) で20万桁なので、ロードは遅いかもね。

f(800万) を大幅に超えていくには、行列の指数計算に落とし込んだりとか、[対数オーダーで計算するアルゴリズム](https://qiita.com/yassu/items/dab30eb2f2070c913451) を使ったりする必要がある。

## 書式通りに出す

まあこれは簡単だよね。 ommit とタイポしないように注意しましょう＞私。

## 一秒で終わる

一秒で処理を終えるのは色々な方法がある。

計算開始前に時計を作って1秒経ったら終える、という処理を書のは一案だけど、1秒経過したかどうかの判断をどこに入れるか迷うところである。

私の場合色々悩んだ上

```ruby:ruby
child=fork do
 # ここで計算
end
sleep(0.85)
Process.kill 15, child
```

としてみた。

手元の環境で実時間で 0.9〜1.0秒ぐらいで終わるようにしたら `sleep(0.85)` になった。

`fork` なので CRuby on Windows や jruby では実行できないけど、macOS でも linux でも OK だと思う。

## 数値の表現

ruby などで素朴に計算すると、100万桁とかもっと大きな整数をあつかうことになるけど、出力すべきは 最初の二桁と最後の二桁だけなので大半は無駄になる。
無駄な部分を端折ると計算量を大幅に減らせる。

桁数が n桁の場合。
普通に多倍長整数を使うと加算の計算量が $O(n)$ になるけど、不要な計算をサボることができれば $O(1)$ になる。 1000桁ぐらいだとそんなに遅くないかもだけど、 100万桁とかになるとてきめん。

### 誤差の蓄積の回避

下の方の桁の計算をサボると、計算があってるかどうか確信が持てなくなる。
なにせ、うまくいったら 10の無量大数乗 をはるかに超える値を扱うので、誤差の蓄積具合は一般市民の直感に従うことは期待できない。

上の二桁が真の値から外れていないことを確認しながら計算する必要がある。

# 実装例

## 素朴な実装

```ruby:ruby
b = (ARGV[0] || 2 ).to_i
c = (ARGV[1] || 1 ).to_i

child=fork do
  def omi(x)
    s = x.to_s
    return s if s.size<6
    s[0,2] + "(omit #{s.size-4} digits)" + s[-2,2]
  end

  r=[0,1]
  n=0
  (0...).each do |ix|
    if ix==b**n+c
      n+=1
      puts( "f(#{ix})=#{omi(r[0])}" )
    end
    r = [r[1], r[0]+r[1]]
  end
end
sleep(0.85)
Process.kill 15, child
```

課題を投稿したマシンとは別の PC だけど、実行するとこう。

```text:実行例
% time ruby simple_int.rb 34 25
f(26)=12(omit 2 digits)93
f(59)=95(omit 8 digits)41
f(1181)=29(omit 243 digits)81
f(39329)=84(omit 8215 digits)29
ruby simple_int.rb 34 25  0.04s user 0.03s system 8% cpu 0.937 total
```

まあこんなもん。

## 対数オーダーのアルゴリズムを入れる

前述の対数オーダーのアルゴリズムを入れるとこんな感じ。

```ruby:ruby
b = (ARGV[0] || 2 ).to_i
c = (ARGV[1] || 1 ).to_i
t0 = Time.now

child=fork do
  def omi(x)
    s = x.to_s
    return s if s.size<6
    s[0,2] + "(omit #{s.size-4} digits)" + s[-2,2]
  end
  
  class Fibo
    def initialize
      @c={}
    end
    def impl(n)
      # see https://qiita.com/yassu/items/dab30eb2f2070c913451
      h=(n/2).floor
      fc = calc(h)
      fp = calc(h-1)
      if n.even?
        fc  *(fc + 2 * fp)
      else
        fc * 2 * (fc + fp) + fp * fp
      end
    end
    def calc(n)
      return n if n<2
      @c[n] ||= impl(n)
    end
  end
  f = Fibo.new
  (0...).each do |x|
    i = b**x+c
    v = omi(f.calc(i))
    puts "f(#{i})=#{v}"
  end
end
sleep(0.85)
Process.kill 15, child
```

実行するとこんな感じ。

```text:実行例
% time ruby rec_int.rb 34 25
f(26)=12(omit 2 digits)93
f(59)=95(omit 8 digits)41
f(1181)=29(omit 243 digits)81
f(39329)=84(omit 8215 digits)29
f(1336361)=38(omit 279279 digits)61
ruby rec_int.rb 34 25  0.04s user 0.04s system 8% cpu 0.938 total
```

一件増えた。

## さらに一部の計算をサボる

```ruby:ruby
b = (ARGV[0] || 2 ).to_i
c = (ARGV[1] || 1 ).to_i

child=fork do
  require "./nnum.rb"

  def omi(x)
    return x.to_i.to_s if x.size<6
    x.head + "(omit #{x.size-4} digits)" + x.tail
  end
  
  class Fibo
    def initialize
      @c={}
    end
    TWO = Nnum.new(2)
    def impl(n)
      # see https://qiita.com/yassu/items/dab30eb2f2070c913451
      h=(n/2).floor
      fc = calc(h)
      fp = calc(h-1)
      if n.even?
        fc  *(fc + TWO * fp)
      else
        fc * TWO * (fc + fp) + fp * fp
      end
    end
    def calc(n)
      return Nnum.new(n) if n<2
      @c[n] ||= impl(n)
    end
  end
  f = Fibo.new
  (0...).each do |x|
    i = b**x+c
    v = omi(f.calc(i))
    puts "f(#{i})=#{v}"
  end
end
sleep(0.85)
Process.kill 15, child
```

整数の代わりに　`Nnum.new(1)` のようにして怪しい数値型にする。
`Nnum` は、 `nnum.rb` で定義されているんだけど、その正体は下記の通り。

```ruby:nnum.rb
class Nnum
  MANLEN = 500
  Val = Struct.new( :man, :exp ) do
    def int(base_exp, method)
      (man * 10**(exp-base_exp)).round
    end

    def head
      man.to_s[0,2]
    end

    def size
      return 1 if man==0
      MANLEN + exp + 1
    end

    def add(x, method)
      return x.dup if self.man==0
      return self.dup if x.man==0

      min_exp = [x.exp, self.exp].min
      delta = {floor:0, ceil:1}[method]
      return Val.create( man: self.man + delta, add_exp:self.exp, method:method ) if min_exp < self.exp-MANLEN*2
      return Val.create( man: x.man + delta, add_exp:x.exp, method:method ) if min_exp < x.exp-MANLEN*2

      m = self.int(min_exp, method) + x.int(min_exp, method)
      return Val.create( man: m, add_exp:min_exp, method:method )
    end

    def mul(x, method)
      return Val.new(0,0) if self.man==0 || x.man==0
      base_exp = x.exp + exp
      m = x.man * man
      return Val.create( man: m, add_exp:base_exp, method:method )
    end
  end

  def Val.create(man:, add_exp:, method:)
    case man
    when 0
      return Val.new( 0, 0 )
    else
      exp = Math.log10(man).floor - MANLEN
      m = (man.to_r / (10r**exp).to_r).send(method)
      return Val.new( m, exp+add_exp)
    end
  end

  def initialize(num, add_exp:0, lsd:nil)
    case num
    when 0
      @vals = Array.new(2){ Val.new(0,0) }
      @lsd = 0
    else
      @vals = [
        Val.create(man:num, add_exp:0, method: :floor),
        Val.create(man:num, add_exp:0, method: :ceil),
      ]
      @lsd = (lsd || num.abs) % 100
    end
  end

  attr_accessor :vals, :lsd

  def to_i
    i = [
      vals[0].int(0, :floor),
      vals[1].int(0, :ceil)
    ].uniq
    return i.first if i.size==1
    raise "no accurate value: #{@vals}"
  end

  def head
    h = @vals.map{ |e| e.head }.uniq
    return h.first if h.size==1
    raise "no accurate value: #{@vals}"
  end

  def tail
    return "%02d" % [@lsd % 100]
  end

  def size
    s = @vals.map{ |e| e.size }.uniq
    return s.first if s.size==1
    raise "no accurate value: #{@vals}"
  end
  
  def +(x)
    r=Nnum.new(0)
    r.vals = [
      vals[0].add( x.vals[0], :floor ),
      vals[1].add( x.vals[1], :ceil )
    ]
    r.lsd = (lsd+x.lsd) % 100
    r
  end

  def *(x)
    r=Nnum.new(0)
    r.vals = [
      vals[0].mul( x.vals[0], :floor ),
      vals[1].mul( x.vals[1], :ceil )
    ]
    r.lsd = (lsd * x.lsd) % 100
    r
  end

  def inspect
    "[#{@vals.map{ |e| e.inspect }.join(", ")}]"
  end
end
```

[@angel_p_57 さんの実装](https://ideone.com/TZQXvv) と違って、計算精度を動的には変えていないので、桁数が少なくて良い時はだいぶ遅い。

遅いとはいっても、全桁真面目に計算するのと比べるとだいぶ速く、実行するとこんな感じ。

```text:実行例
% time ruby rec_nnum.rb 34 25
f(26)=12(omit 2 digits)93
f(59)=95(omit 8 digits)41
f(1181)=29(omit 243 digits)81
f(39329)=84(omit 8215 digits)29
f(1336361)=38(omit 279279 digits)61
f(45435449)=83(omit 9495443 digits)49
f(1544804441)=26(omit 322845031 digits)41
f(52523350169)=31(omit 10976731006 digits)69
f(1785793904921)=76(omit 373208854158 digits)21
f(60716992766489)=16(omit 12689101041340 digits)89
f(2064377754059801)=26(omit 431429435405505 digits)01
f(70188843638032409)=20(omit 14668600803787122 digits)09
f(2386420683693101081)=39(omit 498732427328762096 digits)81
f(81138303245565435929)=34(omit 16956902529177911222 digits)29
f(2758702310349224820761)=23(omit 576534685992048981504 digits)61
ruby rec_nnum.rb 34 25  0.04s user 0.03s system 8% cpu 0.935 total
```

5垓桁（50万ペタ桁、5百エクサ桁）を全桁計算するのは手元のマシンでは絶対に不可能。
ましてや

```text:実行例
% time ruby rec_nnum.rb 2 0
f(1)=1
f(2)=1
︙
︙
f(2707685248164858261307045101702230179137145581421695874189921465443966120903931272499975005961073806735733604454495675614232576)=78(omit 565872750553651788747051202317916355188434282645760345693516708863794108339092679317303021965451138893221013004908066653804022 digits)07
ruby rec_nnum.rb 2 0  0.04s user 0.03s system 7% cpu 0.929 total
```

を全桁計算するなんて、この世のコンピュータ全て使っても無理と思う。

それが、うまくさぼれば手元のパソコンで 1秒足らずで計算できる。

# 感想

多倍長整数を扱ったことはなんどもあるけど、グーゴルプレックス ($10^{10^{100}}$) を超える整数の値をある程度正確に求めたいと思ったのはたぶん初めて。なんか新鮮な体験だった。