ENH: Add smallest_normal and smallest_subnormal attributes to finfo

[steff456]

Given the discussions happened in data-apis/array-api#129 of the Data API consortium we decided to create a PR to add smallest_normal and smallest_subnormal numbers to finfo object. All the values added were checked against IEEE-754 and a test was included.

cc @rgommers @kgryte

[rgommers]

Thanks @steff456! Could you also make the corresponding change to MachAr in numpy/core/machar.py?

To add to the rationale: when looking into which attributes of finfo are used in the wild, we found that tiny was used fairly regularly, but in a good amount of cases not for its intended purpose but rather as "just give me a small number". Which is also what the name suggests. gh-16253 improved the docs, but it's still a poor name.

The smallest_subnormal addition was also made by @seberg in gh-16252 (he just called it tiny_subnormal). I think the tiny is a little odd and misleading, so would be great to leave that as an alias but have clearer names as in this PR.

[rgommers]

Looks good, just needs a few tweaks, in addition to updating machar.py

[rgommers]

There are failures for float128 on ARM64, on Shippable and TravisCI:

Two suggestions:

• don't use assert_equal; this is almost never the right thing to do for floating-point numbers. rather use assert_allclose with a relative tolerance
• in this case, the test is likely fragile and may fail on other platforms like AIX. I suggest to only run it for sys.platform == 'linux'.

[charris]

Two questions:

• float128 is quad precision on aarch64.
• behavior of subnormal numbers is determined by compiler flags.

Are these taken into account?

[seberg]

I just checked and nextafter is C99: https://en.cppreference.com/w/c/numeric/math/nextafter So... as annoying as it may be, it might be best to add the C99 loop for longdouble, and then you can use np.nextafter(0, 1) for the min, and np.nextafter(np.inf, 0) for "huge" and not worry about it.

As Chuck mentioned, it is plausible that this value might end up 0, the reference for nextafter notes:

in any case, the returned value is independent of the current rounding mode

So this seems unproblematic if we use np.nextafter, but problematic if we use log, etc. as done here. Maybe we should note that there is a small chance that subnormals are not available during runtime, i.e. 0. + smallest_subnormal == 0..

We should probably bikeshed the name on the mailing list, its annoying that python uses min for our tiny, while we use min for negative-max :). Apparently, so does C: https://en.cppreference.com/w/cpp/types/numeric_limits/denorm_min; "smallest" is probably the best choice to avoid running into similar names...

[pearu]

Here are other ways to compute the smallest subnormal numbers:

1. Using finfo:

>>> f = numpy.finfo('float32')
>>> numpy.float32(2 ** (f.minexp + f.negep + 1))
1e-45
>>> f = numpy.finfo('float64')
>>> numpy.float64(2 ** (f.minexp + f.negep + 1))
5e-324

2. Using bits representation:

>>> numpy.int32(1).view(numpy.float32)
1e-45
>>> numpy.int64(1).view(numpy.float64)
5e-324

I have not tested if the above is suitable for float128: there is no int128, etc.

[pearu]

Can this be the correct value for float128:

>>> numpy.complex128(0+1j*numpy.int64(1).view(numpy.float64)).view(numpy.float128)
3.3621031431120935063e-4932

?
EDIT: the answer is no as the approach gives incorrect value for float64:

>>> numpy.complex64(0+1j*numpy.int32(1).view(numpy.float32)).view(numpy.float64)
2.121995791e-314

[seberg]

Please, lets not cook our own soup here? Just implement the nextafter loop for long double, its probably as simple as adding an l to some loop defs/repeats. If that is doesn't work, we can still cook our own soup, but I would be very surprised.

[steff456]

Can this be the correct value for float128:

>>> numpy.complex128(0+1j*numpy.int64(1).view(numpy.float64)).view(numpy.float128)
3.3621031431120935063e-4932

?
EDIT: the answer is no as the approach gives incorrect value for float64:

>>> numpy.complex64(0+1j*numpy.int32(1).view(numpy.float32)).view(numpy.float64)
2.121995791e-314

I think in here you are calculating the smallest_normal value instead of the smallest_subnormal. But this are correct as the smallest normal for float 128

[pearu]

Please, lets not cook our own soup here?

I am not sure that this would be completely our own soup here. See for instance https://en.wikipedia.org/wiki/Single-precision_floating-point_format which shows the smallest positive number of float32:

0 00000000 00000000000000000000001_2 = 0000 0001_16

The same pattern applies to float64 and apparently to float128 as well:

>>> numpy.fromstring('\x01' + '\x00'*(4-1), dtype=numpy.float32)[0]
1e-45
>>> numpy.fromstring('\x01' + '\x00'*(8-1), dtype=numpy.float64)[0]
5e-324
>>> numpy.fromstring('\x01' + '\x00'*(16-1), dtype=numpy.float128)[0]
4e-4951

Just implement the nextafter loop for long double, its probably as simple as adding an l to some loop defs/repeats. If that is doesn't work, we can still cook our own soup, but I would be very surprised.

I think it would be interesting to see if the nextafter approach gives the same result as finding the subnormals from bits patterns.

[pearu]

>>> numpy.nextafter(0.0, 1.0, dtype=numpy.float32)
1e-45
>>> numpy.nextafter(0.0, 1.0, dtype=numpy.float64)
5e-324
>>> numpy.nextafter(0.0, 1.0, dtype=numpy.float128)
4e-4951

seems to work fine.

[rgommers]

So it sounds like this is converging. nextafter(0.0, 1.0, np.float128) is implemented already and will work (taking care of the aarch64 platform-dependent behavior). The next step here is:

1. Use nextafter(0.0, 1.0, np.float128) for the float128 smallest_subnormal value.
2. Add a doc disclaimer that this value may be zero on some hardware platforms.
3. Send an email to the mailing list that we plan to make this change.

Note that for (1) we should keep in mind that np.float128 does not exist on all platforms, so it should be implemented as something like:

if hasattr(numpy, `float128')
    tiny_f128 = np.nextafter(0.0, 1.0, np.float128)
else:
    # Likely on a 32-bit platform, where float96 exists instead of float128
    tiny_f128 = 4e-4951   # not used anyway, but set to the normal x86_64 number

We should probably bikeshed the name on the mailing list, its annoying that python uses min for our tiny, while we use min for negative-max :). Apparently, so does C: https://en.cppreference.com/w/cpp/types/numeric_limits/denorm_min; "smallest" is probably the best choice to avoid running into similar names...

I agree with the best choice. Let's just do the mailing list message but avoid the bikeshedding.

[pearu]

Note that for (1) we should keep in mind that np.float128 does not exist on all platforms, so it should be implemented as something like:

if hasattr(numpy, `float128')
    tiny_f128 = np.nextafter(0.0, 1.0, np.float128)
else:
    # Likely on a 32-bit platform, where float96 exists instead of float128
    tiny_f128 = 4e-4951   # not used anyway, but set to the normal x86_64 number

In the else block, 4e-4951 would evaluate as 0.0, and reporting tiny as 0.0 is misleading. I would suggest:

else:
   tiny_f128 = NotImplemented   # or None or anything else that would lead to an exception when some code is trying to use it

[seberg]

Before anyone gets more confused... Its not callednp.float128! Its called np.longdouble and always exists.

[pearu]

Before anyone gets more confused... Its not callednp.float128! Its called np.longdouble and always exists.

long double can be many things and is system and compiler-dependent. See, for instance, https://en.wikipedia.org/wiki/Long_double .

I assumed that we are talking about float128 (or binary128) as described in https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format . It would be strange to call float128 as long double when long double happens to be float64 (as in x86/MVC) or float80 (as in x86) or any other variant.

[seberg]

No @pearu we are most definitely not, there is no such thing in NumPy, float128 is just a badly named alias for longdouble when it happens to be stored in 128 bits.

[rgommers]

@seberg good point, but I'm not sure it applies to finfo. Note that the implementation has explicit f80 and f128 attributes, and a comment for f128 saying # IEEE 754 128-bit binary float. And for f80 it says float80 (Intel 80-bit extended precision).

I don't have a machine at hand where np.longdouble is np.float128 is not true, but if you look at the hardcoded numbers for f128 I believe that comment is correct:

    float128_ma = MachArLike(ld,
                             machep=-112,
                             negep=-113,
                             minexp=-16382,
                             maxexp=16384,
                             it=112,
                             iexp=15,
                             ibeta=2,
                             irnd=5,
                             ngrd=0,
                             eps=exp2(ld(-112)),
                             epsneg=epsneg_f128,
                             huge=huge_f128,

Whether finfo is behaving correctly for all >64-bit dtypes on all platforms is another question, I'm not sure (but it's doubtful, the aarch64 failure suggested it doesn't).

[rgommers]

No @pearu we are most definitely not, there is no such thing in NumPy, float128 is just a badly named alias for longdouble when it happens to be stored in 128 bits.

Also note that this is a somewhat implementation-centric view, that's probably good for downstream library devs who simply want to support all dtypes - but not very useful for end users who need extended precision. The only good reason to be using any of the extended precision dtypes in application code is if 64 bits of precision is not enough. The last thing you'd then want is to use longdouble - your code will run everywhere, but it will also silently give imprecise/incorrect answers - an error is likely better.

[seberg]

Of course it also applies to finfo...

That the code here lists all types of possible longdouble implementations is an implementation detail that is not exposed to the user. The code doesn't even attempt to make sure the the right values are stored (how can it possibly do so if it does calculations in longdouble)...
Hardcoding the values (even if they get calculated and are effectively filled with nonsense) seems very fair to me. At least its easy to understand compared to the np.MaArach actually discovering them.

I tried to nudge this thing in direction of using nextafter since Chuck rightly pointed out that double-double, and flush-to-zero can be tricky and seems to be ignored. I thought I can just mention it and not think about this much more, but more clearly:

• The current PR has line smallest_subnormal=exp2(ld(-16445)) three times in it. It is impossible that this is line is correct more than once (even once is not certain!)
• nextafter should be a solution that just works always. I hoped it would help you to either just use it and be done, or at least use it to check all the arches (or even write a test).
• The other question was what Chuck asked, about flush-to-zero effectively disabling subnormals. My guess is, that flush-to-zero is the "worst" that can happen, which means we still have "smallest_subnormal" defined, but can't really use it. But that also means that we can't actually calculate it, since that would give 0! That is even be problematic for double or float if we care about it, so the current version may be strictly speaking incorrect even for double (nextafter seems to avoid this issue!). Although I am not sure flush-to-zero really matters in NumPy, since it should not be enabled.

The whole discussion mentioning float128 seemed to steer away from those question rather than towards them, so I grumply pointed out to ignore it. Probably that was my bad for not being clear enough.

[leofang]

Although I am not sure flush-to-zero really matters in NumPy, since it should not be enabled.

Just FYI. In CuPy we turn flush-to-zero on by default (via setting a compiler flag -ftz=true), see https://github.com/cupy/cupy/blob/eafd82bfb80b5e7061cfd2eebe56bb2e4de658c5/cupy/_math/floating.py#L56-L58. The reason is there could be significant computational cost on accelerators like GPU. So, I suspect that CuPy is not the only array library that does this. I think I've seen some workarounds done in our test suite to avoid direct comparison with NumPy's outcome that is sensitive to denormals, but I can't find it atm.

[charris]

Needs release note.

[charris]

Needs a release note.

[steff456]

The latest commit includes the change to use nextafter on the finfo class. Before using that value, I'm checking that the data type is not missing on the platform, in the case it is not present I'm returning the value for float64 as the case of huge in the double double type.

As well, I added comments to add a disclaimer that the value of the smallest_subnormal can be zero and all the tests are passing now 🎉.

[seberg]

xref gh-8504 which is the PR adding this originally.

I honestly don't know if this is right or not. @charris or even @matthew-brett do you have any insight for the smallest subnormal? In general I would prefer using nextafter everywhere and unconditionally.
(Mainly, because I don't want to try and understand if flush-to-zero mode might be enabled here! I admit, for all I know, some platforms may not support have subnormals at all and nextafter returns the smallest normal, but that seems less likely/bad then incorrectly returning 0?!).

Note that the hasattr(umath, "nextafter") check is something that snuck in incorrectly in the original PR (gh-8504), it is always True and should be removed! (nextafter may be wrong on some platforms, but it is always defined and there is no way to know.)

There seems to be the problem that nextafterl is not defined on AIX (as Chuck posted on the old PR): https://www.ibm.com/docs/en/aix/7.2?topic=sepl-128-bit-long-double-floating-point-data-type which is not C99 conform :(. That means that we may well be getting bogus values there (we do have a custom nextafter which might patch that, but I don't know and its impossible to know from Python which version we use).

[seberg]

I made some last comments, the rest looks fine. I still expect the NaN return won't be great, but since the old value was probably wrong, might as well.

Maybe one of the other early reviewers wants to finalize?

[mattip]

@steff456 there are merge conflicts.

[seberg]

Ok, give this a shot. Thanks everyone!