In this post I will write a C library to calculate sines and cosines using Taylor polynomials. It is impossible to calculate these directly but they can be approximated to any accuracy using this method, and I will show how a Taylor series converges to ever-increasing accuracy.
Taylor Polynomials
I recently wrote a post on linear regression, the process of creating a formula in the format y = ax + b to describe the line of best fit for a set of data. Such a formula allows the value of y to be calculated for any given value of x using simple arithmetic. Unfortunately such formulae do not exist for certain functions such as sine and cosine, and therefore these need to be calculated using techniques such as Taylor polynomials. You might like to read the Wikipedia article, noting particularly their use for calculating sines and cosines.
The formula for sine looks slightly intimidating in the forms presented by Wikipedia
but in the second form we can see that it consists of a simple repeating term.
There are four points to note about this formula:
- We start off with the angle x (in radians) we want to calculate the sine of
- We then alternately subtract or add the repeating term
- The powers and factorials start at 3 and are incremented by 2 on each iteration
- The more terms we calculate the more accuracy we get
The Project
This project will consist of a small library implementing the following two functions
- double taylor_sin_rad(double radians, int degree)
- double taylor_cos_rad(double radians, int degree)
as well as a couple of functions which will be called from main to test and demonstrate the code. We will also use the Memoization of Factorials project I posted recently to avoid repeatedly calculating factorials. If you are calculating a large number of trigonometric functions, for example when plotting a sine wave, the repetitive overhead of calculating factorials can be significant.
Create a folder somewhere convenient and within it create the following empty files. You can also download the source code as a zip or from Github if you prefer.
- memoizationfactorials.h
- memoizationfactorials.c
- taylorseries.h
- taylorseries.c
- main.c
Source Code Links
I won't repeat the memoizationfactorials code here but it is included with the download and of course in the post linked to above.
Firstly let's get the header file out of the way. Open taylorseries.h and enter/paste the following.
taylorseries.h
#include<stdbool.h> #define PI 3.141592654 #define DEGREES_IN_RADIAN (180.0 / PI) //-------------------------------------------------------- // FUNCTION PROTOTYPES //-------------------------------------------------------- bool taylor_init(int degree); void taylor_free(); double taylor_sin_rad(double radians, int degree); double taylor_cos_rad(double radians, int degree);
Now open taylorseries.c and enter the following.
taylorseries.c (part 1)
#include<stdbool.h> #include"memoizationfactorials.h" #include"taylorseries.h" factorials facs; //-------------------------------------------------------- // FUNCTION taylor_init //-------------------------------------------------------- bool taylor_init(int degree) { return factorials_calculate(&facs, degree); } //-------------------------------------------------------- // FUNCTION taylor_free //-------------------------------------------------------- void taylor_free() { factorials_free(&facs); }
This file starts off with a couple of utility functions, the first being taylor_init. The maximum number in the 3,5,7... sequence of the Taylor series formula shown above is known - rather confusingly - as the degree of the polynomial. It has nothing to do with the meaning of the word degree as a unit of angle! The degree to which we want to calculate the Taylor series corresponds to the maximum factorial which we need, so this value is passed to taylor_init, which calls factorials_calculate (from the memoizationfactorials code) to pre-calculate all the factorials we will need. The second function, factorials_free, simply frees up the memory after we have finished.
Now let's move on to the central function of this whole project, the snappily-named taylor_sin_rad. This function will take as arguments the angle in radians we wish to calculate the sine of, and the degree to which we wish to calculate the Taylor series. The function then returns the sine. Note that I have included _rad in the function name to make it obvious that the function takes radians not degrees as it's angle argument. This is a convention I usually adopt if there could be any doubt about units.
Paste the function's code into taylorseries.c.
taylorseries.c (part 2)
//-------------------------------------------------------- // FUNCTION taylor_sin_rad //-------------------------------------------------------- double taylor_sin_rad(double radians, int degree) { int multiplier = -1.0; double sine = radians; for(int currentdegree = 3; currentdegree <= degree; currentdegree += 2) { sine += ((pow(radians, currentdegree) / facs.calculated[currentdegree]) * multiplier); multiplier *= -1; } return sine; }
As I mentioned above calculating a Taylor series involves alternating between subtracting and adding the next term. To cut down on fiddly code we'll always ADD the next term but first multiply it by a value which will flip between -1 and 1. It is initialized to -1 but is multiplied by -1 on each iteration. The sine variable is initialised to the angle before being used in the loop.
The loop iterates from 3 to the specified degree, incrementing by 2 each time. The current term is then evaluated, using currentdegree as an exponent and to index the pre-calculated factorial. The multiplier sign is then "flipped", and when the loop exits we return the sine.
The last function is taylor_cos_rad. The formula for calculating cosines this is almost identical
but subtly different in that it uses even numbers as exponents and factorials. We could write another near-identical function, or we could write a general-purpose function that calculates either, passing a parameter to specify whether we want the sine or cosine and initialising currentdegree to 3 or 2 respectively. The first option is unnecessarily repetitive, the second is better but might irritate users used to being offered separate functions. We'll therefore do neither but will "cheat" by providing a separate cosine function which will simply call the sine function but with the angle shifted by 90 degrees. We therefore get away with writing just one line of code!
taylorseries.c (part 3)
//-------------------------------------------------------- // FUNCTION taylor_cos_rad //-------------------------------------------------------- double taylor_cos_rad(double radians, int degree) { return taylor_sin_rad(radians - (90 / DEGREES_IN_RADIAN), degree + 1); }
That's the code finished so now we can use it. In main paste the following code.
main.c
#include<stdio.h> #include<math.h> #include<stdlib.h> #include"taylorseries.h" //-------------------------------------------------------- // FUNCTION PROTOTYPES //-------------------------------------------------------- void print_sines(double from_degrees, double to_degrees, double step); //-------------------------------------------------------- // FUNCTION main //-------------------------------------------------------- void main(void) { puts("-------------------------------"); puts("| codedrome.com |"); puts("| Taylor Series Trigonometric |"); puts("-------------------------------\n"); print_sines(0.0, 360.0, 30); } //-------------------------------------------------------- // FUNCTION print_sines //-------------------------------------------------------- void print_sines(double from_degrees, double to_degrees, double step) { if(taylor_init(11)) { printf("-----------------------------------------------------------------------------------------\n"); printf("| Degrees | C Library | Taylor 3 | Taylor 5 | Taylor 7 | Taylor 9 | Taylor 11 |\n"); printf("-----------------------------------------------------------------------------------------\n"); for(double degrees = from_degrees; degrees <= to_degrees; degrees += step) { printf("| %7.0lf | %10.6lf | %10.6lf | %10.6lf | %10.6lf | %10.6lf | %10.6lf |\n", degrees, sin(degrees / DEGREES_IN_RADIAN), taylor_sin_rad(degrees / DEGREES_IN_RADIAN, 3), taylor_sin_rad(degrees / DEGREES_IN_RADIAN, 5), taylor_sin_rad(degrees / DEGREES_IN_RADIAN, 7), taylor_sin_rad(degrees / DEGREES_IN_RADIAN, 9), taylor_sin_rad(degrees / DEGREES_IN_RADIAN, 11)); } printf("-----------------------------------------------------------------------------------------\n"); taylor_free(); } }
The print_sines function starts off by calling taylor_init with the argument 11, this being the highest degree we will use. If successful (taylor_init ultimately calls malloc so success cannot be guaranteed) we then iterate from the specified start angle to the specified end angle (in degrees) in the specified steps. For each angle it prints the sine as calculated by the C library's own sin function, which we'll consider to be the "target" value. The sine is then calculated by our own Taylor function to several different degrees to show the result becoming increasingly accurate.
We can now compile and run the code using these commands in terminal.
Compile and run
gcc main.c taylorseries.c memoizationfactorials.c -std=c11 -lm -o main /.main
Program Output
------------------------------- | codedrome.com | | Taylor Series Trigonometric | ------------------------------- ----------------------------------------------------------------------------------------- | Degrees | C Library | Taylor 3 | Taylor 5 | Taylor 7 | Taylor 9 | Taylor 11 | ----------------------------------------------------------------------------------------- | 0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | | 30 | 0.500000 | 0.499674 | 0.500002 | 0.500000 | 0.500000 | 0.500000 | | 60 | 0.866025 | 0.855801 | 0.866295 | 0.866021 | 0.866025 | 0.866025 | | 90 | 1.000000 | 0.924832 | 1.004525 | 0.999843 | 1.000004 | 1.000000 | | 120 | 0.866025 | 0.563221 | 0.899045 | 0.863971 | 0.866108 | 0.866023 | | 150 | 0.500000 | -0.372581 | 0.652273 | 0.485029 | 0.500950 | 0.499958 | | 180 | -0.000000 | -2.026120 | 0.524044 | -0.075221 | 0.006925 | -0.000445 | | 210 | -0.500000 | -4.540945 | 0.970964 | -0.792011 | -0.463078 | -0.503248 | | 240 | -0.866025 | -8.060603 | 2.685767 | -1.803648 | -0.709604 | -0.884114 | | 270 | -1.000000 | -12.728642 | 6.636667 | -3.602330 | -0.444366 | -1.081890 | | 300 | -0.866025 | -18.688608 | 14.106711 | -7.300487 | 0.850770 | -1.180788 | | 330 | -0.500000 | -26.084051 | 26.733139 | -14.983433 | 4.236804 | -1.559467 | | 360 | 0.000000 | -35.058517 | 46.546732 | -30.159127 | 11.899567 | -3.195076 | -----------------------------------------------------------------------------------------
As I mentioned above, for smaller angles the sines are reasonably accurate compared to the C library sin function values, even with smaller degrees of the Taylor polynomial, and completely accurate to 6 decimal places up to 90 degrees at degree 11.
The values then become increasingly inaccurate. For 360 degrees the Taylor value is about -35 at degree 3 when it should be 0, and even at degree 11 is about -3.
Due to the symmetrical and cyclical nature of the sine function the inaccuracy beyond 90 degrees is not a problem - we can easily use the 0 to 90 degree values for any angle by shifting and/or negating them by the appropriate amount. I will go into this in detail in a future post.
It would also be useful to see a graph of the various sine waves of the above values: this will also form a future post.