计算没有浮点运算或日志的对数表达式

Question

我需要在嵌入式处理器 没有浮点运算，也没有ln函数。结果是一个正整数。我知道极限情况 (p=0)，我稍后会处理它们...

我想解决方案涉及 u 和 p 范围超过 0..UINT16_MAX，并呼吁查找 table 以获得对数，但我无法弄清楚究竟如何：查找 table 映射到什么？

结果无需 100% 准确，近似值即可。

谢谢！

Answer 1

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

我们可以使用定点运算来计算对数。为此，我们将原始输入乘以一系列递减幅度接近 1 的因子。考虑序列中的每个因子，我们仅将输入乘以导致乘积更接近 1 但不超过 1 的那些因子。在这样做的同时，我们对 "fit" 的因素求和 log2()。在此过程结束时，我们最终得到一个非常接近 1 的数字，以及一个表示二进制对数的总和。

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

乙。 de Lugish。 "A Class of Algorithms for Automatic Evaluation of Functions and Computations in a Digital Computer"。博士论文，计算机科学系，伊利诺伊大学厄巴纳分校，1970.

J. E. 梅吉特。 "Pseudo division and pseudo multiplication processes"。 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）

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

不止于此，并且您正在使用足够的内存，您可以认真考虑使用 'soft-float' 库的相关部分。例如，k = 16 会产生 256KiB LUT。

我们正在计算 0 < i < 1024 的值 - log2(i / 1024.0)。由于这些值在开区间(0, k)内，我们只需要4位二进制数来存储整数部分。所以我们以 32 位 [4.28] 定点格式存储 precomputed 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 由 [1, 1023] 中的 [0.10] 定点值表示：

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

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

我们可以根据 'naive' 浮点计算轻松测试所有有效的 (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));

仅生成一个 'mismatch'，但 'bignum' 精度表明它是正确的：

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

计算没有浮点运算或日志的对数表达式

Compute logarithmic expression without floating point arithmetics or log

c

math

embedded

logarithm

lookup-tables