我有一个定点类（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++ 代码，所以我不能使用任何汇编程序优化。
一个非常简单的解决方案是使用一个像样的表驱动近似。如果您正确减少输入，您实际上并不需要大量数据。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- 每次迭代通过 table lookup再产生 3 位。这使您的代码和表格大小保持较小。这通常是嵌入式系统的代码，它们没有大缓存。
[编辑 3] 这对我来说仍然不聪明。log(a) = log(2) + log(a/2). 您可以只存储定点值log2=0.6931471805599，计算前导零的数量，转移a到用于查找表的范围，然后将该转移（整数）乘以定点常数log2。可以低至 3 条指令。
使用e减少步骤只会给你一个“不错”log(e)=1.0的常数，但这是错误的优化。0.6931471805599 和 1.0 一样好；两者都是 10.22 定点的 32 位常量。使用 2 作为范围缩小的常数允许您使用位移位进行除法。
[编辑 5] 由于您将其存储在 Q10.22 中，因此您可以更好地存储 log(65536)=11.09035488。(16 x 对数 (2))。“x16”意味着我们还有 4 位可用的精度。
您仍然可以从编辑 2 中获得诀窍，log(a/2^n) = log((a/2^n)^8)/8. 基本上，这会给你一个结果(a + b/8 + c/64 + d/512) * 0.6931471805599- 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)。But x^1/2is justsqrt(x)和x^1/8is (sqrt(sqrt(sqrt(x)))。保存一个操作，Newton-RaphsonNR(x)给我们1/sqrt(x)，所以我们计算1.0/(NR(x)*NR((NR(NR(x))). 我们只反转最终结果，不要直接使用 sqrt 函数。
下面是 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 ] CS Turner，“快速二进制对数算法”，IEEE 信号处理 Mag。，第 124,140 页，2010 年 9 月。
Jack Crenshaw 的书“实时编程的数学工具包”是一个很好的起点。它对各种超越函数的算法和实现进行了很好的讨论。