arctan(x) 函数给出错误答案答案

【问题标题】：arctan(x) function gives wrong answerarctan(x) 函数给出错误答案
【发布时间】：2012-06-17 01:31:17
【问题描述】：

我将 Maclaurin 系列用于 arctan(x)，但没有得到正确答案。我正在以弧度进行计算。到目前为止的功能如下：

fp32 t32rArcTangent(fp32 number)
{
fp32 a, b, c, d;    /* Temp Variables */
fp32 t;             /* Number Temp */
uint32 i;           /* Loop Counter */

/* Time Savers */
if (b32fpcomp(number, MM_FP8INFINITY)) return((fp32)MM_PI / 2);
if (b32fpcomp(number, -MM_FP8INFINITY)) return(-(fp32)MM_PI / 2);

/* Setup */
a = 0;
b = 0;
c = 1;
d = number;
t = number * number;

/* Calculation Loop */
for (i = 0; i < MMPRVT_FP32_TRIG_LIMIT; i++)
  {
    b += d;
    if (b32fpcomp(a, b)) break;
    a = b;
    c += 2;
    d *= -1 * t / c;
  }
#ifdef DEBUG
printf("Loops: %lu\n", i);
#endif

/* Result */
return(a);

fp32 = typedef'd float

uint32 = typedef'd unsigned long int

MM_FP8INFINITY 是 fp32 数据类型可以包含的最大数。

MM_PI 只是大约 50 位的 PI。

MMPRVT_FP32_TRIG_LIMIT 是可用于计算结果的最大循环数。这是为了防止级数展开进入无限循环，如果由于某种原因导致级数无法收敛。

这些是我得到的结果：

Testing arctangent(x) function.
Loops: 0
arctan(0):      0
Loops: 8
arctan(1):      0.724778414
Loops: 13
arctan(R3):     0.709577262
Loops: 6
arctan(1/R3):   0.517280579

R3 只是 3 的平方根，即 1.732050808....

现在我知道arctan级数的收敛半径是|x|

感谢您指出这一点。问题已得到纠正，我也完成了输入减少。这是已完成且已更正的函数，现在可以给出正确答案：

fp32 t32rArcTangent(fp32 number)
{
fp32 a, b, c, d;    /* Temp Variables */
fp32 t;             /* Number Temp */
uint32 i;           /* Loop Counter */
uint8 fr;           /* Reduction Flag */

/* Time Savers */
if (b32isInf(number) == -1) return(-(fp32)MM_PI / 2);
if (b32isInf(number) == 1) return((fp32)MM_PI / 2);
if (b32isNaN(number)) return(number);
if (b32fpcomp(number, MM_FP8INFINITY)) return((fp32)MM_PI / 2);
if (b32fpcomp(number, -MM_FP8INFINITY)) return(-(fp32)MM_PI / 2);
if (b32fpcomp(number, ONE)) return((fp32)MM_PI / 4);
if (b32fpcomp(number, -ONE)) return(-(fp32)MM_PI / 4);

/* Reduce Input */
if (number > ONE)
    {
      number = 1 / number;
      fr = 1;
    }
  else fr = 0;

/* Setup */
a = 0;
b = 0;
c = 1;
d = number;
t = number * number;

/* Calculation Loop */
for (i = 0; i < MMPRVT_FP32_TRIG_LIMIT; i++)
  {
    b += d / c;
    if (b32fpcomp(a, b)) break;
    a = b;
    c += 2;
    d *= -1 * t;
    #ifdef DEBUG
    printf("a=%g b=%g, c=%g d=%g\n", a, b, c, d);
    #endif
  }
#ifdef DEBUG
printf("Loops: %lu\n", i);
#endif

/* Result */
if (fr != 0) a = ((fp32)MM_PI / 2) - a;
return(a);
}

【问题讨论】：

什么是b32fpcomp？此外，您应该使用doubles 而不是floats。
b32fpcomp 是一个 32 位二进制浮点比较函数。我用它来清除 gcc 给出的警告，说将浮点数与 == 或 != 进行比较是不安全的。我不明白为什么，但我会幽默的编译器。至于使用浮点数，这是多精度，所以当我解决所有问题时，这些例程将用于浮点、双精度和长双精度数据类型。 32 是指浮动。下一个数据类型是 fp64，它是 64 位浮点数或双精度数。 fp96 是 96 位或 128 位浮点数或长双精度数。

标签： c algorithm math trigonometry

【解决方案1】：

想想每个循环中的术语在除以c 后会发生什么变化：

c += 2;
d *= -1 * t / c;

首先，您要除以 1 [在此之前隐式地]，然后除以 3，然后再除以 5，这听起来不错，但是因为您将 d 乘以这个术语，所以您实际上是除以乘积每个除数。 IOW，而不是

x - 1/3*x^3 + 1/5*x^5 - 1/7*x^7 + 1/9*x^9

你想要的，你正在计算

x - 1/(1*3)*x^3 + 1/(1*3*5)*x^5 - 1/(1*3*5*7)*x^7 + 1/(1*3*5*7*9)*x^9

您仍然可以使用d *= -t 技巧，但您应该移动除法。

【讨论】：