C++ Endless loop with function sin x problem

Question

I am practicing for my final exam and professor sent us a last year's copy and in it, there is this assignment:

I made a code for it, but it seems it doesn't work properly, it ends up in an endless loop most of the times. This is my code:

int main(){
    double x;
    double sinx;
    int n=1;
    int sign=-1;
    cin>>x;
    sinx=0;
    for(;;) {
        int fact=1;
        for(int i=1;i<=2*n+1;i++){
            fact*=i;
        }
        double term=pow(x,2*n+1)/fact;
        if(term==0){
            break;
        }
        sinx += term * sign;
        sign *= -1;
        n++;

    }
    cout<<sinx;
}

I really tried everything.

Answer 1

EDIT:

The bug is in integer overflow of fact. There should be a check that integer overflow hasn't happened.

#include<cassert>
#include<iostream>
#include<cmath>

int main(){
    double x;
    double sinx;
    int n=1;
    int sign=-1;
    std::cin>>x;
    sinx=0;
    bool overflow = false;

    for(;;) {
        int64_t fact=1;
        for(int i=1;i<=2*n+1;i++){
            old_fact = fact;
            fact = fact * i;
            // check for overflow
            if (fact / i != old_fact) {
               overflow = true;
               break;
            }
        }
        if (overflow)
           break;


        double term=std::pow(x,2*n+1)/fact;
        if(term==0){
            break;
        }
        std::cout<<"n=" << n << ", fact = " << fact << ", term: " << term <<"\n";
        sinx += term * sign;
        sign *= -1;
        n++;
        std::cout<<"sinx=" << sinx<<"\n\n";
    }
    std::cout<<sinx;
}

godblot link: https://godbolt.org/z/zP8sT6ovq

Old answer:

The problem is with checking if a floating point number is exactly equal to zero.

        if(term==0){

This may not ever evaluate to true. It may keep oscillating around a value close to zero. What you could do is to compare if term is arbitrarily close to zero. You can hard code a value like 0.000001 (which I recommend for a class project) or you can use epsilon if you want to get more technical.

So the code will be like:

if (std::abs(term) < 0.000001) {
  break;

Answer 2

Your basic problem is using an int to calculate fact (factorial), which will overflow pretty quickly, leading to undefined behavior.

Change that to double and things should work, as double has much larger range (won't overflow so soon), and even when it does, the overflow is well-defined, giving ∞

---

You can make it much more efficient by not recomputing n! and xⁿ from scratch each iteration of the loop -- just multiply the term by x*x and divide it by (2*n)*(2*n+1.0)

Another problem you may see is loss of precision when x is large. This is a general problem with using series expansions like this -- for a periodic function like sin for real use, you generally want to do range reduction first.

Answer 3

A few comments about your code:

• I very strongly recommend you put all the sin calculations in a function called sin, and leave only the input and output logic in the main function.
• The formula you have been given is the Taylor expansion of sin(x) at x = 0. It gives a good approximation of sin(x) when x is close to 0, but quickly becomes completely wrong when x is large. In fact I recommend to only use it for x in range [-pi/3, +pi/3]. Fortunately, sinus is periodic, so you can remove 2 pi from x until x is in range [-pi, +pi]; and then you can use the trigonometric formula sin(x) = 2 * sin(x/2) * cos(x/2) and cos(x) = cos(x/2) * cos(x/2) - sin(x/2) * sin(x/2) to divide x by 2 before computing its sine.
• The text you have been given says "The algorithm should run until one of the factors reaches 0 value." This sentence is meaningless. The factors converge towards 0, which means they become closer and closer to 0, but they never reach 0. Fortunately, the Taylor expansion is good enough if you have taken care to reduce x to range [-pi/4, +pi/4]: the error will be below 0.00001 even for small values of n.
• Since x^(2n + 3) == x^(2n+1) * x^2 and fact(2n + 3) == fact(2n + 1) * (2n + 2) * (2n + 3), you don't need to recompute the powers of x and the factorials from scratch at every iteration.
• If you can use a proper for-loop with an explicit counter, upper-bound, and increment, then use it! It makes the code a lot easier to read and debug than using for(;;) and then breaking the loop unannounced with break.
• Since we want to reduce x to range [-pi/4, +pi/4] and then use the Taylor expansion only in that range, I recommend putting the range reduction logic in one function, and the Taylor expansion logic in another function.

With all this in mind:

#include <stdio>

double mysin(double x)
{
    double twopi = 6.283185307179586;
    double    pi = 3.141592653589793;
    while (x > pi)
    {
        x -= twopi;
    }
    while (x < -pi)
    {
        x += twopi;
    }
    double sinx4 = mysin_taylor(x / 4);
    double cosx4 = mycos_taylor(x / 4);
    return (4 * sinx4 * cosx4 * (cosx4 * cosx4 - sinx4 * sinx4));
}

/* Only use mysin_taylor(x) for -pi/3 < x < pi/3 */
double mysin_taylor(double x)
{
    unsigned int n = 11;
    double x2 = x * x;
    double xk = x;
    double ck = 1.0;
    double sinx = 0.0;
    for (unsigned int k = 1; k <= n; k += 2)
    {
        sinx += ck * xk;
        xk *= x2;
        ck *= (-1) / ((k+1)*(k+2));
    }
    return (sinx);

/* Only use mycos_taylor(x) for -pi/3 < x < pi/3 */
double mycos_taylor(double x)
{
    unsigned int n = 10;
    double x2 = x * x;
    double xk = 1.0;
    double ck = 1.0;
    double cosx = 0.0;
    for (unsigned int k = 0; k <= n; k += 2)
    {
        cosx += ck * xk;
        xk *= x2;
        ck *= (-1) / ((k+1)*(k+2));
    }
    return (cosx);

int main(){
    double x;
    double sinx;
    std::cin >> x;
    sinx = mysin(x);
    std::cout << sinx;
}