Issue with implementation of natural logarithm (ln) and exponentiation

Question

I try to follow the topic "Efficient implementation of natural logarithm (ln) and exponentiation" to be able to implement a log function without math.h. The described algorithm works well for values between 1 and 2 (normalized values). However, if the values are not normalized and I follow the normalization instructions, then I get wrong values.

Link: Click here

If I follow the code with an exemplary integer value of 12510, I get the following results:

y = 12510 (0x30DE), log2 = 13, divisor = 26, x = 481,1538

float ln(float y) {
    int log2;
    float divisor, x, result;

    log2 = msb((int)y); // See: https://stackoverflow.com/a/4970859/6630230
    divisor = (float)(1 << log2);
    x = y / divisor;    // normalized value between [1.0, 2.0]

    result = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
    result += ((float)log2) * 0.69314718; // ln(2) = 0.69314718

    return result;
}

The expected result for x should be a normalized value of 1 < x < 2. However, I fail in this calculation, because the received result is 481,1538.

Thanks in advance for any help

Code for the msb() function: unsigned int v; unsigned r = 0; while (v >>= 1) { r++; } return r; — Alexander Almendinger, Oct 28 '19 at 17:19

score 5 · Accepted Answer · answered Oct 26 '19 at 12:50

Out of curiosity, I tried to reproduce:

#include <stdio.h>

int msb(unsigned int v) {
  unsigned int r = 0;
  while (v >>= 1) r++;
  return r;
}

float ln(float y)
{
    int log2;
    float divisor, x, result;

    log2 = msb((int)y); // See: https://stackoverflow.com/a/4970859/6630230
    printf("log2: %d\n", log2);
    divisor = (float)(1 << log2);
    printf("divisor: %f\n", divisor);
    x = y / divisor;    // normalized value between [1.0, 2.0]
    printf("x: %f\n", x);
    result = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
    result += ((float)log2) * 0.69314718; // ln(2) = 0.69314718
    return result;
}

int main()
{
  printf("ln(12510): %f\n", ln(12510));
}

Output:

log2: 13
divisor: 8192.000000
x: 1.527100
ln(12510): 9.434252

Live Demo on coliru

I just tried this in my Windows 7 Pocket calculator and got:

9.434283603460956823997266847405

The first 5 digits are identical. – The rest I would consider as rounding issues whithout knowing which one got closer.

However, there is a typo (or mistake) in the question:

y = 12510 (0x30DE), log2 = 13, divisor = 26, x = 481,1538

divisor = (float)(1 << log2); with log2 = 13 yields 8192.

log2 << 1 would result in 26.

Just for fun, I changed the line into divisor = (float)(log2 << 1); and got the following output:

log2: 13
divisor: 26.000000
x: 481.153839
ln(12510): -2982522368.000000

So, this leaves me a bit puzzled:

The exposed code seems to be correct but the OP seems to have it interpreted (or resembled) wrong.

Thanks so much, I interpreted it in the way log2 <<1. Now its solved, thx! — Alexander Almendinger, Oct 28 '19 at 17:05
it works for a log of a number around 1 but for less than 1, the accuracy decreases dramatically. have you any suggestions for this case? — Ebay, Jan 02 '22 at 12:21

Issue with implementation of natural logarithm (ln) and exponentiation

1 Answers1