C 语言中 CORDIC 对数的问题

问题描述 投票：0回答：1

为了开始使用 CORDIC 

log10
，我实现了此 PDF，第 6 页中派生的算法：

#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;
}

这是论文中代码的直接 C99 实现，我使用 

ldexpf(x,n)
来计算 x·2n。该论文声称该方法对于 
0.4194 < x < 3.4627
是收敛的。该程序使用 
1 ≤ x ≤ 2
。下面打印完整的 C99 代码：

+--+--+--+-+-+++--++ 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

所以它按预期工作，除了

x = 1.125, 1.25
，错误很大并且在计算更多迭代时不会减少。我现在盯着那个代码几个小时，但找不到我缺少的东西......

完整的C99代码

#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;
}

作为参考，这里是论文中提出的算法：

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)
}
c numerical-methods c99
1个回答
0
投票

这不是一个完整的答案。这个问题引起了我的注意，我决定潜水并享受一些乐趣。

我已经使用Python测试了论文中给出的算法，并使用浮点数转换为定点，我的结果与你的结果相同。我还发现了一些其他有问题的点，例如

x = 1.60000
、
x = 2.00625
、
x = 2.03750
和
x = 2.06875
。

我将使用以下约定：

b_i+ = (1 + 2^-i)
b_i- = (1 - 2^-i)

算法可以这样写：

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

我不确定我们可以选择 

b_i
形式的 
(1 ± 2^-i)
的说法本身就是错误的，但我注意到 选择每个 
b_i
 的标准可能是错误的。例如，当
x == 1.125
时，当前序列是

x *= b_1- == 0.562500, z -= log10(b_1-) == 0.301030
x *= b_2+ == 0.703125, z -= log10(b_2+) == 0.204120
x *= b_3+ == 0.791016, z -= log10(b_3+) == 0.152967
x *= b_4+ == 0.840454, z -= log10(b_4+) == 0.126639
x *= b_5+ == 0.866718, z -= log10(b_5+) == 0.113275
x *= b_6+ == 0.880261, z -= log10(b_6+) == 0.106541
x *= b_7+ == 0.887138, z -= log10(b_7+) == 0.103161
x *= b_8+ == 0.890603, z -= log10(b_8+) == 0.101468
                        ...


x
永远不会达到
1
并且
z
的值是错误的。但是，如果我们从 
b_1+
（而不是 
b_1-
）开始，则顺序为：

x *= b_1+ == 1.687500; z -= log10(b_1+) == -0.176091
x *= b_2- == 1.265625; z -= log10(b_2-) == -0.051152
x *= b_3- == 1.107421; z -= log10(b_3-) ==  0.006839
x *= b_4- == 1.038208; z -= log10(b_4-) ==  0.034868
x *= b_5- == 1.005764; z -= log10(b_5-) ==  0.048656
x *= b_6- == 0.990048; z -= log10(b_6+) ==  0.055495
x *= b_7+ == 0.997783; z -= log10(b_7+) ==  0.052116
x *= b_8+ == 1.001681; z -= log10(b_8-) ==  0.050422
                        ...

收敛到正确的结果。

对于 

x == 2.00625
，现在算法使用

b_1- b_2- b_3+ b_4+ b_5+ b_6+ b_7+ b_8+, ...

和 

x
稳定在 
x == 0.957
左右。但如果我们改为 
b_2+
，则序列为 

b_1- b_2+ b_3- b_4- b_5- b_6- b_7- b_8-

这使得 

x == 1.0002
和 
z == 0.323317
（再次正确）。

我的直觉是，我们总是可以选择一系列

b_i
，其形式为
(1 ± 2^-i)
，满足
x * b_1 * ... * b_B == 1
，[1]，但使用

b_i = 1 + 2^-i if x * b1 * b2 * ... * b_i-1 < 1
b_i = 1 - 2^-i if x * b1 * b2 * ... * b_i-1 > 1

不是适当的选择标准。 

(x > 1.0)
的情况应该被一个更聪明的情况所取代，但到目前为止我还没有意识到。

[1]，因此使用查找表和简单的移位加法运算来计算对数，这就是该算法的目的。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.