std.numeric

このモジュールは、Alexander Stepanov による Standard Template Library の numeric ヘッダの発展途上の移植版 + いくつかの追加関数です。

Authors:
Andrei Alexandrescu, Don Clugston

enum CustomFloatFlags;

カスタムの浮動小数点型に関するフラグ。IEEE754 の型に関しては、すべてのフラグが適用されます。

storeNormalized: デフォルトで、値を正規化して保持する。つまり、仮数部は 1.nnnn の形であって、0.nnnn でないとします。正規化は IEE754 の型全てで有効です。浮動小数点型が (-1, 1) 範囲の数のみを格納するという定義だった場合、正規化は特に有用でないかもしれません。
allowDenorm: 非正規化数 (IEEE754 非正規化数のような) を許す。このフラグが ON ならば、指数部のビットが全てゼロの値は非正規化数として扱われます。
infinity: 無限大 (IEEE754の無限大のような) への対応。このフラグが ON ならば、指数部のビットが全て1で仮数部がゼロの値は、 (有効な)無限大値を指しているとして扱われます。符号は、もしあれば、符号ビットによって決まります。
nan: NaN (非数) IEEE754 非数のような) 対応。このフラグがONならば、指数部のビットが全て1で仮数部がゼロでない値は、 NaN 値として扱われます。
all: 上記すべてのオプションを含みます

struct CustomFloat(bool signBit,uint fractionBits,uint exponentBits,uint bias = (1u << exponentBits - 1) - 1,CustomFloatFlags flags = CustomFloatFlags.all);

カスタムの浮動小数点型をユーザー定義することを可能にします。これらの形式は値の記録にのみ使用できます; 演算を行うには、値を float か double に変換する必要があります。演算の後、代入演算子を使って元の形式に再度変換することは可能です。

Example:

// 16-bit 浮動小数点型を定義
alias CustomFloat!(1, 5, 10) HalfFloat;
auto x = HalfFloat(0.125);
assert(halfFloat.get!float == 0.125);

this(float input);: float から初期化
this(double input);: double から初期化
template opAssign(F) if (staticIndexOf!(Unqual!(F),float,double) >= 0): float または double からの代入
F get(F)();: float あるいは double として格納されている値を取り出します

alias FPTemporary(F) if (isFloatingPoint!(F))

最終的に F 型の結果が欲しいときに、途中結果として使用するともっとも高速な浮動小数点型を定義します。（F としては float, double, real を指定できます。）複数ステップの計算をするときには、途中結果を FPTemporary!F で保持するとよいでしょう。

Example:

// 配列の平均値を計算
double avg(in double[] a)
{
    if (a.length == 0) return 0;
    FPTemporary!double result = 0;
    foreach (e; a) result += e;
    return result / a.length;
}

FPTemporary は、実質的にほとんどのプロセッサに存在する、最適化された浮動小数点演算とレジスタの活用のために必要となります。上の例で加算をする場合、実際には、内部的には加算は real の精度で行われることがあります。その場合、中間結果を double format で保持するのは精度を落とすだけでなく、(驚くほどの) 速度低下を招きます。（これは、real から double への変換がループの各ステップで行われてしまうためです。）この lose-lose の状況を解決するために、F 以上の精度で計算できる 最速の 型として FPTemporary!F が定義されています。逆に 最も精度の高い 型をalias定義する必要はありません。それは常に real 型だからです。

最後に注意として、FPTemporary!F を使うのが常に最速とは一概には言えません。浮動小数点計算のスピードは、他の非常に多くの要因に左右されることに留意して下さい。

Num secantMethod(alias fun, Num)(Num xn_1, Num xn);

関数 f の根を点 [xn_1, x_n] (根に近いほどよい) から始めて探索する Secant法の実装です。Num には float, double, real を指定できます。

Example:

float f(float x) {
    return cos(x) - x*x*x;
}
auto x = secantMethod!(f)(0f, 1f);
assert(approxEqual(x, 0.865474));

T findRoot(T, R)(R delegate(T) f, T a, T b);

実数関数 f(x) の実数根を囲い込み法 ( bracketing ) によって探索します。

関数 f と、f(a) と f(b) の符号が異なるような区間 [a..b] が与えられると、その区間内で f(x) の根に最も近い値 x を返します。区間内に f(x) の根が複数ある場合、どれか一つだけが返されます。 f(x) が NaN を返した場合、この関数も NaN を返します。それ以外の場合は、このアルゴリズムが成功することは保証されています。

使用するアルゴリズムは TOMS748 を元にしていて、可能な限り逆三次補完 (inverse cubic interpolation) を使用し、それ以外の場合は放物線補間 (parabolic interpolation) またはセカント法 (secant interpolation) を使用します。 TOMS748 と比較すると、この実装は最悪ケースのパフォーマンスを 100倍以上改善し、典型的なケースのパフォーマンスを2倍改善します。 80-bit 実数では、ほとんどの場合 8 ～ 15 回の f(x) の呼び出しでマシン精度いっぱいの精度で根が得られます。最悪ケースでは、 bit数のおおよそ倍の回数の呼び出しが必要です。

References:
"On Enclosing Simple Roots of Nonlinear Equations", G. Alefeld, F.A. Potra, Yixun Shi, Mathematics of Computation 61, pp733-744 (1993). Fortran code available from www.netlib.org as algorithm TOMS478.

Tuple!(T,T,R,R) findRoot(T, R)(R delegate(T) f, T ax, T bx, R fax, R fbx, bool delegate(T lo, T hi) tolerance);

実数関数 f(x) の実数根を囲い込み法 (bracketing) で探索します。終了条件を指定することができます。

Parameters:

f	解析する関数
ax	f の根を含むことが保証されている初期区間の左端
bx	f の根を含むことが保証されている初期区間の右端
fax	f(ax) の値
fbx	f(bx) の値 (f(ax) と f(bx) は、多くの場合前もって知ることができます)
tolerance	探索打ち切り条件を定義します。現在の根の下限と上限を受け取ります。その区間を根として認めるときは true を返します。可能な最大精度が欲しい場合は、常に false を返すことで実現可能です。

Returns:
二つの区間からなるタプルを返します。最初の２つの要素は、根 x が存在する区間で、次の２つの要素は、その区間に対応する関数の値です。正確な根が見つかった場合、最初の２つの要素のどちらも根そのもので、次の２つの要素は 0 になっています。

CommonType!(ElementType!(Range1),ElementType!(Range2)) euclideanDistance(Range1, Range2)(Range1 a, Range2 b); CommonType!(ElementType!(Range1),ElementType!(Range2)) euclideanDistance(Range1, Range2, F)(Range1 a, Range2 b, F limit);

入力レンジ a と b の間のユークリッド距離を計算します。二つのレンジは同じ長さでなければなりません。三引数バージョンは、答えが指定された距離 limit 以上となることが確定したら、すぐに計算を停止します。（距離の近い場合だけが重要なケースで無駄な計算を省きたい際に有用です）

CommonType!(ElementType!(Range1),ElementType!(Range2)) dotProduct(Range1, Range2)(Range1 a, Range2 b); Unqual!(CommonType!(F1,F2)) dotProduct(F1, F2)(in F1[] avector, in F2[] bvector);

入力レンジ a と b の内積を計算します。二つのレンジは同じ長さでなければなりません。どちらのレンジも length プロパティを提供していれば、このチェックは一度だけ行われます; そうでなければ、要素を一つとりだすたびにチェックが入ります。

CommonType!(ElementType!(Range1),ElementType!(Range2)) cosineSimilarity(Range1, Range2)(Range1 a, Range2 b);

入力レンジ a と b のコサイン類似度を計算します。二つのレンジは同じ長さでなければなりません。どちらのレンジも length プロパティを提供していれば、このチェックは一度だけ行われます; そうでなければ、要素を一つとりだすたびにチェックが入ります。どちらかのレンジの要素が全て0ならば、0を返します。

bool normalize(R)(R range, ElementType!(R) sum = 1);

レンジ range の各要素に同じ数値を掛けて、和が sum となるようにします。range の和がゼロの場合、全体を sum / range.length という値にします。normalize は、range の値が全て正の時に意味を持ちます。特別なチェックは行いませんが、すべて正(か0)であることを仮定して実装されています。

Returns:
正規化が完了した場合 true。false は、 range の全要素がゼロであったか、range が空だった場合の返値です。

ElementType!(Range) entropy(Range)(Range r); ElementType!(Range) entropy(Range, F)(Range r, F max);

入力レンジ r のエントロピーをビット単位で計算します。この関数はr 内の値はすべて [0, 1] にあることを（チェック無しで）仮定します。エントロピーが意味を持つには、多くの場合、r を正規化（総和を1にすること）しておく必要があります。 2引数バージョンは、途中結果が max 以上になると直ちに計算を打ち切ります。

CommonType!(ElementType!(Range1),ElementType!(Range2)) kullbackLeiblerDivergence(Range1, Range2)(Range1 a, Range2 b);

レンジ a と b の間のカルバック・ライブラー情報量、つまり ai * log(ai / bi) の総和を計算します。log の底は 2 です。レンジ中の値は [0, 1] の範囲にあると仮定されます。通常、レンジとしては正規化された確率分布を表現したものを与えますが、これは必ずしも kullbackLeiblerDivergence の要求でもなくチェックもされません。もしどれかの要素 bi がゼロで対応する要素 ai がゼロでないなら、無限大を返します。 (そうでない場合、もし ai == 0 && bi == 0 ならば、項 ai * log(ai / bi) はゼロとして計算されます。）入力が正規化されていれば、結果は正です。

CommonType!(ElementType!(Range1),ElementType!(Range2)) jensenShannonDivergence(Range1, Range2)(Range1 a, Range2 b); CommonType!(ElementType!(Range1),ElementType!(Range2)) jensenShannonDivergence(Range1, Range2, F)(Range1 a, Range2 b, F limit);

レンジ a と b の間のジェンセン・シャノン情報量、つまり (ai * log(2 * ai / (ai + bi)) + bi * log(2 * bi / (ai + bi))) / 2 の総和を計算します。log の底は 2 です。レンジ中の値は [0, 1] の範囲にあると仮定されます。通常、レンジとしては正規化された確率分布を表現したものを与えますが、これは必ずしも jensenShannonDivergence の要求でもなくチェックもされません。入力が正規化されていれば、結果は [0, 1] の範囲に入ります。3引数バージョンは、途中結果が limit 以上になると直ちに計算を打ち切ります。

F gapWeightedSimilarity(alias comp = "a == b", R1, R2, F)(R1 s, R2 t, F lambda);

いわゆる "all-lengths gap-weighted string kernel" を計算します。これは、 s と t の間の、全ての長さの共通部分列に基づいて定義された類似度です。途中ギャップのある部分列も含まれます。

gapWeightedSimilarity(s, t, lambda) が何を計算するのか理解するために、まず lambda = 1、s = ["Hello", "brave", "new", "world"]、t = ["Hello", "new", "world"] の例を考えてみましょう。この場合、gapWeightedSimilarity は次のようなマッチを数えます:

長さ1のマッチが3つ: "Hello", "new", "world";
長さ2のマッチが3つ: ("Hello", "new"), ("Hello", "world"), ("new", "world");
長さ3のマッチが1つ: ("Hello", "new", "world").

gapWeightedSimilarity(s, t, 1) を呼び出すと、このマッチを単純に全て数え、足し上げて 7 を返します。

string[] s = ["Hello", "brave", "new", "world"];
string[] t = ["Hello", "new", "world"];
assert(gapWeightedSimilarity(s, t, 1) == 7);

マッチの際にギャップが無視されているのに気づいたでしょうか。例えば ("Hello", "new") は ("new", "world") と同じくらい良いマッチとして評価されています。これは、アプリケーションによっては過剰評価となることがあります。このようなギャップを完全に排除するには、lambda = 0 とします:

string[] s = ["Hello", "brave", "new", "world"];
string[] t = ["Hello", "new", "world"];
assert(gapWeightedSimilarity(s, t, 0) == 4);

上記の呼び出しでは、ギャップ込みのマッチ ("Hello", "new"), ("Hello", "world"), ("Hello", "new", "world") が除かれ、残りのマッチ数 4 が返っています。

もっとも興味深いのは、ギャップ込みのマッチを結果に含めはするけれど、ギャップ無しほど強くは評価しない、というケースです。結果は、問題の粒度にうまく合わせた文字列間の類似度の尺度として使用することができます。これは、lambda を 0 と 1 の間の値とすることで実現します: ギャップ込みのマッチは、ギャップの数に対して指数的にペナルティを負います （指数の底は lambda）。これはつまり、ギャップ無しのマッチは結果に 1 を加え、ギャップ1つのマッチは lambda を加え、………、 n 個のギャップ込みのマッチは pow(lambda, n) を結果に加えるということです。上の例では、4つのギャップ無しマッチと、2個の1ギャップマッチと、 1個の3ギャップマッチがありました。最後のケースは ("Hello", "world") で、一つ目の列に2つと二つ目の列に1つのギャップがあり、合わせて3ギャップとなります。全ての総和は、4 + 2 * lambda + pow(lambda, 3) です。

string[] s = ["Hello", "brave", "new", "world"];
string[] t = ["Hello", "new", "world"];
assert(gapWeightedSimilarity(s, t, 0.5) == 4 + 0.5 * 2 + 0.125);

gapWeightedSimilarity は、列同士の自然な類似度で、近似的なマッチも認めたいようケースにならいつでも有用です。上記の例では単語の列を使いましたが、要素どうしの同値性が比較できるような列（例として、文字や数値の列）なら適用できます。gapWeightedSimilarity は高度に最適化された動的計画法で実装されており、16 * min(s.length, t.length) バイトのメモリと Ο(s.length * t.length) の実行時間を必要とします。

Select!(isFloatingPoint!(F),F,double) gapWeightedSimilarityNormalized(alias comp = "a == b", R1, R2, F)(R1 s, R2 t, F lambda, F sSelfSim = F.init, F tSelfSim = F.init);

gapWeightedSimilarity による類似度の問題点は、二つの列が長くなるとそれだけで、全然類似していない列でも類似度が高くなってしまう点です。例えば、レンジ ["Hello", "world"] は ["Hello", "world", "world", "world",...] という文字列に "world" を増やせば増やすほど類似度が上がっていきます。これを防ぐために、gapWeightedSimilarityNormalized は類似度を正規化して返します。具体的には、 gapWeightedSimilarity(s, t, lambda) / sqrt(gapWeightedSimilarity(s, t, lambda) * gapWeightedSimilarity(s, t, lambda)) を計算します。関数 gapWeightedSimilarityNormalized (正規化カーネルと呼ばれます) の結果は [0, 1] の範囲に入り、全くどこでもマッチしない列どうしで 0、完全一致するレンジどうしで 1 となります。

Example:

string[] s = ["Hello", "brave", "new", "world"];
string[] t = ["Hello", "new", "world"];
assert(gapWeightedSimilarity(s, s, 1) == 15);
assert(gapWeightedSimilarity(t, t, 1) == 7);
assert(gapWeightedSimilarity(s, t, 1) == 7);
assert(gapWeightedSimilarityNormalized(s, t, 1) == 7. / sqrt(15. * 7));

オプション引数 sSelfSim と tSelfSim は、重複計算を避けるためのものです。多くの応用では、 gapWeightedSimilarity(s, s, lambda) と gapWeightedSimilarity(t, t, lambda) が計算済みであることがあります。この場合は、それらの値を sSelfSim と tSelfSim として渡します。

struct GapWeightedSimilarityIncremental(Range,F = double) if (isRandomAccessRange!(Range) && hasLength!(Range)); GapWeightedSimilarityIncremental!(R,F) gapWeightedSimilarityIncremental(R, F)(R r1, R r2, F penalty);

gapWeightedSimilarity と同じ目的の機能ですが、こちらは、長さ1のマッチ、(ギャップ込みの)長さ2のマッチ、……、とインクリメンタルにマッチを列挙していきます。必要なメモリは Ο(s.length * t.length) で、時間計算量は各ステップ Ο(s.length * t.length) です。先ほどの例を使うと:

string[] s = ["Hello", "brave", "new", "world"];
string[] t = ["Hello", "new", "world"];
auto simIter = gapWeightedSimilarityIncremental(s, t, 1);
assert(simIter.front == 3); // 長さ 1 のマッチが 3 つ
simIter.popFront;
assert(simIter.front == 3); // 長さ 2 のマッチが 3 つ
simIter.popFront;
assert(simIter.front == 1); // 長さ 3 のマッチが 1 つ
simIter.popFront;
assert(simIter.empty);     // もうマッチはない

実装は、Rousuらによる論文 "Efficient Computation of Gapped Substring Kernels on Large Alphabets" の図4の擬似コードに、アルゴリズムとシステム面双方からの最適化を加えたものとなっています。

this(Range s, Range t, F lambda);: 二つのレンジ s と t およびペナルティ lambda でオブジェクトを構築します。コンストラクタは、Ο(s.length * t.length) 時間で、長さ1の全てのマッチの個数を数えます。
GapWeightedSimilarityIncremental opSlice();: this を返します。
void popFront();: 次の長さのマッチの個数を数えます。Ο(s.length * t.length) 時間かかります。
F front();: 現在のマッチ長までの類似度 (初期値 1、popFront を呼ぶたびに増える) を返します。
bool empty();: さらにマッチがあるかどうかを返します。

T gcd(T)(T a, T b);

ユークリッドの互除法により、a と b の最小公倍数を計算します。