CORDIC Arcsine实现失败

Question

最近，我实现了一个CORDIC函数库，以减少所需的计算能力(我的项目基于PowerPC，其执行时间规范非常严格)。语言是ANSI-C。

其他函数(sin/cos/atan)在32位实现和64位实现中的精度限制范围内工作。

不幸的是，asin()函数对于某些输入有系统地失败。

为了测试目的，我实现了一个.h文件，用于simulink S-函数中。(这只是为了方便起见，您可以将以下内容编译为一个独立的.exe，只需进行最小的更改)

注意:我强制执行了32次迭代，因为我的工作是32位精度，并且需要尽可能高的精度。

科迪奇：

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

#define FLOAT32 float
#define INT32 signed long int
#define BIT_XOR ^

#define CORDIC_1K_32 0x26DD3B6A                  
#define MUL_32       1073741824.0F     /*needed to scale float -> int*/            
#define INV_MUL_32   9.313225746E-10F  /*needed to scale int -> float*/

INT32 CORDIC_CTAB_32 [] = {0x3243f6a8, 0x1dac6705, 0x0fadbafc, 0x07f56ea6, 0x03feab76, 0x01ffd55b, 0x00fffaaa, 0x007fff55, 
                           0x003fffea, 0x001ffffd, 0x000fffff, 0x0007ffff, 0x0003ffff, 0x0001ffff, 0x0000ffff, 0x00007fff, 
                           0x00003fff, 0x00001fff, 0x00000fff, 0x000007ff, 0x000003ff, 0x000001ff, 0x000000ff, 0x0000007f, 
                           0x0000003f, 0x0000001f, 0x0000000f, 0x00000008, 0x00000004, 0x00000002, 0x00000001, 0x00000000};

/* CORDIC Arcsine Core: vectoring mode */
INT32 CORDIC_asin(INT32 arc_in)
{
  INT32 k;
  INT32 d;
  INT32 tx;
  INT32 ty;
  INT32 x;
  INT32 y;
  INT32 z;

  x=CORDIC_1K_32;
  y=0;
  z=0;

  for (k=0; k<32; ++k)
  {
    d = (arc_in - y)>>(31);
    tx = x - (((y>>k) BIT_XOR d) - d);
    ty = y + (((x>>k) BIT_XOR d) - d);
    z += ((CORDIC_CTAB_32[k] BIT_XOR d) - d);

    x = tx; 
    y = ty; 
  }  
  return z; 
}

/* Wrapper function for scaling in-out of cordic core*/
FLOAT32 asin_wrap(FLOAT32 arc)
{
  return ((FLOAT32)(CORDIC_asin((INT32)(arc*MUL_32))*INV_MUL_32));
}

这种调用方式类似于：

#include "Cordic.h"
#include "math.h"

void main()
{
    y1 = asin_wrap(value_32); /*my implementation*/
    y2 = asinf(value_32);     /*standard math.h for comparison*/
}

结果如下：

​

左上角显示超过2000步的[-1;1]输入(0.001增量)，左下角显示我函数的输出，右下角显示标准输出，右上显示两个输出的差异。

立即看到错误不小于32位的精度。

我分析了我的代码执行的步骤(和中间结果)，在我看来，在某个点上，y的值与arc_in的初始值“足够接近”，而可能与位移位相关的内容会导致解决方案出现分歧。

我的问题：

• 我很困惑，这是CORDIC实现中固有的错误，还是我在执行过程中犯了错误？我原以为在接近极端的时候精度会下降，但是中间的尖峰是相当出乎意料的。(最值得注意的是+/- 0.6之外的那些，但即使去掉这些，也有更小的值，尽管没有那么明显)
• 如果它是CORDIC实现的一部分，是否有已知的解决办法？

编辑：

• 由于一些评论提到了这一点，是的，我测试了INT32的定义，即使编写#define INT32 int32_T也不会丝毫改变结果。
• 目标硬件上的计算时间通过在有效范围内随机输入的函数块10.000次的数百次重复来测量。观察到的平均结果(对于函数的一个调用)如下：math.h asinf() 100.00 microseconds CORDIC asin() 5.15 microseconds (显然，上一次测试是有问题的，一次新的交叉测试在有效性范围内的平均时间为100微秒。)
• 显然我找到了一个更好的方法。它可以在matlab版本的这里和C 这里中下载。我将更多地分析其内部运作情况，并在稍后进行报告。

Answer 1

回顾一下评论中提到的几点：

• 给定的代码输出与另一个CORDIC执行相同的值。这包括所述的不准确之处。
• 最大的错误是当您接近arcsin(1)时。
• 第二大错误是arcsin(0.60726)到arcsin(0.68514)的值都返回0.754805。
• 对于包括arcsin在内的一些函数，CORDIC方法的不精确性存在一些模糊的参考文献。给出的解决方案是执行“双迭代”，尽管我无法使其工作(所有的值都会产生大量错误)。
• 交替CORDIC实现在arcsin()实现中有一个注释/* |a| < 0.98 */，这似乎加强了已知的不准确接近1。

作为几种不同方法的粗略比较，考虑以下结果(在桌面上执行的所有测试、使用MSVC++ 2010的MSVC++计算机、在arcsin()范围为0-1的范围内使用10M迭代进行的基准测试)：

1. 问题CORDIC代码： 1050 ms，平均误差0.008，最大错误0.173
2. (参考):交替CORDIC码 2600 ms，平均误差0.008，最大误差0.173
3. atan() CORDIC代码： 2900 ms，0.21个平均错误，0.28个最大错误
4. 使用双迭代的CORDIC： 4700 ms，平均误差0.26，最大误差(？？)
5. 数学内置的asin(): 200 ms，0个平均错误，0max错误
6. (参考):有理逼近 250 ms，平均误差0.21，最大误差0.26
7. 线性表查找(见下文) 100 ms，平均误差0.000001，最大误差0.00003
8. Taylor系列(7次方， 300 ms，0.01avg误差，0.16max误差)

这些结果都在桌面上，所以它们与嵌入式系统的相关性是一个很好的问题。如果有疑问，将建议对有关系统进行定性/基准制定。大多数被测试的解决方案在范围(0-1)上都没有很好的精度，而且除了一个之外，所有的解决方案实际上都比内置的asin()函数慢。

线性表查找代码张贴在下面，这是我通常的方法，任何昂贵的数学函数时，希望速度超过精度。它只使用一个1024元素表与线性插值。在测试的所有方法中，它似乎都是最快的，也是最精确的，尽管内置的asin()实际上并不慢(测试它！)。它可以很容易地调整，或多或少的准确性，通过改变表的大小。

// Please test this code before using in anything important!
const size_t ASIN_TABLE_SIZE = 1024;
double asin_table[ASIN_TABLE_SIZE];

int init_asin_table (void)
{
    for (size_t i = 0; i < ASIN_TABLE_SIZE; ++i)
    {
        float f = (float) i / ASIN_TABLE_SIZE;
        asin_table[i] = asin(f);
    }    

    return 0;
}

double asin_table (double a)
{
    static int s_Init = init_asin_table(); // Call automatically the first time or call it manually
    double sign = 1.0;

    if (a < 0) 
    {
        a = -a;
        sign = -1.0;
    }

    if (a > 1) return 0;

    double fi = a * ASIN_TABLE_SIZE;
    double decimal = fi - (int)fi;

    size_t i = fi;
    if (i >= ASIN_TABLE_SIZE-1) return Sign * 3.14159265359/2;

    return Sign * ((1.0 - decimal)*asin_table[i] + decimal*asin_table[i+1]);
}

Answer 2

当参数大于'x‘的初始值时，“单旋转”弧正弦就会出问题，这是神奇的缩放因子- 1/An ~= 0.607252935 ~= 0x26DD3B6A。

这是因为，对于所有参数> 0，第一步总是有y=0< arg，所以d= +1，它设置y= 1/An，并留下x= 1/An。看第二步：

• 如果arg <= 1/An，那么d= -1和遵循的步骤收敛到一个很好的答案。
• 如果arg > 1/An，那么d= +1，并且这个步骤进一步偏离正确的答案，对于一个略大于1/An的值范围，随后的步骤都有d= -1，但无法更正结果:-(

我发现：

 arg = 0.607 (ie 0x26D91687), relative error 7.139E-09 -- OK    
 arg = 0.608 (ie 0x26E978D5), relative error 1.550E-01 -- APALLING !!
 arg = 0.685 (ie 0x2BD70A3D), relative error 2.667E-04 -- BAD !!
 arg = 0.686 (ie 0x2BE76C8B), relative error 1.232E-09 -- OK, again

该方法的描述警告了abs(arg) >= 0.98 (大约)，我发现在0.986之后的某个地方进程无法收敛，相对误差跃升到~5E-02，并达到1E-01 (！)在arg=1 :-(

正如您所做的，我还发现，对于0.303 < arg < 0.313，相对误差跃升到~3E-02，并慢慢减小，直到一切恢复正常。(在本例中，步骤2超调到目前为止，其余步骤无法纠正它。)

所以..。在我看来，用于弧正弦的单一旋转CORDIC看起来很垃圾：

后来又加上..。当我仔细观察单个旋转的CORDIC时，我发现了更多的小区域，这些区域的相对误差很大…

...so，我根本不想把它作为一种方法.这不仅仅是垃圾，它是无用的。

顺便说一句:我彻底推荐了“初级功能软件手册”，威廉·科迪和威廉·怀特，普伦蒂斯-霍尔，1980年。计算函数的方法不再那么有趣了(但是对所需的相关范围缩减进行了彻底、实际的讨论)。但是，对于每个函数，它们都给出了一个很好的测试过程。

Answer 3

问题末尾链接的附加源显然包含了解决方案。

拟议的守则可简化如下：

#define M_PI_2_32    1.57079632F
#define SQRT2_2      7.071067811865476e-001F /* sin(45°) = cos(45°) = sqrt(2)/2 */

FLOAT32 angles[] = {
    7.8539816339744830962E-01F, 4.6364760900080611621E-01F, 2.4497866312686415417E-01F, 1.2435499454676143503E-01F,
    6.2418809995957348474E-02F, 3.1239833430268276254E-02F, 1.5623728620476830803E-02F, 7.8123410601011112965E-03F,
    3.9062301319669718276E-03F, 1.9531225164788186851E-03F, 9.7656218955931943040E-04F, 4.8828121119489827547E-04F,
    2.4414062014936176402E-04F, 1.2207031189367020424E-04F, 6.1035156174208775022E-05F, 3.0517578115526096862E-05F,
    1.5258789061315762107E-05F, 7.6293945311019702634E-06F, 3.8146972656064962829E-06F, 1.9073486328101870354E-06F,
    9.5367431640596087942E-07F, 4.7683715820308885993E-07F, 2.3841857910155798249E-07F, 1.1920928955078068531E-07F,
    5.9604644775390554414E-08F, 2.9802322387695303677E-08F, 1.4901161193847655147E-08F, 7.4505805969238279871E-09F,
    3.7252902984619140453E-09F, 1.8626451492309570291E-09F, 9.3132257461547851536E-10F, 4.6566128730773925778E-10F};

FLOAT32 arcsin_cordic(FLOAT32 t)
{            
    INT32 i;
    INT32 j;
    INT32 flip;
    FLOAT32 poweroftwo;
    FLOAT32 sigma;
    FLOAT32 sign_or;
    FLOAT32 theta;
    FLOAT32 x1;
    FLOAT32 x2;
    FLOAT32 y1;
    FLOAT32 y2;

    flip       = 0; 
    theta      = 0.0F;
    x1         = 1.0F;
    y1         = 0.0F;
    poweroftwo = 1.0F;

    /* If the angle is small, use the small angle approximation */
    if ((t >= -0.002F) && (t <= 0.002F))
    {
        return t;
    }

    if (t >= 0.0F) 
    {
        sign_or = 1.0F;
    }
    else
    {
        sign_or = -1.0F;
    }

    /* The inv_sqrt() is the famous Fast Inverse Square Root from the Quake 3 engine
       here used with 3 (!!) Newton iterations */
    if ((t >= SQRT2_2) || (t <= -SQRT2_2))
    {
        t =  1.0F/inv_sqrt(1-t*t);
        flip = 1;
    }

    if (t>=0.0F) 
    {
        sign_or = 1.0F;
    }
    else
    {
        sign_or = -1.0F;
    }

    for ( j = 0; j < 32; j++ ) 
    {
        if (y1 > t)
        {
            sigma = -1.0F;
        }
        else
        {
            sigma = 1.0F;
        }

        /* Here a double iteration is done */
        x2 =                       x1  - (sigma * poweroftwo * y1);
        y2 = (sigma * poweroftwo * x1) +                       y1;

        x1 =                       x2  - (sigma * poweroftwo * y2);
        y1 = (sigma * poweroftwo * x2) +                       y2;

        theta  += 2.0F * sigma * angles[j];

        t *= (1.0F + poweroftwo * poweroftwo);

        poweroftwo *= 0.5F;
    }

    /* Remove bias */
    theta -= sign_or*4.85E-8F;

    if (flip)
    {
        theta = sign_or*(M_PI_2_32-theta);
    }

    return theta;
}

应指出以下几点：

• 它是一个“双迭代”CORDIC实现。
• 因此，angles表的构造与旧表不同。
• 计算采用浮点表示法，这将大大增加目标硬件上的计算时间。
• 输出中存在一个小的偏置，通过theta -= sign_or*4.85E-8F;通道去除。

下图显示了旧实现(顶部)的绝对(左)和相对错误(右)与这个答案中包含的实现(下)。

只有将CORDIC输出除以内置的数学实现的输出，才能获得相对误差。因为这个原因，它是围绕着1而不是0绘制的。

峰值相对误差(不除以零)是1.0728836e-006。

平均相对误差为2.0253509e-007 (几乎符合32位精度)。

​

​