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 n1 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);
// Switch between adding/subtracting
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);
}
3 answers

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 n1 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); // Switch between adding/subtracting 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/taylorsine.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 x_{0} is:
where R_{n} is the Lagrange Remainder
Note that R_{n} grows fast as soon as x moves away from the center x_{0}.
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 yoursine()
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()
frommath.h
, so you can do:#include <math.h> // Your code 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 ownfmod()
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); }
Try it online!