How approximate sine using Taylor series

I need to approximate the sine function without internal libraries. I used Taylor series in 0 to solve this, but my program works for some values, but for others awful results. The program gets x value, unit (degrees or radians) and how many words we want use for the approximation. For values like

1 rad, 5 words results are ok: sin(1 pi) = 0.833333 Difference: 0.00813765;

sin(0.5 pi) = 0.479167 Difference: 0.000258872

But for values like 125.6 rad it returns 2.60145e+008 and C++ sin() function gives -0.0636631, so the result is completely unreal.

Also it doesn't work for values in degrees although I convert degrees on radians. E.g. sin(4.71239 pi) = 6.63667 and real result should be -1. I completely don't have any idea why it doesn't work properly, could anyone tell me what could cause these problems?

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

using namespace std;

double degConverter(double);
void approximate(double, int, int);
int derivative(int);
double compareWithNorm(double, double);
double inline reductionFormula(double);

double inline degConverter(double x)
{
return x / 180.0; //x * M_PI / 180.0    <-- the expression in comment gives wrong results in all situations?
}

void approximate(double x, int unit, int words)
{
double result = 0;
double xVal = 0;
int factorial = 1;

if (unit == 1)
{
xVal = x;
}
else if (unit == 0)
{
xVal = degConverter(x);
x = xVal;
}

if (words > 1)
{
for (int i = 1; i < words; i++)
{
switch (i % 4)
{
case 1:
result += xVal / (double)factorial;
break;

case 3:
result -= xVal / (double)factorial;
break;

}
xVal *= x;
silnia *= i + 1;
}
}

//Comparing with sin() function from math.h lib
cout << "Difference between sin C++ and Taylor: " << compareWithNorm(x, result) << endl;
}

double compareWithNorm(double x, double result)
{
double diff = 0;
double pattern = sin(x);
diff = pattern - result;
if (diff < 0)
diff *= -1;

return diff;

}

double inline reductionFormula(double x)
{

if (x > 2)
{
x = fmod(x, 2.0);
}
return x;
}

int main()
{
double x;
int unit = 1;
int words = 1;

cout << "Type x [degrees / rad]: ";
cin >> x;
cout << "\nType unit [deg = 0 / rad = 1]: ";
cin >> unit;
cout << "\nHow many words?: ";
cin >> words;
cout << endl;

approximate(x, unit, words);

return 0;
}
• Welcome to Computer Science! Please get rid of the source code and replace it with ideas, pseudo code and arguments of correctness. See here and here for related meta discussions. Also, code debugging is offtopic anywhere on SE; do you have a focussed question about your algorithm? – Raphael Mar 12 '18 at 22:03

Others have talked about argument reduction, but just as a little addendum, designing good argument reduction algorithms is an art in itself. General-purpose library functions have to be "correct" (in the sense that the answer should be within 0.5ulp of the arbitrary-precision correct answer) on the entire range of floating-point numbers, and so will use a lot of special cases for efficiency and accuracy.

For example, one common technique is to let the floating point addition operation round for you. In C:

/* 1/2pi */

const double one_over_period = 1.591549430918953456e-1;

/* This is n = round(x / 2pi). */

const double round_to_integer = 1.5 / DBL_EPSILON;
double n = (x * one_over_period + round_to_integer) - round_to_integer;

/* This is x = x - n * 2*pi, using Cody-Waite reduction.
Don't worry, this will be explained in a moment. The
important part for now is that period_1 + period_2 is
an approximation to 2pi. */

const double period1 = 6.28318524360656738280e+0;
const double period2 = 6.35730190941127867664e-8;
x = x - n * period1 - n * period2;

The argument reduction algorithm works for quite large abs(x); calculating sin(1e8) is no problem. That may be good enough for you, but not for your programming language's standard library, which has to give a 0.5ulp-correct answer for sin(1e22). So in a general-purpose implementation, you'd check the range of the argument and use different argument-reduction code as appropriate.

Now let's talk about that Cody-Waite method.

You want to compute $x - nK$ where $n = \left\lfloor \frac{x}{K}\right \rfloor$. We'll ignore how you get $n$; that's the code above. What we're interested in is computing $x - nK$.

The constant $K$ is known to arbitrary precision, but of course you want to compute this using only floating-point arithmetic, and you want the result to be accurate (within 0.5ulp of the true answer).

The problem with doing this the obvious way is that when $x$ is close to a multiple of $K$, doing the obvious subtraction can cause catastrophic cancellation, since you are subtracting two almost-equal numbers.

The solution is to find two constants, $C_1$ and $C_2$, with the following properties:

1. $C_1$ is very close to $K$.
2. $C_1$ is representable using only a few significant digits.
3. $C_1 + C_2 \approx K$, to far enough beyond the working precision of the floating point representation.

Then calculate $(x - n C_1) - n C_2$.

Property 2 means that $n C_1$ is exactly representable for a large range of integers $n$. In the case above, the first number period_1 is $\mathrm{1.921FB50000000_{16}\times 2^{2}}$. All those zeroes at the end of the mantissa mean that you can multiply this by a very large integer indeed and still get a number that is represented exactly in floating point.

(You will note that the terms "a few significant digits" in property 2 and "far enough" in property 3 are vague. We can make these precise, but this is just to give you the idea. Using roughly half the mantissa bits to represent $C_1$ is a good rule of thumb, but it's not a formal proof.)

The way that argument reduction is typically done for larger ranges (e.g. sin(1e22)) is essentially the same as this, only splitting $K$ into more constants. A lot more constants.

• Many processors have a "fused multiply-add" instruction which calculates a +/- b * c with a single rounding error, which makes it a lot easier. Calculate n, find k0 = best approximation to K, find k1 = best approximation to K - k1, then calculate y = x - nk0 with a single rounding error, and z = y - nk1 with a single rounding error. – gnasher729 Mar 13 '18 at 7:10
• Rewording of my previous comment: This approach to argument reduction relies on the fact that the first FMA has no rounding error in the first step (or all but the last step if you use more than two constants). A single ulp of rounding error before that means you don't get at most 0.5ulp error in the final answer. – Pseudonym Mar 14 '18 at 0:56

The Maclaurin series for $\sin x$ (i.e., the Taylor series centred at $x=0$) does converge for all $x$. However, the speed at which it converges depends on the value of $x$ you use. Looking at it from the other side, the Maclaurin series for $\sin x$ truncated at the term in $x^k$ is a polynomial, so it goes to $+\infty$ or $-\infty$ as $x\to\infty$, whereas $\sin x$ remains in the interval $[-1,1]$. So, if you use a fixed truncation of the series as your approximation, the error will grow without bound as $x$ becomes large.

In more detail, because $\sin x$ is infinitely differentiable and the $k$th derivative is bounded between $+1$ and $-1$, the absolute error from truncating the series at the $x^k$ term is at most $x^{k+1}/(k+1)!$. So, if you want to approximate $\sin x$ for large $x$ you need to either take $k$ large enough for the error to be small or, more sensibly, use the fact that the sine function is periodic, so you only need to consider $-\pi<x\leq \pi$. In that case, you're guaranteed that the absolute error is at most $\pi^{16}/16!\approx4\times10^{-6}$ after the $x^{15}$ term (actually, $\pi^{17}/17!$, since the series to $x^{16}$ is the same). And, using the trigonometric identies pointed out by gnasher729, you can use even smaller ranges of $x$ and get even better error bounds from fewer terms.

You use argument reduction. You know that sin (2kπ + x) = sin (x), so given an argument y, you calculate k = round (y / 2π), y' = y - 2kπ, and sin (y) = sin (y'), with y between -π and π.

You can also use sin (π + x) = - sin (x), sin (π/2 + x) = cos (x) to reduce the argument to sin (y') or cos (y') with -π/4 ≤ x ≤ π/4.

The reason why your Taylor series gives absurd results for large x is that you are adding and subtracting huge terms. Mathematically the series converges just fine, but if you have terms > $10^{20}$, these terms will have rounding errors so large that you cannot expect a good result. Reducing the argument to -π/4 ≤ x ≤ π/4 means you have much smaller terms. (Also see Pseudonym's answer how to get the argument reduction as precise as possible).

You are adding the terms of the Taylor series from the largest. Each addition has a rounding error, proportional in size to the result of the addition. So adding in this order each addition gives a large error. It is better to start with the smallest term. The two smallest terms are tiny, adding them gives a tiny rounding error.

And usually we don't actually use the Taylor series. The Taylor series has its largest error at the end of the range. You would use Chebyshev polynomials, which minimise the largest error over the whole interval. And since you don't want any errors at x = 0, you can write sin(x) = x - x * g(x) and minimise the relative error of g(x).