# Creating a custom sine function

I tried to create a custom `sine` function using `c` and the Taylor Series for calculating `sin` with 10 terms in the series, but I'm getting the wrong results when I try to find the `sine(x)` where `x > 6`.

It works well for `-5 < x < 5`, but anything out of that range isn't producing the correct results.

I expect `sin(10)` to return something close to `-0.5440`, but get `1418.0269775391`

I've put everything in a single file so it's easier.

``````#include <stdio.h>
#include <stdlib.h>

double factorial(double n);
double power(double n, double pow);
double sine(double n);

// This is supposed to all go in a .c file and reference the .h stuff above
// This is the actual implementation of the functions declared above
double factorial(double n) {
// 0! = 1 so just return it
if(n == 0) {
return 1;
}
// Recursively call factorial with n-1 until n == 0
return n * (factorial(n - 1));
}

double power(double n, double power) {
double result = n;
// Loop as many times as the power and just multiply itself power amount of times
for(int i = 1; i < power; i++) {
result = n * result;
}
return result;
}

double sine(double n) {
double result = n;
double coefficent = 3; // Increment this by 2 each loop
for(int i = 0; i < 10; i++) { // Change 10 to go out to more/less terms
double pow = power(n, coefficent);
double frac = factorial(coefficent);
printf("Loop %d:\n%2.3f ^ %2.3f = %2.3f\n", i, n, coefficent, pow);
printf("%2.3f! = %2.3f\n", coefficent, frac);

if(i % 2 == 0) { // If the index of the loop is divided by 2, the index is even, so subtract
result = result - (pow/frac); // x - ((x^3)/(3!)) - ((x^5)/(5!))...
} else {
result = result + (pow/frac); // x - ((x^3)/(3!)) + ((x^5)/(5!))...
}
coefficent = coefficent + 2;
printf("Result = %2.3f\n\n", result);
}
return result;
}

// main starting point. This is suppossed to #include "functions.c" which contain the above functions in it
int main(int argc, char** argv) {

double number = atof(argv[1]); // argv[1] = "6"
double sineResult = sine(number);
printf("%1.10f", sineResult);
return (0);
}
``````

after making the corrections, as listed in my comments to the question, the proposed code looks like:

``````#include <stdio.h>
#include <stdlib.h>

double factorial(double n);
double power(double n, double pow);
double sine(double n);

// This is supposed to all go in a .c file and reference the .h stuff above
// This is the actual implementation of the functions declared above
double factorial(double n) {
// 0! = 1 so just return it
if(n == 0) {
return 1;
}
// Recursively call factorial with n-1 until n == 0
return n * (factorial(n - 1));
}

double power(double n, double power) {
double result = n;
// Loop as many times as the power and just multiply itself power amount of times
for(int i = 1; i < power; i++) {
result = n * result;
}
return result;
}

double sine(double n) {
double result = n;
double coefficent = 3.0; // Increment this by 2 each loop

for(int i = 0; i < 10; i++) { // Change 10 to go out to more/less terms
double pow = power(n, coefficent);
double frac = factorial(coefficent);
printf("Loop %d:\n%2.3f ^ %2.3f = %2.3f\n", i, n, coefficent, pow);
printf("%2.3f! = %2.3f\n", coefficent, frac);

if(i % 2 == 0) { // If the index of the loop is divided by 2, the index is even, so subtract
result = result - (pow/frac); // x - ((x^3)/(3!)) - ((x^5)/(5!))...
} else {
result = result + (pow/frac); // x - ((x^3)/(3!)) + ((x^5)/(5!))...
}
coefficent = coefficent + 2;
printf("Result = %2.3f\n\n", result);
}
return result;
}

// main starting point. This is suppossed to #include "functions.c" which contain the above functions in it
int main( void )
{
double number = atof("6");
double sineResult = sine(number);
printf("%1.10f", sineResult);
return (0);
}
``````

and the resulting output looks like:

``````Loop 0:
6.000 ^ 3.000 = 216.000
3.000! = 6.000
Result = -30.000

Loop 1:
6.000 ^ 5.000 = 7776.000
5.000! = 120.000
Result = 34.800

Loop 2:
6.000 ^ 7.000 = 279936.000
7.000! = 5040.000
Result = -20.743

Loop 3:
6.000 ^ 9.000 = 10077696.000
9.000! = 362880.000
Result = 7.029

Loop 4:
6.000 ^ 11.000 = 362797056.000
11.000! = 39916800.000
Result = -2.060

Loop 5:
6.000 ^ 13.000 = 13060694016.000
13.000! = 6227020800.000
Result = 0.037

Loop 6:
6.000 ^ 15.000 = 470184984576.000
15.000! = 1307674368000.000
Result = -0.322

Loop 7:
6.000 ^ 17.000 = 16926659444736.000
17.000! = 355687428096000.000
Result = -0.275

Loop 8:
6.000 ^ 19.000 = 609359740010496.000
19.000! = 121645100408832000.000
Result = -0.280

Loop 9:
6.000 ^ 21.000 = 21936950640377856.000
21.000! = 51090942171709440000.000
Result = -0.279

-0.2793866930
``````

Taylor expansion has an error that depends on the argument scope as well as the order of the Taylor expansion. I believe that you have overreached the boundaries of the argument. See here for more examples: www.dotancohen.com/eng/taylor-sine.php

As I already said in Python: Calculate sine/cosine with a precision of up to 1 million digits

The real Taylor expansion centered in x0 is:

where Rn is the Lagrange Remainder

Note that Rn grows fast as soon as x moves away from the center x0.

Since you are implementing the Maclaurin series (Taylor series centered in 0) and not the general Taylor series, your function will give really wrong results when trying to calculate sin(x) for big values of x.

So before the `for` loop in your `sine()` function you must reduce the domain to at least [-pi, pi]... better if you reduce it to [0, pi] and take advantage of sine's parity.

To fix your code you'll need `fmod()` from `math.h`, so you can do:

``````#include <math.h>

double sine (double n) {
// Define PI
const double my_pi = 3.14159265358979323846;
// Sine's period is 2*PI
n = fmod(n, 2 * my_pi);
// Any negative angle can be brought back
// to it's equivalent positive angle
if (n < 0) {
n = 2 * my_pi - n;
}
// Sine is an odd function...
// let's take advantage of it.
char sign = 1;
if (n > my_pi) {
n -= my_pi;
sign = -1;
}
// Now n is in range [0, PI].

// The rest of your function is fine

return sign * result;
}
``````

Now if you really hate `math.h` module, you can implement your own `fmod()` like this,

``````double fmod(double a, double b)
{
double frac = a / b;
int floor = frac > 0 ? (int)frac : (int)(frac - 0.9999999999999999);
return (a - b * floor);
}
``````