A while ago I noticed some weird fluctuations in my flux calculations in OpenOffice, and thought to myself that they were caused by errors in the sine and cosine functions. So I decided to check out how the trigonometric functions were implemented in gcc. Despite the fact that I run Windows, and most probably my OpenOffice was compiled in Visual Studio and not gcc, I still learned a lot about the labyrinth that is floating point operations, and the way that the core libraries handle them.

After digging a bit in the code, it can be seen that the glibc implementation of trigonometric functions depends on the processor you are using. The most common case is your basic AMD/Intel x86 architecture – just the ordinary household PC. All such modern hardware is capable of performing floating point calculations within the assembly language. Specifically, there is an opcode mnemonic, FSIN, which takes the topmost number in a special floating point stack, and returns its sine function. Originally these calculations were performed within a special processing chip called the FPU (floating point unit), but today they are incorporated within the CPU itself. So with sine and cosine calculation easy at hand, all glibc does is wrap the calls to the assembly. I do not know how the functions are implemented in the hardware level; it could be either with a Taylor series or a lookup table, or some other unknown method, all coded in hardware.

If however, you run on some other architecture, for example, PowerPC, there is no guarantee that your processor has these opcodes. Hence, glibc must write its own the trigonometric functions. Being open source software which is probably less rigorously tested than an FPU, these can have errors (1 , 2). The functions are implemented by using the Taylor’s series representation. For example, the cosine function declares these constants:

one = 1.0000000000e+00, /* 0x3f800000 */ C1 = 4.1666667908e-02, /* 0x3d2aaaab */ C2 = -1.3888889225e-03, /* 0xbab60b61 */ C3 = 2.4801587642e-05, /* 0x37d00d01 */ C4 = -2.7557314297e-07, /* 0xb493f27c */ C5 = 2.0875723372e-09, /* 0x310f74f6 */ C6 = -1.1359647598e-11; /* 0xad47d74e */

And then uses them in a clear and concise manner:

//... z = x*x; r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); if(ix < 0x3e99999a) /* if |x| < 0.3 */ return one - ((float)0.5*z - (z*r - x*y)); else { //...

It may not be immediately evident, but in these segments of code are the first 8 coefficients of the Taylor’s series (note how the constant “one” is declared, and yet in the actual logic of the code, there exists a “(float)0.5”). The rest of the function and its wrappers is not much more readable. The reason for this is the various optimizations given to such a precious operation. The code above is taken from the core of the cosine function, which is intended to be used for values of x < |pi/4|, where it performs the best in terms of speed and accuracy. For all other legal values, it figures out the remainder of division by pi/2, and uses either sine or cosine Taylor representation to figure out the result, depending on the remainder. This is ok, because if we know the sine or cosine of all x that is within a segment of length pi/2, we can figure out the entire function for any x – it can be built up of mirrors of the values in that segment around the x and y axes.

That is all well and good, and I might have continued to research the behaviour of these hard to read functions. However, as previously said, I run Windows, and do not know if that is indeed that function that will be called; furthermore, I have an Intel processor, so the FSIN will probably be the one called anyway. Reading up on it, I saw that it has accuracy errors:

“For Intel’s hardware routine the typical absolute error in sin(x) increases linearly with the argument, approximately as 4e-22 x. When the result is rounded to double precision, the function suffers partial loss of precision (> 1 ulp absolute) for x >= 2^18, and total loss of precision (50% of full scale) for x >= 2^65.”

Wishing to check this out, I plotted out the result of the following:

import math res = [] for i in range(800): res.append[ -math.sin( i * math.pi )]

You would think that every calculation should give only 0, but in fact, starting from sin(pi), which returned 1.2246467991473532e-16, the return value tended to grow in amplitude. This is most probably because of rounding errors when calculating multiples of pi, but errors in the sine function might also contribute. Plotting the results out against i, I got:

Amazingly, the exact same graph as encountered when calculating the energy flux! Turns out I was right in my guess that the errors were caused by trigonometric functions, although I must say, the verdict of where in the the sine function the error is has still to be reached. The slope which might be constructed from this graph does not indicate a 4e-22 x increase as suggested, but then again, that may be valid only for much larger x values (sin is written to handle very large numbers), while I only showed relatively small values of x.

Over the course of my journey through the land of floating point operations, I discovered some very nasty stuff. There are multiple standards and methods to manipulate floating point numbers, and sometimes the same machine will use different ones in the same calculation request. For example, it might use an 80 bit number while processing it, but round it down to 53 or 32 bit when saving it in the RAM. These roundoffs cause many errors, and you simply cannot treat floats like you treat integers. Refer here for an absolutely superb explanation and solution.

So, all in all – I did nothing new, but learned a lot about floating points and low level calculations. If anyone intends to write a program which might deal with such numbers (for example, fluid simulations or climate models), be sure to remember that using these numbers is not intuitive. In order to get accurate results for models, one must take low level considerations that at times may seem rather ugly. It might be icky, but given the nature of today’s computers, you cannot avoid it. Floats tend to be ok with most “normal numbers”, but very large or very small numbers (in human standards) might cause large errors in calculations. Remember this, and use your arithmetic wisely.