Sine function in C / C++

Question

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);
}

Answer 1

UPDATE: (19-March-2015) This answer sounds like an expert answer.

Original answer:

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);
    // add some more if your accuracy requires them.
    // 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;
}

Answer 2

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.

Answer 3

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;
}

Answer 4

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:

inline double _sin8(double radians)
{ 
    double x = radians;
    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>

double _sin8(double radians);   

inline double _sin(double radians, double epsilon, int &count )
{ 
    double x = radians;
    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

Answer 5

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.

If you want to read more about optimization:
http://www.agner.org/optimize/