Why does complex floating-point division underflow weirdly with NumPy?

Question

Consider this code:

import numpy
numpy.seterr(under='warn')
x1 = 1 + 1j / (1 << 533)
x2 = 1 - 1j / (1 << 533)
y1 = x1 * 1.1
y2 = x2 * 1.1
z1 = x1 / 1.1
z2 = x2 / 1.1
print(numpy.divide(1, x1))  #              1-3.55641399918e-161j  # OK
print(numpy.divide(1, x2))  #              1+3.55641399918e-161j  # OK
print(numpy.divide(1, y1))  # 0.909090909091-3.23310363561e-161j  # underflow
print(numpy.divide(1, y2))  # 0.909090909091+3.23310363561e-161j  # underflow
print(numpy.divide(1, z1))  #            1.1-3.91205539909e-161j  # underflow
print(numpy.divide(1, z2))  #            1.1+3.91205539909e-161j  # underflow

The underflow doesn't seem to make sense no matter how I look at it. Like Wikipedia says,

Underflow is a condition in a computer program where the result of a calculation is a number of smaller absolute value than the computer can actually store in memory on its CPU.

But clearly the computer is capable of storing numbers in the general proximity of the values in question, so the definition doesn't seem at all consistent with the behavior I'm seeing here.

Could someone explain why exactly some of these are giving underflows but others are not?
Is this correct behavior or is it a bug?

Answer 1

To do this math, you need to first "scale", so-to-speak, the 3.23310363561e-161. That is what triggers the underflow

  0.909090909091
- 3.23310363561e-161 =
--------------------

  0.909090909091 x 10^0
- 3.23310363561 x 10^-161 =
-------------------------

  0.909090909091 x 10^0
- 0.000...                 0323310363561 =
    ^^^^^^^^^^^^^^^^^^^^^^^^
    160 0s
  --------------------------------------