Common block mystery

Help
6

As you noticed previously in SLSQP, old foxes were quite self-aware of floating-point spacing

Hence without the necessary tooling back then, they adopted some -isms that you see a lot in old code. Another common one is res = 0.5D0 + (0.5D0 - x) from the 70s.

The case they are trying to avoid is the following (sorry for the NumPy code);

>>> import numpy as np
>>> zzz = np.float64(1234567891011.123456789)

>>> zzz + 1e-6 == zzz
True

>>> zzz + 1e-5 == zzz
True

>>> zzz + 1e-1 == zzz
False

In other words, when is the added value cannot be distinguished in that particular dtype for that magnitude? The switching point from True to False is the zzz*del and del is the distance to the next number in that particular dtype for that magnitude (between two consecutive powers of 2). NumPy also has a convenience function with the same name numpy.spacing.

>>> s = np.spacing(zzz)
>>> s
np.float64(0.000244140625)

>>> s / zzz
np.float64(1.9775390788759812e-16)

>>> np.log2(s / zzz)
np.float64(-52.1671433135915)

The test number is always a power of two to get the relevant spacing in that particular magnitude hence you would see 0.5, 1.0 or other powers of two depending on what they are testing.

2 Likes