Why this sin(x) function in C# return NaN instead of a number

Question

I have this function wrote in C# to calc the sin(x). But when I try with x = 3.14, the printed result of sin X is NaN (not a number), but when debugging, its is very near to 0.001592653 The value is not too big, neither too small. So how could the NaN appear here?

static double pow(double x, int mu)
        {
            if (mu == 0)
                return 1;
            if (mu == 1)
                return x;
            return x * pow(x, mu - 1);
        }

        static double fact(int n)
        {
            if (n == 1 || n == 0)
                return 1;
            return n * fact(n - 1);
        }

        static double sin(double x)
        {
            var s = x;

            for (int i = 1; i < 1000; i++)
            {
                s += pow(-1, i) * pow(x, 2 * i + 1) / fact(2 * i + 1);
            }
            return s;
        }

        public static void Main(String[] param)
        {
            try
            {
                while (true)
                {
                    Console.WriteLine("Enter x value: ");
                    double x = double.Parse(Console.ReadLine());
                    var sinX = sin(x);
                    Console.WriteLine("Sin of {0} is {1}: " , x , sinX);

                    Console.ReadLine();
                }
            }
            catch (Exception ex)
            {
                Console.WriteLine(ex.Message);
            }
        }

What value is X? I used the value of 1 and it comes out as 0.8414709848078965 - explain the steps to reproduce the problem.. — Jeremy Thompson, Jan 15 '16 at 01:35

Rob · Accepted Answer · 2016-01-15T02:22:17.953

3

It fails because both pow(x, 2 * i + 1) and fact(2 * i + 1) eventually return Infinity.

In my case, it's when x = 4, i = 256.

Note that pow(x, 2 * i + 1) = 4 ^ (2 * 257) = 2.8763090157797054523668883052624395737887631663 × 10^309 - a stupidly large number which is just over the max value of a double, which is approximately 1.79769313486232 x 10 ^ 308.

You might be interested in just using Math.Sin(x)

Also note that fact(2 * i + 1) = 513! =an even more ridiculously large number which is more than 10^1000 times larger than the estimated number of atoms in the observable universe.

edited Jan 15 '16 at 02:22

answered Jan 15 '16 at 01:41

Rob

26,989
16
82
98

OK, May be I have to change the algorithm :) – Andiana Jan 15 '16 at 01:45
1

I want to write this function for a benchmark, not for the calculation purpose :) – Andiana Jan 15 '16 at 01:47
Use ILSpy and have a look at the Math.Sin internal implementation. – Jeremy Thompson Jan 15 '16 at 01:48
I think verdana is not using `Math.Sin` because the purpose of the exercise it to calculate the sin function as a series. – Jens Meinecke Jan 15 '16 at 02:52
1

@JeremyThompson All of the trig functions are implemented in the CLR in native code that leverages the platform capabilities, so ILSpy isn't going to tell you anything useful. – Corey Jan 15 '16 at 04:01
Hey @verdana see this thread and all the excellent resources that are referenced: http://stackoverflow.com/questions/25327821/how-can-i-see-the-source-code-of-system-math-sin – Jeremy Thompson Jan 15 '16 at 04:17

score 2 · Answer 2 · answered Jan 15 '16 at 01:43

2

When x == 3.14 and i == 314 then you get Infinity:

?pow(-1, 314)
1.0
?pow(x, 2 * 314 + 1)
Infinity
? fact(2 * 314 + 1)
Infinity

answered Jan 15 '16 at 01:43

Jeremy Thompson

61,933
36
195
321

score 0 · Answer 3 · answered Jan 18 '16 at 02:10

The problem here is an understanding of floating point representation of 'real' numbers.

Double numbers while allowing a large range of values only has a precision of 15 to 17 decimal digits.

In this example we are calculating a value between -1 and 1.

We calculate the value of the sin function by using the series expansion of it which is basically a the sum of terms. In that expansion the terms become smaller and smaller as we go along.

When the terms have reached a value less than 1e-17 adding them to what is already there will not make any difference. This is so because we only have 52 bit of precision which are used up by the time we get to a term of less than 1e-17.

So instead of doing a constant 1000 loops you should do something like this:

 static double sin(double x)
    {
        var s = x;

        for (int i = 1; i < 1000; i++)
        {
            var term = pow(x, 2 * i + 1) / fact(2 * i + 1);

            if (term < 1e-17)
               break;

            s += pow(-1, i) * term;
        }
        return s;
    }

Why this sin(x) function in C# return NaN instead of a number

3 Answers3