Problems with CORDIC for Logarithm in C

Question

In order to get started with CORDIC for log10, I implemented the algorithm derived in this PDF, pp6:

#include <stdio.h>
#include <math.h>

// https://www.mikrocontroller.net/attachment/31117/cordic1.pdf
float logf_cordic (float x, int B)
{
    float z = 0.0f;

    for (int i = 1; i <= B; i++)
    {
        if (x > 1.0f)
        {
            printf ("-");
            x = x - ldexpf (x, -i);
            z = z - log10f (1.0f - ldexpf (1.0f, -i));
        }
        else
        {
            printf ("+");
            x = x + ldexpf (x, -i);
            z = z - log10f (1.0f + ldexpf (1.0f, -i));
        }
    }
    return z;
}

This is a straight forward C99 implementation of the code from the paper, where I used ldexpf(x,n) to compute x·2ⁿ. The paper claims that the method converges for 0.4194 < x < 3.4627. The program uses 1 ≤ x ≤ 2. The complete C99 code below prints:

+--+--+--+-+-+++--++ x = 1.00000: y = -0.000000, dy = -1.487272e-07
-+++++++++++++++++++ x = 1.12500: y = +0.099773, dy = +4.862081e-02
-+++++++++++++++++++ x = 1.25000: y = +0.099773, dy = +2.863325e-03
-+++-+--+-+--+-++--+ x = 1.37500: y = +0.138302, dy = -4.023314e-07
-++-+--++----+++--+- x = 1.50000: y = +0.176091, dy = -2.831221e-07
-+-+++++-++-++-++++- x = 1.62500: y = +0.210854, dy = +2.831221e-07
-+-+-+-+++++--+---++ x = 1.75000: y = +0.243038, dy = +2.235174e-07
-+--++-+--+---+---+- x = 1.87500: y = +0.273001, dy = +0.000000e+00
-+---+--+++--------- x = 2.00000: y = +0.301030, dy = -5.960464e-08

So it works as expected, except for x = 1.125, 1.25 where the error is big and doesn't decrease when computing with more iterations. I am now staring at that code for hours, but cannot find what I am missing...

---

• Code in Compiler Explorer (gcc v11)

Complete C99 code

#include <stdio.h>
#include <math.h>

float logf_cordic (float x, int B)
{
    float z = 0.0f;

    for (int i = 1; i <= B; i++)
    {
        if (x > 1.0f)
        {
            printf ("-");
            x = x - ldexpf (x, -i);
            z = z - log10f (1.0f - ldexpf (1.0f, -i));
        }
        else
        {
            printf ("+");
            x = x + ldexpf (x, -i);
            z = z - log10f (1.0f + ldexpf (1.0f, -i));
        }
    }
    return z;
}

int main (int argc, char *argv[])
{
    int ex = 3;
    int B = 20;

    if (argc >= 2)
        sscanf (argv[1], "%i", &ex);

    if (argc >= 3)
        sscanf (argv[2], "%i", &B);

    if (ex < 0) ex = 0;
    if (ex > 16) ex = 16;
    if (B > 100) B = 100;

    int n = 1 << ex;
    float dx = 1.0f / n;
    for (int i = 0; i <= n; ++i)
    {
        float x = 1.0f + i * dx;
        float y = logf_cordic (x, B);
        float dy = y - log10f (x);
        printf (" x = %.5f: y = %+f, dy = %+e\n",
                (double) x, (double) y, (double) dy);
    }
    return 0;
}

For reference, here is the algo as presented in the paper:

log10(x){
    z = 0;
    for ( i=1;i=<B;i++ ){
        if (x > 1)
            x = x - x*2^(-i);
            z = z - log10(1-2^(-i));
        else
            x = x + x*2^(-i);
            z = z - log10(1+2^(-i));
    }
return(z)
}

Answer 1

The algorithm is bogus. The problem is that not all numbers can be presented in the form

x = Π (1 ± 2^−k)

See this discussion on MSE. There doesn't even exist an interval [a,b] such that all x's in that interval can be represented in that form.

---

What does exist though is a representation of the form

x = Π (1 + a_k·2^−k)

with a_k in { 0, 1 } for all 1 ≤ x ≤ 2.38 = x₀, and this is enough to fix the algorithm:

To compute log_B(x), first normalize the argument x such that 1 / x₀ ≤ x ≤ 1.

Then apply the following algorithm to the reduced x:

logB_cordic (x, N)
    z = 0
    xmax = 1 + 2^{−N−1}
    FOR k = 1 ... N
        xk = x + x·2^{-k}
        IF xk ≤ xmax
            x = xk;
            z = z - logB (1 + 2^{−k})
    return z

In practice, N is fixed and the N logB values come from a lookup-table. Apart from that, the method requires additions, comparisions, and shifts by variable offsets.

The comparison IF xk ≤ xmax makes the final absolute error (almost) symmetric around 0. With IF xk ≤ 1 the absolute error would be asymmetric and twice as high.

Results

For reference, below is an implementation in C99 that adds some confidence that the fixed algorithm is actually working. It's B = 2, N = 10, and the plot is for 16385 values of x evenly spaced over [0.5, 1].
(click to enlarge)

C99 Code

#include <math.h>

double log2_cordic (double x, int N)
{
    double xmax = 1.0 + ldexp (1.0, -N - 1);
    double z = 0.0;

    for (int k = 1; k <= N; k++)
    {
        double x_bk = x + ldexp (x, -k);

        if (x_bk <= xmax)
        {
            x = x_bk;
            z = z - log2 (1.0 + ldexp (1.0, -k));
        }
    }
    return z;
}

---

Edit:

Just found out that I reinvented the wheel... According to Wikipedia, it's a method Feynman already used during the Manhattan Project.