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...
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?
The result needs not be 100% exact, approximations are OK.
Thanks!
Since the logarithm is used in both dividend and divisor, there is no need to use
log()
; we can uselog2()
instead. Due to the restrictions on the inputsu
andp
the logarithms are known to be both negative, so we can restrict ourselves to compute the positive quantity-log2()
.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.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. de Lugish. "A Class of Algorithms for Automatic Evaluation of Functions and Computations in a Digital Computer". PhD thesis, Dept. of Computer Science, University of Illinois, Urbana, 1970.
J. E. Meggitt. "Pseudo division and pseudo multiplication processes". IBM Journal of Research and Development, Vol. 6, No. 2, April 1962, pp. 210-226
As the chosen set of factors comprises 2i 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.
Since the inputs
u
andp
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 afloor()
operation at the same time, because for positive numbers,floor(x)
is identical totrunc(x)
and integer division is truncating.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
.The maximum value of this function basically depends on the precision limit; that is, how arbitrarily close to the limits
(u -> 0)
or(1 - p -> 1)
the fixed point values can be.If we assume
(k)
fractional bits, e.g., with the limits:u = (2^-k)
and1 - p = 1 - (2^-k)
,then the maximum value is:
k / (k - log2(2^k - 1))
(As the ratio of natural logarithms, we are free to use any base e.g.,
lb(x)
orlog2
)Unlike njuffa's answer, I went with a lookup table approach, settling on
k = 10
fractional bits to represent0 < frac(u) < 1024
and0 < frac(p) < 1024
. This requires a log table with2^k
entries. Using 32-bit table values, we're only looking at a4KiB
table.Any more than that, and you are using enough memory that you could seriously consider using the relevant parts of a 'soft-float' library. e.g.,
k = 16
would yield a256KiB
LUT.We're computing the values
- log2(i / 1024.0)
for0 < i < 1024
. Since these values are in the open interval(0, k)
, we only need 4 binary digits to store the integral part. So we store the precomputed LUT in 32-bit[4.28]
fixed-point format:Given:
u, p
represented by[0.10]
fixed-point values in[1, 1023]
:We can easily test all valid
(u, p)
pairs against the 'naive' floating-point evaluation:floor(log(u / 1024.0) / log(1.0 - p / 1024.0))
and only get a mismatch (+1 too high) on the following cases:
Finally, it turns out that using the natural logarithm in a
[3.29]
fixed-point format gives us even higher precision, where:only yields a single 'mismatch', though 'bignum' precision suggests it's correct: