問題2改

Project Euler

ぎゃー、毎日更新するとか云っておきながら、さっそく三日坊主……。
いや、風邪が……。体調が………。
はい、すみません、四の五の云わずに問題2です。

前回: http://d.hatena.ne.jp/cocoatomo/20090419/1240148001

数学的考察

数列 a_0 = 1, a_1 = 2, a_n = a_{n-1} + a_{n-2} (n \geq 2) の一般項は a_n = \frac{1}{\sqrt{5}}\{\phi^{n+2}-(-\frac{1}{\phi})^{n+2}\}.
(式が煩雑になるので, 以下 \phi = \frac{1+\sqrt{5}}{2} と置く. 因みに, \phiは黄金比. )

Fibonacci 数列を \bmod 2 で見たときの周期は 3 で 101101101... \pmod 2 となっているので, 偶数の項は一般に, a_{3m-2} = \frac{1}{\sqrt{5}}\{\phi^{3m}-(-\frac{1}{\phi})^{3m}\} (m = 1, 2,...) と書ける.

|\phi| \approx 1.6... より (-\frac{1}{\phi})^{3m} \to 0 (m \to \infty) であることに着目すると, M = 4000000 を越えない偶数項のインデックスは以下の式を満たす.

\begin{eqnarray} \frac{1}{\sqrt{5}}\{\phi^{3m}-(-\frac{1}{\phi})^{3m}\} &\leq& M, \\ \phi^{3m} &\leq& (\sqrt{5}M + (-\frac{1}{\phi})^{3m})\\ &\leq& (\sqrt{5}M + 1), \\ m &\leq& \frac{1}{3}\log_\phi{(\sqrt{5}M + 1)}. \end{eqnarray}

これを満たす m の中で最大のものを m_0 と置くと, 求める値は

\begin{eqnarray} S &=& \sum_{m=1}^{m_0} \frac{1}{\sqrt{5}}\{\phi^{3m}-(-\frac{1}{\phi})^{3m}\}\\ &=& \frac{1}{\sqrt{5}}\{\frac{\phi^{3(m_0 + 1)} - 1}{\phi^3 - 1} - \frac{(-\phi^{-1})^{3(m_0 + 1)} - 1}{-\phi^{-3} - 1}\}\\ &=& [\frac{1}{\sqrt{5}}\frac{\phi^{3(m_0 + 1)} - 1}{\phi^3 - 1}].\end{eqnarray}

ソースコード

上の考察をそのままコードにしたような, 誤差計算も好い加減なコード.
しかし, 意外とすんなり答えが出てビックリ. 数学の偉大さを改めて感じる. (まぁ, 俺が数学大好き, ってのがあるが.)

解法バージョンの切り替えには, そろそろ getopt でも使おうかな.
http://www.gnu.org/software/hello/manual/libc/Using-Getopt.html#Using-Getopt

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

long solveVer00(long);
long solveVer01(long);

int main(int argc, char **argv) {
	int version;
	long upperBound;

	upperBound = 4000000;
	if (argc == 1) {
		version = 0;
	} else if (argc == 2) {
		version = atoi(argv[1]);
	} else {
		version = atoi(argv[1]);
		upperBound = atol(argv[2]);
	}

	switch (version) {
	case 0:
		printf("Solution: %ld\n", solveVer00(upperBound));
		break;
	case 1:
		printf("Solution: %ld\n", solveVer01(upperBound));
		break;
	default:
		break;
	}

	return 0;
}

long solveVer00(long upperBound) {
	long firstTerm, secondTerm, temp, sum;

	sum = 0;
	firstTerm = 1;
	secondTerm = 2;
	while (secondTerm <= upperBound) {
#if DEBUG
		printf("%ld, %ld\n", firstTerm, secondTerm);
#endif
		if ((secondTerm % 2) == 0) {
			sum += secondTerm;
		}
		temp = firstTerm + secondTerm;
		firstTerm = secondTerm;
		secondTerm = temp;
	}
	return sum;
}

long solveVer01(long upperBound) {
	int maxIndex;
	double phi; // golden ratio

	phi = (1 + sqrt(5)) / 2;

	maxIndex = floor(
		log(sqrt(5) * upperBound + 1)
		/
		(3 * log(phi))
		);
	return floor(
		(pow(phi, 3 * (maxIndex + 1)) - 1)
		/
		(sqrt(5) * (pow(phi, 3) - 1))
		);
}