-nan return value / e (euler) raised to a power calculation loop

Question

I'm learning C programming and made the algorithm below to solve this problem:

Link to a pic of the problem. The code actually works, but initially the loop was with only 10 repetitions (rep <= 10), and the anwer for p = 3 was almost correct, so I changed rep <= 20. And It gave me just the exact answer from my calculator. And then I tried with a higher number, 12, and the output again was inaccurate. So I ended raising rep <= 35. If I get the loop for higher repetitions I get "-nan", and if the input for p is too high it will be the same. So just have to see the pattern to know that the problem of inaccuracy will get back as I input higher numbers which is not the case because the output will be NaN if I input a high value.

Is it possible to solve it without higher level functions? just want to know if my program is ok for the level in which I am now...

#include <stdio.h>

int main()
{
    float p; //the power for e
    float power; //the copy of p for the loop
    float e = 1; //the e number I wanna raise to the power of p
    int x = 1; //the starting number for each factorial generation
    float factorial = 1;
    int rep = 1; //the repeater for the loop

    printf( "Enter the power you want to raise: " );
    scanf( "%f", &p );

    power = p;

    while ( rep <= 35) {
        while ( x > 1) {
            factorial *= x;
            x--;
        }
        e += p / factorial;

        //printf("\nthe value of p: %f", p); (TESTER)
        //printf("\nthe value of factorial: %f", factorial); (TESTER)

        p *= power; //the new value for p
        rep++;
        factorial = 1;
        x = rep; //the new value for the next factorial to be generated

        //printf("\n%f", e); (TESTER)
    }
    printf("%.3f", e);

    return 0;
}

Sorry if I had syntax/orthography errors, I'm still learning the language.

Hint: What's the range (minimum/maximum value) of a `float`? When will it overflow? — user202729, Mar 10 '18 at 16:15
Please do not take the `!` from the picture if you edit. Having to click the link to get the picture is a hassle. With the picture in the question it's much easier to read. — Pablo, Mar 10 '18 at 16:21
@Pablo I doubt the OP, with 1 rep, can add image to the post. Also, we don't have Latex enabled, so... — user202729, Mar 10 '18 at 16:21
@user202729 the OP uploaded the picture when writing the question, Stack Overflow provides `[enter picture description][1]` but most people don't write any description nor add the `!` before the link to embed the picture. I edited to show the picture and a few minutes later the OP edited the question and removed the `!`. That was what I was pointing out. — Pablo, Mar 10 '18 at 16:23
A common error when working with floating points is that floating point errors accumulate over time. See https://stackoverflow.com/questions/588004/is-floating-point-math-broken — Pablo, Mar 10 '18 at 16:28
1) you can compute each term iteratively *together* instead of computing the power and factorial separately. 2) can use [Kahan/Neumaier summation](https://en.wikipedia.org/wiki/Kahan_summation_algorithm) to reduce error somewhat. — meowgoesthedog, Mar 10 '18 at 16:32
The factorial calculation can be done by doing `factorial *= rep; ` in the loop, then you can remove `x` and the whole `while(x > 1)` loop. You get the NaN because `float` cannot hold the value if `e` and the factorial. In the 30th iteration, `e` is already `inf` and in the 34th iteration `factorial` is already NaN. Use `double` instead of `float` for more precision. — Pablo, Mar 10 '18 at 16:38

score 1 · Accepted Answer · answered Mar 10 '18 at 17:11

Before we begin, let's write your original code as a function, with some clean-ups:

float exp_original(float x, int rep = 35)
{
   float sum = 1.0f;
   float power = 1.0f;
   for (int i = 1; i <= rep; i++)
   {
      float factorial = 1.0f;
      for (int j = 2; j <= i; j++)
         factorial *= j;
      power *= x;
      sum += power / factorial;
   }
   return sum;
}

There were some unnecessary variables you used which were removed, but otherwise the procedure is the same: compute the factorial from scratch.

Let's look at the ratio between successive terms in the series:

We can thus simply multiply the current term by this expression to get the next term:

float exp_iterative(float x, int rep = 35)
{
   float sum = 1.0f;
   float term = 1.0f;
   for (int i = 1; i <= rep; i++)
   {
      term *= x / i;
      sum += term;
   }
   return sum;
}

Seems much simpler, but is it better? Comparison against the C-library exp function (which we assume to be maximally precise):

x   exp (C)       exp_orig    exp_iter
-------------------------------------------
1   2.7182817     2.718282    2.718282
2   7.3890562     7.3890567   7.3890567
3   20.085537     20.085539   20.085539
4   54.598148     54.598152   54.598152
5   148.41316     148.41318   148.41316
6   403.4288      403.42871   403.42877
7   1096.6332     1096.6334   1096.6334
8   2980.958      2980.9583   2980.9587
9   8103.084      8103.083    8103.083
10  22026.465     22026.467   22026.465
11  59874.141     59874.148   59874.152
12  162754.8      162754.77   162754.78
13  442413.41     -nan(ind)   442413.38
14  1202604.3     -nan(ind)   1202603.5
15  3269017.3     -nan(ind)   3269007.3
16  8886111       -nan(ind)   8886009
17  24154952      -nan(ind)   24153986
18  65659968      -nan(ind)   65652048
19  1.784823e+08   -nan(ind)  1.7842389e+08
20  4.8516518e+08  -nan(ind)  4.8477536e+08

The two custom implementations are neck-and-neck in-terms of precision, until x = 13 where the original gives NaN. This is because the highest power term 13^35 = 9.7278604e+38 exceeds the maximum value FLT_MAX = 3.40282e+38. The accumulated term in the iterative version never reaches anywhere near the limit.

-nan return value / e (euler) raised to a power calculation loop

1 Answers1