Source code for trigonometric functions calculations

Question

For program that needs to be deterministic and provide the same result on different platforms (compilers), the built-in trigonometric functions can't be used, since the algorithm to compute it is different on different systems. It was tested, that the result values are different.

(Edit: the results need to be exactly the same to the last bit as it is used in game simulation that is ran on all the clients. These clients need to have the state of the simulation exactly the same to make it work. Any small error could result in bigger and bigger error over time and also the crc of the game state is used as check of synchronisation).

So the only solution that I came up with was to use our own custom code to calculate these values, the problem is, that (surprisingly) it is very hard to find any easy to use source code for all the set of the trigonometric functions.

This is my modification of the code I got (https://codereview.stackexchange.com/questions/5211/sine-function-in-c-c) for the sin function. It is deterministic on all platforms and the value is almost the same as the value of standard sin (both tested).

#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;
}

But I didn't find anything suitable for other functions, like asin, atan, tan (other than the sin/cos) etc.

These functions doesn't have be as precise as the standard ones, but at least 8 figures would be nice.

Answer 1

"It was tested, that the result values are different."

How different is different enough to matter? You claim to want 8 significant (decimal?) digits of agreement. I don't believe that you've found less than that in any implementation that conforms to ISO/IEC 10967-3:2006 §5.3.2.

Do you understand how trivial a trigonometric error of one part per billion represents? It would be under 3 kilometers on a circle the size of the earth's orbit. Unless you are planning voyages to Mars, and using sub-standard implementation, your claimed "different" ain't going to matter.

added in response to comment:

What Every Programmer Should Know About Floating-Point Arithmetic. Read it. Seriously.

Since you claim that:

1. precision isn't as important as bit for bit equality
2. you need only 8 significant digits

then you should truncate your values to 8 significant digits.