Skip to content

Instantly share code, notes, and snippets.

@giangnguyen2412
Created May 16, 2018 09:50
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save giangnguyen2412/bcab883b5a53b437b980d7be9745beaf to your computer and use it in GitHub Desktop.
Save giangnguyen2412/bcab883b5a53b437b980d7be9745beaf to your computer and use it in GitHub Desktop.
sine and cosine function implementation in C using Taylor series
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
int compare_float(double f1, double f2)
{
double precision = 0.00000000000000000001;
if ((f1 - precision) < f2)
{
return -1;
}
else if ((f1 + precision) > f2)
{
return 1;
}
else
{
return 0;
}
}
double cos(double x){
if( x < 0.0f )
x = -x;
if (0 <= compare_float(x,M_PI_M_2))
{
do {
x -= M_PI_M_2;
}while(0 <= compare_float(x,M_PI_M_2));
}
if ((0 <= compare_float(x, M_PI)) && (-1 == compare_float(x, M_PI_M_2)))
{
x -= M_PI;
return ((-1)*(1.0f - (x*x/2.0f)*( 1.0f - (x*x/12.0f) * ( 1.0f - (x*x/30.0f) * (1.0f - (x*x/56.0f )*(1.0f - (x*x/90.0f)*(1.0f - (x*x/132.0f)*(1.0f - (x*x/182.0f)))))))));
}
return 1.0f - (x*x/2.0f)*( 1.0f - (x*x/12.0f) * ( 1.0f - (x*x/30.0f) * (1.0f - (x*x/56.0f )*(1.0f - (x*x/90.0f)*(1.0f - (x*x/132.0f)*(1.0f - (x*x/182.0f)))))));
}
double sin(double x){return cos(x-M_PI_2);}
@giangnguyen2412
Copy link
Author

As we have known that, in embedded, math library is sometimes unavailable and unreasonable. Hence, the necessity of diagonal functions such as sine or cosine which are portable, fast and accurate is quite high.

So I refer from this to get the formula:
http://www.wolframalpha.com/widgets/view.jsp?id=f9476968629e1163bd4a3ba839d60925

In C, subsequently, I implemented these functions which may help you.

@Lewiscowles1986
Copy link

Lewiscowles1986 commented Mar 3, 2020

What is the value and type of M_PI_M_2 without endian-ness (literal value), and which compiler did you use to compile this?

https://pubs.opengroup.org/onlinepubs/7908799/xsh/math.h.html suggests it's just M_PI_2 and represents a double-precision of Value of M_PI/2

Although you also use that, so could it be M_2_PI which is 2/M_PI

@Hachem-H
Copy link

Sorry I know this is quite late but I am going to guess it might be pi/2, as cos(x) = 0; with x = n*pi/2, and since it is a periodic function, I am assuming he is just decreasing all factors of n in order for x to be between 0-1 (he turns the negative branch positives).

@Lewiscowles1986
Copy link

Sorry what is a periodic function? M_PI_M_2 reads to be a constant

@robertdetian
Copy link

robertdetian commented Sep 1, 2022

#define M_PI (3.1415927f)
#define M_PI_2 (M_PI/2.0f)
#define M_PI_M_2 (M_PI*2.0f)

@Lewiscowles1986
Copy link

Thanks very nice. Was this taken from a specific libraries headers?

@giangnguyen2412
Copy link
Author

@Lewiscowles1986 No i wrote it myself based on Taylor's series.

@Lewiscowles1986
Copy link

Lewiscowles1986 commented Sep 4, 2022

@giangnguyen2412 nice to see you finally. Apologies, the question was assumed as directed at @robertdetian unless these are both your github usernames, or you have an answer with regards the defines in https://gist.github.com/giangnguyen2412/bcab883b5a53b437b980d7be9745beaf?permalink_comment_id=4286644#gistcomment-4286644

@giangnguyen2412
Copy link
Author

@Lewiscowles1986 Sorry for my (super) late reply. No we are not :D But I think that definition works unless you want the results more accurate (i.e. extending the numbers after the floating point).

@robertdetian
Copy link

robertdetian commented Sep 5, 2022

Thanks @giangnguyen2412 for this example.

I managed to implement the float port cosine, and rewrote above code as following, which can make Taylor's series more scalable, although 7 times works well for float.

float c_cos(float x)
{
const float pi = 3.1415927;

x = std::copysign(x, 1.0f);

auto pi_count = (uint32_t)(x / pi);
x -= (float)pi_count * pi;

float taylor_series = 1.0;
for (int i = 14; i > 0; i-=2)
{
    taylor_series = 1.0f - (x*x)/(float)(i*(i-1)) * taylor_series;
}

uint32_t n_one_f = 0xBF800000;
return (pi_count & 1) ? taylor_series * *(float*)&n_one_f : taylor_series;

}

Ref: https://en.wikipedia.org/wiki/Sine_and_cosine

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment