UPDATE: (19-March-2015) This answer sounds like an expert answer.
There are books full on how to calculate those functions efficiently and accurately. So it's too much for a short answer here.
I'm hardly an expert, but I can easily spot that your loop starts at the large terms and ends with the small terms. This has a larger accumulated error than going from small to large.
Secondly, subtracting is a very good way to lose a lot of accuracy. It might be better to partially unroll the loop, combining each positive term with the next negative term, and only add the difference between those two. Although that still might result in a large loss of accuracy.
As none of the other answers even touches on those point, you have not been answered by an expert yet. I only know enough to know that I don't know enough.
So clearly work for specialists. And indeed, solved a long time ago.
A quick google found What Every Computer Scientist Should Know About Floating-Point Arithmic.
Note – This appendix is an edited reprint of the paper What Every
Computer Scientist Should Know About Floating-Point Arithmetic, by
David Goldberg, published in the March, 1991 issue of Computing
Surveys. Copyright 1991, Association for Computing Machinery, Inc.,
reprinted by permission.
I just googled for that title, betting that someone would have used that title at some point in time. It explains the problem I mentioned above and much more, but I'm not really qualified to judge it. As it is SUN's SPARC documentation, I expect it to be among the best you can find online, although it could be too technical for some.
It's the appendix of some part of SUN's documentation for SPARC. Scanning the full document, In the references is this paper:
Tang, Peter Ping Tak, Some Software Implementations of the Functions
Sin and Cos, Technical Report ANL-90/3, Mathematics and Computer
Science Division, Argonne National Laboratory, Argonne, Illinois,
I can't find it online using google, so good luck hunting that one!
Drawback of the backwards loop is that one can't use the result of previous step. An even better approach might be a forward loop to calculate the terms. But don't add them, store them in an array, then use a backwards loop to sum them.
That might sound a lot of work, but when the size of the array is fixed at compile time, one can unroll the loop (see Artur Mustafin's answer), and just use 4 local doubles (each one holds the diference between two terms, so that's the same 8 terms he calculates).
Based on my remarks, here is my version. Inlined, loop-unrolled, and parallizable. I usually don't bother with re-using local variables: the compilers nowadays are smart enough to do that for me. I also assume that e.g.
/ (10*11) is optimized into a
* some_constant where
1.0 / 110, assuming the CPU is faster on multiplications than on divisions.
But I bet that the version used by real libraries is even smarter.
double sine_taylor(double 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;