切り捨てられたテイラー級数のパデ近似に基づく近似自然対数関数を実装しました。精度は許容範囲内 (±0.000025) ですが、数回の最適化にもかかわらず、実行時間は標準ライブラリln関数の約 2.5 倍です! 速くなく、正確でなければ意味がありません。それにもかかわらず、Rust コードを最適化する方法を学ぶ方法としてこれを使用しています。criterion(私のタイミングはクレートの使用から来ています。私はブラックボックスを使用し、ループ内の値を合計し、結果から文字列を作成してオプティマイザーを無効にしました。)
Rust Playground では、私のコードは次のとおりです。
アルゴリズム
符号なし整数の比率で機能する私のアルゴリズムの概要:
- 値を超えない最大の 2 の累乗で除算することにより、区間 [1, 2] に範囲を縮小します。
- 分子の表記を変更 →
2ⁿ·N where 1 ≤ N ≤ 2 - 分母の表記変更→
2ᵈ·D where 1 ≤ D ≤ 2
- 分子の表記を変更 →
- これにより結果が得られます
log(numerator/denominator) = log(2ⁿ·N / 2ᵈ·D) = (n-d)·log(2) + log(N) - log(D) - log(N) を実行するために、テイラー級数は 0 付近では収束しませんが、1 付近では収束します...
- ... N は 1 に近いので、x = N - 1 に置き換えて、log(1 + x) を評価する必要があります。
- の置換を実行します
y = x/(2+x) - 関連する機能を検討する
f(y) = Log((1+y)/(1-y))= Log((1 + x/(2+x)) / (1 - x/(2+x)))= Log( (2+2x) / 2)= Log(1 + x)
- f(y) には、Log(1+x) の展開よりも速く収束する必要があるテイラー展開があります ...
- Log(1+x) の場合 →
x - x²/2 + x³/3 - y⁴/4 + ... - Log((1+y)/(1-y)) の場合 →
y + y³/3 + y⁵/5 + ...
- Log(1+x) の場合 →
- 切り捨てられた系列にパデ近似を使用する
y + y³/3 + y⁵/5 ... - ...どちらが
2y·(15 - 4y²)/(15 - 9y²) - 分母について繰り返し、結果を結合します。
パデ近似
コードのパデ近似部分は次のとおりです。
/// Approximate the natural logarithm of one plus a number in the range (0..1).
///
/// Use a Padé Approximation for the truncated Taylor series for Log((1+y)/(1-y)).
///
/// - x - must be a value between zero and one, inclusive.
#[inline]
fn log_1_plus_x(x : f64) -> f64 {
// This is private and its caller already checks for negatives, so no need to check again here.
// Also, though ln(1 + 0) == 0 is an easy case, it is not so much more likely to be the argument
// than other values, so no need for a special test.
let y = x / (2.0 + x);
let y_squared = y * y;
// Original Formula is this: 2y·(15 - 4y²)/(15 - 9y²)
// 2.0 * y * (15.0 - 4.0 * y_squared) / (15.0 - 9.0 * y_squared)
// Reduce multiplications: (8/9)y·(3.75 - y²)/((5/3) - y²)
0.8888888888888889 * y * (3.75 - y_squared) / (1.6666666666666667 - y_squared)
}
明らかに、そこまで高速化する必要はありません。
上位ビット
これまでで最も大きな影響を与えた変更は、最上位ビットの位置を取得する計算を最適化したことです。範囲を縮小するにはそれが必要です。
これが私のmsb機能です：
/// Provide `msb` method for numeric types to obtain the zero-based
/// position of the most significant bit set.
///
/// Algorithms used based on this article:
/// https://prismoskills.appspot.com/lessons/Bitwise_Operators/Find_position_of_MSB.jsp
pub trait MostSignificantBit {
/// Get the zero-based position of the most significant bit of an integer type.
/// If the number is zero, return zero.
///
/// ## Examples:
///
/// ```
/// use clusterphobia::clustering::msb::MostSignificantBit;
///
/// assert!(0_u64.msb() == 0);
/// assert!(1_u64.msb() == 0);
/// assert!(2_u64.msb() == 1);
/// assert!(3_u64.msb() == 1);
/// assert!(4_u64.msb() == 2);
/// assert!(255_u64.msb() == 7);
/// assert!(1023_u64.msb() == 9);
/// ```
fn msb(self) -> usize;
}
#[inline]
/// Return whether floor(log2(x))!=floor(log2(y))
/// with zero for false and 1 for true, because this came from C!
fn ld_neq(x : u64, y : u64) -> u64 {
let neq = (x^y) > (x&y);
if neq { 1 } else { 0 }
}
impl MostSignificantBit for u64 {
#[inline]
fn msb(self) -> usize {
/*
// SLOWER CODE THAT I REPLACED:
// Bisection guarantees performance of O(Log B) where B is number of bits in integer.
let mut high = 63_usize;
let mut low = 0_usize;
while (high - low) > 1
{
let mid = (high+low)/2;
let mask_high = (1 << high) - (1 << mid);
if (mask_high & self) != 0 { low = mid; }
else { high = mid; }
}
low
*/
// This algorithm found on pg 16 of "Matters Computational" at https://www.jjj.de/fxt/fxtbook.pdf
// It avoids most if-branches and has no looping.
// Using this instead of Bisection and looping shaved off 1/3 of the time.
const MU0 : u64 = 0x5555555555555555; // MU0 == ((-1UL)/3UL) == ...01010101_2
const MU1 : u64 = 0x3333333333333333; // MU1 == ((-1UL)/5UL) == ...00110011_2
const MU2 : u64 = 0x0f0f0f0f0f0f0f0f; // MU2 == ((-1UL)/17UL) == ...00001111_2
const MU3 : u64 = 0x00ff00ff00ff00ff; // MU3 == ((-1UL)/257UL) == (8 ones)
const MU4 : u64 = 0x0000ffff0000ffff; // MU4 == ((-1UL)/65537UL) == (16 ones)
const MU5 : u64 = 0x00000000ffffffff; // MU5 == ((-1UL)/4294967297UL) == (32 ones)
let r : u64 = ld_neq(self, self & MU0)
+ (ld_neq(self, self & MU1) << 1)
+ (ld_neq(self, self & MU2) << 2)
+ (ld_neq(self, self & MU3) << 3)
+ (ld_neq(self, self & MU4) << 4)
+ (ld_neq(self, self & MU5) << 5);
r as usize
}
}
さびた u64::next_power_of_two、安全でないコードと組み込み関数
これで、Rust には数値以上の最小の 2 乗を見つけるための高速な方法があることがわかりました。これが必要ですが、ビット位置も必要です。これは、数値の基数 2 の対数に相当するからです。(例: next_power_of_two(255) は 256 を生成しますが、8 番目のビットが設定されているため、8 が必要です。) のソース コードを見ると、next_power_of_twoというプライベート ヘルパー メソッド内に次の行がありますfn one_less_than_next_power_of_two。
let z = unsafe { intrinsics::ctlz_nonzero(p) };
同じ方法でビット位置を取得するために使用できる組み込み関数はありますか? 私がアクセスできるパブリックメソッドで使用されていますか? または、私が知らない組み込み関数を呼び出すための安全でないコードを作成する方法はありますか (これがほとんどです)。
私が呼び出すことができるそのようなメソッドまたは組み込み関数があれば、それが私のプログラムを大幅に高速化すると思いますが、おそらく他にも役立つことがあります。
アップデート：
頭叩き！63 - x.leading_zeros()最上位ビットの位置を見つけるために使用できます! 反対側から来るとは思いもしませんでした。これを試して、速度が上がるかどうかを確認します...