optimization - 如何加快我的近似自然对数函数？

Question

我已经实现了一个基于 的近似自然对数函数。帕德近似 截断的泰勒级数。准确度可以接受(±0.000025)，但尽管经过多轮优化，它的执行时间仍然是标准库ln的2.5倍左右功能!如果它不更快并且不那么准确，那将毫无值(value)!尽管如此，我还是用它来学习如何优化我的 Rust 代码。 (我的时间来自于使用 criterion crate 。我使用了黑盒，对循环中的值求和，并从结果中创建了一个字符串来击败优化器。)

在 Rust Playground 上，我的代码是:

https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=94246553cd7cc0c7a540dcbeff3667b9

算法

我的算法概述，它适用于无符号整数的比率:

通过除以不超过值的 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)，泰勒级数不会在零附近收敛，但会在 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 + ...

对截断序列使用 Padé 近似 y + y³/3 + y⁵/5 ...

...这是2y·(15 - 4y²)/(15 - 9y²)

重复分母并合并结果。

Padé 逼近

这是代码的 Padé Approximation 部分:

/// 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
    }
}

Rust 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 的私有(private)辅助方法中看到这一行:

    let z = unsafe { intrinsics::ctlz_nonzero(p) };

那么有没有 内在 我可以用同样的方式获得位位置吗？它是否用于我可以访问的公共(public)方法中？或者有没有办法编写不安全的代码来调用一些我不知道的内在代码(其中大部分)？

如果有这样的方法或内在我可以调用，我怀疑这会大大加快我的程序，但也许还有其他的东西也会有所帮助。

更新:

打头!我可以使用 63 - x.leading_zeros()找到最高位的位置!我只是没想到从另一端过来。我会试试这个，看看它是否会加快速度......

Answer 1

一次优化可将时间缩短一半!
我重写了我的 msb(最高有效位)函数以使用库函数 u64::leading_zeroes内部使用内在函数:

    fn msb(self) -> usize {
        // THIRD ATTEMPT
        let z = self.leading_zeros();
        if z == 64 { 0 }
        else { 63 - z as usize }
    }

现在我的对数近似值只比固有的 ln 函数长 6%。我不太可能做得更好。
经验教训:使用内置日志!
更新:
以上适用于我的用例，但如果您想要一个在此范围内工作的通用 log(x) 函数:[1,e=2.718281828]，我为以后的项目想出了这个:
我最终决定使用三阶 Padé Approximant，它对于我的 x ∈ [1,2] 范围来说足够准确。在 Release模式下，它的速度大约是系统日志的两倍。

    fn pade_approximant_log(x: f64, log_base : f64) -> f64 {
        let x2 = x * x;
        let x3 = x2 * x;

        // 2nd order: 
        // let numerator = 3.0 * (x2 - 1.0);
        // let denominator = log_base * (x2 + 4.0*x + 1.0);

        // 3rd order: 
        let numerator = (11.0/27.0) * x3 + x2 - x - (11.0/27.0);
        let denominator = log_base * (x3/9.0 + x2 + x + (1.0/9.0));

        let ln_x = numerator / denominator;
        ln_x
    }

为了得到 Pade Approximant，我去了 Wolfram Alpha 并运行了这个:

PadeApproximant(ln[x], {x, 1, {3, 3}})

然后我重写了多项式，它给了我使用更少的乘法。
在上面，请注意，我让调用者将其数字基数的自然对数作为预先计算的值传递。我在所有调用中都使用相同的基数，因此这消除了在我的主循环中进行日志调用的需要。因此，以上适用于任何数字基数。