C++(和数学): fast approximation of a trigonometric function

Question

我知道这是一个反复出现的问题，但我还没有真正找到有用的答案。我基本上是在寻找 C++ 中函数 acos 的快速近似值，我想知道我是否可以显着击败标准函数。

但是你们中的一些人可能对我的具体问题有见解:我正在编写一个科学程序，我需要非常快。主要算法的复杂性归结为计算以下表达式(多次使用不同的参数):

sin( acos(t_1) + acos(t_2) + ... + acos(t_n) )

其中 t_i 是已知的实数( double )，而 n 非常小(例如小于 6)。我需要至少 1e-10 的精度。我目前正在使用标准的 sin 和 acos C++ 函数。

你认为我能以某种方式显着提高速度吗？对于那些知道一些数学的人，你认为扩展该正弦以获得根据 t_i 的代数表达式(仅涉及平方根)是否明智？

感谢您的回答。

Answer 1

下面的代码提供了 sin() 和 acos() 的简单实现，它们应该满足您的准确性要求并且您可能想尝试一下。请注意，您平台上的数学库实现很可能针对该平台的特定硬件功能进行了高度调整，并且可能还以汇编形式编码以获得最大效率，因此不太可能提供不满足硬件细节的简单编译 C 代码更高的性能，即使精度要求从全 double 有所放宽。正如 Viktor Latypov 指出的那样，寻找不需要昂贵调用超越数学函数的算法替代方案也可能是值得的。

在下面的代码中，我尝试使用简单、可移植的结构。如果您的编译器支持 rint() 函数 [由 C99 和 C++11 指定]，您可能希望使用它而不是 my_rint()。在某些平台上，对 floor() 的调用可能很昂贵，因为它需要动态更改机器状态。函数 my_rint()、sin_core()、cos_core() 和 asin_core() 想要成为内联以获得最佳性能。您的编译器可能会在高优化级别自动执行此操作(例如，使用 -O3 编译时)，或者您可以为这些函数添加适当的内联属性，例如inline 或 __inline 取决于您的工具链。

对您的平台一无所知，我选择了简单的多项式近似，使用 Estrin 方案和 Horner 方案进行评估。有关这些评估方案的描述，请参见 Wikipedia:

http://en.wikipedia.org/wiki/Estrin%27s_scheme , http://en.wikipedia.org/wiki/Horner_scheme

近似值本身属于极大极小值类型，并且是使用 Remez 算法为此答案自定义生成的:

http://en.wikipedia.org/wiki/Minimax_approximation_algorithm , http://en.wikipedia.org/wiki/Remez_algorithm

acos() 的参数缩减中使用的标识在注释中注明，对于 sin() 我使用了 Cody/Waite 风格的参数缩减，如如下书中所述:

W。 J. Cody, W. Waite，基本功能软件手册。普伦蒂斯·霍尔，1980 年

评论中提到的误差范围是近似值，未经严格测试或证明。

/* not quite rint(), i.e. results not properly rounded to nearest-or-even */
double my_rint (double x)
{
  double t = floor (fabs(x) + 0.5);
  return (x < 0.0) ? -t : t;
}

/* minimax approximation to cos on [-pi/4, pi/4] with rel. err. ~= 7.5e-13 */
double cos_core (double x)
{
  double x8, x4, x2;
  x2 = x * x;
  x4 = x2 * x2;
  x8 = x4 * x4;
  /* evaluate polynomial using Estrin's scheme */
  return (-2.7236370439787708e-7 * x2 + 2.4799852696610628e-5) * x8 +
         (-1.3888885054799695e-3 * x2 + 4.1666666636943683e-2) * x4 +
         (-4.9999999999963024e-1 * x2 + 1.0000000000000000e+0);
}

/* minimax approximation to sin on [-pi/4, pi/4] with rel. err. ~= 5.5e-12 */
double sin_core (double x)
{
  double x4, x2, t;
  x2 = x * x;
  x4 = x2 * x2;
  /* evaluate polynomial using a mix of Estrin's and Horner's scheme */
  return ((2.7181216275479732e-6 * x2 - 1.9839312269456257e-4) * x4 + 
          (8.3333293048425631e-3 * x2 - 1.6666666640797048e-1)) * x2 * x + x;
}

/* minimax approximation to arcsin on [0, 0.5625] with rel. err. ~= 1.5e-11 */
double asin_core (double x)
{
  double x8, x4, x2;
  x2 = x * x;
  x4 = x2 * x2;
  x8 = x4 * x4;
  /* evaluate polynomial using a mix of Estrin's and Horner's scheme */
  return (((4.5334220547132049e-2 * x2 - 1.1226216762576600e-2) * x4 +
           (2.6334281471361822e-2 * x2 + 2.0596336163223834e-2)) * x8 +
          (3.0582043602875735e-2 * x2 + 4.4630538556294605e-2) * x4 +
          (7.5000364034134126e-2 * x2 + 1.6666666300567365e-1)) * x2 * x + x; 
}

/* relative error < 7e-12 on [-50000, 50000] */
double my_sin (double x)
{
  double q, t;
  int quadrant;
  /* Cody-Waite style argument reduction */
  q = my_rint (x * 6.3661977236758138e-1);
  quadrant = (int)q;
  t = x - q * 1.5707963267923333e+00;
  t = t - q * 2.5633441515945189e-12;
  if (quadrant & 1) {
    t = cos_core(t);
  } else {
    t = sin_core(t);
  }
  return (quadrant & 2) ? -t : t;
}

/* relative error < 2e-11 on [-1, 1] */
double my_acos (double x)
{
  double xa, t;
  xa = fabs (x);
  /* arcsin(x) = pi/2 - 2 * arcsin (sqrt ((1-x) / 2)) 
   * arccos(x) = pi/2 - arcsin(x)
   * arccos(x) = 2 * arcsin (sqrt ((1-x) / 2))
   */
  if (xa > 0.5625) {
    t = 2.0 * asin_core (sqrt (0.5 * (1.0 - xa)));
  } else {
    t = 1.5707963267948966 - asin_core (xa);
  }
  /* arccos (-x) = pi - arccos(x) */
  return (x < 0.0) ? (3.1415926535897932 - t) : t;
}