#1877 Tiny double numbers arithmetics

[Anonymous]

I've tried to plot simple graph with the following command:

gnuplot -p -e 'plot [-2:2][-2:2] x**1000'

You can see the result in attachment (no solid line approximately between -0.5 and 0.5)

It turns out that gnuplot handles tiny numbers very specific:

$ gnuplot -e 'print 0.47**1000'
print 0.47**1000
      ^
line 0: undefined value

This behaviour differs from IEEE-754

System information:

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 16.04 LTS
Release:        16.04
Codename:       xenial
$ gnuplot --version
gnuplot 5.0 patchlevel 3

[Ethan Merritt]

"This behaviour differs from IEEE-754"
How so? This is a floating underflow error.

gnuplot uses the standard math library, in this case the pow() function.
man pow() says:

   Range error: the result underflows
          errno is set to  ERANGE.   An  underflow  floating-point  exception
          (FE_UNDERFLOW) is raised.

Gnuplot traps the exception and sets the returned variable to UNDEFINED.

[Sergei]

Yes, you are right. I just was confused with C++ standard (not C) which states for pow function

If a range error occurs due to underflow, the correct result (after rounding) is returned.

It would be nice to have some wrapper or flag to plot functions with underflow. Is there any way to check for underflow in gnuplot?

[Ethan Merritt]

I've looked harder at the underflow code path and will amend my earlier statement. In the case at hand the underflow occurs from gnuplot's expression evaluation interpreter, evaluate_at().
1) On entry the interpreter disables FPE trapping and clears the global errno.
2) It runs through the sequence of operations in the current expression
3) before returning a final value to the caller it checks to see if errno is set. If so then it returns UNDEFINED rather than a numeric result.

Unfortunately errno is too generic to distinguish underflow from overflow or other range errors and for that matter it might have been set multiple times during evaluation of a single expression. So we couldn't get what you want by modifying just this one place in the code.

C used to have a mechanism matherr(), or math_error(), but if I'm interpreting the man pages correctly these have been deprecated in favor of fpclassify().

One possibility is to add explicit tests at the individual math library call sites, in this case pow(a,b) using (fpclassify(result) == FP_ZERO). Something like

    temp = pow(a,b);
    if (errno == ERANGE  &&  fpclassify(temp) == FP_ZERO) {
        temp = 0.0;
        errno = 0;
    }
    Gcomplex(&result, temp, 0.0); /* push result onto stack */

The disadvantages I see are that equivalent code would have to be added in a lot of places, some other cases are not handled (e.g. FP_SUBNORMAL), and I do not know if the same code would work on all platforms.

Another possibility that I have not investigated is to put equivalent code in the routines that push an intermediate result onto the stack, Gcomplex() in this case. The problem I see there is that these routines are called from other places besides the expression parser and I don't know whether the errno/fpclassify tests would be valid in those other contexts.

[Karl Ratzsch]

Anybody know why I don't see this? The mingw64 windows build (from the sf download) returns zero for 0.47^1000. It's also compiled with gcc and the same standard libraries, isn't it?

[Ethan Merritt]

internal.c:50 bypasses the matherr code for MINGW64. I don't know why.

I'm less clear about the fpe_env FPE signal trap code. So far as I can see it is enabled on all platforms, but the details vary in syscfg.h so I suppose it is possible that the call to fpe() via fpe_env may not get set up on MINGW64. Those are the two places that the code sets the "undefined" flag if a floating point error occurs.

On the other hand, both of those paths may be relevant only to older platforms. So far as I can tell on my current linux systems the matherr and the fpe paths are set up but are never executed. The floating point error test that actually triggers is the test in evaluate_at() made after walking through the actions in the expression evaluation table.

    if (errno == EDOM || errno == ERANGE) {
        undefined = TRUE;

[Hans-Bernhard Broeker]

internal.c:50 bypasses the matherr code for MINGW64. I don't know why.

MinGW64 brings its own implementation of matherr, which collides with our own. It's unclear why this is being pulled in forcibly, instead of allowing itself to be overridden by user code, but that's how it is.

[Ethan Merritt]

f_power() modifed to test for underflow in 5.1 CVS.
I do not see many other places that reasonably require it, but the 3-4 that we do have could be consolidated at some point to use the same test.