C语言中对数的CORDIC问题

91zkwejq  于 2023-10-16  发布在  其他
关注(0)|答案(3)|浏览(106)

为了开始使用log10的CORDIC,我实现了在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;
}

这是论文中代码的直接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)
}
23c0lvtd

23c0lvtd1#

这不是一个完整的答案。这个问题吸引了我的眼球,我决定潜水,找点乐子。
我已经在Python中使用浮点数和转换为定点测试了论文中给出的算法,我的结果和你的一样。我还发现了其他一些问题点,例如x = 1.60000x = 2.00625x = 2.03750x = 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)
}

我不确定我们可以选择(1 ± 2^-i)形式的b_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永远不会到达1z的值是错误的。然而,如果我们从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.0002z == 0.323317(再次正确)。
我的直觉[1]是,我们总是可以选择一系列(1 ± 2^-i)b_i s,满足x * b_1 * ... * b_B == 1,[2]但使用

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]结果我的直觉又错了!)。见OP的回答。
[2]因此,使用查找表和简单的移位和加法运算来计算最小值,这是该算法的目的。

oogrdqng

oogrdqng2#

算法是假的。问题是不是所有的数字都可以用
x =(1 ± 2−k)
看看这个discussion on MSE。甚至不存在一个区间[a,B],使得该区间内的所有x都可以用这种形式表示。
但 * 确实存在 * 的是形式的表示
x =(1 + ak·2−k)
对于所有1 ≤ x ≤ 2.38 = x 0,ak在{ 0,1 }中,这足以修复算法:
为了计算logB(x),首先将参数x归一化,使得1 /x 0 ≤ x ≤ 1。
然后将以下算法应用于约化的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

实际上,N是固定的,N个logB值来自查找表。除此之外,该方法还需要加法、比较和可变偏移量的移位。
比较IF xk ≤ xmax使得最终的绝对误差(几乎)关于0对称。对于IF xk ≤ 1,绝对误差将是不对称的,并且是两倍高。

结果

作为参考,下面是C99中的一个实现,它增加了固定算法实际工作的信心。它是B = 2,N = 10,图中是16385个均匀分布在[0.5,1]上的x值。
(点击放大)

C99编码

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

编辑:

刚发现我发明了轮子...根据维基百科,这是费曼在曼哈顿计划中已经使用的方法。

yzuktlbb

yzuktlbb3#

我不太愿意称之为CORDIC计算。CORDIC在双曲矢量化模式下可以计算tanh-1,并且可以通过一些预处理和后处理来计算双曲余弦。这种对数计算与二进制手写除法紧密相关,因此该算法被称为伪除法。一般的方法可以追溯到亨利布里格斯如何计算他的表的。
在二进制除法中,可以使用各种数字集,例如在恢复或“非执行”变体中的数字集{0,1},或者如果使用非恢复变体,则为{-1,1}(我在这里使用-1代替带上划线的1)。问题中链接的文档的作者似乎错误地假设非恢复除法算法将直接转换为对应的伪除法,在乘法归一化过程中使用数字集{-1,1}来表示因子(1-2-i)和(1+2-i)。然而,正如其他答案所指出的那样,这并不起作用,因为不是所有的论点都可以通过这种方式规范化为统一。
通过伪除法计算除法的最简单方法是在“非执行”算法中使用数字集{0,1},正如我在对this question的回答中所演示的那样。也可以使用对应于归一化因子(1+2-i),1,(1-2-i)的数字集{-1,0,1},如De Lugish在他的博士论文中所示。他为乘法、除法、平方根和各种超越函数的按位计算设计了共享硬件:
布鲁斯基因德Lugish,* 一类算法自动评估某些基本功能的二进制计算机 *,报告编号。399,计算机科学系,伊利诺伊大学香槟分校,1970年6月。
下面是直接基于本报告描述的log()的ISO-C99实现。由于我知道De Lugish的工作,但直到现在才使用它,所以我编写代码是为了自己的利益,而不是试图最大化性能。

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

// ln(1+2*(-i))
float tabp [24];
// ln(1-2*(-i))
float tabm [24];

float logf_delugish (float a)
{
    const float three_eigth = 0.375f;
    const float ln2 = 0.69314718056f;
    int sk, expo;
    float sum = tabm[0], x = frexpf (a, &expo), ex2 = 1.0f;
    x = 2.0f * x;
    for (int k = 0; k < 24; k++) {
        sk = 0;
        if ((x - 1.0f) <  (-three_eigth * ex2)) sk = +1;
        if ((x - 1.0f) >= (+three_eigth * ex2)) sk = -1;
        ex2 *= 0.5f;
        if (sk == 1) {
            x = x + x * ex2;
            sum = sum - tabp [k];
        } 
        if (sk == -1) {
            x = x - x * ex2;
            sum = sum - tabm [k];
        }
    }
    return expo * ln2 + sum;
}

float tabp [24] =
{ 
    0.40546510810816f,
    0.22314355131421f,
    0.11778303565638f,
    0.06062462181643f,
    0.03077165866675f,
    0.01550418653597f,
    0.00778214044205f,
    0.00389864041566f,
    0.00195122013126f,
    0.00097608597306f,
    0.00048816207950f,
    0.00024411082753f,
    0.00012206286253f,
    0.00006103329368f,
    0.00003051711247f,
    0.00001525867265f,
    0.00000762936543f,
    0.00000381468999f,
    0.00000190734681f,
    0.00000095367386f,
    0.00000047683704f,
    0.00000023841855f,
    0.00000011920928f,
    0.00000005960464f,
};

// ln(1-2*(-i))
float tabm [24] =
{
    -0.69314718055995f,
    -0.28768207245178f,
    -0.13353139262452f,
    -0.06453852113757f,
    -0.03174869831458f,
    -0.01574835696814f,
    -0.00784317746103f,
    -0.00391389932114f,
    -0.00195503483580f,
    -0.00097703964783f,
    -0.00048840049811f,
    -0.00024417043217f,
    -0.00012207776369f,
    -0.00006103701897f,
    -0.00003051804380f,
    -0.00001525890548f,
    -0.00000762942364f,
    -0.00000381470454f,
    -0.00000190735045f,
    -0.00000095367477f,
    -0.00000047683727f,
    -0.00000023841861f,
    -0.00000011920930f,
    -0.00000005960465f,
};

#define LIMIT (256)

int main (void)
{
    unsigned int count = 0;
    float sumerrsq = 0;
    float maxerr = 0, maxerr_loc = INFINITY;
    float a = 1.0f / LIMIT;
    do {
        float res = logf_delugish (a);
        double ref = log ((double)a);
        float err = fabsf (res - ref);
        sumerrsq += err * err;
        if (err > maxerr) {
            maxerr = err;
            maxerr_loc = a;
        }
        a = nextafterf (a, INFINITY);
        count++;
    } while (a < LIMIT);
    printf ("maxerr = %15.8e @ %15.8e  RMS err = %15.8e\n", 
            maxerr, maxerr_loc, sqrt(sumerrsq / count));
    return EXIT_SUCCESS;
}

相关问题