无需浮点运算或对数计算对数表达式

Question

我需要在 没有浮点运算且没有 @987654324 的嵌入式处理器上计算 C 中 0 < u < 1 和 0 < p < 1 的数学表达式 floor(ln(u)/ln(1-p)) @ 功能。结果是一个正整数。我知道极限情况（p=0），稍后我会处理它们......

我想该解决方案涉及使u 和p 的范围超过0..UINT16_MAX，并使用查找表来查找对数，但我无法弄清楚究竟如何：查找表映射到什么？

结果不必是 100% 准确的，近似值即可。

谢谢！

Answer 1

由于除数和除数都用对数，所以不用log()；我们可以改用log2()。由于对输入 u 和 p 的限制，已知对数都是负数，因此我们可以限制自己计算正数 -log2()。

我们可以使用定点算法来计算对数。我们通过将原始输入乘以接近 1 的一系列递减因子来做到这一点。考虑到序列中的每个因子，我们仅将输入乘以那些导致乘积接近 1 但不超过它的因子。这样做时，我们将“适合”的因素的log2() 相加。在这个过程的最后，我们得到一个非常接近 1 的数字作为我们的最终乘积，以及一个表示二进制对数的和。

这个过程在文献中被称为乘法归一化或伪除法，描述它的一些早期出版物是 De Lugish 和 Meggitt 的作品。后者表明原点基本上是亨利布里格斯计算常用对数的方法。

B.德卢吉什。 “一类用于自动评估数字计算机中的功能和计算的算法”。博士论文，计算机科学系，伊利诺伊大学厄巴纳分校，1970 年。

J. E. 美吉特。 “伪除法和伪乘法过程”。 IBM 研究与开发杂志，卷。 6，第 2 期，1962 年 4 月，第 210-226 页

由于所选的因子集包括 2ⁱ 和 (1+2^-i)，因此可以在不需要乘法指令的情况下执行必要的乘法：乘积可以通过移位或移位加加来计算。

由于输入 u 和 p 是 16 位的纯小数，我们可能希望为对数选择 5.16 定点结果。通过简单地将两个对数值相除，我们去除了定点比例因子，同时应用了floor() 运算，因为对于正数，floor(x) 与trunc(x) 相同，并且整数除法被截断。

请注意，对数的定点计算会导致输入接近 1 时产生较大的相对误差。这反过来意味着，如果 p 很小，使用定点算法计算的整个函数可能会提供与参考显着不同的结果.下面的测试用例就是一个例子：u=55af p=0052 res=848 ref=874。

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>

/* input x is a 0.16 fixed-point number in [0,1)
   function returns -log2(x) as a 5.16 fixed-point number in (0, 16]
*/   
uint32_t nlog2_16 (uint16_t x)
{
    uint32_t r = 0;
    uint32_t t, a = x;

    /* try factors 2**i with i = 8, 4, 2, 1 */
    if ((t = a << 8       ) < 0x10000) { a = t; r += 0x80000; }
    if ((t = a << 4       ) < 0x10000) { a = t; r += 0x40000; }
    if ((t = a << 2       ) < 0x10000) { a = t; r += 0x20000; }
    if ((t = a << 1       ) < 0x10000) { a = t; r += 0x10000; }
    /* try factors (1+2**(-i)) with i = 1, .., 16 */
    if ((t = a + (a >>  1)) < 0x10000) { a = t; r += 0x095c0; }
    if ((t = a + (a >>  2)) < 0x10000) { a = t; r += 0x0526a; }
    if ((t = a + (a >>  3)) < 0x10000) { a = t; r += 0x02b80; }
    if ((t = a + (a >>  4)) < 0x10000) { a = t; r += 0x01664; }
    if ((t = a + (a >>  5)) < 0x10000) { a = t; r += 0x00b5d; }
    if ((t = a + (a >>  6)) < 0x10000) { a = t; r += 0x005ba; }
    if ((t = a + (a >>  7)) < 0x10000) { a = t; r += 0x002e0; }
    if ((t = a + (a >>  8)) < 0x10000) { a = t; r += 0x00171; }
    if ((t = a + (a >>  9)) < 0x10000) { a = t; r += 0x000b8; }
    if ((t = a + (a >> 10)) < 0x10000) { a = t; r += 0x0005c; }
    if ((t = a + (a >> 11)) < 0x10000) { a = t; r += 0x0002e; }
    if ((t = a + (a >> 12)) < 0x10000) { a = t; r += 0x00017; }
    if ((t = a + (a >> 13)) < 0x10000) { a = t; r += 0x0000c; }
    if ((t = a + (a >> 14)) < 0x10000) { a = t; r += 0x00006; }
    if ((t = a + (a >> 15)) < 0x10000) { a = t; r += 0x00003; }
    if ((t = a + (a >> 16)) < 0x10000) { a = t; r += 0x00001; }
    return r;
}

/* Compute floor(log(u)/log(1-p)) for 0 < u < 1 and 0 < p < 1,
   where 'u' and 'p' are represented as 0.16 fixed-point numbers
   Result is an integer in range [0, 1048676]
*/
uint32_t func (uint16_t u, uint16_t p)
{
    uint16_t one_minus_p = 0x10000 - p; // 1.0 - p
    uint32_t log_u = nlog2_16 (u);
    uint32_t log_p = nlog2_16 (one_minus_p);
    uint32_t res = log_u / log_p;  // divide and floor in one go
    return res;
}

Answer 2

这个函数的最大值基本上取决于精度限制；也就是说，定点值可以任意接近极限(u -> 0) 或(1 - p -> 1)。

如果我们假设 (k) 小数位，例如，有以下限制：u = (2^-k) 和 1 - p = 1 - (2^-k)，
那么最大值为：k / (k - log2(2^k - 1))

（作为自然对数的比率，我们可以随意使用任何底，例如lb(x) 或log2）

---

与njuffa's 答案不同，我采用了查找表方法，使用k = 10 小数位来表示0 < frac(u) < 1024 和0 < frac(p) < 1024。这需要一个包含2^k 条目的日志表。使用 32 位表值，我们只查看 4KiB 表。

不仅如此，而且您正在使用足够的内存，您可以认真考虑使用“软浮点”库的相关部分。例如，k = 16 将产生 256KiB LUT。

我们正在为0 < i < 1024 计算值- log2(i / 1024.0)。由于这些值在开区间(0, k) 中，我们只需要 4 个二进制数字来存储整数部分。所以我们以 32 位 [4.28] 定点格式存储 预计算 LUT：

uint32_t lut[1024]; /* never use lut[0] */

for (uint32_t i = 1; i < 1024; i++)
    lut[i] = (uint32_t) (- (log2(i / 1024.0) * (268435456.0));

给定：u, p 由[0.10] 表示[1, 1023] 中的定点值：

uint32_t func (uint16_t u, uint16_t p)
{
    /* assert: 0 < u, p < 1024 */

    return lut[u] / lut[1024 - p];
}

我们可以根据“天真”浮点计算轻松测试所有有效的 (u, p) 对：

floor(log(u / 1024.0) / log(1.0 - p / 1024.0))

并且仅在以下情况下得到不匹配（+1 太高）：

u =  193, p =    1 :  1708 vs  1707 (1.7079978488147417e+03)
u =  250, p =  384 :     3 vs     2 (2.9999999999999996e+00)
u =  413, p =    4 :   232 vs   231 (2.3199989016957960e+02)
u =  603, p =    1 :   542 vs   541 (5.4199909906444600e+02)
u =  680, p =    1 :   419 vs   418 (4.1899938077226307e+02)

---

最后，事实证明，在[3.29] 定点格式中使用自然对数可以为我们提供更高的精度，其中：

lut[i] = (uint32_t) (- (log(i / 1024.0) * (536870912.0));

只产生一个“不匹配”，尽管'bignum'precision 表明它是正确的：

u =  250, p =  384 :     3 vs     2 (2.9999999999999996e+00)