wrong values coming from cos() function

[tquiva]

I need to create my own function to compute the cos of an angle (given in radians) using the Taylor series. The task is to create a function in which a float value is passed to it, and a double gets returned. My program compiled, but I am getting wrong values.

If my input were 2, then I would get 20.987860
The answer I need is cos(2pi) = 2

Here's my code:

#include <stdio.h>
#define DEBUG

float cosi(float x);
int factorial(int x);

main()
{
 /* Declarations */
 double pi = 3.14159;
 float input;
 double rad;
 double cosine;

       /* Ask user for an angle (in radians) */
       printf("Enter an angle in rads (don't include pi): ");

       /* While user has more input */
       while ( scanf("%f", &input) != EOF)
       {
        /* Convert input into its decimal form */
        rad = input * pi;
        #ifdef DEBUG
        printf("debug: rad = %f\n", rad);
        #endif

        /* Call cosine function */
        cosine = cosi(rad);

        /* Print result */
        printf("The cosine value is: %f\n", cosine);

        /* Update loop */
        printf("Enter an angle (in radians): ");

       }

}

float cosi(float x)
/* Given  : an angle in radians (float)
 * Returns: the cosine value of the angle
 */

{
 /*  Declare variables */
 int i,fact=1;
 float sum,term;

 /* Initialize at 1 */
 sum=1.0;
 term=1.0;

       for(i=1;i<5;i++)
       {
        term *= -x * x;
        #ifdef DEBUG
        printf("term = %f\n", term);
        #endif

        fact=factorial(2*i);
        #ifdef DEBUG
        printf("factorial = %d\n", fact);
        #endif

        sum=sum+(term/fact);
        #ifdef DEBUG
        printf("sum = %f\n", sum);
        #endif

       }

 return(sum);
}



int factorial(int x)
/* Given   : a multiple of 2
 * Returns :the value of */

{
 int fact;
       if(x==1)
               return(1);

       else
               fact=x*factorial(x-1);
               return(fact);
}

Does anyone see where I went wrong?

[jonsca]

Line 22 should be rad = input*pi/180; (since pi rad = 180 degrees)

(also cos(2*pi) = 1 not 2)

I didn't check your series approximation but see if that fixes it first.

[jephthah]

Line 22 should be rad = input*pi/180; (since pi rad = 180 degrees)

(also cos(2*pi) = 1 not 2)

pi / 180 is what you'd use if you were making the conversion from degrees. i don't think he is doing that, his input is just in radians, actually the input is the multiplier, like 1/4 * PI or 2 * PI, etc.

and yes, cos(2 * PI) = 1

using 6 taylor series terms will give you good values for all angles between -PI and +PI radians. however if you try to do 2*PI, you will need more series terms and the factorial will be too large for a 32-bit integer. any angles greater magnitude than +/- PI will have your series summation become increasingly inaccurate.

here's your corrected code. i think your problem was that you weren't alternately subtracting and adding the terms correctly. i cleaned up the indentation and teh debug output so it's more clear what's going on.

#include <stdio.h>

#define PI                3.14159   // make PI a #defined value, not a variable
#define TOTAL_SERIES      6         // the number of taylor series terms to use

#define DEBUG

float cosi(float x);
int factorial(int x);

main()
{
    /* Declarations */
    float input;
    double rad;
    double cosine;

    /* Ask user for an angle (in radians) */
    printf("\nEnter the angle multiplier in rads (don't include pi): ");

    /* While user has more input */
    while ( scanf("%f", &input) > 0)
    {
        /* Convert input into its decimal form */
        rad = input * PI;
        
        #ifdef DEBUG
        printf("\ndebug: radians = %06.4f\n", rad);
        #endif

        /* Call cosine function */
        cosine = cosi(rad);

        /* Print result */
        printf("\nThe cosine value is: %06.4f\n", cosine);

        /* Update loop */
        printf("\nEnter the angle multiplier in rads (don't include pi): ");
    }
    
    return 0;
}

float cosi(float x)
/* Given  : an angle in radians (float)
 * Returns: the cosine value of the angle */
{
    /*  Declare variables */
    int     series, fact=1;
    float   sum, term;
    
    sum = 1.0;
    term = 1.0;

    for(series = 1; series <= TOTAL_SERIES; series++)
    {
        #ifdef DEBUG
        printf("debug: sum = %05.3f\n", sum);
        #endif

        term *= x * x;    // multiply x^2 each iteration  (x^2, x^4, x^6, etc.)
        
        fact = factorial(2 * series);
        
        if (series % 2 == 1)         // odd sequences (first, third, fifth, etc...) 
            sum -= term / fact;      // subtract the quotient
 
        else if (series % 2 == 0)    // even sequences (second, fourth, sixth, etc...) 
            sum += term / fact;      // add the quotient
        
        #ifdef DEBUG
        printf("debug: %c %05.3f <- %05.3f / %d\n", 
                        (series % 2) ? '-' : '+', term/fact, term, fact);
        #endif
    }
    return(sum);
}


int factorial(int x)
/* Given   : a multiple of 2
 * Returns :the value of */
{
    int fact;
    
    if(x==1)
        fact = 1;

    else
        fact=x*factorial(x-1);
    
    return(fact);
}

[tquiva]

thank you both for the input. and thanks jephthah for helping me clean it up a bit. I have one more question. It seems to be working for values of 0.5, 0.24, 0.166667, and 1. But why am I getting a cosine value of 2.4653 for an input of 2?

[jonsca]

pi / 180 is what you'd use if you were making the conversion from degrees. i don't think he is doing that, his input is just in radians, actually the input is the multiplier, like 1/4 * PI or 2 * PI, etc.

I didn't read the OP's prompts, so that's my bad.

@OP: To illustrate the problem of approximations that jephthah was explaining look at this graph of sin x versus it's series approximation http://en.wikipedia.org/wiki/File:Taylorsine.svg. It only holds over a limited interval. I think you can in fact improve this a bit by specifying a point about which you are approximating the function (see http://en.wikipedia.org/wiki/Taylor_series#Definition where they have the (x-a) instead of x) but there's only so much you can do.

[jephthah]

thank you both for the input. and thanks jephthah for helping me clean it up a bit. I have one more question. It seems to be working for values of 0.5, 0.24, 0.166667, and 1. But why am I getting a cosine value of 2.4653 for an input of 2?

taylor series is an approximation. at n=0, the series is perfectly accurate. afterwards, everything else is an approximation.

the approximation will get better as you use more and more terms.

with 5 or 6 terms, the accuracy is "good enough" up to values of n=PI. the accuracy quickly degrades beyond 1*PI, and is totally unacceptable by the time you get to 2*PI.

if you increase the number of terms, the accuracy will hold out longer. however, as you increase the number of terms, the factorial in the denominator starts to get ridiculously huge. a 32-bit int will overflow after 6 terms. even a 64-bit int will overflow after 10 terms.

[0x69]

But why am I getting a cosine value of 2.4653 for an input of 2?

As jephthah explained - that is a property of Taylor series approximation. If you want some continuous approximation (which doesn't have increasing error) then you must use some other polynomial approximations - for example by using quadratic curve approximation or etc. =>
http://www.devmaster.net/forums/showthread.php?t=5784
http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/

Good luck!

[nezachem]

One thing to add. A straightforward implementation of Taylor series, as in the post #3, quickly leads to numerical problems. A 32-bit integer cannot hold a factorial of numbers beyond 12, which essentially limits the series to the first six members (hence the line 54 of the original post).

A proper way to calculate series is a Horner schema, which allows for literally thousands of member, and never forces a calculation of unreasonable large numbers. Here's how it works for a cosine:

cos x = 
1 - x^2/2! + x^4/4! - x^6/6! + ... =
1 - (x^2/(1*2)) * (1 - x^2/(3*4) + x^4/(3*4*5*6) - ...) =
1 - (x^2/(1*2)) * ( 1 - x^2/(3*4) * (1 - x^2/(5*6) + ...))

etc., which leads to the code

double cos(double x, int n)
{
    int i;
    double rc = 1.0;
    double xs = x*x;

    for(i = n*2; i > 0; i -= 2)
    {
        rc = 1.0 - xs / (i * (i - 1)) * rc;
    }

    return rc;
}

[tquiva]

I see it now. Thanks everyone. I have one more question. My assignment asks me to create another function called close_enough() in which it is given 2 double values, and returns TRUE if both values are within 0.00005 of each other. I really don't understand how I can implement this function into cosi(). How would it be related to finding the cos of an angle? My program seems alright now and I don't know what values the assignment is asking me to pass into this function. This is the function:

close_enough(double x, double y);

[0x69]

It seems that 0.00005 number should be used to accumulate required number of Taylor series terms. That is - you must continually add series terms, until you find that min(real_sin(),your_sin())/max(real_sin(),your_sin()) <= 0.00005 Or something like that.

[tquiva]

It seems that 0.00005 number should be used to accumulate required number of Taylor series terms. That is - you must continually add series terms, until you find that min(real_sin(),your_sin())/max(real_sin(),your_sin()) <= 0.00005 Or something like that.

I'm sorry, but I don't quite understand. So I would have to find the sin also? Could you please elaborate on this?

[0x69]

Oops, made mistake, I mean-

[TEX]\frac{1}{N} \sum_{x=0}^{2\pi} \frac{min(cos_{real}(x), cos_{approx}(x))}{max(cos_{real}(x), cos_{approx}(x))} <= 0.00005[/TEX]
So that averaged error should be not greater than 0.00005. So you keep adding Taylor terms, and after each additional series term re-compute average error over chosen set of angle values,- if it reaches 0.00005 limit - stop iteration through new Taylor series terms. This procedure helps to define How much terms you need to accumulate with series. It is useful when you don't have direct requirement how much terms you need to compute.

[tquiva]

Oops, made mistake, I mean-

[TEX]\frac{1}{N} \sum_{x=0}^{2\pi} \frac{min(cos_{real}(x), cos_{approx}(x))}{max(cos_{real}(x), cos_{approx}(x))} <= 0.00005[/TEX]
So that averaged error should be not greater than 0.00005. So you keep adding Taylor terms, and after each additional series term re-compute average error over chosen set of angle values,- if it reaches 0.00005 limit - stop iteration through new Taylor series terms. This procedure helps to define How much terms you need to accumulate with series. It is useful when you don't have direct requirement how much terms you need to compute.

So for the function:
close_enough(double x, double y);
The two values I would pass to it are term and..... what would the other value be?

[0x69]

Hm..., i guess second parameter should be 0.00005. And you sum series while NOT close_enough(new_series_term, 0.00005); But this is a bit different approach than I said before, because in this case we just verifying how much new series term adds precision. But maybe you can do that also.

[tquiva]

would it look something like this?

float total_term = 0;

total_term = total_term + term;

[0x69]

sort of. But i'm not sure how this integrates with Horner schema. Maybe in that case total_term = term * total_term , that is term = 1.0 - xs / (i * (i - 1)) . Try this out.

[tquiva]

I have honestly never heard of or worked with horner schema before, so could you please help me get the code for this function started?

So in my cosi() function, I would include this value:

new_series_term = 1.0 - xs / (i * (i - 1))

and I would then pass this value into the function

double NOT close_enough(double new_series_term, 0.00005)

Is this correct?

[jephthah]

if this is a school assignment, i think you should stick with the Taylor Series.

6 Taylor Series terms can still fit the factorial part in a 32-bit int, and this is a "close enough" approximation for all Theta values from -PI to +PI, therefore covering all possible output values of COS()

if someone sticks in a theta multiplier greater than +1 (or less than -1), just subtract (or add) enough multiples of 2 from that multiplier to get the total theta angle back between -PI and +PI radians.

for instance

cos( 5.0 * PI) = cos( 1.0 * PI)  -->   3.0 - ( 2 * 2) =  1.0
cos( 3.0 * PI) = cos( 1.0 * PI)  -->   3.0 - ( 1 * 2) =  1.0
cos( 2.4 * PI) = cos( 0.4 * PI)  -->   2.4 - ( 1 * 2) =  0.4
cos( 1.5 * PI) = cos(-0.5 * PI)  -->   1.5 - ( 1 * 2) = -0.5
cos(-1.5 * PI) = cos( 0.5 * PI)  -->  -1.5 - (-1 * 2) =  0.5
cos(-3.9 * PI) = cos( 0.1 * PI)  -->  -3.9 - (-2 * 2) =  0.1

cos(77.5 * PI) = cos( 1.5 * PI)  -->  77.5 - (38 * 2) =  1.5

of course there is benefit in learning new ways, so i'm not saying don't try it. but you it might be "easier" and you get the appropriate answer for your assignment to stick to the well-understood method.

.

[tquiva]

if this is a school assignment, i think you should stick with the Taylor Series.

6 Taylor Series terms can still fit the factorial part in a 32-bit int, and this is a "close enough" approximation for all Theta values from -PI to +PI, therefore covering all possible output values of COS()

if someone sticks in a theta multiplier greater than +1 (or less than -1), just subtract (or add) enough multiples of 2 from that multiplier to get the total theta angle back between -PI and +PI radians.

for instance

cos( 5.0 * PI) = cos( 1.0 * PI)  -->   3.0 - ( 2 * 2) =  1.0
cos( 3.0 * PI) = cos( 1.0 * PI)  -->   3.0 - ( 1 * 2) =  1.0
cos( 2.4 * PI) = cos( 0.4 * PI)  -->   2.4 - ( 1 * 2) =  0.4
cos( 1.5 * PI) = cos(-0.5 * PI)  -->   1.5 - ( 1 * 2) = -0.5
cos(-1.5 * PI) = cos( 0.5 * PI)  -->  -1.5 - (-1 * 2) =  0.5
cos(-3.9 * PI) = cos( 0.1 * PI)  -->  -3.9 - (-2 * 2) =  0.1

cos(77.5 * PI) = cos( 1.5 * PI)  -->  77.5 - (38 * 2) =  1.5

of course there is benefit in learning new ways, so i'm not saying don't try it. but you it might be "easier" and you get the appropriate answer for your assignment to stick to the well-understood method.

.

Yeah I think I will just stick to the Taylor series. But in one of my other posts above, I posted another question and was wondering if you knew the answer.

My assignment asks me to create another function called close_enough() in which it is given 2 double values, and returns TRUE if both values are within 0.00005 of each other. I really don't understand how I can implement this function into cosi(). How would it be related to finding the cos of an angle? My program seems alright now and I don't know what values the assignment is asking me to pass into this function. This is the function:

close_enough(double x, double y);

[nezachem]

I really don't understand how I can implement this function into cosi().

I don't think you should.
In my understanding (from the teacher's point of view), you have to prove somehow that your cosi indeed calculates cosine; in other words, that the values returned by cosi are close enough to the "golden standard" ones returned by cos from a math library.
That is, for each angle you have to calculate both cosi and cos , and feed the results to close_enough .
Yet again, it is my understanding; better check with your teacher.

[jephthah]

sorry i misunderstood the nature of the discussion, i thought you were trying to implement a more clever method to approximate COS.

...

as for "close_enough" don't make it harder than it needs to be.

say for instance, for 6 Taylor terms has your cosi() function returning the value of cosi( 1 * PI ) = -0.999971 (i don't know what it really returns)

and the *real* value of cos(), as implemented by the <math.h> library returns the value of cos ( 1 * PI ) = -1.000000 .

the close_enough() function just subtracts the two results and sees if the absolute value (magnitude) is less than 0.00005. if so, then it's "close enough". in the above made-up example, -0.999971 - (-1.000000) = 0.000029 so it is close enough.

this function is NOT part of cosi(). it's a separate function. the two arguments are the return values from cosi() and the math library's cos()

#include <math.h>

#define NUM_TAYLOR_TERMS                  6
#define PI                                3.14159

double approxVal, realVal, multiplier, theta;
int    isClose;

int main(void)
{
    //... get mulitplier input

    theta = PI * multiplier;

    approxVal = cosi(theta);
    realVal = cos(theta);

    //... print results

    isClose = close_enough(approxVal, realVal);

    printf("Using %d Taylor Series terms to compute cos(%06.4f) %s \"close enough\".\n",
                      NUM_TAYLOR_TERMS, theta, (isClose == 1) ? "is" : "is NOT");

    //... etc.