In a42 the almost_equal test is probably wrong

[ikirill]

I mean this test:

cc03ba6#diff-9471aab4b037f803f311ca8d2d35fe8fR558

It potentially uses the wrong type (Float64, which is not necessarily the type it's working with), and also it seems to use realmin in place of eps, where eps is the appropriate value.

I want to change it to

abs(x - y) < 32 * max(eps(x), eps(y))

which is probably what was intended anyway.

[ikirill]

Okay, I now found it in the original paper, the condition there is specified as x-y == 0. The test abs(x-y) <= 32 * 1e-308 essentially accomplishes the same thing, because either x-y will be on the order of max(abs(x),abs(y))*eps(), or exactly zero, so I think this is identical, slightly roundabout way of checking x-y==0.

See http://www.netlib.org/toms-2014-06-10/748 at the lines starting with PROF=(FA-FB)...

Perhaps the reason for writing it this way was the infamous maxim "don't test floating point numbers for equality", so it was a way of "fixing" x==y?