如何在c中实现定点数的取幂

编程问答 2022-05-07

问题描述

我需要将一个定点数提高到三次方和五次方，但标准的 pow 方法不起作用。这种情况怎么办

解决方法

（注意，当我回答这个问题时，它被标记为 C++。基本技术仍然适用于 C，但库解决方案不起作用。）

给定一个足够宽的整数类型，将底层整数提高必要的幂，并将其作为小数位数乘以相同幂的定点类型返回。

例如对于 x^3，如果输入 x 是一个 std::int32_t 整数，表示 15.16 定点数（即有 16 个小数位），那么结果将是值为 x*x*x 的 94 位宽整数，表示 45.48 定点数。而对于 x^5，结果将是一个 156 位宽的整数，值为 x*x*x*x*x，代表一个 75.80 的定点数。这些结果将是无损的。

不幸的是，即使使用扩展的基本整数类型 __int128_t，您也无法表示 156 位值。相反，您可以通过缩放转换回输入类型来跟踪每个乘法。这可能会失去精度或超出类型的范围。但是对于定点算法，用户有责任确保他们的类型足够宽。 （毕竟，定点只是整数算术的一般情况。）

以下是使用乘法/缩放交替 using raw integers 的 pow 函数示例：

// integer representing a signed 15.16 fixed-point number
using fp15_16 = std::int32_t;

constexpr auto to_fp{65536};
constexpr auto from_fp{1./65536};

fp15_16 pow_fixed(fp15_16 x,int y)
{
    assert(y>=0);
    return (y == 0) ? to_fp : ((std::int64_t{x} * pow_fixed(x,y-1)) * from_fp);
}

int main()
{
    // value 3.5 represented manually as fixed-point
    auto x{fp15_16(3.5 * to_fp)};
    std::cout << (pow_fixed(x,3) * from_fp) << '\n'; 
    std::cout << (pow_fixed(x,5) * from_fp) << '\n'; 
}

定点库可以管理扩展和转换 (example) 并留意溢出 (example)。

,

a^b = a * a ... * a

b 次。结果，做类似的事情

double prod = 1;
for (int power = 0; power < b; power++) prod *= a;

现在，您可以减少步骤数。如果 b 恰好是以 2 为底的对数的整数结果，那么您可以：

double prod = a;
for (int power = 1; power < log2b; power++) prod *= prod;

然而，这是一个不可靠的假设，所以假设您有 lo2b 作为小于 b 的 2 的最大整数幂。在这种情况下，您可以这样做：

double prod = a;
for (int power = 1; power < log2b; power++) prod *= prod;
for (int furtherPower = power; furtherPower < b; furtherPower++) prod *= a;

最后，你可以先计算 b 的最大整数对数的乘方，然后相应地减小 b 并对新的 b 重复这个过程，直到你完成操作。

有人可能会说计算对数的成本非常高。真的。然而，我们可以预先计算一组整数对数，并一次又一次地重复使用它来优化性能。

最后，数据类型具有的边界以及要计算的值高于给定边界存在问题。在这种情况下，需要使用一些支持更大值的数据结构。由于这不是本次讨论的一部分，我不会详细介绍。

编辑：尾随零的问题

如果您对小数的限制有疑问，那么您可以将您的数字相乘，得到整数，然后最后除以您需要的数字，以获得正确的数字。

,

对于整数指数，您可以通过平方来计算幂。但是，如果您的指数也是固定点，您必须使用不同的方法，例如 log,exp。

有关两者的更多信息，请参阅（这里有你所需要的）：

现在有了定点，您必须记住，您需要在每次操作后将数字校正回其比例：

n = bits_after_decimal_point
m = 1<<n
Integer(a) = Real(x*m)
Real(x) = Integer(a)/m
x*m + y*m = (x+y)*m         -> (a+b)
x*m - y*m = (x-y)*m         -> (a-b)
x*m * y*m = (x*y)*m*m       -> (a*b)>>n 
x*m / y*m = (x/y) + (x%y)/y -> (a/b)<<n + ((a%b)<<n)/b

所以 +,- 是相同的，但是 *,/ 需要将结果移位 n 位，并且还需要更多位来存储子结果，因此您需要

16bit*16bit = 32bit
32bit/16bit = 16bit
32bit%16bit = 16bit

解决此问题的最简单方法是在 *,/ 操作期间临时使用 32 位变量，或者使用通常本机涵盖此内容的 CPU 指令。

如果这不是一个选项，如您在上面看到的那样，可以部分地使用模直接在 16 位上完成除法（和/或迭代或重写到您自己的除法器），但是对于乘法，您可以使用朴素或 Karatsuba 乘法。在日志中，链接答案中的 exp 子链接是带有调试信息的 C++ 实现。

也看看这个：

它是我的 C++ 32 位 ALU 实现，我用作 bignum 算术的构建块。如果您注意 mul,div 操作，它们有两种变体，一种使用 x86 CPU 指令集，另一种使用本机 C++ 来完全满足您的需要...

如果你想使用数学中的浮点 pow（它会很慢）用 z=x^y 和 Q1.15 x,z 计算 int y 然后试试这个：

z = float(pow(float(x)/32768.0,y)*32768.0);

然而，这否定了使用定点数学的所有优点！！！

cfixed-pointmath