Why my natural log function is so imprecise?

Question

Firstly, I'm using this approximation of a natural log. Or look here (4.1.27) for a better representation of formula.

Here's my implementation:

constexpr double eps = 1e-12;

constexpr double my_exp(const double& power)
{
    double numerator = 1;
    ull denominator = 1;
    size_t count = 1;
    double term = numerator / denominator;
    double sum = 0;
    while (count < 20)
    {
        sum += term;
        numerator *= power;
        #ifdef _DEBUG
            if (denominator > std::numeric_limits<ull>::max() / count)
                throw std::overflow_error("Denominator has overflown at count " + std::to_string(count));
        #endif // _DEBUG
        denominator *= count++;
        term = numerator / denominator;
    }
    return sum;
}

constexpr double E = my_exp(1);

constexpr double my_log(const double& num)
{
    if (num < 1)
        return my_log(num * E) - 1;
    else if (num > E)
        return my_log(num / E) + 1;
    else
    {
        double s = 0;
        size_t tmp_odd = 1;
        double tmp = (num - 1) / (num + 1);
        double mul = tmp * tmp;
        while (tmp >= eps)
        {
            s += tmp;
            tmp_odd += 2;
            tmp *= mul / tmp_odd;
        }
        return 2 * s;
    }
}

You probably can see why I want to implement these functions. Basically, I want to implement a pow function. But still my approach gives very imprecise answers, for example my_log(10) = 2.30256, but according to google (ln 10 ~ 2.30259).

my_exp() is very precise since it's taylor expansion is highly convergant. my_exp(1) = 2.718281828459, meanwhile e^1 = 2.71828182846 according to google. But unfortunately it's not the same case for natural log, and I don't even know how is this series for a natural log derived (I mean from the links above). And I couldn't find any source about this series.

Where's the precision errors coming from?

Answer 1

if (num < 1) return my_log(num * E) - 1; has an imprecision in the multiplication. Multiplying by 2 is more accurate. Of course, my_log(num) = my_log(2*num) - ln(2) so you'll need to change the 1.0 constant.

Yes, now you'll have a rounding error in -ln(2) instead of a rounding error in *E. That's typically less bad.

Also, you can save repeated rounding errors by first checking if (num<1/16) and then use my_log(num) = my_log(16*num) - ln(16). That's only a single rounding error.

As for the error in your core loop, I suspect the culprit is s += tmp;. This is a repeated addition. You can use Kahan summation there.

Answer 2

The line tmp *= mul / tmp_odd; means that each term is also being divided by the denominators of all previous terms, i.e. 1, 1*3, 1*3*5, 1*3*5*7, ... rather than 1, 3, 5, 7, ... as the formula states.

The numerator and denominator should therefore be computed independently:

double sum = 0;
double value = (num - 1) / (num + 1);
double mul = value * value;
size_t denom = 1;
double power = value;
double term = value;
while (term > eps)
{
    sum += term;
    power *= mul;
    denom += 2;
    term = power / denom;
}
return 2 * sum;

...

// Output for num = 1.5, eps = 1e-12
My func:   0.405465108108004513
Cmath log: 0.405465108108164385
             ------------

Much better!

Reducing the epsilon to 1e-18, we hit the accuracy limits of naïve summation:

// Output for num = 1.5, eps = 1e-18
My func:   0.40546510810816444
Cmath log: 0.405465108108164385
             ---------------

Kahan-Neumaier to the rescue:

double sum = 0;
double error = 0;
double value = (num - 1) / (num + 1);
double mul = value * value;
size_t denom = 1;
double power = value;
double term = value;
while (term > eps)
{
    double temp = sum + term;
    if (abs(sum) >= abs(term))
        error += (sum - temp) + term;
    else
        error += (term - temp) + sum;
    sum = temp;
    power *= mul;
    denom += 2;
    term = power / denom;
}
return 2 * (sum + error);

...

// Output for num = 1.5, eps = 1e-18
My func:   0.405465108108164385
Cmath log: 0.405465108108164385