三角参数归约（归约模 2π）

Question

我正在尝试创建一个在技术上仅在 0-pi/2 范围内运行的正弦和余弦计算器。现在这可能看起来很傻，但稍后会使用它，所以我可以使用泰勒级数。

我有一个大部分工作的实现，但是当 theta 采用 x * (pi/2) 的形式时，我遇到了一个严重的问题，其中 x 是任意整数。似乎在这些值上，有时它们会被推入它们不属于的附近象限。还有一些我无法解释的偶尔的彻底错误。

我怎样才能使它更有效和正确？

这是执行此操作的代码。

#define T_PI (2.0 * M_PI)
#define H_PI (0.5 * M_PI)
void sincos(float theta, float* cosine, float* cosine) {
  int mode;
  prepareForRange(&theta, cosine, sine);
  Assert(!(f < 0.0 || f > H_PI));
  *cosine = cos(theta);
  *sine = sin(theta);
  range_output(mode, cosine, sine);
}
void prepareForRange(float* theta, int* mode, float *cosine, float* sine) {
  if (*theta < 0.0) *theta += ceil(-*theta / T_PI) * T_PI;
  *mode = (int)floor(*theta / H_PI) % 4 + 1;
  *theta = fmodf(*theta, H_PI);
}
void range_output(int mode, float *cos, float *sin) {
  float temp;
  switch (mode) {
    case 1:
      break;
    case 2:
      temp = *cos;
      *cos = -*sin;
      *sin = temp;
      break;
    case 3:
      *cos = -*cos;
      *sin = -*sin;
      break;
    case 4:
      temp = *cos;
      *cos = *sin;
      *sin = -temp;
      break;
    default:
      break;
    }
}

Answer 1

您遇到了一个长期存在且经常未被识别的问题区域，称为range reduction。基本问题是浮点 PI 常数仅定义为 8 位精度，因此当您尝试计算 (x - n * PI) 为 n ~ 10^4 时，结果中的精度已丢失一半，并且随着 n 变大而变得更糟。这个问题没有简单的软件解决方案。为了在我自己的numerical library 中真正解决它，我必须有效地实现任意精度的浮点运算并存储一个 320 位（1078 位）的 PI​​ 常数。 libc 的一些实现有效地为您做到了这一点，但不是全部，因此您不能安全地假设它。