【高速】フィボナッチ数列をビネの公式を使いGoで実装する

結論

O(log (n))となり

1000項目めを求めるには0.000199s

10000項目めは0.000405sで求められました。

完成品

https://github.com/kooooohe/fibonacci/blob/master/main.go

通常の再帰を使った簡単な解き方

func fib(num int) int {
	if num <= 1 {
		return num
	}
	return fib(num-2) + fib(num-1)
}

これだとO(2^n)となるためかなり遅い。

※メモ化については今回は割愛メモ化すれば O(n)となる

一般項を求める公式

証明は漸化式を使うのだが、今回の主題ではないので割愛。

興味ある方は: https://mathtrain.jp/fibonacci

これをGoで実装していく前に考えることがある。

`√5`は無理数であるため、プログラムでは表現できない。

丸め誤差が起きてしまい、100項目までいくと全然違う数値になってしまう。ただ今回は結果的に出てくる値が自然数になるため、それをうまく使うことにしよう。

いわゆる複素数的なものを作れば良い。

ということで今回の無理数√5を表現するための構造体からまずは作成する。

type fibNum struct {
	num               *big.Rat
	numMultipliedBySR *big.Rat
}

となるためこの表現には

a := big.NewRat(1, 2)
b := big.NewRat(1, 2)
fibNum{a, b}

これを使う。

同様にしてこのビネの公式をこの構造体を使い表すと

a := big.NewRat(1, 2)
b := big.NewRat(1, 2)
fibNum{a, b}

このようになる。

各関数の実装

乗算実装

上記の式変形を利用して、 fibNum構造体の乗算を実装する

a,c= fibNum.num

b,d = fibNum.numMultipliedBySR と置き換えてもらえれば良い。


//(a+b√5)(c+d√5) = (ac+5bd)+(ad+bc)√5
func fibMul(t1, t2 fibNum) (r fibNum) {
	//(ac+5bd)
	var r1, r2 *big.Rat
	tmp1 := new(big.Rat).Mul(t1.num, t2.num)
	tmp2 := new(big.Rat).Mul(
		new(big.Rat).Mul(t1.numMultipliedBySR, t2.numMultipliedBySR),
		big.NewRat(5, 1),
	)
	r1 = new(big.Rat).Add(tmp1, tmp2)

	//(ad+bc)√5
	tmp3 := new(big.Rat).Mul(t1.num, t2.numMultipliedBySR)
	tmp4 := new(big.Rat).Mul(t1.numMultipliedBySR, t2.num)
	r2 = new(big.Rat).Add(tmp3, tmp4)

	//(ac+5bd)+(ad+bc)√5
	//new(big.Rat).Add(r3, rr3)

	r.num = r1
	r.numMultipliedBySR = r2
	return
}

累乗実装

これは作成済みのfibMul関数を流用すれば簡単である。

func fibExp(t fibNum, n int) (r fibNum) {
	r := t
	for i := 1; i < n; i++ {
		r = fibMul(t, r)
	}
	return
}

簡単な実装はこうだが、このアルゴリズムでは

O(n)の時間がかかってしまう。

高速指数計算アルゴリズム

そこで高速指数計算アルゴリズムを実装していく


func fibExp(t fibNum, n int) fibNum {

	r := t

	// 1
	a := fibNum{num: big.NewRat(1, 1), numMultipliedBySR: big.NewRat(0, 999)}

	for n > 1 {
		if n%2 == 1 {
			a = fibMul(a, r)
		}
		r = fibMul(r, r)
		n = n / 2
	}

	return fibMul(r, a)
}

すでに計算したものの2乗が使えるものは使っていこうという作戦

例: 2¹⁶ = (((2²)²)²)² といった具合である。

もちろんこれのみでは表現できない場合もあるが、それはaとして計算をしておいて最後に掛け合わせている。

これにより O(log(n)) まで高速化

6行目は、初期値である1を表現している。

減算実装

func fibMin(t1, t2 fibNum) (r fibNum) {
	r1 := new(big.Rat).Add(
		t1.num,
		new(big.Rat).Mul(t2.num, big.NewRat(-1, 1)),
	)

	r2 := new(big.Rat).Add(
		t1.numMultipliedBySR,
		new(big.Rat).Mul(t2.numMultipliedBySR, big.NewRat(-1, 1)),
	)

	r.num = r1
	r.numMultipliedBySR = r2
	return
}

これは特に難しいことはしていないのだが、一点、

math/big pkgには減算関数はなかったため、-1倍したものを足す形で実装した。

除算実装

最後に√5で割る処理が入っているが、1/√5 = 1/5 * √5 なため。

これをfibMulに渡せば同等の処理をできるため、特に実装することはなし。

計測

1000項目目を求めるには0.000199s

10000項目目は0.000405s

とかなり高速に求められる。

全体のコード


package main

import (
	"fmt"
	"math/big"
)

type fibNum struct {
	num               *big.Rat
	numMultipliedBySR *big.Rat
}

func main() {
	//println(fib(10000))
	fmt.Println(calcFib(10000))
}

// (((1+√5)/2)^n - ((1-√5)/2)^n) / √5
func calcFib(n int) *big.Int {
	// (1+√5)/2
	a := big.NewRat(1, 2)
	b := big.NewRat(1, 2)
	x := fibNum{a, b}

	// (1-√5)/2
	c := big.NewRat(1, 2)
	d := big.NewRat(-1, 2)
	y := fibNum{c, d}

	// (1/√5)
	e := big.NewRat(0, 999)
	f := big.NewRat(1, 5)
	z := fibNum{e, f}

	tmp := fibMin(fibExp(x, n), fibExp(y, n))

	// *(1/√5)
	r := fibMul(tmp, z)
	return r.num.Num()

}

func fibExp(t fibNum, n int) fibNum {

	r := t

	// 1
	a := fibNum{num: big.NewRat(1, 1), numMultipliedBySR: big.NewRat(0, 999)}

	for n > 1 {
		if n%2 == 1 {
			a = fibMul(a, r)
		}
		r = fibMul(r, r)
		n = n / 2
	}

	return fibMul(r, a)
}

/*
func fibExp(t fibNum, n int) (r fibNum) {
	tmp := t
	for i := 1; i < n; i++ {
		tmp = fibMul(t, tmp)
	}
	r = tmp
	return
}
*/

func fibMin(t1, t2 fibNum) (r fibNum) {
	r1 := new(big.Rat).Add(
		t1.num,
		new(big.Rat).Mul(t2.num, big.NewRat(-1, 1)),
	)

	r2 := new(big.Rat).Add(
		t1.numMultipliedBySR,
		new(big.Rat).Mul(t2.numMultipliedBySR, big.NewRat(-1, 1)),
	)

	r.num = r1
	r.numMultipliedBySR = r2
	return
}

//(a+b√5)(c+d√5) = (ac+5bd)+(ad+bc)√5
func fibMul(t1, t2 fibNum) (r fibNum) {
	//(ac+5bd)
	var r1, r2 *big.Rat
	tmp1 := new(big.Rat).Mul(t1.num, t2.num)
	tmp2 := new(big.Rat).Mul(
		new(big.Rat).Mul(t1.numMultipliedBySR, t2.numMultipliedBySR),
		big.NewRat(5, 1),
	)
	r1 = new(big.Rat).Add(tmp1, tmp2)

	//(ad+bc)√5
	tmp3 := new(big.Rat).Mul(t1.num, t2.numMultipliedBySR)
	tmp4 := new(big.Rat).Mul(t1.numMultipliedBySR, t2.num)
	r2 = new(big.Rat).Add(tmp3, tmp4)

	//(ac+5bd)+(ad+bc)√5
	//rrr := new(big.Rat).Add(r3, rr3)

	r.num = r1
	r.numMultipliedBySR = r2
	return
}

func fib(num int) int {
	if num <= 1 {
		return num
	}
	return fib(num-2) + fib(num-1)
}