トップ > math

14問目

Project Euler math

正の整数の数列が以下の反復規則で定められている.
\begin{align}\begin{cases}n\quad\to n/2 & (n \text{ is even})\\n\to 3n+1 & (n \text{ is odd})\end{cases}\end{align}.
上記のルールを使って13から始めると, 次の数列が得られる.
\begin{align}13\to40\to20\to10\to5\to16\to8\to4\to2\to1\end{align}.
この数列には(13から1まで)10個の数字が含まれている. これはまだ未解決問題(Collatz 問題)なのだが, どの数字から始めても1に辿り着くと考えられている.
1000000 未満のどの数字で始めた場合に, 数列が一番長くなるか?

http://projecteuler.net/index.php?section=problems&id=14

題材が有名な未解決問題だということで, 非常にやる気が上がりました. (変な云い方だが, 私はこういう問題に萌えます.)

しっかしこの問題も糸口ってどこにあるんでしょうね.
まず第一感として, 「保存量を見付ける」か「1回の反復で1だけ減る指標を見付ける」という方針が思い付くのですが, 反復ルールが mod しか使ってない上に, 半分にしたり3倍して1を足したりされてしまうと, 保存する量というのはなかなか考え付かないですね. 2つ目の方針も行き当たりばったりでは見付からなさそうだし. さすが未解決問題!

さてプログラムの方ですが, 如何に計算量を節約できるかをテーマに組んでみました.
問題文の例で云うと, 13から始めたときの数列の長さを調べたら, 同時に1, 2, 4, 5, 8, 10, 16, 20, 40から始めたときの数列の長さも分かっています. なので二度手間を省くため, 上の反復規則を逆にして,
\begin{align}\begin{cases}n\to 2n, \frac{n-1}{3} & (\text{if } n = 4 \pmod{6})\\ n\to 2n & (\text{otherwise})\end{cases}\end{align}
1から始めて長さを計算しています.
かと云ってこれで全てを計算すると関数呼び出しが多くなってしまうのでほどほどで止めて, 計算し切っていない部分は最後に個々に計算してあります. gdb 使ってバグを修正しつつ進め, 無事クリア.

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

#define MAX_NUMBER 1000000
#define MAGNIF 30

int numbers[MAGNIF * MAX_NUMBER];

int solveVer00();
int solveVer01();
void initArray();
void CollatzCountUp(int, int);
long long CollatzCountDown(long long);

int main() {
	fprintf(stdout, "Solution: %d\n", solveVer01());

	return EXIT_SUCCESS;
}

int solveVer00() {
	int i;
	int maxCount;
	int maxNum;
	int count;

	numbers[1] = 1;
	for (i = 1, maxCount = 0; i < MAX_NUMBER; i++) {
		count = CollatzCountDown(i);
		if (count > maxCount) {
			maxCount = count;
			maxNum = i;
		}
	}

	return maxNum;
}

int solveVer01() {
	int i;
	int max;
	int maxNum;

	initArray();
	CollatzCountUp(1, 1);

	for (i = 1, max = 0; i < MAX_NUMBER; i++) {
		if (numbers[i] == 0) {
			numbers[i] = CollatzCountDown(i);
		}
		if (numbers[i] > max) {
			max = numbers[i];
			maxNum = i;
		}
	}

	return maxNum;
}

void initArray() {
	int i;
	for (i = 0; i < MAGNIF * MAX_NUMBER; i++) {
		numbers[i] = 0;
	}
}

void CollatzCountUp(int number, int count) {
	numbers[number] = count;
	if ((number % 6) == 4
		&& numbers[(number-1)/3] == 0) {
		CollatzCountUp((number - 1) / 3, count + 1);
	}

	if (number > MAGNIF * MAX_NUMBER / 2) {
		return;
	} else {
		CollatzCountUp(number * 2, count + 1);
	}
}

long long CollatzCountDown(long long number) {
	if (number < MAGNIF * MAX_NUMBER
		&& numbers[number] > 0) {
		return numbers[number];
	} else if ((number % 2) == 0) {
		return CollatzCountDown(number / 2) + 1;
	} else {
		return CollatzCountDown(3 * number + 1) + 1;
	}
}

48問目

Project Euler math

1^1 + 2^2 + 3^3 + \dots + 10^{10} = 10405071317 である.
1^1 + 2^2 + 3^3 + \dots + 1000^{1000}の下10桁を求めよ.

http://projecteuler.net/index.php?section=problems&id=48

実装自体は簡単だと高を括って powMod 関数を作ってみたが, 桁あふれに上手く対処できず, 結局 gmp さんの力を借りることにしました.

#include <stdio.h>
#include <stdlib.h>
#include "powMod.h"
#include <gmp.h>

const long long MODULO = 10000000000LL;

long long solveVer00(int);
char *solveVer01(int);

int main(int argc, char **argv) {
	int number;

	if (argc < 1 + 1) {
		fprintf(stdout, "Error: too few argument.\n");
		return EXIT_FAILURE;
	}

	number = atoi(argv[1]);
//	fprintf(stdout, "Solution: %lld.\n", solveVer00(number));
	fprintf(stdout, "Solution: %s.\n", solveVer01(number));

	return EXIT_SUCCESS;
}

long long solveVer00(int number) {
	int i;
	long long result;

	for (i = 1, result = 0; i <= number; i++) {
		result += powMod(i, i, MODULO);
		result %= MODULO;
	}

	return result;
}

char *solveVer01(int number) {
	mpz_t i;
	mpz_t result;
	mpz_t term;
	mpz_t modulo;
	char *resultString;

	mpz_init(i);
	mpz_init(result);
	mpz_init(term);
	mpz_init(modulo);

	mpz_set_si(result, 0);
	mpz_set_str(modulo, "10000000000", 10);
	for (mpz_set_si(i, 1); mpz_cmp_si(i, number) < 0; mpz_add_ui(i, i, 1)) {
		mpz_powm(term, i, i, modulo);
		mpz_add(result, result, term);
	}

	mpz_clear(i);
	mpz_clear(term);
	mpz_clear(result);
	mpz_clear(modulo);

	mpz_get_str(resultString, 10, result);

	return resultString;
}

18問目

Project Euler math

下図の三角形に並べられた数字の一番上からスタートして, 1つ下の行の隣接する数字に移動していったとき, 通った数の合計の最大値は23.
\begin{array}&&&3&&& \\ &&7&&5&& \\ &2&&4&&6& \\ 8&&5&&9&&3 \end{array}
つまり, 3 + 7 + 4 + 9 = 23. (訳注. 実は, 3+5+6+9 = 23 というルートでも良い.)

では, 次の三角形に並べられた数字の一番上から一番下までの合計の最大値を求めよ.
(図は省略.)

メモ: 取り得るルートが16384 通りしかないので, この問題はしらみつぶしでも解ける. しかし, 問題67では, 100行ある三角形に対し同じ問題を解くことになる. これは力尽くでは解けず, 賢いやり方が求められる ;o)

http://projecteuler.net/index.php?section=problems&id=18

さて, まずはやはり数学的考察から.

この三角形のサイズを n としたときに, 力尽くの方法だと O(2^n) の計算量が必要だがこれは明かに二度手間の無駄がある. これは文字通り「トップダウン」で計算を行っているからで, 「ボトムアップ」で計算を行うと O(n^2) の計算量で済む. (やり方はシンプルなものなので, 詳細はソースコードを参照してください. )
それを実装したのが以下のソース. 計算機科学の分野のアルゴリズムでは良く使われる手ですね.

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

#define SIZE 15

int triangle[][SIZE] =
{{75},
 {95, 64},
 {17, 47, 82},
 {18, 35, 87, 10},
 {20, 4, 82, 47, 65},
 {19, 1, 23, 75, 3, 34},
 {88, 2, 77, 73, 7, 63, 67},
 {99, 65, 4, 28, 6, 16, 70, 92},
 {41, 41, 26, 56, 83, 40, 80, 70, 33},
 {41, 48, 72, 33, 47, 32, 37, 16, 94, 29},
 {53, 71, 44, 65, 25, 43, 91, 52, 97, 51, 14},
 {70, 11, 33, 28, 77, 73, 17, 78, 39, 68, 17, 57},
 {91, 71, 52, 38, 17, 14, 91, 43, 58, 50, 27, 29, 48},
 {63, 66, 4, 68, 89, 53, 67, 30, 73, 16, 69, 87, 40, 31},
 {4, 62, 98, 27, 23, 9, 70, 98, 73, 93, 38, 53, 60, 4, 23}};

int solveVer00();
int max(int, int);

int main() {
	fprintf(stdout, "Solution: %d\n", solveVer00());
	return EXIT_SUCCESS;
}

int solveVer00() {
	int row;
	int column;

	for (row = SIZE - 2; row >= 0; row--) {
		for (column = 0; column <= row; column++) {
			triangle[row][column] 
				+= max(triangle[row+1][column], triangle[row+1][column+1]);
		}
	}
	return triangle[0][0];
}

int max(int num1, int num2) {
	return (num1 > num2 ? num1 : num2);
}

12問目

Project Euler math

三角数からなる数列は, 自然数を次々と足していくことで得られる. 7番目の三角数は 1+2+3+4+5+6+7 = 28 になる. 数列の最初の10項は,
1,\quad 3,\quad 6,\quad 10,\quad 15,\quad 21,\quad 28,\quad 36,\quad 45,\quad 55,\quad...
となる.

三角数の最初の7項について約数を列挙すると,
\begin{eqnarray} 1 &:& 1 \\ 3 &:& 3 \\ 6 &:& 1,\quad 2,\quad 3,\quad 6\\ 10 &:& 1,\quad 2,\quad 5,\quad 10 \\ 15 &:& 1,\quad 3,\quad 5,\quad 15 \\ 21 &:& 1,\quad 3,\quad 7,\quad 21 \\ 28 &:& 1,\quad 2,\quad 4,\quad 7,\quad 14,\quad 28 \end{eqnarray}
となる.
これから, 28 が約数の数が5を越える最小の三角数だと分かる.

では, 約数の数が 50 を越える最小の三角数はいくつか?

http://projecteuler.net/index.php?section=problems&id=12

最初は数学を使って考察. 計算量を減らしていきます.

一般に三角数は \frac{n(n+1)}{2} と表せるが,
n \equiv 0 \quad \pmod 2 のとき, \gcd\(\frac{n}{2},\quad n+1\) = 1,
n \equiv 1 \quad \pmod 2 のとき, \gcd\(n,\quad \frac{n+1}{2}\) = 1
となる. よって,
\begin{equation}d\(\frac{n(n+1)}{2}\) = \begin{cases} d\(\frac{n}{2}\)d\(n+1\)\quad \text{if}\quad n \equiv 0 \quad \pmod 2, \\ d\(n\)d\(\frac{n+1}{2}\)\quad \text{if}\quad n \equiv 1 \quad \pmod 2 \end{cases}\end{equation}
となる. ( d(n) は, n の正の約数の個数. )

これを踏まえてコーディング.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include "prime.h"

const int LIMIT = 500;
const int NUM_LIMIT = LONG_MAX;

long long solveVer00();
long numDivisors(long);

int main(int argc, char **argv) {
	initPrime();
	fprintf(stdout, "Solution: %lld\n", solveVer00());

	return 0;
}

long long solveVer00() {
	long n;


	for (n = 1; n < NUM_LIMIT; n++) {
		if ((n % 2) == 0) {
			if (numDivisors(n/2) * numDivisors(n+1) > LIMIT) {
				return (long long) (n / 2) * (n + 1);
			}
		} else {
			if (numDivisors(n) * numDivisors((n+1)/2) > LIMIT) {
				return (long long) n * (n + 1) / 2;
			}
		}
	}

	return 0;
}

long numDivisors(long number) {
	long long p;
	int power;
	long result;

	prime(-1);
	result = 1;
	while (number > 1
		   && (p = prime(0)) != -1) {
		power = 0;
		while ((number % p) == 0) {
			number /= p;
			power++;
		}
		result *= power + 1;
	}

	return result;
}

この問題を機に書き直した素数関連の関数たちはこちら.

#include "prime-2.0.h"

#define SIZE 1000000

static int numbers[SIZE];
static int primes[SIZE];
const int NOT_PRIME = 0;
const int PRIME = 1;

/**
 * numbers, primes の初期化
 */
void initPrime() {
	long i;
	long numberIndex;
	long primeIndex;

	numbers[0] = NOT_PRIME;
	numbers[1] = NOT_PRIME;
	for (numberIndex = 2; numberIndex < SIZE; numberIndex++) {
		numbers[numberIndex] = PRIME;
	}
	for (primeIndex = 0; primeIndex < SIZE; primeIndex++) {
		primes[primeIndex] = NOT_PRIME;
	}

	for (numberIndex = 2, primeIndex = 0; numberIndex < SIZE; numberIndex++) {
		if (numbers[numberIndex] == PRIME) {
			for (i = 2; i * numberIndex < SIZE; i++) {
				numbers[i * numberIndex] = NOT_PRIME;
			}
			primes[primeIndex] = numberIndex;
			primeIndex++;
		}
	}
}

/**
 * 小さい方から value 番目の素数を返す.
 * 0 を入れると小さい方から順々に素数を返す.
 * -1 を入れることでリセットでき, 再び2から返していく.
 * インデックスが大きすぎる場合は -1 を返す.
 */
long long prime(long long value) {
	static long long primeIndex = 0;

	if (value > SIZE) {
		return -1;
	}

	if (value == -1) {
		primeIndex = 0;
		return 0;
	} else if (value == 0) {
		return primes[primeIndex++];
	} else if (value > 0) {
		return primes[value];
	} else {
		return -1;
	}
}

int isPrime(long long value) {
	return numbers[value] == PRIME;
}

問題34

Project Euler math

145 は興味深い数字である. なぜなら, 1! + 4! + 5! = 1 + 24 + 120 = 145 となるからである.
各桁の数字の階乗の合計が自身と等しくなる数の合計を求めよ.
ただし, 1! = 1 と 2! = 2 は含まない.

http://projecteuler.net/index.php?section=problems&id=34

「問題30と同じじゃん. プログラムがそのまま使えるぞ! 上限は…(計算中)…7桁までの数字でいいんだな.」と意気込んでやってみた結果は……まず失敗……. 145 以外の答えが見付からない……orz
「なんでー?」と思いながら, アルゴリズムの改良なんかをやってみたりしてたら, ふと気が付いた. そうだ 0! = 1 じゃん, 0! = 0 で計算してたら合うわけないよ.
最終学歴を抹消した方が良いでしょうか? いやー, めっちゃ凹んだ. やっぱ四六時中数学に触れてないと, そういう細かいところが鈍るのね…….

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

const long long RADIX = 10;
const long long UPPER_BOUND = 10000000LL;
const int fact[10] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};

long long solveVer00();
long long solveVer01();
long long sumOfFactorials(long long);
int curious(long long);

int main() {
	fprintf(stdout, "Solution: %lld\n", solveVer01());

	return EXIT_SUCCESS;
}

long long solveVer00() {
	long long number;
	long long result;

	for (number = 3, result = 0; number < UPPER_BOUND; number++) {
		// fprintf(stdout, "%lld, %lld\n", number, sumOfFactorials(number));
		if (curious(number)) {
			fprintf(stdout, "%lld\n", number);
			result += number;
		}
	}
	return result;
}

long long solveVer01() {
	long long number;
	long long result;
	int lastDigit;
	long long upperPart;
	long long sum;

	for (upperPart = 0, result = 0; upperPart < UPPER_BOUND / RADIX; upperPart++) {
		for (lastDigit = 0; lastDigit < RADIX; lastDigit++) {
			number = upperPart * RADIX + lastDigit;
			sum = sumOfFactorials(number);
			if (sum < 3) {
				continue;
			}
			if (sum > number) {
				break;
			} else if (sum == number) {
				fprintf(stdout, "%lld\n", number);
				result += number;
			}
		}
	}
	return result;
}

int curious(long long number) {
	int digit;
	long long sum;
	long long copy;

	sum = 0;
	copy = number;
	while (copy) {
		digit = copy % RADIX;
		copy /= RADIX;
		sum += fact[digit];
		if (number < sum) {
			return 0;
		}
	}
	return number == sum;
}

long long sumOfFactorials(long long number) {
	int digit;
	long long result;

	result = 0;
	while (number) {
		digit = number % RADIX;
		number /= RADIX;
		result += fact[digit];
	}

	return result;
}

問題30

Project Euler math

驚くべきことに, 各桁の数字を4乗して合計すると元の数字になる数字は次の3つだけである.
\begin{eqnarray} 1634 &=& 1^4 + 6^4 + 3^4 + 4^4 \\ 8208 &=& 8^4 + 2^4 + 0^4 + 8^4 \\ 9474 &=& 9^4 + 4^4 + 7^4 + 4^4 \end{eqnarray}
ただし, 1 = 1^4 は含まない.
それらの数の合計は, 1634 + 8208 + 9474 = 19316.

では, 各桁の数字を5乗して合計すると元の数字になる数字の合計を求めよ.

http://projecteuler.net/index.php?section=problems&id=30

小さい方から探索していけば良いのですが, 上限を決めなければいけません. そこで数学の登場です!

n 桁の数字の各桁の4乗の合計は, 高々 n*9^5 = 59049n. よって, 7桁以上の数になると絶対に問題の条件を満たさない. てなわけでコーディング.

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

const int RADIX = 10;
const int UPPER_BOUND = 1000000;

int sumOfFifthPowers(int);

int main() {
	int number;
	int result;

	for (number = 2, result = 0; number < UPPER_BOUND; number++) {
		if (number == sumOfFifthPowers(number)) {
			fprintf(stdout, "%d\n", number);
			result += number;
		}
	}
	fprintf(stdout, "Solution: %d.\n", result);

	return EXIT_SUCCESS;
}

int sumOfFifthPowers(int number) {
	int digit;
	int sum;

	sum = 0;
	while (number) {
		digit = number % RADIX;
		number /= RADIX;
		sum += pow(digit, 5);
	}

	return sum;
}