Skip to main content

How to implement Sine and Cosine

· 10 min read
Strider

Hi, after a long break spent doing math, I would like to show how to tinker with sine and cosine from Scratch.

Why these two? Well, it is simply because everyone uses the functions but only few know how it looks inside. There are many implementations for these functions but still it is good to know how to implement the function itself for e.g. embedded systems or even your own libc. But to implement sine and cosine, we first have to understand what the functions are and how they work.

dia1.png

In the picture above I have shown sine and cosine. Both functions look similar. But there is one difference, and that is that with cosine you have an offset of π2\frac{\pi}{2}.

What you can also see quite quickly is that both functions have a target range of [−1;1][-1; 1], and are periodic. This means no matter what X-values you give the function you will always get values in the range [−1;1][-1; 1]. Both functions can also be represented in the unit circle where you can see what exactly both represent. We all used these two functions once in school in trigonometry, where we calculated triangles back then. Well back to the topic.

dia2.png

Sine and cosine in the unit circle are to be seen as XX, YY. The sine returns the Y-value and the cosine the X-value. This is interesting e.g. when programming analog clocks on the screen. If you take a closer look at the functions, you will see the following distinctive points, which you can take into account later in the test.

00π2\frac{\pi}{2}π\pi32π\frac{3}{2}\pi2π2\piπ4\frac{\pi}{4}π+π4\pi + \frac{\pi}{4}
sinsin001100−1-10012\frac{1}{\sqrt{2}}−12-\frac{1}{\sqrt{2}}
coscos1100−1-1001112\frac{1}{\sqrt{2}}−12-\frac{1}{\sqrt{2}}

These values will help us later when we implement both functions. What we should also know to implement sine and cosine are the power series. Because both functions can be approximated as power series. I simply write down the power series for both functions.

sin(x)≈x−13!x3+15!x5−17!x7+19!x9−111!x11±...sin(x) \approx x - \frac{1}{3!}x^{3} + \frac{1}{5!}x^{5} - \frac{1}{7!}x^{7} + \frac{1}{9!}x^{9} - \frac{1}{11!}x^{11} \pm ... cos(x)≈1−12!x2+14!x4−16!x6+18!x8−110!x10±...cos(x) \approx 1 - \frac{1}{2!}x^{2} + \frac{1}{4!}x^{4} - \frac{1}{6!}x^{6} + \frac{1}{8!}x^{8} - \frac{1}{10!}x^{10} \pm ...

Just let me show you why this represents sine and cosine. What you can see here is that sine has only odd powers and cosine only even powers. If you enter the power series, member by member, you can see that the series represents a function, which slowly becomes more and more similar to sine and cosine. This can be extended as often as you like.

dia3.png Power series approximates sine

dia4.png Power series approximates cosine

One can see in the two pictures above that in each case the middle range is approximated or approximated very well. At the ends of the pictures, however, the power series deviate further and further, which is not bad, because the complete definition range of [−π;π][-\pi; \pi] is completely covered. But if we now extend both power series to n=10n = 10, we can map the definition range of [−π;2π][-\pi; 2\pi]. As formulas or as a series, it then looks like this.

sin(x)≈∑n=010((−1)n(2n+1)!x2n+1)sin(x) \approx \sum_{n=0}^{10}(\frac{(-1)^n}{(2n+1)!}x^{2n+1}) cos(x)≈∑n=010((−1)n(2n)!x2n)cos(x) \approx \sum_{n=0}^{10}(\frac{(-1)^n}{(2n)!}x^{2n})

The two compact notations can now be implemented easily. Important !! is factorial and means the multiplication from 1 to counted up nn. Or simpler:

n=3→3!=1⋅2⋅3=6n = 3 \rightarrow 3! = 1\cdot2\cdot3 = 6

The sigma, or simpler sum symbol, can be built as a research loop that goes from n=0n=0 to 1010. As an example:

∑n=010a⇔\sum_{n=0}^{10} a \Leftrightarrow for(int n = 0; n < 10; n++){ ... }

And we can start already 😄. First we should implement the factorial to calculate the denominator later in the power series. Here one should take a long long int since a factorial grows very fast and otherwise it comes to overflows. Nevertheless we have to think about what 0!0! and 1!1! should be. The best would be that both deliver the value 1, because 0 makes no sense.

long long int factorial(unsigned n)
{
    long long int f = 1;
    for(int i = 1; i <= n; i++)
        f *= i;

    return f;
}

Here in the snippet we have implemented the factorial and can test this e.g. with an online calculator, whether the right comes out.

After we have implemented the factorial, we can actually start with the implementation of sine and cosine. The implementation itself is actually just a copy of the power series. You can see how the sigma or sum symbol is implemented, as a for-loop. At the end, to get a total sum, a variable y was created, which is summed up further and further. The expression in the sum symbol can be written as a one liner.

double my_sin(double x)
{
    double y = 0.0;
    for(int n = 0; n < 10; n++)
    {
        y += (pow(-1, n) / factorial((2*n)+1)) * pow(x, (2*n)+1);
    }

    return y;
}

From the principle quite simple. Now we can implement cosine.

double my_cos(double x)
{
    double y = 0.0; 
    for(int n = 0; n < 10; n++)
    {
        y += (pow(-1, n) / factorial(2*n))*pow(x, 2*n);
    }

    return y;
}

That was not really difficult either. And now it's time to test.

int main()
{
     for(int i = -20; i < 20; i++)
     {
         std::cout << "sin(" << i << ") = " << my_sin(double(i)) << 
                      "\t\tcos(" << i << ") = " << my_cos(double(i)) 
                   << std::endl;
     }

     return 0;
}

If you compile the whole thing and start it, you can see that it runs from [−20;20][-20; 20] along the X-axis and at each position xx the corresponding yy is output from the function.

sin(-20) = 2.21005e+07		cos(-20) = -2.2134e+07
sin(-19) = 7.89887e+06		cos(-19) = -8.35082e+06
sin(-18) = 2.66196e+06		cos(-18) = -2.97936e+06
sin(-17) = 840111		cos(-17) = -998616
sin(-16) = 246299		cos(-16) = -312036
sin(-15) = 66435.5		cos(-15) = -90063.8
sin(-14) = 16297		cos(-14) = -23747.7
sin(-13) = 3584.88		cos(-13) = -5643.05
sin(-12) = 695.555		cos(-12) = -1188.17
sin(-11) = 117.147		cos(-11) = -217.429
sin(-10) = 16.8119		cos(-10) = -34.4386
sin(-9) = 1.42813		cos(-9) = -5.1463
sin(-8) = -0.829433		cos(-8) = -0.560662
sin(-7) = -0.647032		cos(-7) = 0.724297
sin(-6) = 0.279816		cos(-6) = 0.958777
sin(-5) = 0.958933		cos(-5) = 0.283625
sin(-4) = 0.756803		cos(-4) = -0.653644
sin(-3) = -0.14112		cos(-3) = -0.989992
sin(-2) = -0.909297		cos(-2) = -0.416147
sin(-1) = -0.841471		cos(-1) = 0.540302
sin(0) = 0		cos(0) = 1
sin(1) = 0.841471		cos(1) = 0.540302
sin(2) = 0.909297		cos(2) = -0.416147
sin(3) = 0.14112		cos(3) = -0.989992
sin(4) = -0.756803		cos(4) = -0.653644
sin(5) = -0.958933		cos(5) = 0.283625
sin(6) = -0.279816		cos(6) = 0.958777
sin(7) = 0.647032		cos(7) = 0.724297
sin(8) = 0.829433		cos(8) = -0.560662
sin(9) = -1.42813		cos(9) = -5.1463
sin(10) = -16.8119		cos(10) = -34.4386
sin(11) = -117.147		cos(11) = -217.429
sin(12) = -695.555		cos(12) = -1188.17
sin(13) = -3584.88		cos(13) = -5643.05
sin(14) = -16297		cos(14) = -23747.7
sin(15) = -66435.5		cos(15) = -90063.8
sin(16) = -246299		cos(16) = -312036
sin(17) = -840111		cos(17) = -998616
sin(18) = -2.66196e+06		cos(18) = -2.97936e+06
sin(19) = -7.89887e+06		cos(19) = -8.35082e+06

Here you can see that the first 3 negative/positive values are correct but after that, everything is neither sine nor cosine.

I said to begin that sine and cosine are periodic, but we haven't implemented that yet. A period goes over 2π2\pi which is about 6.3...6.3... is. Here we can say that everything greater than 2π2\pi must be corrected with −2π-2\pi. But what about negative numbers? Likewise here, anything less than −2π-2\pi must be incremented by 2π. This can then be implemented as a While loop.

while(x > 2*M_PI)
     x -=2*M_PI;

while(x < -2*M_PI)
     x +=2*M_PI;

We can insert these two loops before the summation and should now arrive at the correct values. But here it is better, since the same loops come in with both functions anyway, to lay out these loops in a separate function.

double correct_value(double x)
{
    while(x > 2*M_PI)
        x -=2*M_PI;

    while(x < -2*M_PI)
        x +=2*M_PI;

    return x;
}

Now we just have to check and adjust the values before summing. If we have implemented this now, the code should look like this.

#include <iostream>
#include <math.h>

long long int factorial(unsigned n)
{
    long long int f = 1;
    for(int i = 1; i <= n; i++)
        f *= i;

    return f;
}

double correct_value(double x)
{
    while(x > 2*M_PI)
        x -=2*M_PI;

    while(x < -2*M_PI)
        x +=2*M_PI;

    return x;
}

double my_sin(double x)
{
    double y = 0.0;
    x = correct_value(x);
    for(int n = 0; n < 10; n++)
    {
        y += (pow(-1, n) / factorial((2*n)+1)) * pow(x, (2*n)+1);
    }

    return y;
}

double my_cos(double x)
{
    double y = 0.0;
    x = correct_value(x);
    for(int n = 0; n < 10; n++)
    {
        y += (pow(-1, n) / factorial(2*n))*pow(x, 2*n);
    }

    return y;
}

int main()
{
    for(int i = -20; i < 20; i++)
    {
        std::cout << "sin(" << i << ") = " << my_sin(double(i)) << 
              "\t\tcos(" << i << ") = " << my_cos(double(i)) 
           << std::endl;
    }

    return 0;
}

If you run the program, our self-implemented sine and cosine functions should now display the correct values.

sin(-20) = -0.912945		cos(-20) = 0.408082
sin(-19) = -0.149877		cos(-19) = 0.988705
sin(-18) = 0.751038		cos(-18) = 0.660122
sin(-17) = 0.961398		cos(-17) = -0.275167
sin(-16) = 0.287903		cos(-16) = -0.95766
sin(-15) = -0.650288		cos(-15) = -0.759688
sin(-14) = -0.990607		cos(-14) = 0.136737
sin(-13) = -0.420167		cos(-13) = 0.907447
sin(-12) = 0.536719		cos(-12) = 0.843321
sin(-11) = 0.999993		cos(-11) = 0.00441405
sin(-10) = 0.544021		cos(-10) = -0.839072
sin(-9) = -0.412118		cos(-9) = -0.91113
sin(-8) = -0.989358		cos(-8) = -0.1455
sin(-7) = -0.656987		cos(-7) = 0.753902
sin(-6) = 0.279816		cos(-6) = 0.958777
sin(-5) = 0.958933		cos(-5) = 0.283625
sin(-4) = 0.756803		cos(-4) = -0.653644
sin(-3) = -0.14112		cos(-3) = -0.989992
sin(-2) = -0.909297		cos(-2) = -0.416147
sin(-1) = -0.841471		cos(-1) = 0.540302
sin(0) = 0		cos(0) = 1
sin(1) = 0.841471		cos(1) = 0.540302
sin(2) = 0.909297		cos(2) = -0.416147
sin(3) = 0.14112		cos(3) = -0.989992
sin(4) = -0.756803		cos(4) = -0.653644
sin(5) = -0.958933		cos(5) = 0.283625
sin(6) = -0.279816		cos(6) = 0.958777
sin(7) = 0.656987		cos(7) = 0.753902
sin(8) = 0.989358		cos(8) = -0.1455
sin(9) = 0.412118		cos(9) = -0.91113
sin(10) = -0.544021		cos(10) = -0.839072
sin(11) = -0.999993		cos(11) = 0.00441405
sin(12) = -0.536719		cos(12) = 0.843321
sin(13) = 0.420167		cos(13) = 0.907447
sin(14) = 0.990607		cos(14) = 0.136737
sin(15) = 0.650288		cos(15) = -0.759688
sin(16) = -0.287903		cos(16) = -0.95766
sin(17) = -0.961398		cos(17) = -0.275167
sin(18) = -0.751038		cos(18) = 0.660122
sin(19) = 0.149877		cos(19) = 0.988705

This already looks right, but we should query a few distinctive points. Here the points from the value table from the beginning offer themselves.

 sin(0) = 0                      cos(0) = 1
 sin(1.5708) = 1                 cos(1.5708) = -3.43243e-15
 sin(3.14159) = -5.28918e-10     cos(3.14159) = -1
 sin(4.71239) = -1               cos(4.71239) = -1.14329e-05
 sin(6.28319) = -0.00104818      cos(6.28319) = 0.996521
 sin(0.785398) = 0.707107        cos(0.785398) = 0.707107
 sin(3.92699) = -0.707107        cos(3.92699) = -0.707107

So we see, sine and cosine can be easily implemented using power series. I have uploaded the code on Github.

I hope you enjoyed it, see you later 😄