Sine function in C / C++

Due to software constraints, I cannot use the standard libraries, cmath, algorithm, templates, inline, or boost. I am also using standard C (ISO C99) such that array is not a reserved keyword like it is in Visual Studio.

Here is the custom sine, pow and factorial functions I wrote (assume they are all part of the same namespace and not using cmath). I looking for suggestions to improve robustness and efficiency. For example, both the pow and factorial are constrained by the int i in the for loops.

double factorial(double x)
{
double result = 1;
for (int i=1;i<=x;i++)
result *= i;
return result;
}

double pow(double x,double y)
{
double result = 1;
for (int i=0;i<y;i++)
result *= x;
return result;
}

double sine_taylor(double x,double n)
{
double sine = x;
bool isNeg;
isNeg = true;
for (double i=3;i<=n;i+=2)
{
if(isNeg)
{
sine -= pow(x,i)/factorial(i);
isNeg = false;
}
else
{
sine += pow(x,i)/factorial(i);
isNeg = true;
}
}
return sine;
}

int main()
{
const double pi = 3.14159265;
double sineTaylor = sine_taylor(pi/2,7);
}
• Explanation for -1? I don't see a duplicate question, code formatted correctly, and problem statement clearly defined. Oct 6 '11 at 18:58
• @Ant, I didn't think I had to request for code review at Code Review, anyway I updated looking to increase robustness and efficiency. Oct 6 '11 at 19:31
• These are solved problems. Crack open the standard library and copy and paste the code you need. Oct 7 '11 at 16:49
• There is no point to the question (because it is a well know SOLVED problem). Asking it is silly. There is a difference between not being able to use the standard library and copying code (which just happens to come from the standard library presumably you can still get a copy to look at). Oct 7 '11 at 16:55
• @Sjoerd: If you have an FPU. see: stackoverflow.com/questions/2284860/… Oct 25 '11 at 5:34

There are books full on how to calculate those functions efficiently and accurately. So it's too much for a short answer here.

I'm hardly an expert, but I can easily spot that your loop starts at the large terms and ends with the small terms. This has a larger accumulated error than going from small to large.

Secondly, subtracting is a very good way to lose a lot of accuracy. It might be better to partially unroll the loop, combining each positive term with the next negative term, and only add the difference between those two. Although that still might result in a large loss of accuracy.

As none of the other answers even touches on those point, you have not been answered by an expert yet. I only know enough to know that I don't know enough.

So clearly work for specialists. And indeed, solved a long time ago.

A quick google found What Every Computer Scientist Should Know About Floating-Point Arithmic.

Note – This appendix is an edited reprint of the paper What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, published in the March, 1991 issue of Computing Surveys. Copyright 1991, Association for Computing Machinery, Inc., reprinted by permission.

I just googled for that title, betting that someone would have used that title at some point in time. It explains the problem I mentioned above and much more, but I'm not really qualified to judge it. As it is SUN's SPARC documentation, I expect it to be among the best you can find online, although it could be too technical for some.

It's the appendix of some part of SUN's documentation for SPARC. Scanning the full document, In the references is this paper:

Tang, Peter Ping Tak, Some Software Implementations of the Functions Sin and Cos, Technical Report ANL-90/3, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Illinois, February 1990.

I can't find it online using google, so good luck hunting that one!

Drawback of the backwards loop is that one can't use the result of previous step. An even better approach might be a forward loop to calculate the terms. But don't add them, store them in an array, then use a backwards loop to sum them.

That might sound a lot of work, but when the size of the array is fixed at compile time, one can unroll the loop (see Artur Mustafin's answer), and just use 4 local doubles (each one holds the diference between two terms, so that's the same 8 terms he calculates).

Based on my remarks, here is my version. Inlined, loop-unrolled, and parallizable. I usually don't bother with re-using local variables: the compilers nowadays are smart enough to do that for me. I also assume that e.g. / (10*11) is optimized into a * some_constant where some_constant equals 1.0 / 110, assuming the CPU is faster on multiplications than on divisions.

But I bet that the version used by real libraries is even smarter.

double sine_taylor(double x)
{
// useful to pre-calculate
double x2 = x*x;
double x4 = x2*x2;

// Calculate the terms
// As long as abs(x) < sqrt(6), which is 2.45, all terms will be positive.
// Values outside this range should be reduced to [-pi/2, pi/2] anyway for accuracy.
// Some care has to be given to the factorials.
// They can be pre-calculated by the compiler,
// but the value for the higher ones will exceed the storage capacity of int.
// so force the compiler to use unsigned long longs (if available) or doubles.
double t1 = x * (1.0 - x2 / (2*3));
double x5 = x * x4;
double t2 = x5 * (1.0 - x2 / (6*7)) / (1.0* 2*3*4*5);
double x9 = x5 * x4;
double t3 = x9 * (1.0 - x2 / (10*11)) / (1.0* 2*3*4*5*6*7*8*9);
double x13 = x9 * x4;
double t4 = x13 * (1.0 - x2 / (14*15)) / (1.0* 2*3*4*5*6*7*8*9*10*11*12*13);
// But remember that x is smaller than 2, and the factorial grows very fast
// so I doubt that 2^17 / 17! will add anything.
// Even t4 might already be too small to matter when compared with t1.

// Sum backwards
double result = t4;
result += t3;
result += t2;
result += t1;

return result;
}

Not being able to use the standard libraries suggests you are doing embedded work. A more typical approach for calculating sine in that situation would be to:

• Translate angles from a 360 degree system to a 256 degree system;
• Store sine values for all 256 possible degrees in a lookup table.

You lose a lot of accuracy, but you gain a lot of speed.

• Accuracy is very important and this is not necessarily on an embedded system. Oct 6 '11 at 19:32
• @Elpezmuerto: Using lookup table for 256 values you can interpolate Sin functions with very high level of accuracy using well known triginometrical equations and, if you accuracy level is really tremendous you can just use the biiger lookup table (65536 values) Oct 7 '11 at 4:04

The factorial() function is an integral function. It's arguments and return type should be integers, not doubles. By using floating-point arithmetic here, it is also probably not as fast as it should be. Also, your loop may start at 2, save yourself the iteration. :)

There's no sense in making y in your pow() function a double. They should be typed based on how you use it in this case so it is absolutely clear. It would also be a good idea to rename it as it is not a standard pow() function (that only works for positive integer exponents). That assumes you're not intending to expose it for use outside your module. If so, you should probably make it static as well (same for factorial()).

Also, the similar arguments with your sine_taylor() function. There's no sense in making your loop variable a double, it should be an integer anyway. With the above two fixes in place, you can easily fix this.

Rather than looping through all the terms, alternating between positive and negative terms, it might be a good idea to help the compiler and processor out here and split these up as separate loops to remove the branches. That way the compiler could possibly unroll this better and you're not potentially confusing the branch predictors.

Try it out like this:

static int factorial(int n)
{
int i;
int result = 1.0;
for (i = 2; i <= n; i++)
result *= i;
return result;
}

static double pow_int(double x, int y)
{
int i;
double result = 1.0;
for (i = 0; i < y; i++)
result *= x;
return result;
}

double sine_taylor(double x, int n)
{
int i;
double result = x;
for (i = 3; i <= n; i += 4)
result += pow_int(x, i) / factorial(i);
for (i = 5; i <= n; i += 4)
result -= pow_int(x, i) / factorial(i);
return result;
}
• both the power and the factorial are needlessly computed from scratch each iteration; you should get significant performance gains if you inline both Oct 7 '11 at 0:31
• There's plenty of room for more optimizations. Memoization, alternative algorithms, lookup tables, etc. My main focus here was to address other aspects of the existing code without changing the overall structure too much. Most of those changes are not exactly trivial to implement given the restrictions, may or may not be worth it. Oct 7 '11 at 0:46
• What about the inline version? I suppose if you can't optimize something you'd better not doing it, because you will not teach something really interesting, for example, you need to analyze input algorithms and data before you can have a chance to refactor the code. I think, every time you post a message, you have to provide really good examples, not duplicate common student's faults and bad ideas, bla-bla-something about it. If you know somethig, you can show, how-to, otherwize you'd better keep it within yourself. Oct 7 '11 at 7:50
• Inline is not allowed in my application (software certification issue) Oct 7 '11 at 14:14
• @JeffMercado, I agree with @ArturMustafin, feel free to change other aspects of the code. Its better to post the optimized version within the constraints of the problem (no inline) and teach why it is better. Oct 7 '11 at 16:54

This is the correct answer. By providing the Pi number with accuracy 1e-8, you give me an opportunity to optimize you function to inline:

{
double x2 = x*x;
double f = 1;
double s = 0;
int i = 1;

s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;

return s;
}

I hope you have done you homework well, as I did in my past:

#include <math.h>

inline double _sin(double radians, double epsilon, int &count )
{
double x2 = x*x;
double f = 1;
double s = 0;
int i = 1;

while (x/f >= epsilon)
{
s += x/f; x *= x2; f *= ++i; f *= ++i;
s -= x/f; x *= x2; f *= ++i; f *= ++i;
count++;
}

return s;
}

int main()
{
double e_x = 0.00000001; // 1e-8, because pi differs from real pi at 8 symbol after decimal point

double e_pi = 3.14159265;
double m_pi = M_PI;

int n = 10;

double e_dx = e_pi/(2.0*n);
double m_dx = m_pi/(2.0*n);

for (double e_i=e_pi/2.0, m_i=m_pi; n>=0; e_i-=e_dx, m_i-=m_dx, n--)
{
int c = 0;
double e_f = _sin(e_i, e_x, c);
double m_f = sin(m_i);
printf("e_pi=%.16e e_f=%.16e c2=%d\n", e_pi, e_f, c);
printf("m_pi=%.16e m_f=%.16e\n", m_pi, m_f);
}
// by the experiment, for the given epsion, we can easily take double _sin8(double) function, because:
// e(f(x)) >= e(x), for any given f(x) = a0*x^0 + a1*x^1 + ... that's called "accumuilation of errors"
// in my opinion, you'd better go to scool library, i do not think it is good for you to
// ask community to do your homework. as you see, i got bachelor degree in CS years ago,
// so we know something about it, and you'd better to discover this knowledge yourself
}

Some mathiematical constants, to remember:

/* Definitions of useful mathematical constants
* M_E        - e
* M_LOG2E    - log2(e)
* M_LOG10E   - log10(e)
* M_LN2      - ln(2)
* M_LN10     - ln(10)
* M_PI       - pi
* M_PI_2     - pi/2
* M_PI_4     - pi/4
* M_1_PI     - 1/pi
* M_2_PI     - 2/pi
* M_2_SQRTPI - 2/sqrt(pi)
* M_SQRT2    - sqrt(2)
* M_SQRT1_2  - 1/sqrt(2)
*/

#define M_E        2.71828182845904523536
#define M_LOG2E    1.44269504088896340736
#define M_LOG10E   0.434294481903251827651
#define M_LN2      0.693147180559945309417
#define M_LN10     2.30258509299404568402
#define M_PI       3.14159265358979323846
#define M_PI_2     1.57079632679489661923
#define M_PI_4     0.785398163397448309616
#define M_1_PI     0.318309886183790671538
#define M_2_PI     0.636619772367581343076
#define M_2_SQRTPI 1.12837916709551257390
#define M_SQRT2    1.41421356237309504880
#define M_SQRT1_2  0.707106781186547524401
• I cannot use inline unfortunately Oct 7 '11 at 14:16

factorial begs to be an array (1,2,6,24,120...) pow can be optimized a bit... 3^101= 3^64*3^32*3^4*3^1= 3^1*3^32*...* by using exponentiation by squaring.

Problem is that you don't have int as power but still you can to the "int specialization" that way. Regarding sin you could manually unroll the for loop and make sure that you have some parallel processing (not multiple threads but that maximum amount of instruction level parallelism is possible). This is highly speculative but let's say that compiler might not be smart enough to optimize

x=(a+b)+(c  +d)+(e+f);`

in a way that is is faster than

x=a+b+c+d+e+f;

so you might want to that explicitly for the taylor sum.

But I would guess that compiler will unroll and do the transformation automatically. Also it might be the case that floating point multiplication is faster than division so you might want to do something like having 1/1,1/2,1/6... array for taylor calculation.