快速定点 pow、log、exp 和 sqrt

Posted 2023-02-18

【中文标题】快速定点 pow、log、exp 和 sqrt

【英文标题】：Fast fixed point pow, log, exp and sqrt

【发布时间】：2011-06-07 03:31:35

【问题描述】：

我有一个定点类 (10.22)，我需要一个 pow、一个 sqrt、一个 exp 和一个 log 函数。

唉，我什至不知道从哪里开始。谁能给我一些有用文章的链接，或者更好的是，给我一些代码？

我假设一旦我有了一个 exp 函数，那么实现 pow 和 sqrt 就变得相对容易了。

pow( x, y ) => exp( y * log( x ) )
sqrt( x )   => pow( x, 0.5 )

我觉得困难的只是那些 exp 和 log 函数（好像我记得我的一些日志规则，但我记不起关于它们的更多内容）。

据推测，sqrt 和 pow 也会有一种更快的方法，因此即使它只是说使用我上面概述的方法，也可以理解这方面的任何指针。

请注意：这必须是跨平台和纯 C/C++ 代码，所以我不能使用任何汇编程序优化。

【问题讨论】：

如果你想要 fast 功能，那些 exp( y * log( x ) ) 实现不会削减它。

@MSalters：可能是真的，尤其是 sqrt ......但是，真的可以在 pow 上做得更好吗？

是的，正如我在回答中指出的那样。分解 y 的整数部分。

【参考方案1】：

下面是 Clay S. Turner 的定点对数基 2 算法 [1] 的示例 C 实现。该算法不需要任何类型的查找表。这在内存限制严格且处理器缺少 FPU 的系统上很有用，例如许多微控制器的情况。以 e 为底的对数和以 10 为底的对数也可以通过使用对数的属性来支持，对于任何以 n 为底的：

          logₘ(x)
logₙ(x) = ───────
          logₘ(n)

其中，对于该算法，m 等于 2。

这个实现的一个很好的特性是它支持可变精度：精度可以在运行时确定，但以范围为代价。按照我实现它的方式，处理器（或编译器）必须能够进行 64 位数学运算以保存一些中间结果。它可以轻松适应，不需要 64 位支持，但范围会缩小。

使用这些函数时，x 预计是一个定点值，根据 指定precision。例如，如果 precision 是 16，那么 x 应该按 2^16 (65536) 缩放。结果是一个与输入具有相同比例因子的定点值。返回值INT32_MIN 表示负无穷大。返回值INT32_MAX表示错误，errno将设置为EINVAL，表示输入精度无效。

#include <errno.h>
#include <stddef.h>

#include "log2fix.h"

#define INV_LOG2_E_Q1DOT31  UINT64_C(0x58b90bfc) // Inverse log base 2 of e
#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10

int32_t log2fix (uint32_t x, size_t precision)

    int32_t b = 1U << (precision - 1);
    int32_t y = 0;

    if (precision < 1 || precision > 31) 
        errno = EINVAL;
        return INT32_MAX; // indicates an error
    

    if (x == 0) 
        return INT32_MIN; // represents negative infinity
    

    while (x < 1U << precision) 
        x <<= 1;
        y -= 1U << precision;
    

    while (x >= 2U << precision) 
        x >>= 1;
        y += 1U << precision;
    

    uint64_t z = x;

    for (size_t i = 0; i < precision; i++) 
        z = z * z >> precision;
        if (z >= 2U << (uint64_t)precision) 
            z >>= 1;
            y += b;
        
        b >>= 1;
    

    return y;


int32_t logfix (uint32_t x, size_t precision)

    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31;

    return t >> 31;


int32_t log10fix (uint32_t x, size_t precision)

    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31;

    return t >> 31;

此实现的代码也位于Github，以及一个示例/测试程序，该程序说明如何使用此函数计算和显示从标准输入读取的数字的对数。

[1] C. S. Turner，"A Fast Binary Logarithm Algorithm"，IEEE 信号处理杂志，第 124,140 页，2010 年 9 月。

【讨论】：

“精度”到底是什么意思？这是用于小数部分的位数吗？ IE。精度 = 10 意味着 int32_t 变量被解释为具有 1 个符号位、21 位整数部分和 10 位小数部分的浮点数。对吗？

@Joerg 是的，除了没有符号位（输入值 x 是无符号的，因为对于负值未定义实数对数）。所以对于精度 10，有 22 个整数位和 10 个小数位。

@DanMoulding 是否可以使用这种技术来计算具有固定点的 2 的幂？我对此提出了另一个问题：***.com/questions/61471447/…

感谢您的参考。这是一个非常漂亮的算法，并且由于其简单性而易于移植。

【参考方案2】：

一个非常简单的解决方案是使用合适的表驱动近似值。如果您正确减少输入，您实际上并不需要大量数据。 exp(a)==exp(a/2)*exp(a/2)，这意味着你真的只需要为1 < x < 2计算exp(x)。在该范围内，runga-kutta 近似将给出合理的结果，大约 16 个条目 IIRC。

同样，sqrt(a) == 2 * sqrt(a/4) == sqrt(4*a) / 2 这意味着您只需要1 < a < 4 的表条目。 Log(a) 有点难：log(a) == 1 + log(a/e)。这是一个相当慢的迭代，但 log(1024) 只有 6.9，所以你不会有很多迭代。

您可以对 pow 使用类似的“整数优先”算法：pow(x,y)==pow(x, floor(y)) * pow(x, frac(y))。这是因为pow(double, int) 是微不足道的（分而治之）。

[编辑] 对于log(a) 的整数部分，存储一个表1, e, e^2, e^3, e^4, e^5, e^6, e^7 可能很有用，因此您可以通过在该表中对a 进行简单的硬编码二进制搜索来减少log(a) == n + log(a/e^n)。从 7 步到 3 步的改进不是很大，但这意味着您只需将 e^n 除以一次，而不是 n 乘以 e。

[编辑 2] 对于最后一个log(a/e^n) 术语，您可以使用log(a/e^n) = log((a/e^n)^8)/8 - 每次迭代按表查找 产生3 个更多位。这使您的代码和表格大小保持较小。这通常是嵌入式系统的代码，它们没有大缓存。

[编辑 3] 这对我来说仍然不聪明。 log(a) = log(2) + log(a/2)。您可以只存储定点值log2=0.30102999566，计算前导零的数量，将a 移动到用于查找表的范围内，然后将该移动（整数）乘以定点常量log2。可以低至 3 条指令。

使用e 进行缩减步骤只会给你一个“不错的”log(e)=1.0 常量，但这是错误的优化。 0.30102999566 和 1.0 一样好；两者都是 10.22 定点的 32 位常量。使用 2 作为范围缩小的常数允许您使用位移位进行除法。

您仍然可以从编辑 2 中获得诀窍，log(a/2^n) = log((a/2^n)^8)/8。基本上，这会给你一个结果(a + b/8 + c/64 + d/512) * 0.30102999566 - b,c,d 在 [0,7] 范围内。 a.bcd 真的是一个八进制数。毫不奇怪，因为我们使用 8 作为电源。 （这个技巧同样适用于 2、4 或 16 次方。）

[编辑 4] 仍然有一个开放的结局。 pow(x, frac(y) 只是 pow(sqrt(x), 2 * frac(y))，我们有一个不错的 1/sqrt(x)。这为我们提供了更有效的方法。说frac(y)=0.101二进制，即1/2加1/8。那么这意味着x^0.101 是(x^1/2 * x^1/8)。但是x^1/2 只是sqrt(x) 而x^1/8 是(sqrt(sqrt(sqrt(x)))。再保存一个操作，Newton-Raphson NR(x) 给我们1/sqrt(x) 所以我们计算1.0/(NR(x)*NR((NR(NR(x)))。我们只反转最终结果，不要直接使用 sqrt 函数。

【讨论】：

对于 exp 和 log，您的方法是可以的（除了我会使用 Taylor 或 Pade 扩展在 1 附近，并使用 -0.5 和 0.5 之间的参数作为 exp，以及 1 和 2 作为 log）。对于 sqrt，这可能是矫枉过正：牛顿法似乎相当合适（你必须用牛顿法计算 1 / sqrt(x)，只有乘法）

顺便说一句，我已经将 sqrt 实现为 newton raphson 迭代。性能很好，只需要几个步骤就可以比我的 10.22 固定的更精确......

你是怎么做 pow(x, frac(y)) 的？

@AdamTegen：可能是exp(frac(y)*log(x))，使用上面的优化。由于frac(y) < 1 和log(x) 无论如何都不会很大，因此您不需要多次迭代exp(a)==exp(a/2)*exp(a/2)。我也可以考虑=pow(sqrt(x), 2*frac(y)。

【参考方案3】：

一个好的起点是Jack Crenshaw's book, "Math Toolkit for Real-Time Programming"。它对各种超越函数的算法和实现进行了很好的讨论。

【讨论】：

链接似乎已损坏。

@DougMcClean：谢谢 - Kindle 格式的它似乎不再可用 - 我现在更新了指向平装本的链接。

【参考方案4】：

仅使用整数运算检查我的定点 sqrt 实现。 发明很有趣。现在已经很老了。

https://groups.google.com/forum/?hl=fr%05aacf5997b615c37&fromgroups#!topic/comp.lang.c/IpwKbw0MAxw/discussion

否则检查CORDIC 算法集。这就是实现您列出的所有函数和三角函数的方法。

编辑：我在 GitHub here

上发布了经过审核的源代码