计算对数表达式而无需浮点算术或对数 [英] Compute logarithmic expression without floating point arithmetics or log

问题描述

我需要计算 0<的数学表达式 floor（ln（u）/ ln（1-p）） u< 1 和 0< p ＜ 1 在具有无浮点算术且没有 ln 函数的嵌入式处理器上的 C 中。结果是一个正整数。我知道极限情况（p = 0），稍后再处理...

I need to compute the mathematical expression floor(ln(u)/ln(1-p)) for 0 < u < 1 and 0 < p < 1 in C on an embedded processor with no floating point arithmetics and no ln function. The result is a positive integer. I know about the limit cases (p=0), I'll deal with them later...

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

I imagine that the solution involves having u and p range over 0..UINT16_MAX, and appeal to a lookup table for the logarithm, but I cannot figure out how exactly: what does the lookup table map to?

结果不必是100％精确的，近似值就可以。

The result needs not be 100% exact, approximations are OK.

谢谢！

推荐答案

由于对数用于除数和除数，无需使用 log（）;我们可以改用 log2（）。由于输入 u 和 p 的限制，对数都是负数，因此我们可以限制自己计算正数 -log2（）。

Since the logarithm is used in both dividend and divisor, there is no need to use log(); we can use log2() instead. Due to the restrictions on the inputs u and p the logarithms are known to be both negative, so we can restrict ourselves to compute the positive quantity -log2().

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

We can use fixed-point arithmetic to compute the logarithm. We do so by multiplying the original input by a sequence of factors of decreasing magnitude that approach 1. Considering each of the factor in sequence, we multiply the input only by those factors that result in a product closer to 1, but without exceeding it. While doing so, we sum the log2() of the factors that "fit". At the end of this procedure we wind up with a number very close to 1 as our final product, and a sum that represents the binary logarithm.

文献将其称为乘法归一化或伪除法，一些早期的出版物将其描述为De Lugish和Meggitt的作品。后者表明原点基本上是亨利·布里格斯（Henry Briggs）计算常用对数的方法。

This process is known in the literature as multiplicative normalization or pseudo division, and some early publications describing it are the works by De Lugish and Meggitt. The latter indicates that the origin is basically Henry Briggs's method for computing common logarithms.

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

J。 E.梅吉特。 伪除法和伪乘法过程。 IBM研究与发展杂志，第1卷。 1962年4月，第2号，第6卷，第210-226页

由于所选因子集包括2 ⁱ 和（1 + 2 ^-i ）可以执行必要的乘法运算，而无需乘法指令：可以通过shift或shift plus add来计算乘积。

As the chosen set of factors comprises 2ⁱ and (1+2^-i) the necessary multiplications can be performed without the need for a multiplication instruction: the products can be computed by either shift or shift plus add.

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

Since the inputs u and p are purely fractional numbers with 16 bits, we may want to chose a 5.16 fixed-point result for the logarithm. By simply dividing the two logarithm values, we remove the fixed-point scale factor, and apply a floor() operation at the same time, because for positive numbers, floor(x) is identical to trunc(x) and integer division is truncating.

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

Note that the fixed-point computation of the logarithm results in large relative error for inputs near 1. This in turn means the entire function computed using fixed-point arithmetic may deliver results significantly different from the reference if p is small. An example of this is the following test case: 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;
}

这篇关于计算对数表达式而无需浮点算术或对数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！