トップページ > プログラミング研究室 > Algorithm!

Algorithm!

アルゴリズム (algorithm) とは「計算可能」な数学的問題を解くための計算手順や処理手順のことです。ここでは、簡単な数値計算に関するアルゴリズムを集めてみました。ソースコードでは関数リテラルとしていますが、関数宣言(function文)でも同様に動作します。

数値計算以外の小技については "Tips!" を参照してください。

• 累乗
• 冪剰余
• n乗根の整数部
• ホーナー法
• 漸化式による数列
• ハミング数
• ハンブル数
• ズッカーマン数
• 分割数
• 約数の総和
• 利用上の注意
• 参考文献

[ 戻る ]

累乗

実数 (real number) xの冪 (power) を求めるときはMath.pow関数を使います。ところが、Math.pow関数は小数の冪乗 (exponentiation) もサポートしているため演算速度が遅い関数 (function) です。

そこで、底 (base) がx、冪指数 (exponent) が自然数 (natural number) nである累乗 (power) 計算xⁿに繰り返し二乗法 (repeated squaring / exponentiation by squaring) を適用することで高速化を図ることができます。例えば、x¹¹は x⁸⁺²⁺¹=x⁸⋅x²⋅x¹ と表されるため、x²=x⋅x、x⁴=x²⋅x²、x⁸=x⁴⋅x⁴ を順に求めて必要なものだけを乗じます。

var pow = function(x, n) {
/* xの累乗(x^n) */
	var v = 1;
	while (n > 0) {
		if ((n & 1) == 1) {
			v *= x;	// nのLSBが1ならば、x^(2^i)をかける
		}
		x *= x;
		n >>= 1;	// 1bit右シフト
	}
	return v;
};

底に0〜94906266、冪指数に2〜53の値を入力すると、結果を即時に出力します。

[ 戻る ]

冪剰余

冪剰余 (modular exponentiation) とは冪乗の剰余 (remainder) のことで、底が自然数x、冪指数が自然数nのとき、自然数mを法 (modulus) とするxⁿの剰余演算 (modulo operation) xⁿ mod m になります。

xⁿを直接計算してからmで割った剰余を求めると、冪乗時にすぐに桁あふれ (overflow) を起こしてしまいます。そこで、前述の累乗計算の pow(x, n) 関数に、(a×b) mod m=((a mod m)×(b mod m)) mod m の関係を適用し、計算の各段階で剰余を求めることで効率よく計算ができるようになります。

var mod_pow = function(x, n, m) {
/* xの冪累乗(x^n mod m) */
	var v = 1;
	x %= m;
	while (n > 0) {
		if (n & 1) {
			v = v * x % m;	// nのLSBが1ならば、x^(2^i)をかけ、mの剰余演算
		}
		x = x * x % m;
		n >>= 1;	// 1bit右シフト
	}
	return v;
};

底に0〜94906266、冪指数に1〜2147483647、法に2〜1000000009の値を入力すると、結果を即時に出力します。

[ 戻る ]

n乗根の整数部

正数 (positive number) xを被開平数 (radicand) としたn乗根 (n-th root, 2≤n≤20) の整数部 (integer part) ⌊ⁿ√x⌋ をニュートン法 (Newton's method) で求めます。初期値 (initial value) はxを2進数 (binary number) で表したときの桁数 (number of digits) を1/n倍した桁を1とした値にしています。累乗計算は前述のpow関数を使用して高速化を図りました。

var int_root = function(x, n) {
/* n乗根の整数部 */
	var s = 1, t = x
	var m = n - 1;
	var b = pow(2, m);
	while (s < t) {
	// 初期値の決定(xの2進桁数の1/nの桁数の値→初期値)
		s *= 2;					// s <<= 1
		t = Math.floor(t / b);	// t >>= m
	}
	do {
		t = s;
		s = Math.floor((x / pow(s, m) + m * s) / n);
	} while(s < t)
	return t;
};

被開平数に0〜9007199254740992、指数に2〜20の値を入力すると、結果を即時に出力します。

[ 戻る ]

ホーナー法

ホーナー法 (Horner's method) とは、最も少ない加算・乗算回数でn次多項式 (n-th degree polynomial) を評価するアルゴリズムです。

係数 (coefficient) $a_n,a_{n-1},\; \cdots ,a_1,a_0$からなる一変数xのn次多項式

において、$x=x_0$ のときの値 $f\left(x_0\right)$ を求めるとき、直接計算する方法では n(n+1)/2 回の乗算が必要になります。ところが、多項式を以下のように変形すると乗算を n 回で済ますことができます。

この多項式の評価方法がホーナー法です。

配列aに降冪の順 (descending order of powers) で設定された係数を持つ多項式で、xのときの式の値を求めます。引数dには式全体を除する値を設定します。係数a[]が昇冪の順 (ascending order of powers) のときは、4行目を "var p = a[n];" に、5行目のfor文の条件を "var i = n - 1; i >= 0; i--" にそれぞれ変更します。

var horner = function(x, a, d) {
/* ホーナー法による多項式の値の計算 */
	var n = a.length;	// 多項式の次数
	var p = a[0];
	for (var i = 1; i < n; i++) {
		p = p * x + a[i];
	}
	return p / d;
};

係数は降冪の順にカンマ若しくはスペースで区切って入力してください。Sampleを選択すると選択した多項式の値が設定されます。内部演算で9007199254740992=2⁵³以上になると下位桁に丸め誤差 (round-off error / rounding error) が生じ、結果に影響します。

[ 戻る ]

漸化式による数列

d-階定数係数線形漸化式 (d-order linear recurrence relations with constant coefficients) とは、

の形で表されるd項+定数項c₀の漸化式 (recurrence relations) で、係数c_i は全ての i について定数 (constant) となっています。c₀=0のときは斉次 (homogeneous) 、c₀≠0のときは非斉次 (non-homogeneous / inhomogeneous) と呼ばれます。

係数c₀≠0, c₁=1と初期値a₀で初項 (first term) a₀、公差 (common difference) c₀の等差数列 (arithmetic sequence)、係数c₀=0, c₁≠0と初期値a₀≠0で初項a₀、公比 (common ratio) c₁の等比数列 (geometric sequence) になります。定数係数線形漸化式で表される数列 (sequence) には以下のものがあります。

配列cに定数係数、配列aに数列の初期値を設定します。cには定数項c₀が含まれるので要素数はd+1、aの要素数はdである必要があります。計算アルゴリズムはすでに計算された値を用いて計算する動的計画法 (dynamic programming: DP) を用いており、a[i−d]から a[i−1]の値より i 番目の値 a[i]をn項目まで順次求めています。

var recurrence = function(n, c, a) {
/* 定数係数線形漸化式 */
	var ia, ic, s;
	var d = c.length - 1;	// 漸化式の階数
	for (var i = d; i <= n; i++) {
		s = c[0];
		for (ia = i - 1, ic = 1; ic <= d; ia--, ic++) {
			s += c[ic] * a[ia];
		}
		a[i] = s;
	}
	return a;
};

係数はc₀, c₁, c₂, …, c_d 、初期値はa₀, a₁, …, a_d−1の順にそれぞれカンマ若しくはスペースで区切って入力してください。Sampleで数列を選択すると係数と初期値が設定されます。n項目までの数列を表示しますが、nが漸化式の階数 (order) d以下のときは設定した初期値を全て表示します。内部演算で9007199254740992=2⁵³以上になると下位桁に丸め誤差が生じ、結果に影響します。

[ 戻る ]

ハミング数

ハミング数 (Hamming number) とは60の累乗を割り切ることができる自然数、即ち素因数 (prime factor) として2、3、5しか持たない数 =2ⁱ×3^j×5^k (i, j, k≥0) のことです。正則数 (regular number)、5-平滑数 (5-smooth number) とも呼ばれています。

配列qに1番目 (q₀=1) からn番目 (q_n−1) までのハミング数を代入していきます。i番目までの任意のハミング数をq_j (0≤j≤i) とすると、i+1番目q_i+1は2q_j、3q_j、5q_jのうち、q_iより大きい最小のものになります。これをi=n−1になるまで繰り返し、戻り値として配列qを戻します。

n番目のみのハミング数の取得は、関数呼び出しを hamming(n).slice(-1)[0] とします。

var hamming = function(n) {
/* n番目までのハミング数リストの生成 */
	var q = [];
	var j2 = j3 = j5 = 0;
	var x2 = x3 = x5 = 1;
	var min;
	for (var i = 0; i < n; i++) {
		min = x2;
		if (x3 < min) min = x3;
		if (x5 < min) min = x5;
		q[i] = min;
		while (x2 <= min) x2 = 2 * q[j2++];
		while (x3 <= min) x3 = 3 * q[j3++];
		while (x5 <= min) x5 = 5 * q[j5++];
	}
	return q;
};

n番目 (1≤n≤7716) までのハミング数を全て列挙しています。

[ 戻る ]

ハンブル数

ハンブル数 (humble number) とは素因数として2、3、5、7しか持たない数 =2ⁱ×3^j×5^k×7^l (i, j, k, l≥0) のことで、7-平滑数 (7-smooth number) とも呼ばれています。

ハミング数と同じ考え方でハンブル数を列挙することができます。n番目のみのハンブル数の取得は、関数呼び出しを humble(n).slice(-1)[0] とします。

var humble = function(n) {
/* n番目までのハンブル数リストの生成 */
	var q = [];
	var j2 = j3 = j5 = j7 = 0;
	var x2 = x3 = x5 = x7 = 1;
	var min;
	for (var i = 0; i < n; i++) {
		min = x2;
		if (x3 < min) min = x3;
		if (x5 < min) min = x5;
		if (x7 < min) min = x7;
		q[i] = min;
		while (x2 <= min) x2 = 2 * q[j2++];
		while (x3 <= min) x3 = 3 * q[j3++];
		while (x5 <= min) x5 = 5 * q[j5++];
		while (x7 <= min) x7 = 7 * q[j7++];
	}
	return q;
};

n番目 (1≤n≤42038) までのハンブル数を全て列挙しています。

[ 戻る ]

ズッカーマン数

ズッカーマン数 (Zuckerman number) とは、各桁の数字を総乗した値である数字積 (digit product) が元の数の約数 (divisor) に含まれている自然数のことです。

数字積は最初の4つの素数2、3、5、7のみを素因数に持つハンブル数ですので、ズッカーマン数の素因数にもこれら4つの素数が含まれています。逆にズッカーマン数にならない自然数Nの特徴 (property) は下記の通りで、これに当てはまらない自然数がズッカーマン数になる可能性があります。

• 1個以上の0を含む（∵数字積が0）
• 10の倍数≡素因数として2と5を含む（∵数字積が0）
• 1桁の数を除くレピュニット (repunit) でない素数
• 数字積が10の倍数≡素因数として2と5を含む（∵数字積で割り切れるNは10の倍数）
• 十の位が偶数で一の位が奇数（∵Nは奇数だが数字積が偶数）
• 十の位が奇数で一の位が4の倍数（∵Nは2の奇数倍だが数字積が4の倍数）

n1からn2までの自然数でズッカーマン数になるものを列挙します。まずは上記の1番目の条件に当てはまる数を除外します。そして10以上では100単位でまとめ、配列ldに設定した下2桁 (last two digits, lower two digits) を持つ数に限定することで、5番目、6番目の条件に当てはまる数を除外しています。これらの枝刈り (pruning) 後に残った数のうち、数字積で割り切れた数をズッカーマン数として配列zに格納し、関数の戻り値としています。

数字積を求めるdigit_product関数については、PrimeFactor!を参照してください。

var zuckerman =  function(n1, n2) {
/* ズッカーマン数の計算 */
	var z = [];		// Zuckerman数
	var ld = [11, 12, 13, 15, 16, 17, 19,
		22, 24, 26, 28, 31, 32, 33, 35, 36, 37, 39,
		42, 44, 46, 48, 51, 52, 53, 55, 56, 57, 59,
		62, 64, 66, 68, 71, 72, 73, 75, 76, 77, 79,
		82, 84, 86, 88, 91, 92, 93, 95, 96, 97, 99];
	var i, j, k = 0, n, dp;

	if (n1 > n2) {
		n = n1;
		n1 = n2;
		n2 = n1;
	}
	for (i = n1; i < 10; i++) {
		z[k++] = i;	// 1〜9は全てZuckerman数
	}
	var i1 = Math.floor(n1 / 100) * 100;	// n1を100単位で切り捨てた数
	var i2 = Math.ceil(n2 / 100) * 100;		// n2を100単位で切り上げた数
	for (i = i1; i <= i2; i += 100) {
		if (digit_product(Math.floor(i / 100)) == 0) continue;
		for (j = 0; j <= 50; j++) {
			n = i + ld[j];
			dp = digit_product(n);
			if (dp == 0) continue;
			if (n % dp == 0 && n >= n1 && n <= n2) {
				z[k++] = n;
			}
		}
	}
	return z;
};

計算例を下表に示します。計算範囲が広いとその分時間がかかります。途中で打ち切るときは、ブラウザのタスクマネージャーから終了させてください。

[ 戻る ]

分割数

分割数 (partition function, OEIS A000041) p(n)は自然数nを1つ以上の自然数の和として表すときの項の組み合わせの総数を表す数論的関数です。p(n)はp(0)=1、任意の負の整数nに対してp(n)=0と定められています。

オイラーの五角数定理 (Euler's pentagonal number theorem) により、p(i)の漸化式 (recursion formula) が得られます。

p(i) = p(i−1) + p(i−2) − p(i−5) − p(i−7) + p(i−12) + p(i−15) − p(i−22) −⋯

ここで、pの添字の数字は一般五角数 (generalized pentagonal number, OEIS A001318) gpn=m(3m−1)/2 になっています。mを1, −1, 2, −2, 3, −3, 4, −4, …とすると、gpn=1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, …の数列が得られます。項の符号は交互に+, +, −, −, +, +, −, −, …と続きます。

p(k)を格納するための配列pを作成します。初期値 (initial value) としてp[0]に1を代入してから、1からnまでの分割数を順に求めていき、戻り値として配列pを戻します。sgnは項の符号 (sign) で、mが奇数 (odd number) のときは1、偶数 (even number) のときは−1になります。

n番目のみの分割数の取得は、関数呼び出しを partition(n).slice(-1)[0] とします。

var partition = function(n) {
/* 分割数(partition function) */
	var p = [];	// 分割数
	var i, m;
	var gpn;	// 一般五角数(generalized pentagonal number)
	var sgn;	// 和の符号+/-
	p[0] = 1;
	for (i = 1; i <= n; i++){
		p[i] = 0;
		for (m = 1, sgn = 1; ; m++, sgn *= -1) {
			gpn = m * (3 * m - 1) / 2;	// m>0
			if (i < gpn) break;
			p[i] += p[i - gpn] * sgn;
			gpn = m * (3 * m + 1) / 2;	// m<0
			if (i < gpn) break;
			p[i] += p[i - gpn] * sgn;
		}
	}
	return p;
};

n番目 (1≤n≤299) までの分割数を全て列挙しています。

[ 戻る ]

約数の総和

自然数nの約数の総和 (summation of divisor, OEIS A000203) は1次の約数関数 (divisor function) σ₁(n)になります。nを素因数分解 (prime factorization) して求めることが多いです。

約数の総和の漸化式は前述の分割数の漸化式と同じですが、初期値σ₁(0)がnになる点が異なります。そこで、σ₁(i)を計算するときの初期値σ₁(0)をiとすることで、分割数の計算と同様な手順で求めることもできます。しかし、時間計算量 (time complexity) がΟ(n^1.5)であり、n=65536ではn=1024のときの約500倍、n=4194304では約26万倍、n=4294967295では約86億倍の時間がかかることになります。このため、nが大きくなるに連れて実用的ではなくなってしまいます。

戻り値は1からnまでの自然数の約数の総和を格納した配列sigmaです。先頭の要素sigma[0]は戻す直前に削除し、1の約数の総和が先頭の要素になるようにしています。

var divisor_sum = function(n) {
/* 1次の約数関数(divisor function) */
	var sigma = [];	// 約数関数
	var i, m;
	var gpn;		// 一般五角数(generalized pentagonal number)
	var sgn;		// 和の符号+/-
	for (i = 1; i <= n; i++) {
		sigma[0] = i;
		sigma[i] = 0;
		for (m = 1, sgn = 1; ; m++, sgn *= -1) {
			gpn = m * (3 * m - 1) / 2;	// m>0
			if (i < gpn) break;
			sigma[i] += sigma[i - gpn] * sgn;
			gpn = m * (3 * m + 1) / 2;	// m<0
			if (i < gpn) break;
			sigma[i] += sigma[i - gpn] * sgn;
		}
	}
	sigma.shift();	// 先頭の要素を削除
	return sigma;
};

4194304以下の自然数nに対応しており、指定されたnまでの各自然数の約数の総和を列挙しています。

[ 戻る ]

利用上の注意

• Microsoft Edge、Internet Exploler 11、FireFox 45.0以降での動作を確認しています。
• その他のブラウザでの動作は保証いたしません。
• サポートは一切行いません。
• 気が向いたときに改良や追加をします。

[ 戻る ]

参考文献

• 奥村晴彦 『C言語による最新アルゴリズム事典』 技術評論社、平成6年9月25日初版第12刷（平成3年2月25日初版第1刷）
• Paula Catarino, Helena Campos and Paulo Vasco: "On the Mersenne sequence", Annales Mathematicae et Informaticae, 46, pp.37-53, 2016

[ 戻る ]

• トップページ
• 京王
• 鉄道
• 短波
• プログラミング
• 路上観察
• 備忘録
• 談話室
• Annex
• 研究所概要
• サイトマップ
• 更新履歴