三角变元约简（约简模2π)

Question

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

我有一个大部分工作的实现，但是当θ是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

您正遇到一个长期存在且经常无法识别的问题领域，称为范围缩减。基本问题是，浮点PI常数的精度只有8位数，所以当你试图计算(x-n*PI)n~10^4时，你的结果中损失了一半的精度，而且随着n变大，情况会变得更糟。这个问题没有简单的软件解决方案。为了在我自己的数值库中真正解决这个问题，我必须有效地实现任意精度浮点算术，并存储一个320位（1078位）的PI常数。libc的一些实现有效地为您实现了这一点，但不是所有的实现，所以您不能放心地假设它。