アルゴリズム (algorithm) とは「計算可能」な数学的問題を解くための計算手順や処理手順のことです。ここでは、簡単な数値計算に関するアルゴリズムを集めてみました。ソースコードでは関数リテラルとしていますが、関数宣言(function文)でも同様に動作します。
数値計算以外の小技については "Tips!" を参照してください。
[ 戻る ]
実数 (real number) xの冪 (power) を求めるときはMath.pow関数を使います。ところが、Math.pow関数は小数の冪乗 (exponentiation) もサポートしているため演算速度が遅い関数 (function) です。
そこで、底 (base) がx、冪指数 (exponent) が自然数 (natural number) nである累乗 (power) 計算xnに繰り返し二乗法 (repeated squaring / exponentiation by squaring) を適用することで高速化を図ることができます。例えば、x11は x8+2+1=x8⋅x2⋅x1 と表されるため、x2=x⋅x、x4=x2⋅x2、x8=x4⋅x4 を順に求めて必要なものだけを乗じます。
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) とするxnの剰余演算 (modulo operation) xn mod m になります。
xnを直接計算してから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の値を入力すると、結果を即時に出力します。
[ 戻る ]
正数 (positive number) xを被開平数 (radicand) としたn乗根 (n-th root, 2≤n≤20) の整数部 (integer part) ⌊n√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) からなる一変数xのn次多項式
において、 のときの値
を求めるとき、直接計算する方法では 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=253以上になると下位桁に丸め誤差 (round-off error / rounding error) が生じ、結果に影響します。
[ 戻る ]
d-階定数係数線形漸化式 (d-order linear recurrence relations with constant coefficients) とは、
の形で表されるd項+定数項c0の漸化式 (recurrence relations) で、係数ci は全ての i について定数 (constant) となっています。c0=0のときは斉次 (homogeneous) 、c0≠0のときは非斉次 (non-homogeneous / inhomogeneous) と呼ばれます。
係数c0≠0, c1=1と初期値a0で初項 (first term) a0、公差 (common difference) c0の等差数列 (arithmetic sequence)、係数c0=0, c1≠0と初期値a0≠0で初項a0、公比 (common ratio) c1の等比数列 (geometric sequence) になります。定数係数線形漸化式で表される数列 (sequence) には以下のものがあります。
| 数列 | 漸化式 | 初期値 |
|---|---|---|
| フィボナッチ数 | Fn=Fn−1+Fn−2 | F0=0, F1=1 |
| リュカ数 | Ln=Ln−1+Ln−2 | L0=2, L1=1 |
| ペル数 | Pn=2Pn−1+Pn−2 | P0=0, P1=1 |
| ペル-リュカ数 | Qn=2Qn−1+Qn−2 | Q0=2, Q1=2 |
| NSW数 | Sn=2Sn−1+Sn−2 | S0=1, S1=1 |
| ヤコブスタール数 | Jn=Jn−1+2Jn−2 | J0=0, J1=1 |
| ヤコブスタール-リュカ数 | jn=jn−1+2jn−2 | j0=2, j1=1 |
| ペラン数 | Pn=Pn−2+Pn−3 | P0=3, P1=0, P2=2 |
| パドヴァン数 | Pn=Pn−2+Pn−3 | P0=1, P1=1, P2=1 |
| レオナルド数 | Ln=Ln−1+Ln−2+1 | L0=1, L1=1 |
| トリボナッチ数 | Tn=Tn−1+Tn−2+Tn−3 | T0=T1=0, T2=1 |
| テトラナッチ数 | Tn=Tn−1+Tn−2+Tn−3+Tn−4 | T0=T1=T2=0, T3=1 |
| 偶数 / 奇数 | Nn=Nn−1+2 | N0=0 / 1 |
| 6n±1 | Nn=Nn−2+6 | N0=1, N1=5 |
| ヒルベルト数 | Hn=Hn−1+4 Hn=2Hn−1−Hn−2 |
H0=1, H1=5 |
| mの累乗数 | pn=mpn−1 pn=(m+1)pn−1−mpn−2 |
p0=1, p0=m |
| メルセンヌ数 | Mn=2Mn−1+1 Mn=3Mn−1−2Mn−2 |
M0=0, M1=1 |
| 2n−1Mn | pn=6pn−1−8pn−2 | p0=0, p1=1 |
| カレン数 | Cn=4Cn−1−4Cn−2+1 Cn=5Cn−1−8Cn−2+4Cn−3 |
C0=1, C1=3, C2=9 |
| ウッダル数 | Wn=4Wn−1−4Wn−2−1 Wn=5Wn−1−8Wn−2+4Wn−3 |
W0=−1, W1=1, W2=7 |
| キャロル数 | Cn=6Cn−1−7Cn−2−6Cn−3+8Cn−4 | C0=−2, C1=−1, C2=7, C3=47 |
| キニア数 | Kn=6Kn−1−7Kn−2−6Kn−3+8Kn−4 | K0=2, K1=7, K2=23, K3=79 |
| 三角数 | Tn=2Tn−1−Tn−2+1 Tn=3Tn−1−3Tn−2+Tn−3 |
T0=0, T1=1, T2=3 |
| 怠け仕出し屋数 | ln=2ln−1−ln−2+1 ln=3ln−1−3ln−2+ln−3 |
l0=1, l1=2, l2=4 |
| 平方数 | sn=2sn−1−sn−2+2 sn=3sn−1−3sn−2+sn−3 |
s0=0, s1=1, s2=4 |
| 矩形数 | sn=2sn−1−sn−2+2 sn=3sn−1−3sn−2+sn−3 |
s0=0, s1=2, s2=6 |
| m角数 | Pn=2Pn−1−Pn−2+(m−2) Pn=3Pn−1−3Pn−2+Pn−3 |
P0=0, P1=1, P2=m |
| 中心つきm角数 | Pn=2Pn−1−Pn−2+m Pn=3Pn−1−3Pn−2+Pn−3 |
P0=1, P1=m+1, P2=3m+1 |
| 立方数 | cn=3cn−1−3cn−2+cn−3+6 cn=4cn−1−6cn−2+4cn−3−cn−4 |
c0=0, c1=1, c2=8, c3=27 |
| 三連続積数 | cn=3cn−1−3cn−2+cn−3+6 cn=4cn−1−6cn−2+4cn−3−cn−4 |
c0=0, c1=6, c2=24, c3=60 |
| ケーキ数 | cn=3cn−1−3cn−2+cn−3+1 cn=4cn−1−6cn−2+4cn−3−cn−4 |
c0=1, c1=2, c2=4, c3=8 |
| m角錐数 | Pn=3Pn−1−3Pn−2+Pn−3+(m−2) Pn=4Pn−1−6Pn−2+4Pn−3−Pn−4 |
P0=0, P1=1, P2=m+1, P3=4m−2 |
| 一般五角数 | pn=2pn−2−pn−4+3 pn=pn−1+2pn−2−2pn−3−pn−4+pn−5 |
p0=0, p1=1, p2=2, p3=5, p3=7 |
配列cに定数係数、配列aに数列の初期値を設定します。cには定数項c0が含まれるので要素数は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;
};係数はc0, c1, c2, …, cd 、初期値はa0, a1, …, ad−1の順にそれぞれカンマ若しくはスペースで区切って入力してください。Sampleで数列を選択すると係数と初期値が設定されます。n項目までの数列を表示しますが、nが漸化式の階数 (order) d以下のときは設定した初期値を全て表示します。内部演算で9007199254740992=253以上になると下位桁に丸め誤差が生じ、結果に影響します。
[ 戻る ]
ハミング数 (Hamming number) とは60の累乗を割り切ることができる自然数、即ち素因数 (prime factor) として2、3、5しか持たない数 =2i×3j×5k (i, j, k≥0) のことです。正則数 (regular number)、5-平滑数 (5-smooth number) とも呼ばれています。
配列qに1番目 (q0=1) からn番目 (qn−1) までのハミング数を代入していきます。i番目までの任意のハミング数をqj (0≤j≤i) とすると、i+1番目qi+1は2qj、3qj、5qjのうち、qiより大きい最小のものになります。これを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しか持たない数 =2i×3j×5k×7l (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) は下記の通りで、これに当てはまらない自然数がズッカーマン数になる可能性があります。
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;
};計算例を下表に示します。計算範囲が広いとその分時間がかかります。途中で打ち切るときは、ブラウザのタスクマネージャーから終了させてください。
| 計算範囲 | ズッカーマン数 | 数字積 |
|---|---|---|
| 8999111111111100〜8999261111111100 | 8999111281773312 8999112113183232 8999114219421696 8999118331121664 8999121233498112 8999126272346112 8999131122237312 8999132111111232 |
82301184 10077696 1088391168 120932352 120932352 282175488 17635968 419904 |
| 8999300000000000〜8999319787593872 | 8999311943817216 | 1269789696 |
| 8999400000000000〜8999479999999999 | 8999411241222144 8999413182111744 |
23887872 125411328 |
| 8999500000000000〜9007199254740992 | 8999621812813824 8999641242611712 8999723911211712 8999741132734464 8999761131922176 8999821431816192 |
1719926784 188116992 61725888 7900913664 1111065984 967458816 |
[ 戻る ]
分割数 (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) −⋯
ここで、和は一般五角数 (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) σ1(n)になります。nを素因数分解 (prime factorization) して求めることが多いです。
約数の総和の漸化式は前述の分割数の漸化式と同じですが、初期値σ1(0)がnになる点が異なります。そこで、σ1(i)を計算するときの初期値σ1(0)をiとすることで、分割数の計算と同様な手順で求めることもできます。しかし、時間計算量 (time complexity) がΟ(n1.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までの各自然数の約数の総和を列挙しています。
[ 戻る ]
[ 戻る ]
[ 戻る ]
