C++(和数学):三角函数的快速逼近

Question

我知道这是一个反复出现的问题,但我还没有真正找到有用的答案.我基本上是在寻找acosC++ 函数的快速近似,我想知道我是否能够大大超过标准函数.

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

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

其中t_i已知的真实(双)数字,并且n非常小(如小于6).我需要至少1e-10的精度.我目前正在使用标准sin和acosC++函数.

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

谢谢你的回答.

Answer 1

下面的代码提供了简单的实现,sin()并且acos()应该满足您的准确性要求并且您可能想要尝试.请注意,您平台上的数学库实现很可能针对该平台的特定硬件功能进行了高度调整,并且可能还在汇编中进行编码以实现最高效率,因此不能满足硬件细节的简单编译C代码不太可能提供更高的性能,即使精度要求从完全双倍精度稍微放松.正如Viktor Latypov指出的那样,搜索不需要昂贵调用超越数学函数的算法替代方案也是值得的.

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

对平台一无所知,我选择了简单的多项式近似,使用Estrin方案和Horner方案进行评估.有关这些评估方案的说明,请参阅Wikipedia:

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

近似值本身是minimax类型,并使用Remez算法为此答案自定义生成:

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

参数减少中使用的身份在acos()注释中注明,因为sin()我使用了Cody/Waite风格的参数减少,如下面的书中所述:

WJ Cody,W.Waite,基本功能软件手册.Prentice-Hall,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;
}

Answer 2

正如 Jonas Wielicki 在评论中提到的那样，您无法做出太多精确的权衡。

最好的选择是尝试使用函数的处理器内在函数（如果您的编译器尚未执行此操作）并使用一些数学来减少必要的计算量。

同样非常重要的是，将所有内容保持为 CPU 友好的格式，确保很少有缓存未命中等。

如果您正在计算大量函数，例如acos迁移到 GPU 可能是您的一个选择？