CORDIC Arcsine实现失败

Question

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

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

不幸的是,asin()函数系统地失败了某些输入.

出于测试目的,我已经实现了一个.h用于simulink S-Function的文件.(这只是为了方便起见,您可以将以下内容编译为独立版.exe,只需进行少量更改)

注意:我强制进行32次迭代,因为我工作在32位精度,并且需要最大可能的精度.

Cordic.h:

#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*/
}

结果如下所示:

左上角显示[-1;1]输入超过2000步(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)to arcsin(0.68514)all return的值0.754805。
• 对于包括 arcsin 在内的某些函数，对 CORDIC 方法中的不准确性有一些模糊的参考。给定的解决方案是执行“双迭代”，尽管我一直无法使其工作（所有值都会产生大量错误）。
• 该备用CORDIC FPGA实现有评论/* |a| < 0.98 */在这似乎强化了公知有不准确之处接近1的反正弦（）实现。

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

1. 问题 CORDIC 代码： 1050 毫秒，0.008 平均误差，0.173 最大误差
2. 替代 CORDIC 代码 ( ref )： 2600 毫秒，0.008 平均误差，0.173 最大误差
3. atan() CORDIC 代码： 2900 毫秒，0.21 平均误差，0.28 最大误差
4. CORDIC 使用双迭代： 4700 毫秒，0.26 平均误差，0.917 最大误差 (???)
5. Math 内置 asin()： 200 毫秒，0 平均误差，0 最大误差
6. 有理近似 ( ref )： 250 毫秒，0.21 平均误差，0.26 最大误差
7. 线性表查找（见下文） 100 毫秒，0.000001 平均误差，0.00003 最大误差
8. 泰勒级数（7 次方，参考）： 300 毫秒，0.01 平均误差，0.16 最大误差

这些结果在桌面上，所以它们与嵌入式系统的相关性是一个很好的问题。如果有疑问，建议在相关系统上进行分析/基准测试。测试的大多数解决方案在范围 (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]);
}