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

You can use Taylor series (actually it seems that it is what you are using, maybe without knowing)

Take a look on wikipedia (or everywhere else): https://en.wikipedia.org/wiki/Taylor_series

You have here the list for the most common functions (exp, log, cos, sin etc ...) https://en.wikipedia.org/wiki/Taylor_series#List_of_Maclaurin_series_of_some_common_functions but with some mathematical knowledge you can find/calculate quite everything (ok clearly not everything but ...)

Some examples (there are many others)

Notes:

1. the more terms you add the more precision you have.
2. I don't think it's the most efficient way to calculate what you need but it's a quite "simple" one (the idea I mean)
3. A factorial(n) function could be really useful if you decide to use that

I hope it will help.

Answer 2

I'd suggest looking into using lookup tables and linear/bicubic interpolation.

that way you control exactly the values at each point, and you don't have to perform a awful lot of multiplications.

Taylor expansions for sin/cos functions sucks anyway

spring rts fought ages against this kind of desync error: try posting on their forum, not many old developers remain but those that do should still remember the issues and the fixes.

in this thread http://springrts.com/phpbb/viewtopic.php?f=1&t=8265 they talk specifically on libm determinism (but different os might have different libc with subtle optimization differences, so you need to take the approach and throw the library)

Answer 3

"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.

Answer 4

I guess the easiest would be to pick a liberal runtime library which implements the required math functions:

• FreeBSD
• go, will need transliteration, but I think all functions have a non-assembly implementation
• MinGW-w64
• ...

And just use their implementations. Note the ones listed above are either public domain or BSD licensed or some other liberal license. Make sure to abide by the licenses if you use the code.