Machine limit functions

[ben-albrecht]

Feature request for something like numpy.finfo and numpy.iinfo, which return the machine epsilon value for float and int types, respectively.

I believe one can work around this by doing something like:

/* Machine epsilon for real(64) */
private proc epsilon(type t: real(64)) : real {
  extern const DBL_EPSILON: real;
  return DBL_EPSILON;
}

[bradcray]

Random question: Is there any chance of these values varying from one compute node to the next? (on a user cluster, for example). If so, should these be made methods on the locale type?

[ben-albrecht]

I am not sure if there are cases where machine epsilon can depend on factors beyond the type itself. I am also not sure when it is required to computing machine epsilon vs. looking it up a pre-computed value, like my example that uses the float.h header.

Looking into the numpy implementations could be a good first step to answering these questions.

[damianmoz]

I shall endeavour to put a draft of paper on my implementation of such limits before Christmas. Sorry, been distracted for a while from my R+D doing billable work and generally keeping customers happy.

[RedHatTurtle]

I think a more interesting implementation of this might be the Fortran EPSILON(X) procedure defined as

EPSILON(X) returns the smallest number E of the same kind as X such that 1 + E > 1.

I'm looking if I can find the GCC/gfortran implementation of this procedure but I've never messed with the gfortran source code before.

Edit: Found some stuff but it's way over my head to understand it

[damianmoz]

I respond very belatedly to the question by @bradcray which slipped my attention at the time. As long as the hardware complies with IEEE 754, and computations are being done at the same precision on both nodes, say real(w), then the value of epsilon is the same.

If the supposition that those computations are being done at the same precision is violated, then the effective value of epsilon is different. But that should be obvious from your Chapel code, if Chapel is running on both nodes. But, if you are downloading something to a GPU and asking it to do things in real(32) while you CPU is working at real(64)(, then off course that epsilon will be different of those two locales.

The above is different to rounding error which is a dynamic quantity, and for a simple operation such as the addition of two floating point quantities, may have a absolute value that is as large as epsilon Also, if for some strange reason, there are different rounding modes in effect on those two locales the same point in time, the rounding error inherent in the same expression on those two machine could be different. But that has nothing to do with machine limits. The user should ensure that rounding modes are the same for different locales, although some exotic CPUs such as the Hitachi SH-4s may not allow such a change .

I hope I have answered that question.

[damianmoz]

In response to the question about the Fortran Intrinsic, this is never expand as a function, more as a constant.

I will look to see what is required in terms of test case requirements to get our ieee754 module in the normal release. It contains these Machine Limits, as well as some routines to manipulate the same.

This module provides constants such as epsilon, its reciprocal, the radix, precision, the largest and small floating point numbers, infinity, NaN and several other useful numbers which characterise the floating point model under the hood of Chapel's real(w) type. They are param methods so there is no run-time overhead in obtaining them within a program. They are available as methods against an identifier, say x of some floating point type or the type name itself, say T, as in

x.epsilon
T.epsilon

This provides the constant in question. Where a return value is non-integral or non-boolean, it has the same type as x or T. The use of a method avoids namespace issues. Constants which are integral have type int(16), which may change to int(32) on production release.

The module has guaranteed inline versions of some of the standard library functions like

isinf(x)
isnan(x)

It includes some use-with-care routines to provide introspection of floating point values to yield things like the negative bit, the raw biased or unbiased exponent, facilitate manipulating the sign-magnitude representation of those same floating point values, and convert that representation back into a floating point value.

It is not-portable to a Chapel implementation which does not support a uint(w) type for every real(w) type.

[RedHatTurtle]

I took a refresher on IEEE754 and I guess the simple solution is just to take the exponent from a floating point and set the fraction to 0...01 to get epsilon, or something like that?

Anyway it would be great to have access to the ieee754 module in some circumstances.

[damianmoz]

What you are explaining is not what Fortran calls EPSILON(X). Maybe you are asking for a way to compute the rounding error in X or what is called ULP(x) or Unit in Last Place of x. That is not a machine limit function but is something totally different to the machine epsilon. In that case, if this was C, there is routine nextafter(x, y)that is available in the Math library for this task. Assuming that say x is a double, the computation is just

ulpx = nextafter(x, (double) INFINITY) - x;

Unless I have totally misunderstood your explanation, your issue has nothing to do with Chapel nor even machine limit functions.

That nextafter(x, t) routine is not exposed to Chapel directly with a wrapper. You can access it with a simple Chapel wrapper interface to C's Math library, Others are better qualified to give you advice on how to do this optimally.