快速定点pow,log,exp和sqrt

Question

我有一个固定点类(10.22),我需要一个pow,一个sqrt,一个exp和一个日志函数.

唉,我不知道从哪里开始.任何人都可以提供一些有用的文章的链接,或者,还没有给我提供一些代码？

我认为,一旦我有了exp函数,那么实现pow和sqrt变得相对容易.

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

它只是那些我发现很难的exp和日志函数(好像我记得我的一些日志规则,我记不起其他的很多了).

据推测,顺便说一句,对于sqrt和pow也会有一个更快的方法,所以即使只是说使用上面概述的方法:)前面的任何指针都会受到赞赏:)

请注意:这是跨平台和纯C/C++代码,所以我不能使用任何汇编程序优化.

Answer 1

一个非常简单的解决方案是使用适当的表驱动近似.如果正确减少输入,实际上并不需要大量数据.exp(a)==exp(a/2)*exp(a/2),这意味着你真的只需要计算exp(x)的1 < x < 2.在该范围内,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这样您就可以log(a) == n + log(a/e^n)通过对该表中的a进行简单的硬编码二进制搜索来减少.从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函数.

Answer 2

下面是Clay S. Turner的定点对数库2算法的示例C实现[ 1 ].该算法不需要任何类型的查找表.这对于内存限制紧张且处理器缺少FPU的系统非常有用,例如许多微控制器就是这种情况.然后还使用对数属性支持日志库e和日志库10,对于任何基数n:

          log?(x)
log?(x) = ???????
          log?(n)

其中,对于此算法,y等于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 ] CS Turner,"快速二进制对数算法",IEEE信号处理Mag.,第124,140页,2010年9月.

Answer 3

杰克·克伦肖(Jack Crenshaw)的着作"实时编程的数学工具包"(Math Toolkit for Real-Time Programming)是一个很好的起点.它对各种超越函数的算法和实现进行了很好的讨论.