fast_polynomial

该crate实现了一种混合的Estrin/Horner方法，适合通过利用指令级并行性快速评估多项式。

关于融合乘加的重要说明

只有在您的二进制文件使用适当的Rust标志编译时，Rust才会使用FMA

RUSTFLAGS="-C target-feature=+fma"

或者

# .cargo/config.toml
[build]
rustflags = ["-C", "target-feature=+fma"]

否则将使用单独的乘法和加法操作。

动机

考虑以下简单的多项式评估函数

fn horners_method(x: f32, coefficients: &[f32]) -> f32 {
    let mut sum = 0.0;
    for coeff in coefficients.iter().rev().copied() {
        sum = x * sum + coeff;
    }
    sum
}

assert_eq!(horners_method(0.5, &[1.0, 0.3, 0.4, 1.6]), 1.45);

简单且干净，这是Horner方法。然而，请注意，每次迭代都依赖于前一次的结果，从而创建了一个无法并行化的依赖链，必须按顺序执行

vxorps      %xmm1,    %xmm1, %xmm1
vfmadd213ss 12(%rdx), %xmm0, %xmm1 /* Note the reuse of xmm1 for all vfmadd213ss */
vfmadd213ss 8(%rdx),  %xmm0, %xmm1
vfmadd213ss 4(%rdx),  %xmm0, %xmm1
vfmadd213ss (%rdx),   %xmm1, %xmm0

Estrin方案是一种组织多项式计算的方式，使得它们可以利用指令级并行性并行计算多项式的一部分。ILP是现代CPU可以同时排队多个计算，只要它们不相互依赖。

例如，(a + b) + (c + d)可能将此表达式中的每个括号内的部分分别使用不同的寄存器计算，同时进行。

此crate利用这一点，对所有至多15次幂的多项式进行操作，在这一点上，它切换到可以一次处理任意高次多项式（至多15个系数）的混合方法。

以下示例中，使用4个系数，通过poly_array将生成以下汇编代码

vmovss      4(%rdx),  %xmm1
vmovss      12(%rdx), %xmm3
vmulss      %xmm0,    %xmm0, %xmm2
vfmadd213ss 8(%rdx),  %xmm0, %xmm3
vfmadd213ss (%rdx),   %xmm0, %xmm1
vfmadd231ss %xmm3,    %xmm2, %xmm1

注意，它使用了多个xmm寄存器，因为前两个vfmadd213ss指令将并行运行。vmulss指令也可能与这些FMAs并行运行。尽管它们是更独立的指令，但由于它们在硬件上并行运行，这将显著提高速度。

有理多项式

fast_polynomial支持评估有理多项式，例如在Padé近似中找到的多项式，但有一个重要的注意事项：为了避免输入x的幂爆炸，我们执行了一种技术，将x替换为z = 1/x并在反向中有效地评估多项式

\begin{align}
    \frac{a_0 + a_1 x + a_2 x^2}{b_0 + b_1 x + b_2 x^2} &= \frac{a_0 + a_1 z^{-1} + a_2 z^{-2}}{b_0 + b_1 z^{-1} + b_2 z^{-2}} \\
        &= \frac{a_0 z^2 + a_1 z + a_2}{b_0 z^2 + b_1 z + b_2} \\
        &= \frac{a_2 + a_1 z + a_0 z^2}{b_2 + b_1 z + b_0 z^2} \\
\end{align}

然而，如果分子和分母的次数不同，则需要额外的校正步骤来调整次数以匹配，这可能会降低性能和准确性，因此应避免。填充多项式到相同的次数可能真正更快，尤其是在使用rational_array以避免过多的代码生成时。

其他缺点

Estrin的方案在非常高次的多项式中稍微数值上不太稳定。但是，使用FMA和提供的有理多项式评估例程都有可能在可能的情况下提高数值稳定性。

附加说明

对于固定次数的多项式，使用poly_array可以显著提高性能。在优化构建中，单态化的代码生成将几乎是理想的，并且可以避免不必要的分支。

但是，如果您需要评估具有相同X值的多个多项式，则存在polynomials模块，它提供直接固定次数的函数，允许重用X的次数直到15次方。

Cargo功能

默认的std和libm包功能会传递给num-traits。