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

Question

我需要计算的数学表达式floor(ln(u)/ln(1-p))为0 < u < 1和0 < p < 1在Ç具有嵌入式处理器上没有浮点算术和无ln功能.结果是正整数.我知道极限情况(p = 0),我稍后会处理它们......

我想这个解决方案包括拥有u和p覆盖范围0..UINT16_MAX,并且在查找表中寻找对数,但我无法弄清楚究竟是什么:查找表映射到了什么？

结果不一定是100%精确,近似值可以.

谢谢!

Answer 1

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

我们可以使用定点算法来计算对数.我们通过将原始输入乘以接近1的递减幅度因子序列来实现这一点.考虑到序列中的每个因素,我们仅将输入乘以导致乘积接近1的那些因子,但不超过它.在这样做的同时,我们总结了log2()"适合"的因素.在这个过程的最后,我们得到一个非常接近1的数字作为我们的最终产品,以及代表二进制对数的总和.

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

B. de Lugish."一种自动评估数字计算机功能和计算的算法".博士论文,伊利诺伊大学计算机科学系,Urbana,1970.

JE Meggitt."伪划分和伪乘法过程".IBM Journal of Research and Development,Vol.1962年4月6日第2期,第210-226页

由于所选择的一组因子包括2 ⁱ和(1 + 2 ^-i),所以可以在不需要乘法指令的情况下执行必要的乘法:可以通过shift或shift加add来计算乘积.

由于输入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;
}