Explain this code in K&R 2-1

I'm trying to determine range of the various floating-point types. When I read this code:

#include <stdio.h>

main()
{
    float fl, fltest, last;
    double dbl, dbltest, dblast;

    fl = 0.0;
    fltest = 0.0;
    while (fl == 0.0) {
        last = fltest;
        fltest = fltest + 1111e28;
        fl = (fl + fltest) - fltest;
    }
    printf("Maximum range of float variable: %e\n", last);

    dbl = 0.0;
    dbltest = 0.0;
    while (dbl == 0.0) {
        dblast = dbltest;
        dbltest = dbltest + 1111e297;
        dbl = (dbl + dbltest) - dbltest;
    }
    printf("Maximum range of double variable: %e\n", dblast);
    return 0;
}

I don't understand why author added 1111e28 at fltest variable ?

回答1:

The loop terminates when fltest reaches +Inf, as at that point fl = (fl + fltest) - fltest becomes NaN, which is unequal to 0.0. last contains a value which when added to 1111e28 produces +Inf and so is close to the upper limit of float.

1111e28 is chosen to reach +Inf reasonably quickly; it also needs to be large enough that when added to large values the loop continues to progress i.e. it is at least as large as the gap between the largest and second-largest non-infinite float values.

回答2:

OP: ... why author added 1111e28 at fltest variable ?
A: [Edit] For code to work using float, 1111e28, or 1.111e31 this delta value needs careful selection. It should be big enough such that if fltest was FLT_MAX, the sum of fltest + delta would overflow and become float.infinity. With round to nearest mode, this is FLT_MAX*FLT_EPSILON/4. On my machine:

min_delta           1.014120601e+31 1/2 step between 2nd largest and FLT_MAX
FLT_MAX             3.402823466e+38
        FLT_EPSILON 8.388608000e+06
FLT_MAX*FLT_EPSILON 4.056481679e+31

delta needs to be small enough so if f1test is the 2nd largest number, adding delta, would not sum right up to float.infinity and skip FLT_MAX. This is 3x min_delta

max_delta           3.042361441e+31

So 1.014120601e+31 <= 1111e28 < 3.042361441e+31.

@david.pfx Yes. 1111e28 is a cute number and it is in range.

Note: Complications occur when the math and its intermediate values, even though the variables are float may calcuate at higher precsison like double. This is allowed in C and control by FLT_EVAL_METHOD or very careful coding.

1111e28 is a curious value that makes sense if the author all ready knew the general range ofFLT_MAX.

The below code is expected to loop many times (24946069 on one test platform). Hopefully, the value fltest eventually becomes "infinite". Then f1 will becomes NaN as the difference of Infinity - Infinity. The the while loop ends as Nan != 0.0. @ecatmur

while (fl == 0.0) {
    last = fltest;
    fltest = fltest + 1111e28;
    fl = (fl + fltest) - fltest;
}

The looping, if done in small enough increments, will arrive at a precise answer. Prior knowledge of FLT_MAX and FLT_EPSILON are needed to insure this.

The problem with this is that C does not define the range FLT_MAX and DBL_MAX other than they must be at least 1E+37. So if the maximum value was quite large, the increment value of 1111e28 or 1111e297 would have no effect. Example: dbltest = dbltest + 1111e297;, for dbltest = 1e400 would certainly not increase 1e400 unless dbltest a hundred decimal digits of precision.

If DBL_MAX was smaller than 1111e297, the method fails too. Note: On simple platforms in 2014, it is not surprising to find double and float to be the same 4-byte IEEE binary32 ) The first time though the loop, dbltest becomes infinity and the loop stops, reporting "Maximum range of double variable: 0.000000e+00".

There are many ways to efficiently derive the maximum float point value. A sample follows that uses a random initial value to help show its resilience to potential variant FLT_MAX.

float float_max(void) {
  float nextx = 1.0 + rand()/RAND_MAX;
  float x;
  do {
    x = nextx;
    nextx *= 2;
  } while (!isinf(nextx));
  float delta = x;
  do {
    nextx = x + delta/2;
    if (!isinf(nextx)) {
      x = nextx;
    }
    delta /= 2;
  } while (delta >= 1.0);
  return x;
}

isinf() is a new-ish C function. Simple enough to roll your own if needed.

In re: @didierc comment

[Edit]
The precision of a float and double is implied with "epsilon": "the difference between 1 and the least value greater than 1 that is representable in the given floating point type ...". The maximum values follow

FLT_EPSILON 1E-5
DBL_EPSILON 1E-9

Per @Pascal Cuoq comment. "... 1111e28 being chosen larger than FLT_MAX*FLT_EPSILON.", 1111e28 needs to be at least FLT_MAX*FLT_EPSILON to impact the loop's addition, yet small enough to precisely reach the number before infinity. Again, prior knowledge of FLT_MAX and FLT_EPSILON are needed to make this determination. If these values are known ahead of time, then the code simple could have been:

printf("Maximum range of float variable: %e\n", FLT_MAX);

回答3:

The largest value representable in a float is 3.40282e+38. The constant 1111e28 is chosen such that adding that constant to a number in the range of 10^38 still produces a different floating point value, so that the value of fltest will continue to increase as the function runs. It needs to be large enough that it will still be significant at the 10^38 range, and small enough that the result will be accurate.