4

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.

Community
  • 1
  • 1
kovarex
  • 1,766
  • 18
  • 34
  • 1
    Don't forget floating point math is implementation defined and usually architecture and compiler option dependent! If all you need is 8 digits of precision, use the standard library. I doubt you're saving time doing this. – rubenvb May 26 '14 at 20:04
  • 2
    I know, that the standard library is precise enough, but the difference on the 20th figure can still change the value on the 8th figure when rounded (0.99999999999999999 versus 1.0 for example). As I stated at the beginning of my post, the calculations have to be deterministic to the last bit to achieve determinism. This is simulation of networked computer game, where the game state is crc-ed and checked to be the same on all clients. – kovarex May 27 '14 at 06:35
  • then you must (in the strict sense of the word) use some form of integer-based fixed-point math. No way you can use floating point to achieve this unless your compiler can guarantee strict IEEE 754 conformance. – rubenvb May 27 '14 at 07:14
  • 2
    The compilers we use can guarantee strict IEEE 754, we tested that. Everything is strict and precise, unless we use the custom-implemented trigonometric functions. Having to rewrite the whole program to integer-based fixed-point math would be really A LOT of work (months of work I guess), with lot of trouble. – kovarex May 27 '14 at 07:29
  • Switching to a custom numerical class should only take one rename throughout the whole codebase (or a typedef change if you already are using a typedef). Not what I would call a lot of work. Performance though, will suffer. Anyways, I suggest using an existing implementation (see my answer). – rubenvb May 27 '14 at 07:45

4 Answers4

4

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

msw
  • 42,753
  • 9
  • 87
  • 112
  • 3
    Any difference is too much. The result has to be exactly the same to the last bit. The value doesn't need to be super precise, but it needs to be always the same. When I want to have deterministic simulator, nothing less is sufficient. – kovarex May 27 '14 at 06:33
  • @kovarex see "added" above. – msw May 27 '14 at 17:44
  • 1
    Truncating doesn't solve the problem, there will be always two "neighbour" double numbers, that will go to different truncated interval. So even small difference on 15th digit can affect (although unprobably) the value of number truncated to 8th digit. – kovarex May 27 '14 at 19:08
3

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.

rubenvb
  • 74,642
  • 33
  • 187
  • 332
2

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)

Taylor1

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.

smagnan
  • 1,197
  • 15
  • 29
  • Yes, the function is using taylor. This is something I will have to do probably, but I supposed, that this was transfered to C++ code with the proper count of steps and preparations many times already, I wanted to avoid doing the same job again. – kovarex May 27 '14 at 06:38
0

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)

Lorenzo Boccaccia
  • 6,041
  • 2
  • 19
  • 29