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

Question

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

我想解决方案涉及 u和 p范围超过 0..UINT16_MAX ，并求助于对数的查找表，但我无法弄清楚到底是怎样的:查找表映射到什么？

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

谢谢!

Answer 1

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

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

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

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

J。 E. 梅 git 。 “伪除法和伪乘法过程”。 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;
}