What should I do with my functions that return wrong answers?

Question

We could not use math.h

I get wrong answers in cases with inputs greater than 54 in my sine function and inputs greater than 37 in my exp functions; I guess it overflows, what should I do? I want to write sine function and exp function on my own and I use Taylor expansion. And I just need 6 digits in the fraction;

    //Start of exp function
float exp(float x)
{
    double expreturn=0;
    long long int fctrl=1;
    for(int n=0;n<=12;n++){

        fctrl=1;
        for(int i=2;i<=n;i++)
            fctrl *= i;
        expreturn += (pow(x,n)/fctrl);

    }

    return (float)expreturn;
}//End of exp function



//Start of sin function
float sin(float x)
{
    double sinus=0;
    long long int fctrl=1;
   while(!(0<=x&&x<6.3))
    {
        if(x<0)
            x += PI*2;
        else
            x -= PI*2;

    }
    for(int n=0;n<=8;n++)
    {
        fctrl=1;
        for(int i=2;i<=(2*n+1);i++)
            fctrl *= i;
       sinus += ((pow(-1,n)*pow(x,2*n+1))/fctrl);
    }
    return (float)sinus;
}//End of sin function


//Start of pow function
float pow(float x,int y)
{
    double pow = 1;
    for(int i=0;i<y;i++)
        pow *= x;
    return (float)pow;
}//End of pow Function

Here are some examples:

Input

sin(200)

Desired output

-0.873297

My function output

-0.872985

But it works properly with small values

I use this:

What Could I do now?

Answer 1

The way you calculate 'x modulus 2 PI' in sin() doesn't produce accurate results, e.g. with x = 200 you get 5.221206 whilst fmod() returns 5.221255.

Making x a double instead of float solved that inaccuracy for me.

That's the code I've tested with. It has produced the same results wihtout optimization or with -O3 on gcc (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI  3.14159265358979323846264338

int main( int argc, char *argv[] )
{
    double v = strtod( argv[1], 0 );
    float  f;
    double d;

    f = v;
    while( f > 6.3 ) {
        f -= 2*PI;
    }

    d = v;
    while( d > 6.3 ) {
        d -= 2*PI;
    }

    printf( "%f\n", f );
    printf( "%lf\n", d );
    printf( "%lf\n", fmod( v, 2*PI ) );
    return( 0 );
}

Edit

As @Eric Postpischil@ has pointed out, that still doesn't explain the false results, so 'Ive staretd playing around with sin() itself, and it seems increasing the number of iterations can increase accuracy. Making fctrl an unsigned long long you can iterate up to n=10 without overflowing and for x=200 or x=5.221255 you will now get the result -0.873260 which is at least closer.

That's the actual code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI  3.14159265358979323846264338

int iter = 8;
//Start of sin function
float mysin(double x)
{
    double sinus=0;
    unsigned long long int fctrl=1;
    while(!(0<=x&&x<6.3))
    {
        if(x<0)
            x += PI*2;
        else
            x -= PI*2;

    }
    for(int n=0;n<=iter;n++)
    {
        fctrl=1;
        for(int i=2;i<=(2*n+1);i++)
            fctrl *= i;

        double pr = (pow(-1,n)*pow(x,2*n+1));
        printf( "%lf / %llu\n", pr, fctrl );
        sinus += pr / fctrl;
    }
    return (float)sinus;
}//End of sin function

int main( int argc, char *argv[] )
{
    double v = strtod( argv[1], 0 );
    if( argc > 2 ) {
        iter = strtol( argv[2], 0, 10 );
    }
    float  f;
    double d;

    f = v;
    while( f > 6.3 ) {
        f -= 2*PI;
    }

    d = v;
    while( d > 6.3 ) {
        d -= 2*PI;
    }

/*
    printf( "%f\n", f );
    printf( "%lf\n", d );
    printf( "%lf\n", fmod( v, 2*PI ) );
*/

    printf( "%f\n", mysin( v ) );
    return( 0 );
}

Edit 2

Another improvement is to avoid calculating x ^ (2n+1) and n! first (which will produce huge values) and then divide but use a loop like that:

for(int n=0;n<=iter;n++)
{
    double v = n % 2 == 0 ? 1.0 : -1.0;
    for( int j=1; j <= 2*n+1; j++ ) {
            v *= x;
            v /= j;
    }
    sinus += v;
}

With n=10, mysin(200) returns -0.873296 now

Answer 2

There are several ways to increase the accuracy of the posted code.

• Using double as preferred type, instead on int and float, for variables storing intermediates, arguments and result.
• Take advantage of the symmetry of the function to limit the calculations in the range where the numerical approximation is more accurate.
  The sin function, for example, could be implemented using two function, an entry point which transforms the angle passed as argument and corrects the result of the other function implementing the series computation:

const double PI = 3.14159265358979323846264338; // It will be rounded, anyway.
const double TAU = PI * 2;
const double HALF_PI = PI / 2;
const double THREE_HALF_PI = PI * 1.5;

// This is called only for x in [-Pi/2, Pi/2]
double my_sin_impl(double x);

double my_sin(double x)
{
    // I want to transform all the angles in the range [-pi/2, 3pi/2]
    if ( x < -HALF_PI )
    {
        long long y = (THREE_HALF_PI - x) / TAU;
        x += y * TAU;
    }
    if ( x > THREE_HALF_PI )
    {
        long long y = (x + HALF_PI) / TAU;
        x -= y * TAU;
    }

    if ( x < HALF_PI )
        return my_sin_impl(x);
    else                              //    When x is in [pi/2, 3pi/2],
        return -my_sin_impl(x - PI);  // <- the result should be transformed
}

• Unless it's mandatory for the assignment, avoid using pow directly and calculating the factorial at every iteration. It's inefficient and could lead to overflow with integer types. See, e.g. the following implementation of a different approach.
• In the posted code, the number of iterations to calculate sin is fixed to 8, but we can stop when the numerical accuracy of the type is reached. In the following snippet I've also splitted the sum of terms with alternating signs into two partial sums, to limit (slightly) the loss of significance.

double my_sin_impl(double x)
{
    double old_positive_sum, positive_sum = x;
    double old_negative_sum, negative_sum = 0.0;
    double x2 = x * x;
    double term = x;
    int n = 2; 
    do
    {
        term *= x2 / (n * (n + 1));
        old_negative_sum = negative_sum;
        negative_sum += term;
        term *= x2 / ((n + 2) * (n + 3));
        old_positive_sum = positive_sum;
        positive_sum += term;
        n += 4;
    }
    while (positive_sum != old_positive_sum  &&  negative_sum != old_negative_sum);
    //                  ^^ The two values are still numerically different 
    return positive_sum - negative_sum;
}

Otherwise, noticing that the loop in the previous snippet is executed at most 5 times (give or take), we can decide to keep the number of titerations fixed and precalculate the coefficients. Something like this (testable here).

double my_sin_impl(double x)
{
    static const double k_neg[] = {1.0/ 6.0, 1.0/42.0, 1.0/110.0, 1.0/210.0, 1.0/342.0};
    static const double k_pos[] = {1.0/20.0, 1.0/72.0, 1.0/156.0, 1.0/272.0, 1.0/420.0};

    double positive_sum = 0.0;
    double negative_sum = 0.0;
    double x2 = x * x;
    double term = 1.0;
    for (int i = 0; i < 5; ++i) 
    {
        term *= x2 * k_neg[i];
        negative_sum += term;
        term *= x2 * k_pos[i];
        positive_sum += term;
    }
    return x * (1.0 + positive_sum - negative_sum);
}