A couple weeks ago, I was bored and decided to try and make a sine function. I had recently learned about the Taylor Series approximations of sine. And if you don’t know, this is how it looks :

Obviously, you can’t do an infinite sum on a computer, because it will run forever; which is why it is only an approximation. But with the floating-point nature of computers, any non-perfect calculation is an approximation anyways, so it doesn’t matter.

So, first, I did what I do when trying to solve any problem, I plugged it into Desmos.

And, immediately, you can see a very obvious problem.

It has a very dramatic falloff. After some experimentation, the more iterations you have in the sum, the less the falloff is, and the more accurate the entire approximation is.

Very quickly, I found a fix for the fallout. I simply modulus-ed x by 2π. And, it results in a very nice sine function.

 

Are we done? Nope, not even close, The next step is to code it.

And after very simply coding it with 20 iterations, it runs much, much slower than Math.sin – so this will not work.

The way to reduce the time it takes is to reduce the iterations. The problem is, that the more you reduce the iterations, the lower the accuracy of the approximation, and the more it falls off.

So, playing with a couple numbers, the lowest we can hope to get it, is around 14 iterations.

This is where my next big realization came, the way the Taylor Series works, is that it falls off the further you are from 0, so it’s much more accurate from, let’s say [0, π/2] then [3π/2, 2π].

The reason this is important is that we can use the values from [0, π/2] for the other 3 quarters of the function. Hopefully, this gif will explain what I mean.

With this realization, it means that we can lower the iterations until it stops being accurate from [0, π/2], I used 7 as a very good mix between accuracy and speed (accurate within 10 decimal places).

 

So our task is to have θ be between 0 and π/2 and still return the value of its quadrant.

Let’s approach this quadrant by quadrant.

In the first quadrant, θ is already between 0 and π/2, so we don’t have to do anything there.

In the second quadrant, θ is between π/2 and π, so the simplest way to get it between 0 and π/2 would be to just subtract π/2. The obvious problem with this is that all the values are flipped, as we’re essentially rotating the unit circle π/2 radians clockwise.

So to fix this problem, we can just apply a property of sines that sin(π – θ) = sin(θ)

Which is perfect for us because, in the 2nd quadrant, π – θ puts θ between 0 and π/2 and returns the correct value. As illustrated in this diagram :

 

Now, for the third quadrant, θ is between π and 3π/2, so we can get it into 0 and π/2 by simply subtracting π. The only problem here is that it’ll be positive, so we just have to make negative.

And this works as -sin(θ – π) = sin(θ)

 

Finally, for the fourth quadrant, θ is between 3π/2 and 2π so we can subtract 3π/2 to get it between 0 and 2π, but there are two problems. First, it has the same problem as quadrant 2 where the values are inverted, and secondly, the value is negative and not positive.

To fix this, we can simply subtract 2π by θ and then make it negative afterwards.

So this works as -sin(2π – θ) = sin(θ)

 

Now with all of that logic out of the way, there’s one more thing I have to discuss, then we’ll be ready to code.

While experimenting, I found that Java’s modulus was incredibly slow, so in the final code, you won’t see modulus anywhere, but you will see this :

theta -= ((int)(theta * INV_TWO_PI) * TWO_PI);

This is my “modulus”, it works for any value and is probably faster than Java’s because we don’t have to divide but just multiply.

And finally, this is the code

public class FastSine {
    private static final double PI = Math.PI;
    private static final double TWO_PI = 2*PI;
    private static final double INV_TWO_PI = 1.0 / TWO_PI;

    private static final double[] INV_FACTORIALS = { 1/6.0, 1/120.0, 1/5040.0, 1/362880.0, 1/39916800.0, 6227020800.0 };

    public static double sin(double theta) {
        theta -= ((int)(theta * INV_TWO_PI) * TWO_PI);
        if (theta < 0) theta += TWO_PI;
        
        if (theta <= HALF_PI)
            return taylorSeriesSin(theta); // Just plug in theta if it's in the 1st quadrant
        else if (theta <= PI)
            return taylorSeriesSin(PI - theta); // In the 2nd quadrant, sin(PI - theta) = sin(theta)
        else if (theta <= ONE_HALF_PI)
            return -taylorSeriesSin(theta - PI); // In the 3rd quadrant, sin(theta - PI) = -sin(theta)
        else
            return -taylorSeriesSin(TWO_PI - theta); // In the 4th quadrant, sin(2PI - theta) = -sin(theta)
    }

    private static double taylorSeriesSin(double x) {
        double x2 = x*x;
        return x*(1 - x2*(INV_FACTORIALS[0] - x2*(INV_FACTORIALS[1] - x2*(INV_FACTORIALS[2]
                    - x2*(INV_FACTORIALS[3] - x2*(INV_FACTORIALS[4] - x2*(INV_FACTORIALS[5])))))));
    }
}