計算對數表達而不浮點算術或日誌

Question

我需要計算用於0 < u < 1和0 < p < 1數學表達式floor(ln(u)/ln(1-p))在Ç的嵌入式處理器上與沒有浮點算術和沒有ln功能。結果是一個正整數。我知道了極限的情況下（P = 0），我會處理他們後來......

我想象的解決方案包括具有u和p射程超過0..UINT16_MAX，並呼籲查找表的對數，但我無法弄清楚：查找表映射到什麼地方？

結果不需要100％準確，近似值是可以的。

謝謝！

Answer 1

由於對數用於除數和除數，因此不需要使用log();我們可以用log2()代替。由於對輸入u和p的限制，已知對數都是負數，所以我們可以限制自己計算正數量-log2()。

我們可以使用定點算法來計算對數。我們這樣做是通過將原始輸入乘以一系列減小幅度的因子的方法來實現的。考慮到序列中的每個因子，我們只將輸入乘以那些導致產品接近1但不超過它的因子。在這樣做的時候，我們總結了「合適」因素的log2()。在這個過程結束時，我們得出一個非常接近1的數字作爲我們的最終產品，以及代表二進制對數的總和。

這個過程在文獻中被稱爲乘法規範化或僞劃分，描述它的一些早期出版物是De Lugish和Meggitt的着作。後者表明起源基本上是亨利布里格斯計算常用對數的方法。

B. de Lugish。 「一種數字計算機功能和計算自動評估算法」。博士論文，伊利諾伊大學計算機科學系，Urbana，1970.

J. E. Meggitt。 「僞劃分和僞乘法過程」。 IBM Journal of Research and Development，Vol。 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。

我們計算值- log2(i/1024.0)爲0 < i < 1024。由於這些值在開放區間(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)); 

考慮：通過[0.10]定點值表示u, p在[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'精密表明它是正確的：

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