我怎样才能加快我的近似自然对数函数？
How can I speed up my approximate natural log function?
我已经实现了基于截断泰勒级数的 Padé 近似 的近似自然对数函数。精度是可以接受的（±0.000025），但尽管经过几轮优化，它的执行时间仍然是标准库 ln 函数的 2.5 倍左右！如果它不是更快并且不是那么准确，那它就毫无价值！尽管如此，我还是用它来学习如何优化我的 Rust 代码。 （我的计时来自使用 criterion 箱子。我使用黑盒，对循环中的值求和并根据结果创建一个字符串来击败优化器。）
在 Rust Playground 上，我的代码是：
算法
我的算法概述，它适用于无符号整数的比率：
- 范围缩小到区间[1, 2]除以不超过值的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)，泰勒级数不收敛于零附近，但收敛于一附近...
- ...因为 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é 近似部分：
/// 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。 （例如：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) };
那么是否有一个 intrinsic 可以用来以相同的方式获取位位置？它是否用于我有权访问的 public 方法？或者有没有一种方法可以编写不安全的代码来调用一些我不知道的内在函数（这是其中的大部分）？
如果有这样的方法或内在函数我可以调用，我怀疑这会大大加快我的程序，但也许还有其他的东西也会有帮助。
更新：
砸脑袋！我可以使用 63 - x.leading_zeros() 找到最高有效位的位置！我只是没想到来自另一端。我会试试这个，看看它是否加快速度...
一次优化将时间减半！
我重写了我的 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) 函数来处理这个 运行ge: [1,e=2.718281828]，我想出了这个一个后来的项目：
我最终决定使用三阶 Padé Approximant，它对我的 运行ge of 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}})
然后我重写了它给我的多项式以使用更少的乘法。
在上面，请注意，我让调用者将其基数的自然对数作为预先计算的值传递。我在我的所有调用中都使用相同的基数，因此这消除了在我的主循环中进行日志调用的需要。因此，以上适用于任何数字基数。
我已经实现了基于截断泰勒级数的 Padé 近似 的近似自然对数函数。精度是可以接受的（±0.000025），但尽管经过几轮优化，它的执行时间仍然是标准库 ln 函数的 2.5 倍左右！如果它不是更快并且不是那么准确，那它就毫无价值！尽管如此，我还是用它来学习如何优化我的 Rust 代码。 （我的计时来自使用 criterion 箱子。我使用黑盒，对循环中的值求和并根据结果创建一个字符串来击败优化器。）
在 Rust Playground 上，我的代码是：
算法
我的算法概述，它适用于无符号整数的比率：
- 范围缩小到区间[1, 2]除以不超过值的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)，泰勒级数不收敛于零附近，但收敛于一附近...
- ...因为 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) →
- 对截断序列使用 Padé 近似
y + y³/3 + y⁵/5 ... - ... 即
2y·(15 - 4y²)/(15 - 9y²) - 重复计算分母并合并结果。
帕德近似
这是代码的 Padé 近似部分：
/// 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。 （例如：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) };
那么是否有一个 intrinsic 可以用来以相同的方式获取位位置？它是否用于我有权访问的 public 方法？或者有没有一种方法可以编写不安全的代码来调用一些我不知道的内在函数（这是其中的大部分）？
如果有这样的方法或内在函数我可以调用，我怀疑这会大大加快我的程序，但也许还有其他的东西也会有帮助。
更新：
砸脑袋！我可以使用 63 - x.leading_zeros() 找到最高有效位的位置！我只是没想到来自另一端。我会试试这个，看看它是否加快速度...
一次优化将时间减半！
我重写了我的 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) 函数来处理这个 运行ge: [1,e=2.718281828]，我想出了这个一个后来的项目：
我最终决定使用三阶 Padé Approximant，它对我的 运行ge of 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}})
然后我重写了它给我的多项式以使用更少的乘法。 在上面，请注意，我让调用者将其基数的自然对数作为预先计算的值传递。我在我的所有调用中都使用相同的基数，因此这消除了在我的主循环中进行日志调用的需要。因此，以上适用于任何数字基数。