Inconsistent underflow warnings

[WarrenWeckesser]

This inconsistent warning of underflow was reported on stackoverflow: https://stackoverflow.com/questions/45703529/why-does-complex-floating-point-division-underflow-weirdly-with-numpy

I think the following are cases of the same inconsistency (using numpy 1.13.1):

In [302]: np.seterr(under='warn')
Out[302]: {'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'warn'}

In [303]: r = 1/(1<<533)

In [304]: r
Out[304]: 3.556413999176124e-161

The following multiplications do not trigger an underflow warning:

In [305]: np.multiply(r, r)
Out[305]: 1.2648080533535912e-321

In [306]: np.multiply(r, 1.5*r)
Out[306]: 1.8972120800303867e-321

In [307]: np.multiply(r, 2.0*r)
Out[307]: 2.5296161067071823e-321

But this does:

In [308]: np.multiply(r, 1.6*r)
/Users/warren/miniconda3scipy/bin/ipython:1: RuntimeWarning: underflow encountered in multiply
  #!/Users/warren/miniconda3scipy/bin/python
Out[308]: 2.0256691479491108e-321

In all the above cases, the result is a denormal number. Shouldn't they all trigger the underflow warning?

Does this have something to do with the underlying bit patterns of the multiplicands?

In [314]: bin(np.array(r).view(np.uint64))
Out[314]: '0b1111010100000000000000000000000000000000000000000000000000000'

In [315]: bin(np.array(1.5*r).view(np.uint64))
Out[315]: '0b1111010101000000000000000000000000000000000000000000000000000'

In [316]: bin(np.array(1.6*r).view(np.uint64))
Out[316]: '0b1111010101001100110011001100110011001100110011001100110011010'

[mdickinson]

Under IEEE 754 default exception handling (IEEE 754-2008 section 7.5), the underflow flag is raised only if the exact result of the arithmetic operation was not representable. As the standard says

This is the only case in this standard of an exception signal receiving default handling that does not raise the corresponding flag.

I suspect that that's the cause of the apparent inconsistency here.

[WarrenWeckesser]

Thanks Mark, that does look like what is happening.

In the stackoverflow question, the warnings actually arise during a division of complex numbers.

In [7]: np.seterr(under='warn')
Out[7]: {'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'ignore'}

In [8]: x1 = 1 + 1j / (1 << 533)

In [9]: y1 = 1.1*x1

In [10]: np.divide(1, x1)  # No underflow warning.
Out[10]: (1-3.556413999176124e-161j)

In [11]: np.divide(1, y1)
/Users/warren/miniconda3scipy/bin/ipython:1: RuntimeWarning: underflow encountered in true_divide
  #!/Users/warren/miniconda3scipy/bin/python
Out[11]: (0.90909090909090906-3.2331036356146581e-161j)

I think this is because the multiplication underflow described earlier occurs in an intermediate calculation in https://github.com/numpy/numpy/blob/master/numpy/core/src/umath/loops.c.src#L2517. In this case, the underflow warning seems spurious--it is the result of an implementation detail, and probably of no use to a user. Does anyone agree, and if so, is it worth the effort to fix it?

[seberg]

Closing, seems correct as per IEEE (as per comments). And frankly, I don't think we should mess with it if was incorrect... making sure that underflow flags are correct across compilers and math libs sounds like the stuff nightmares are made from.

EDIT: Hmmpf, I suppose it is correct that possible the division shouldn't do this (i.e. it drops out of divmod or so). But, I still don't think we should spend any cycles on it, so leaving it closed.