Mark Dickinson <dicki...@gmail.com> added the comment:
> am still studying the new one
Sorry - I should have added at least _some_ comments to it.
Here's the underlying principle, which I think of as "The magic of
round-to-odd":
Suppose x is a positive real number and e is an integer satisfying 2^e <= x.
Then assuming IEEE 754 binary64 floating-point format, the quantity:
2^(e-54) * to-odd(x / 2^(e-54))
rounds to the same value as x does under _any_ of the standard IEEE 754
rounding modes, including the round-ties-to-even rounding mode that we care
about here.
Here, to-odd is the function R -> Z that maps each integer to itself, but maps
each non-integral real number x to the *odd* integer nearest x. (So for example
all of 2.1, 2.3, 2.9, 3.0, 3.1, 3.9 map to 3.)
This works because (x / 2^(e-54)) gives us an integer with at least 55 bits:
the 53 bits we'll need in the eventual significand, a rounding bit, and then
the to-odd supplies a "sticky" bit that records inexactness. Note that the
principle works in the subnormal range too - no additional tricks are needed
for that. In that case we just end up wastefully computing a few more bits than
we actually _need_ to determine the correctly-rounded value.
The code applies this principle to the case x = sqrt(n/m) and
e = (n.bit_length() - m.bit_length() - 1)//2. The condition 2^e <= x holds
because:
2^(n.bit_length() - 1) <= n
m < 2^m.bit_length()
so
2^(n.bit_length() - 1 - m.bit_length()) < n/m
and taking square roots gives
2^((n.bit_length() - 1 - m.bit_length())/2) < √(n/m)
so taking e = (n.bit_length() - 1 - m.bit_length())//2, we have
2^e <= 2^((n.bit_length() - 1 - m.bit_length())/2) < √(n/m)
Now putting q = e - 54, we need to compute
2^q * round-to-odd(√(n/m) / 2^q)
rounded to a float. The two branches both do this computation, by moving 2^q
into either the numerator or denominator of the fraction as appropriate
depending on the sign of q.
----------
_______________________________________
Python tracker <rep...@bugs.python.org>
<https://bugs.python.org/issue45876>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe:
https://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
