三角函数计算的源代码

问题描述 投票:1回答:4

对于需要确定性并在不同平台(编译器)上提供相同结果的程序,不能使用内置的三角函数,因为计算它的算法在不同系统上是不同的。经过测试,结果值不同。

(编辑:结果必须与在所有客户端上运行的游戏模拟中使用的最后一位完全相同。这些客户端需要具有完全相同的模拟状态才能进行模拟任何小错误都会随着时间的流逝而产生越来越大的错误,并且游戏状态的crc也将用作同步检查)。

因此,我想到的唯一解决方案是使用我们自己的自定义代码来计算这些值,问题是,(令人惊讶地)很难为所有三角函数集找到任何易于使用的源代码功能。

这是我对sin函数获得的代码(https://codereview.stackexchange.com/questions/5211/sine-function-in-c-c)的修改。它在所有平台上都是确定性的,其值几乎与标准sin(均已测试)的值相同。

#define M_1_2_PI 0.159154943091895335769 // 1 / (2 * pi)

double Math::sin(double x)
{
  // Normalize the x to be in [-pi, pi]
  x += M_PI;
  x *= M_1_2_PI;
  double notUsed;
  x = modf(modf(x, &notUsed) + 1, &notUsed);
  x *= M_PI * 2;
  x -= M_PI;

  // the algorithm works for [-pi/2, pi/2], so we change the values of x, to fit in the interval,
  // while having the same value of sin(x)
  if (x < -M_PI_2)
    x = -M_PI - x;
  else if (x > M_PI_2)
    x = M_PI - x;
  // useful to pre-calculate
  double x2 = x*x;
  double x4 = x2*x2;

  // Calculate the terms
  // As long as abs(x) < sqrt(6), which is 2.45, all terms will be positive.
  // Values outside this range should be reduced to [-pi/2, pi/2] anyway for accuracy.
  // Some care has to be given to the factorials.
  // They can be pre-calculated by the compiler,
  // but the value for the higher ones will exceed the storage capacity of int.
  // so force the compiler to use unsigned long longs (if available) or doubles.
  double t1 = x * (1.0 - x2 / (2*3));
  double x5 = x * x4;
  double t2 = x5 * (1.0 - x2 / (6*7)) / (1.0* 2*3*4*5);
  double x9 = x5 * x4;
  double t3 = x9 * (1.0 - x2 / (10*11)) / (1.0* 2*3*4*5*6*7*8*9);
  double x13 = x9 * x4;
  double t4 = x13 * (1.0 - x2 / (14*15)) / (1.0* 2*3*4*5*6*7*8*9*10*11*12*13);
  // add some more if your accuracy requires them.
  // But remember that x is smaller than 2, and the factorial grows very fast
  // so I doubt that 2^17 / 17! will add anything.
  // Even t4 might already be too small to matter when compared with t1.

  // Sum backwards
  double result = t4;
  result += t3;
  result += t2;
  result += t1;

  return result;
}

但是我没有找到适合其他功能的任何东西,例如asin,atan,tan(除了sin / cos等等)

这些功能没有标准功能那么精确,但是至少要有8位数字才可以。

c++ trigonometry numerical-stability
4个回答
2
投票
我想最容易的是选择一个实现所需数学功能的自由运行时库:

    FreeBSD
  • go,将需要音译,但我认为所有功能都具有非汇编实现
  • MinGW-w64
  • ...
  • 然后使用它们的实现。请注意,上面列出的是公共领域或BSD许可或其他一些自由许可。如果使用代码,请确保遵守许可证。

  • 4
    投票
    “经过测试,结果值不同。”

    区别有多大足以重要?您声称需要8个有效(十进制)数字。我相信您在符合ISO/IEC 10967-3:2006§5.3.2。的任何实现中所发现的内容都不会少于

    您了解十亿分之一的三角误差代表什么呢?绕地球轨道一圈的距离不到3公里。除非您计划前往火星并使用不合标准的实施方式,否则您声称的“不同”将无关紧要。

    为回应评论添加:

    What Every Programmer Should Know About Floating-Point Arithmetic。阅读。认真。

    由于您声称:

      对于位平等,精度不如位重要
    1. 您只需要8位有效数字
    2. 然后,您应该将值截断为8个有效数字。

    2
    投票
    您可以使用泰勒级数(实际上似乎是您正在使用的,也许不知道)

    0
    投票
    我建议您使用查找表和线性/双曲线插值法。
    © www.soinside.com 2019 - 2024. All rights reserved.