std.numeric
このモジュールは、Alexander Stepanov による Standard Template Library の numeric ヘッダの発展途上の移植版 + いくつかの追加関数です。 Source:std/numeric.d License:
Boost License 1.0. Authors:
Andrei Alexandrescu, Don Clugston, Robert Jacques
- CustomFloatに関するフラグ
- 符号付き数を格納可能とするために符号ビットを追加する。
- デフォルトで、値を正規化して保持する。つまり、仮数部は 1.nnnn の形であって、0.nnnn でないとします。 正規化は IEE754 の型全てで有効です。
- 指数部のビットが全てゼロの時 IEEE754 非正規化数 を格納可能とする。0 を表現するには必要です。
- IEEE754の無限大 を格納可能とする。
- IEEE754 非数 を格納可能とする。
- これをセットすると、指数の格納範囲を、最大値が max_exp = 1 となるようにします。 つまり、格納できる最大の値が 1.0 以上 2.0 未満となります。 指数の格納範囲が手動指定された場合は無視されます。
- セットすると、符号無し値は負の値として取り扱われます。
- セットすると IEEE754 非正規化数 として 0 だけを格納可能にします。 allowDenorm と storeNormalized と同時に指定する必要があります。
- IEEE754 オプションを全て有効にします。
- 以上の全てのオプションを指定しません。
- template CustomFloat(uint bits) if (bits == 8 || bits == 16 || bits == 32 || bits == 64 || bits == 80)
template CustomFloat(uint precision,uint exponentWidth,CustomFloatFlags flags = CustomFloatFlags.ieee) if (((flags & flags.signed) + precision + exponentWidth) % 8 == 0 && precision + exponentWidth > 0)
struct CustomFloat(uint precision,uint exponentWidth,CustomFloatFlags flags,uint bias) if (((flags & flags.signed) + precision + exponentWidth) % 8 == 0 && precision + exponentWidth > 0); - カスタムの浮動小数点型をユーザー定義することを可能にします。
これらの形式は値の記録にのみ使用できます; 演算を行うには、
値を real に暗黙変換する必要があります。
演算の後、代入演算子を使って元の形式に再度変換することは可能です。
Example:
// 16-bit 浮動小数点型を定義 CustomFloat!16 x; // ビット数指定 CustomFloat!(10, 5) y; // 仮数部と指数部のビット数指定 CustomFloat!(10, 5,CustomFloatFlags.ieee) z; // 仮数部と指数部のビット数指定に加えフォーマットフラグ CustomFloat!(10, 5,CustomFloatFlags.ieee, 15) w; // U数部と指数部のビット数指定、フォーマットフラグに加え、指数オフセット指定 // 16-bit 浮動小数点数をほぼ通常の数のように扱う w = x*y - 1; writeln(w); // 関数呼び出しには変換が必要 z = sin(+x) + cos(+y); // 単項プラス演算子はrealへの変換を行う簡単な方法です z = sin(x.re) + cos(y.re); // .re プロパティを使う方法もあります z = sin(x.get!float) + cos(y.get!float); // Or use get!T z = sin(cast(float)x) + cos(cast(float)y); // あるいは cast(T) で明示キャスト // 8-bit 浮動小数点型を確率保持用に定義 alias CustomFloat!(4, 4, CustomFloatFlags.ieee^CustomFloatFlags.probability^CustomFloatFlags.signed ) Probability; auto p = Probability(0.5);
- this(float input);
- float から初期化
- this(double input);
- double から初期化
- float または double からの代入
- float あるいは double として格納されている値を取り出します
- 最終的に 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 を使うのが常に最速とは一概には言えません。浮動小数点計算のスピードは、 他の非常に多くの要因に左右されることに留意して下さい。 - 関数 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));
- 実数関数 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. - 実数関数 f(x) の実数根を囲い込み法 (bracketing) で探索します。
終了条件を指定することができます。
Parameters:
Returns:f 解析する関数 ax f の根を含むことが保証されている初期区間 の左端 bx f の根を含むことが保証されている初期区間 の右端 fax f(ax) の値 fbx f(bx) の値 (f(ax) と f(bx) は、多くの場合 前もって知ることができます) tolerance 探索打ち切り条件を定義します。 現在の根の下限と上限を受け取ります。 その区間を根として認めるときは true を返します。 可能な最大精度が欲しい場合は、常に false を返すことで実現可能です。
二つの区間からなるタプルを返します。最初の２つの要素は、根 x が存在する区間で、 次の２つの要素は、その区間に対応する関数の値です。 正確な根が見つかった場合、最初の２つの要素のどちらも根そのもので、 次の２つの要素は 0 になっています。 - 入力レンジ a と b の間の ユークリッド距離 を計算します。二つのレンジは同じ長さでなければなりません。 三引数バージョンは、答えが指定された距離 limit 以上となることが確定したら、すぐに計算を停止します。 （距離の近い場合だけが重要なケースで無駄な計算を省きたい際に有用です）
- 入力レンジ a と b の 内積 を計算します。二つのレンジは同じ長さでなければなりません。 どちらのレンジも length プロパティを提供していれば、このチェックは一度だけ行われます; そうでなければ、要素を一つとりだすたびにチェックが入ります。
- 入力レンジ a と b の コサイン類似度 を計算します。二つのレンジは同じ長さでなければなりません。 どちらのレンジも length プロパティを提供していれば、このチェックは一度だけ行われます; そうでなければ、要素を一つとりだすたびにチェックが入ります。どちらかのレンジの要素が全て0ならば、0を返します。
- レンジ range の各要素に同じ数値を掛けて、
和が sum となるようにします。range の和がゼロの場合、全体を sum / range.length
という値にします。normalize は、range
の値が全て正の時に意味を持ちます。特別なチェックは行いませんが、すべて正(か0)であることを仮定して実装されています。
Returns:
正規化が完了した場合 true。false は、 range の全要素がゼロであったか、range が空だった場合の返値です。 - 入力レンジ r の エントロピー をビット単位で計算します。 この関数はr 内の値はすべて [0, 1] にあることを（チェック無しで）仮定します。エントロピーが意味を持つには、多くの場合、r を正規化（総和を1にすること）しておく必要があります。 2引数バージョンは、 途中結果が max 以上になると直ちに計算を打ち切ります。
- レンジ a と b の間の カルバック・ライブラー情報量、 つまり ai * log(ai / bi) の総和を計算します。log の底は 2 です。 レンジ中の値は [0, 1] の範囲にあると仮定されます。通常、レンジとしては正規化された確率分布を表現したものを与えますが、 これは必ずしも kullbackLeiblerDivergence の要求でもなくチェックもされません。もしどれかの要素 bi がゼロで対応する要素 ai がゼロでないなら、無限大を返します。 (そうでない場合、もし ai == 0 && bi == 0 ならば、項 ai * log(ai / bi) はゼロとして計算されます。） 入力が正規化されていれば、結果は正です。
- レンジ a と b の間の ジェンセン・シャノン情報量 、つまり (ai * log(2 * ai / (ai + bi)) + bi * log(2 * bi / (ai + bi))) / 2 の総和を計算します。log の底は 2 です。 レンジ中の値は [0, 1] の範囲にあると仮定されます。通常、レンジとしては正規化された確率分布を表現したものを与えますが、 これは必ずしも jensenShannonDivergence の要求でもなくチェックもされません。入力が正規化されていれば、 結果は [0, 1] の範囲に入ります。3引数バージョンは、 途中結果が limit 以上になると直ちに計算を打ち切ります。
- いわゆる "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").
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) の実行時間を必要とします。 - 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 として渡します。 - 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の全てのマッチの個数を数えます。
- this を返します。
- 次の長さのマッチの個数を数えます。Ο(s.length * t.length) 時間かかります。
- 現在のマッチ長までの類似度 (初期値 1、popFront を呼ぶたびに増える) を返します。
- さらにマッチがあるかどうかを返します。
- ユークリッドの互除法により、a と b の最小公倍数を計算します。
- 2のベキ乗サイズでフーリエ変換を行うクラスです。
このクラスは、同一サイズのFFTを行う際に再利用できる大量のデータをカプセル化しています。
この効果で、
同じサイズのFFTを何度も繰り返すときはフリー関数APIで行うより効率があがります。
しかし、１回のFFTをするためてであれば、
フリー関数APIの方が便利に使うことができます。
References:
en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm- this(size_t size);
- size サイズの高速フーリエ変換を行うための Fft オブジェクトを構築します。 size は2のベキ乗である必要があります。
- Ο(N log N)
の Cooley-Tukey アルゴリズムでフーリエ変換を計算します。range
には、スライシング可能なランダムアクセスレンジで、長さがコンストラクタに渡した size
と等しいものを指定できます。レンジの要素には、
数値型（実数値と解釈される）か、
プロパティ .re と .im を持つ複素数型のいずれかを使用できます。
Returns:
データを周波数領域へと変換した結果の複素数の配列 - 上の同名の関数と同じ機能ですが、結果はユーザの指定したバッファに格納されます。 バッファは入力レンジと同じ長さで、ランダムアクセス可能で、 スライシング可能で、 複素数ライクな型を要素に持たなければなりません。複素数ライクとは、浮動小数点数を .re と .im というメンバ/プロパティに読み書きできることをいいます。
- レンジを逆フーリエ変換した結果を計算します。
レンジはスライシング可能なランダムアクセスレンジで、長さがコンストラクタに渡した size
と等しいものを指定できます。レンジの要素には、
std.complex.Complex 型か、
それと本質的に同じコンパイル時インターフェイスをもった型が使えます。
Returns:
時間領域の信号 - ユーザの指定したバッファに逆FFTの結果を格納します。 スライシング可能なランダムアクセスレンジで、 複素数ライクな型を格納できるレンジが指定できます。
- Fft オブジェクトを作り、FFT や逆FFTを実行し、
結果を返す便利関数です。単発でFFTを実行するのに便利です。
Note:
便利さのだけでなく、１回だけの仕様であれば、 手動でFftオブジェクトを作るのに比べてわずかに効率的でもあります。 理由は、この関数の終了時に決定的にFftオブジェクトを破棄できるためです。