众所周知,斐波那契数列的定义是:f(0) = 0, f(1) = 1, 且当 n > 1 时 f(n) = f(n-1) + f(n-2)。今天讨论个相对简单的问题:斐波那契数列第一千万项怎么求?

暴力算法(V0)

最容易想到的是直接翻译递推式,二话不说撸出代码:

func fibV0(n int) Number {
	defer logElapsed("fibV0")()
	a, b := Num1, Num0
	for i := 0; i < n; i++ {
		a, b = NumAdd(a, b), a
	}
	return b
}

很可惜当输入为 n = 10000000 时这段代码在我的电脑上跑了足足有 15 分钟,风扇呼呼地转。

矩阵快速幂优化(V1)

因为从定义我们很容易得到:

\[R = \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} ^ n = T^n\]

问题转成了如何快速求出 R,这里对 n 分奇偶讨论:

\[R = T^n = \begin{cases} T^k \cdot T^k , & \text{if $n = 2k$} \\ T^k \cdot T^k \cdot T , & \text{if $n = 2k+1$} \\ \end{cases}\]

根据这一版的递推公式来实现的代码:

func fibV1Recursive(n int) Number {
	defer logElapsed("fibV1Recursive")()
	return fibV1Helper(n)[1][0]
}

func fibV1Helper(n int) Matrix {
	if n <= 0 {
		return Mat0
	}
	r := fibV1Helper(n / 2)
	if n % 2 != 0 {
		r = MatMul(MatMul(r, r), Mat1)
	} else {
		r = MatMul(r, r)
	}
	return r
}

当输入为 n = 10000000 时大约 2.435 秒出结果,当然我们也可以把递归改成迭代(速度没有明显差距):

func fibV1Iterative(n int) Number {
	defer logElapsed("fibV1Iterative")()
	r := Mat0
	bs := strconv.FormatInt(int64(n), 2)
	bs = strings.TrimLeft(bs, "0")
	for _, b := range bs {
		if b != '0' {
			r = MatMul(MatMul(r, r), Mat1)
		} else {
			r = MatMul(r, r)
		}
	}
	return r[1][0]
}

矩阵对称性优化(V2)

之所以写这篇文章,主要就是因为网上大部分的资料都止步于矩阵快速幂,其实在这基础上还有可以秀微操的空间。我们再来仔细看看 R 的定义:

\[R = \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} = \begin{pmatrix} x & z \\ y & t \end{pmatrix}\]

你会发现之前的运算其实做了很多无用功:

  • z 和 y 必然相等,z 没必要再计算一遍。
  • t = x - y,因此 t 也没必要再计算一遍。

最后你会发现,原来只需要计算矩阵第一列的那两个元素即可:

\[x_n = \begin{cases} x_k \cdot x_k + y_k \cdot y_k , & \text{if $n = 2k$} \\ (x_k + y_k) \cdot x_k + x_k \cdot y_k , & \text{if $n = 2k+1$} \\ \end{cases} \\ y_n = \begin{cases} x_k \cdot y_k + y_k \cdot (x_k - y_k) , & \text{if $n = 2k$} \\ x_k \cdot x_k + y_k \cdot y_k , & \text{if $n = 2k+1$} \\ \end{cases}\]

根据这一版的公式再次实现代码:

func fibV2Recursive(n int) Number {
	defer logElapsed("fibV2Recursive")()
	_, y := fibV2Helper(n)
	return y
}

func fibV2Helper(n int) (Number, Number) {
	if n <= 0 {
		return Num1, Num0
	}
	x, y := fibV2Helper(n / 2)
	xx := NumMul(x, x)
	xy := NumMul(x, y)
	yy := NumMul(y, y)
	if n % 2 != 0 {
		x = NumAdd(NumAdd(xx, xy), xy)
		y = NumAdd(xx, yy)
	} else {
		x = NumAdd(xx, yy)
		y = NumAdd(xy, NumSub(xy, yy))
	}
	return x, y
}

这次只需要 0.843 秒了,同样我们可以把递归改成迭代(速度没有明显差距):

func fibV2Iterative(n int) Number {
	defer logElapsed("fibV2Iterative")()
	x, y := Num1, Num0
	bs := strconv.FormatInt(int64(n), 2)
	bs = strings.TrimLeft(bs, "0")
	for _, b := range bs {
		xx := NumMul(x, x)
		xy := NumMul(x, y)
		yy := NumMul(y, y)
		if b != '0' {
			x = NumAdd(NumAdd(xx, xy), xy)
			y = NumAdd(xx, yy)
		} else {
			x = NumAdd(xx, yy)
			y = NumAdd(xy, NumSub(xy, yy))
		}
	}
	return y
}

一种更慢的写法(Slow)

其实不管是V1还是V2,除了上面的递归和迭代两种写法,网上还能见到另一种写法:

func fibV1Slow(n int) Number {
	defer logElapsed("fibV1Slow")()
	r, a := Mat0, Mat1
	for m := n; m > 0; m /= 2 {
		if m != n {
			a = MatMul(a, a)
		}
		if m%2 != 0 {
			r = MatMul(r, a)
		}
	}
	return r[1][0]
}

func fibV2Slow(n int) Number {
	defer logElapsed("fibV2Slow")()
	x, y, a, b := Num1, Num0, Num1, Num1
	for m := n; m > 0; m /= 2 {
		if m != n {
			aa := NumMul(a, a)
			ab := NumMul(a, b)
			bb := NumMul(b, b)
			a = NumAdd(aa, bb)
			b = NumAdd(ab, NumSub(ab, bb))
		}
		if m%2 != 0 {
			xa := NumMul(x, a)
			xb := NumMul(x, b)
			ya := NumMul(y, a)
			yb := NumMul(y, b)
			x = NumAdd(xa, yb)
			y = NumAdd(ya, NumSub(xb, yb))
		}
	}
	return y
}

很不幸这种写法虽然更容易实现,但在绝大多数情况下要更慢一些,下表是各个函数在我电脑上的运行时间汇总:

n V1Rec V1Iter V1Slow V2Rec V2Iter V2Slow
10000000 2.435s 2.423s 3.254s 0.843s 0.837s 1.374s
16777215 5.353s 5.329s 7.785s 1.909s 1.887s 3.558s
16777216 5.314s 5.315s 5.460s 2.023s 1.969s 1.875s

想知道为什么可以查一查秦九韶算法(Horner’s method),我们原始的写法实际上是利用了这个方法进行优化的。

补充版本:Cassini’s formula(V3)

V2版本看起来似乎已经做到极致了,因为对于一个 n 它只需要执行 3*log(n) 次乘法运算(数据大到一定程度时乘法运算是主要开销)。

还能再优化吗?可以的。

再来分析V2版本里面的 3 次核心乘法,其本质上就是:

\[\begin{cases} xx = f_{n+1} \cdot f_{n+1} \\ xy = f_{n+1} \cdot f_{n} \\ yy = f_{n} \cdot f_{n} \\ \end{cases}\]

这时候我们需要用到一个叫做 Cassini’s formula 的公式了:

\[f_{n+1} \cdot f_{n-1} - f_{n} \cdot f_{n} = (-1)^{n}\]

证明很简单,考虑 R 矩阵的行列式就可以了,这里不多赘述。利用这个公式,我们可以进行以下推导:

\[f_{n+1} \cdot f_{n} = f_{n+1} \cdot (f_{n+1} - f_{n-1}) = f_{n+1} \cdot f_{n+1} - f_{n} \cdot f_{n} - (-1)^{n}\]

也就是说,原来 xy 还可以由 xx 和 yy 得到,免去了一次乘法,基于此优化的代码如下:

func fibV3Recursive(n int) Number {
	defer logElapsed("fibV3Recursive")()
	_, y := fibV3Helper(n)
	return y
}

func fibV3Helper(n int) (Number, Number) {
	if n <= 0 {
		return Num1, Num0
	}
	x, y := fibV3Helper(n / 2)
	xx := NumMul(x, x)
	yy := NumMul(y, y)
	xy := NumSub(xx, yy)
	if n / 2 % 2 != 0 {
		xy = NumAdd(xy, Num1)
	} else {
		xy = NumSub(xy, Num1)
	}
	if n % 2 != 0 {
		x = NumAdd(NumAdd(xx, xy), xy)
		y = NumAdd(xx, yy)
	} else {
		x = NumAdd(xx, yy)
		y = NumAdd(xy, NumSub(xy, yy))
	}
	return x, y
}

这段代码在 n = 10000000 时运行了 0.522 秒,比起V2版本的 0.843 秒还要好不少,当然我们也可以把递归改成迭代(速度没有明显差距):

func fibV3Iterative(n int) Number {
	defer logElapsed("fibV3Iterative")()
	x, y := Num1, Num0
	bs := strconv.FormatInt(int64(n), 2)
	bs = strings.TrimLeft(bs, "0")
	prev := '0'
	for _, curr := range bs {
		xx := NumMul(x, x)
		yy := NumMul(y, y)
		xy := NumSub(xx, yy)
		if prev != '0' {
			xy = NumAdd(xy, Num1)
		} else {
			xy = NumSub(xy, Num1)
		}
		if curr != '0' {
			x = NumAdd(NumAdd(xx, xy), xy)
			y = NumAdd(xx, yy)
		} else {
			x = NumAdd(xx, yy)
			y = NumAdd(xy, NumSub(xy, yy))
		}
		prev = curr
	}
	return y
}

最后附上相关的辅助代码:

type Number = *big.Int

var (
	Num0 = big.NewInt(0)
	Num1 = big.NewInt(1)
)

func NumAdd(x Number, y Number) Number {
	return big.NewInt(0).Add(x, y)
}

func NumSub(x Number, y Number) Number {
	return big.NewInt(0).Sub(x, y)
}

func NumMul(x Number, y Number) Number {
	return big.NewInt(0).Mul(x, y)
}

type Matrix = [2][2]Number

var (
	Mat0 = Matrix{ {Num1, Num0}, {Num0, Num1} }
	Mat1 = Matrix{ {Num1, Num1}, {Num1, Num0} }
)

func MatAdd(x Matrix, y Matrix) (z Matrix) {
	z[0][0] = NumAdd(x[0][0], y[0][0])
	z[0][1] = NumAdd(x[0][1], y[0][1])
	z[1][0] = NumAdd(x[1][0], y[1][0])
	z[1][1] = NumAdd(x[1][1], y[1][1])
	return
}

func MatSub(x Matrix, y Matrix) (z Matrix) {
	z[0][0] = NumSub(x[0][0], y[0][0])
	z[0][1] = NumSub(x[0][1], y[0][1])
	z[1][0] = NumSub(x[1][0], y[1][0])
	z[1][1] = NumSub(x[1][1], y[1][1])
	return
}

func MatMul(x Matrix, y Matrix) (z Matrix) {
	z[0][0] = NumAdd(NumMul(x[0][0], y[0][0]), NumMul(x[0][1], y[1][0]))
	z[0][1] = NumAdd(NumMul(x[0][0], y[0][1]), NumMul(x[0][1], y[1][1]))
	z[1][0] = NumAdd(NumMul(x[1][0], y[0][0]), NumMul(x[1][1], y[1][0]))
	z[1][1] = NumAdd(NumMul(x[1][0], y[0][1]), NumMul(x[1][1], y[1][1]))
	return
}

func logElapsed(name string) func() {
	start := time.Now()
	return func() {
		log.Printf("%s: %v", name, time.Since(start).Seconds())
	}
}