This post is about how sin and cos are calculated in the c language standard math library on a typical Linux system. Because Python floats are really c doubles, what I write here will apply to Python too. For background on how floating point arithmetic works, What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg is useful as is the wikipedia page on double precision floating point format.

To find out exactly which code is used on your system you can write a piece of c code with a sin or cos call that gcc can’t optimize out and then examine it with the debugger gdb. For example:

#include <math.h>
#include <stdlib.h>
#include <time.h>

int main() {
  srand(time(0)); // now gcc can't know the random seed
  int r = rand();
  double s = sin(r);
  return 0;
}

Compile with gcc -g filename.c -lm -o sin (the g option lets you use the debugger, lm links the math library). Start the debugger with gdb sin and set a breakpoint at the sin call with break 8, 8 being the line number on which sin is used. Run the program with run and when it breaks, hit s then enter to step through the code. You’ll see a path like ../sysdeps/ieee754/dbl-64/s_sin.c indicating which source file is used. The ieee754/dbl-64 means that my system is using 64 bit double precision IEEE754 floating point numbers which I will refer to as doubles. The source file specified won’t normally be present on your computer - under Debian you can do apt-get source libc6 to install it (in the current directory).

If you step through further with gdb you probably see a lot of file not found errors. This unix and linux SE question has answers explaining how to fix that, if you want to get rid of them.

Instead of installing the source package, you can look at the official glibc repository on sourceware.org for the most up-to-date code. The file mentioned above is here. There’s a github repository here which says it is a mirror of the sourceware repository updated daily; it might be easier to browse.

Outline of __sin(double x)

Roughly speaking, to compute sin of a double x, a multiple of \(\pi/2\) is subtracted so that the result has a small enough magnitude. The sine of x will be equal to plus or minus sin or cos of this small number. Polynomial approximations (not Taylor polynomials!) to sin and cos, possibly also lookup tables, are then used to compute the output value.

The main sin function __sin(double x) begins on line 193. It calls two auxiliary functions:

  • do_sin(x, dx) line 124 computes an approximation to \(\sin(x+dx)\) when dx is small relative to x. Depending on the magnitude of x it either uses polynomial approximations of degree 6 and a lookup table or a higher degree polynomial approximation computed by a macro.
  • do_cos(x, dx) line 100 which always uses a degree 6 polynomial and a lookup table to approximate \(\cos(x+dx)\).

__sin(double x) is computed using different techniques according to the magnitude of x.

  • For small enough |x| it just returns x.
  • For |x| up to about 0.86, it returns the value of do_sin(x, 0).
  • For |x| up to about 2.4, the do_cos function is used to compute an approximation to \(\cos(\pi/2 - \vert x\vert )\).
  • For |x| up to about 100 million, a multiple of \(\pi/2\) is subtracted using reduce_sincos and either do_sin or do_cos is used on the result.
  • For larger finite |x| a different reduction function is used, followed by either do_sin or do_cos.
  • Finally, infinities and NaNs are dealt with.

The main conditional in __sin(double x) starts on line 212, but it isn’t a simple branch on the absolute value of x. Before the conditional we have, omitting a few lines,

mynumber u;
int4 k, m, n;

u.x = x;
m = u.i[HIGH_HALF];
k = 0x7fffffff & m;    /* no sign */

An int4 is just an int (perhaps 4 refers to the 4 bytes in a 32b int?), and a mynumber is a c union

typedef union { int4 i[2]; double x; double d; } mynumber;

defined in the mydefs.h header file. It lets us access the two halves of the bit representation of a double x interpreted as two 32b ints stored in an array i. An IEEE754 64b double is stored in memory as a sign bit, 11 bits representing the exponent, then 52 bits of significand. If we split one of these into two 32b pieces, one contains the sign bit, exponent, and the most significant 20 bits of the significand and the other has the least significant 32 bits of the significand. Which of these two pieces ends up as i[0] and which as i[1] (when a mynumber is created from a double) depends on whether your system is big or little endian . For little-endian systems like mine, the sign bit, exponent, and most significant part of the significand will be in i[1]. I don’t know where HIGH_HALF is defined, but from context it must be 1 on little-endian systems so that m is a 32b int formed from this part of the bit representation of x.

(To check the endianness of your system, in Python do from sys import byteorder then print(byteorder), or in a terminal lscpu | grep "Byte Order".)

The integer constant 0x7fffffff has a bit representation of 0 followed by 31 1s, so the result of anding it with m is to set the sign bit of m to 0 if it wasn’t already and not change anything else. So k will be the absolute value of m. The branching on the size of |x| is done by branching on k.

Very very small x

The first condition (line 216) is

if (k < 0x3e500000)    /* if x->0 =>sin(x)=x */

To see how this translates into an inequality on |x|, note that the 32b int 0x3e500000 is

s eee eeee eeee ffff ... ffff
0 011 1110 0101 0000     0000

The labelling on the top row shows how the bits would be interpreted if this were the first part of a double: s means the sign bit, the 11 es are exponent bits, and the 20 fs are the 20 most significant bits of the significand or fraction part. The bits of k are the first 32 bits of |x|. Thus the inequality is equivalent to |x| being less than the double with bit representation the same as that of k, followed by 32 zeros. If we treat the bits of the exponent part of a double as an unsigned integer \(e\) then (assuming the double isn’t denormal) the number it represents is \((-1)^s 1.\mathbf{f} \times 2^{e-1023}\) where \(s\) is the sign bit and \(\mathbf{f}\) the sequence of significand bits. 0x3e5 is 997 in decimal and 997-1023=-26, so the inequality is equivalent to |x| being less than \(2^{-26}\).

This is sufficient for x to be the nearest double to \(\sin(x)\). Using Taylor’s theorem with the Lagrange remainder, \(\vert\sin(x) - x\vert \leqslant \vert x\vert^3 / 3!\). If \(\vert x\vert < 2^{-26}\) then an ULP for \(x\) is at most \(2^{-26-52} = 2^{-78}\) so half an ULP is at most \(2^{-79}\), while \(\vert x\vert ^3/6! < 2^{-78 -2} = 2^{-80}\).

x of magnitude up to 0.855469

For |x| up to 0.85469, __sin returns the value of do_sin(x, 0). This function is described in another post.

x of magnitude between 0.855469 and 2.426265

This time plus or minus do_cos(hp0 - fabs(x), hp1) is returned, where \(hp0 + hp1\) is a double-double approximation of \(\pi/2\). There is a separate post on do_cos.

x of magnitude between 2.426265 and 105414350

This time x is first range reduced using reduce_sincos, described in a separate post. Afterwards do_sincos is called, which calls do_cos or do_sin as appropriate.

finite x of magnitude larger than 105414350

This time a different range reduction function called branred is used. It is defined in a separate file.

NaNs and infinities

The code signals a domain error and returns x/x. This is a NaN in both cases (according to What Every Computer Scientist Should Know About Floating-Point Arithmetic sections 2.2.1 and 2.2.2).