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;
}