Egyptian!

単位分数展開とは、1つの有理数 (分数) を分子が1で分母が自然数である単位分数 (unit fraction) の和として表すことです。エジプト式分数 (Egyptian fraction)、単位分数分解とも呼ばれます。単位分数展開は一意に決まるものではなく、あらゆる分数は無数の単位分数展開を持っています。

単位分数展開

項	単位分数の分母	部分和

[ 戻る ]

使用方法

分子 (numerator)・分母 (denominator) ともに1以上9007199254740991 (=2⁵³−1) 以下の自然数 (natural number) を扱うことができます。
[計算]ボタンで計算します。
[クリア]ボタンで入力値と結果を消去できます。
算出する項目の意味は次の通りです。

単位分数の分母

単位分数展開したときの各々の単位分数の分母の値

部分和

最大の単位分数からその項の分数までの総和

[ 戻る ]

利用上の注意

Microsoft Edge、Internet Exploler 11、FireFox 45.0以降での動作を確認しています。
その他のブラウザでの動作は保証いたしません。
内部演算で9007199254740992 (=2⁵³) 以上になったときは不正確な値になることがあります。
サポートは一切行いません。
気が向いたときに改良や修正をします。

[ 戻る ]

貪欲法

貪欲法 (greedy algorithm) は逐次差分アルゴリズムの一種で、フィボナッチ＝シルベスター法 (Fibonacci-Sylvester method) とも呼ばれます。与えられた分数 a/b を超えない最大の単位分数 1/⌈b/a⌉ を引き算 (差分を計算) し、残った残差に対して同じ操作を繰り返す手法で、次式の2項目である残差の分子が0になった時点で終了します。

\frac{a}{b} = \frac{1}{⌈ \frac{b}{a} ⌉} + \frac{a - (b mod a)}{b \times ⌈ \frac{b}{a} ⌉} = \frac{1}{⌈ \frac{b}{a} ⌉} + \frac{a \times ⌈ \frac{b}{a} ⌉ - b}{b \times ⌈ \frac{b}{a} ⌉}

ここで、⌈b/a⌉ の括弧は切り上げを意味する天井関数 (ceiling function) です。この方法の特徴は、必ず有限回数で終了しますが、分母が爆発的に巨大化しやすいことです。

var greedy_algorithm = function(a, b) {
/* 貪欲法 (greedy algorithm) */
	var d = [];
	while (a > 0) {
		// 既約分数化
		var g = extended_euclidean(a, b).gcd;
		a /= g;
		b /= g;

		// a/bを越えない最大の単位分数の分母b_n
		var b_n = Math.ceil(b / a);
		d.push(b_n);

		// 残差 (a/b - 1/b_n) の計算
		a = a * b_n - b;
		b *= b_n;
	}
	return d;
};

[ 戻る ]

Golombの方法

Golombの方法 (Golomb's method) は拡張ユークリッドの互除法 (Extended Euclidean algorithm) を適用して a・a^-1≡1 (mod b) を満たす逆元 (modular multiplicative inverse) a^-1 (コード内の mod_inv_a) を見つけ、等式 a・a^-1 = b×r+1 を作ることで展開する手法です。この等式の両辺を b×a^-1 で割ると次式が得られます。

\frac{a}{b} = \frac{1}{b \times a^{−1}} + \frac{r}{a^{−1}}

1項目が単位分数として分離されますので、2項目を新たに a/b と置いて、2項目の分子が r=1 となるまで繰り返して計算します。この方法の最大のメリットは展開された分母が最大でも b×(b−1) を超えないことが保証されている点です。

var golomb_method = function(a, b) {
/* Golombの方法 (Golomb's method) */
	var d = [];

	// 分子が 0 になるまで非再帰ループ
	while (a > 0) {
		// 既約分数化
		var g = extended_euclidean(a, b).gcd;
		a /= g;
		b /= g;

		// 分子が 1 になったら終了
		if (a === 1) {
			d.sort(function(a, b) {
				return a - b;	// 昇順に並び替え
			});
			d.push(b);
			break;
		}

		// b*x - num*y = gcd(a, b) = 1 を満たす x, y を求める
		var ext_gcd = extended_euclidean(a, b);
		var mod_inv_a = ext_gcd.x;	// a の b を法とする逆元(候補)

		// 逆元 mod_inv_a を 0 < mod_inv_a < b に調整
		mod_inv_a = mod_inv_a % b;
		if (mod_inv_a <= 0) {
			mod_inv_a += b;
		}

		// 公式: a * mod_inv_a = b * next_a + 1 から、次の分子 next_a を一発計算
		var next_a = Math.round((a * mod_inv_a - 1) / b);

		// Golombの単位分数の分母 = b * mod_inv_a
		var next_b = b * mod_inv_a;
		d.push(next_b);
		a = next_a;
		b = mod_inv_a;
	};
	d.sort(function(a, b) {
		return a - b;	// 昇順に並び替え
	});
	return d;
};

[ 戻る ]

連分数法

分数 a/b を連分数 (continued fraction) に展開 [cf₀ ; cf₁, cf₂, …] したときに得られる近似分数 (approximate fraction) 列 p_n/q_n の階差 (difference) p_n/q_n − p_n-1/q_n-1は正または負の単位分数になります。このことを利用すると交替符号型単位分数展開（Alternating Egyptian Fraction Expansion）が実現します。

\frac{p_{n}}{q_{n}} - \frac{p_{n - 1}}{q_{n - 1}} = \frac{p_{n} q_{n - 1} - p_{n - 1} q_{n}}{q_{n} q_{n - 1}} = {(−1)}^{n - 1} \frac{1}{q_{n} q_{n - 1}}

通常の単位分数展開 (エジプト分数展開) は全て足し算で表すのに対し、交替符号型は「足す、引く、足す、引く…」を繰り返しています。なお、交替符号型単位分数展開には、貪欲法を適用したパース展開 (Pierce expansion) でも実現できます。

var alternate_fraction = function(a, b) {
/* 交替型連分数法 (alternating continued fraction method) */
	var i, j;

	// 連分数展開の係数計算
	var cf = [];
	cf = (function(a, b) {
		var q = [], r;
		var i = 0;
		while (b !== 0) {
			q.push(Math.floor(a / b));
			r = a % b;
			a = b;
			b = r;
		}
		return q;
	})(a, b);

	// 近似分数の分母数列 Q_k の計算
	var Q = [];
	var q_minus2 = 1; // Q_{-2} = 1
	var q_minus1 = 0; // Q_{-1} = 0
	for (i = 0; i < cf.length; i++) {
		var nextQ = q_minus1 * cf[i] + q_minus2;
		Q.push(nextQ);
		q_minus2 = q_minus1;
		q_minus1 = nextQ;
	}

	// 隣接する Q の項の積を分母とした単位分数リストの作成
	var denominators = [];
	var sign = 1;
	for (i = 0; i < Q.length - 1; i++) {
		denominators.push(sign * Q[i] * Q[i + 1]);
		sign *= -1;
	}

	return denominators;
};

連分数を使った通常の単位分数展開は、数学者 David Eppstein 教授が提唱した「近似分数を奇数の項数に調整し、中間近似分数 (Secondary Convergent) の分母の積を用いたアルゴリズム」で実現できます。このアルゴリズムによる展開結果は、Golombの方法の計算結果と完全に一致します。

var continued_fraction = function(a, b) {
/* 連分数法 (continued fraction method) */
	var i, j;

	// 連分数展開の係数計算
	var cf = [];
	cf = (function(a, b) {
		var q = [], r;
		var i = 0;
		while (b !== 0) {
			q.push(Math.floor(a / b));
			r = a % b;
			a = b;
			b = r;
		}
		return q;
	})(a, b);

	// アルゴリズムの性質上、連分数の長さは必ず「奇数」
	// 長さが偶数の場合は最後の項 a_n を (a_n - 1) と 1 に分割
	if (cf.length % 2 === 0) {
		cf[cf.length - 1] = cf[cf.length - 1] - 1;
		cf.push(1);
	}

	// 近似分数の分母数列 Q_k の計算
	var Q = [];
	var q_minus2 = 1; // Q_{-2} = 1
	var q_minus1 = 0; // Q_{-1} = 0
	for (i = 0; i < cf.length; i++) {
		var nextQ = q_minus1 * cf[i] + q_minus2;
		Q.push(nextQ);
		q_minus2 = q_minus1;
		q_minus1 = nextQ;
	}

	// 中間近似分数の分母数列 M の生成
	var M = [];
	i = 0;
	while (i < Q.length) {
		if (i + 2 < Q.length) {
			var q0 = Q[i];
			var q1 = Q[i + 1];
			var q2 = Q[i + 2];
			var limit = Math.floor((q2 - q0) / q1);
			for (j = 0; j < limit; j++) {
				M.push(q0 + j * q1);
			}
			i += 2;
		} else {
			M.push(Q[i]);
			i += 1;
		}
	}

	// 隣り合う M の項の積を分母とした単位分数リストの作成
	var denominators = [];
	for (i = 0; i < M.length - 1; i++) {
		denominators.push(M[i] * M[i + 1]);
	}

	return denominators;
};

[ 戻る ]

約数法

約数法 (divisor method) は、元の分数の分母 b をその約数の和で分子を表現できるような都合のよい倍数 (multiple) に拡張し、分子を約数の和に分解した後に個別に約分して単位分数を得る手法です。貪欲法とは異なり項数が少なく綺麗に分解しやすい特徴があります。

var divisor_method = function(a, b) {
/* 約数法 (divisor method) */
	// 既約分数化
	var g = extended_euclidean(a, b).gcd;
	a /= g;
	b /= g;

	// 分子が1なら終了
	var d = [];
	if (a === 1) {
		d.push(b);
		return d;
	}

	// 分母のk倍の数 b_k の約数の組み合わせを探索
	for (var k = 1; k <= 1000; k++) {	// 無限ループ防止のため上限を設定
		var b_k = b * k;
		var a_k = a * k;

		// b_k の全約数を求める
		var divisors = [];		// b_kの約数を格納
		for (var i = 1; i * i <= b_k; i++) {
			if (b_k % i == 0) {
				divisors.push(i);
				if (i * i != b_k) {
					divisors.push(b_k / i);
				}
			}
		}
		divisors.sort(function(a, b) {
			return b - a;	// 降順に並び替え
		});

		// 分子が分母の約数のときは単位分数になるので、和が a_k になる約数の組み合わせを選ぶ
		//（貪欲法的に大きい約数から比較）
		var current_sum = 0;
		var chosen_divisors = [];

		for (var j= 0; j < divisors.length; j++) {
			var div = divisors[j];
			if (current_sum + div <= a_k) {
				current_sum += div;
				chosen_divisors.push(div);
			}
			if (current_sum === a_k) {
				break;
			}
		}

		// 選択完了
		if (current_sum === a_k) {
			for (var i = 0; i < chosen_divisors.length; i++) {
				var d = b_k / chosen_divisors[i];
				d.push(d);
			}
			return d;
 		}
	}
	alert('単位分数に分解できません。');
	return [];
};

[ 戻る ]

バイナリ剰余法

バイナリ剰余法 (binary remainder method) は分母を2のべき乗に揃えるように処理を施し、分子を2進数 (binary number) のビット表現 (bit representation / binary representation) して展開・処理していく手法です。

分数 a/bに対し、b<2^kとなる2のべき乗P=2^kを求めます。a×Pをbで割ったときの商をq、剰余をrとすると (a×P=b×q+r)、以下の関係が得られます。

\frac{a}{b} = \frac{q}{P} + \frac{r}{P \times b}

ここで、分母が2のべき乗であるq/Pとr/Pは、分子を2進数 (2のべき乗和) 展開するだけで一発で単位分数に分解できます。後半の項はそれに1/bを掛けるだけとなっています。

var binary_method = function(a, b) {
/* バイナリ剰余法 (binary remainder method) */
	// 既約分数化
	var g = extended_euclidean(a, b).gcd;
	a /= g;
	b /= g;

	// 分子がのときは分母 b を追加し終了
	if (a == 1) {
		return [b];
	}

	// b < 2^k となる最小の k を見つける (P=2^k)
	var k = 0;
	var P = 1;
	while (P <= b) {
		P *= 2;
		k++;
	}

	// 分子に 2^k を掛けた値を分母で割り、商と剰余を求める
	// a * 2^k = q * b + r;
	var dividend = a * P;				// 割られる数 (dividend)
	var q = Math.floor(dividend / b);	// 商 (quotient)
	var r = dividend % b;				// 剰余 (remainder)

	var bs = [];		// 2のべき乗和

	// 商 q のバイナリ (2進数) 展開 (q/Pの部分) 
	var temp_q = q;
	for (var i = 0; i <= k; i++) {
		if (temp_q % 2 == 1) {
		// 2^i / 2^k = 1 / 2^(k-i)
			bs.push(pow(2, k - i));
		}
		temp_q = Math.floor(temp_q / 2);
	}

	// 剰余の2進数展開 (r / (b * P) の部分)
	var temp_r =r;
	for (var j = 0; j <= k; j++) {
		if (temp_r % 2 == 1) {
		// 1 / (2^(k-j) * b)
			bs.push(pow(2, k - j) * b);
		}
		temp_r = Math.floor(temp_r / 2);

	}

	// 配列内の重複要素を削除
	var d = bs.filter(function(item, pos) {
		return bs.indexOf(item) == pos;
	});

	// 昇順に並べ替え
	return d.sort(function(a, b) {
		return a - b;
	});
};

[ 戻る ]

エンゲル展開

エンゲル展開 (Engel expansion) は一意の単位分数展開が得られ、分母が直前項の倍数となる特長を持っています。まずは初期値 x₀ = a/b とします。そして、c_n=⌈1/x_n−1⌉ (1/x_n−1 の天井関数・切り上げ) の計算と、次の残差 x_n=c_nx_n-1−1 の更新を x_n=0になるまで繰り返します。これにより得られた c_n を使って単位分数展開をします。

\frac{a}{b} = \frac{1}{c_{1}} + \frac{1}{c_{1} c_{2}} + \frac{1}{c_{1} c_{2} c_{3}} + \dots + \frac{1}{c_{1} c_{2} \dots c_{n}} + \dots

ここで、c₁, c₂, c₃, … は c₁<c₂<c₃<… を満たす厳密に増加する自然数列です。

var engel_expansion = function(a, b) {
/* エンゲル展開 (Engel expansion) */
	var denominators = [];
	var c = 1;
	while (num > 0) {
		ci = Math.ceil(b / a);		
		a = a * ci - b;
		c *= ci;
		denominators.push(c);
	}
	return denominators
};

[ 戻る ]

パース展開

パース展開 (Pierce expansion) はエンゲル展開に類似した展開法で、一意の交替符号型単位分数展開が得られます。交互エジプト積 (alternated Egyptian product) とも呼ばれます。まずは初期値 x₀ = a/b とします。そして、c_n=⌊1/x_n−1⌋ (1/x_n−1 の床関数・整数部分) の計算と、次の残差 x_n=1−c_nx_n-1 の更新を x_n=0になるまで繰り返します。これにより得られた c_n を使って単位分数展開をします。

\frac{a}{b} = \frac{1}{c_{1}} - \frac{1}{c_{1} c_{2}} + \frac{1}{c_{1} c_{2} c_{3}} - \dots + \frac{{(−1)}^{n - 1}}{c_{1} c_{2} \dots c_{n}} + \dots

ここで、c₁, c₂, c₃, … は c₁<c₂<c₃<… を満たす厳密に増加する自然数列です。

var pierce_expansion = function(a, b) {
/* パース展開 (Pierce expansion) */
	var denominators = [];
	var c = 1;
	var sign = 1;
	while (a > 0) {
		ci = Math.floor(b / a);		
		a = b - a * ci;
		c *= ci;
		denominators.push(sign * c);
		sign *= -1;
	}
	return denominators
};

[ 戻る ]

アルゴリズム比較

各アルゴリズムの大まかな特徴を比較すると以下のようになります。

アルゴリズム名	最大分母	項数	計算の難易度	実装の難易度
貪欲法	爆発的に巨大化	比較的少ない	非常に簡単	極めて容易
Golombの方法	`b`(`b`−1) 以下	中程度	複雑	やや複雑
交替符号型連分数法	`b`(`b`−1) 以下	比較的少ない	比較的容易	比較的容易
連分数法	`b`(`b`−1) 以下	中程度	比較的容易	比較的容易
約数法	比較的抑えられる	比較的コンパクト	やや複雑	中程度
バイナリ剰余法	2のべき乗依存	ビット数に比例	機械的で簡単	やや複雑
エンゲル展開	直前項の倍数	比較的少ない	非常に簡単	極めて容易
パース展開	直前項の倍数	比較的少ない	非常に簡単	極めて容易

[ 戻る ]

謝辞?

アルゴリズムなどの検証でGoogle検索のAIチャット (Gemini) には大変お世話になりました…

[ 戻る ]