Author mark.dickinson
Recipients mark.dickinson, rhettinger, steven.daprano
Date 2021-11-23.17:47:27
SpamBayes Score -1.0
Marked as misclassified Yes
Message-id <1637689647.84.0.163294872354.issue45876@roundup.psfhosted.org>
In-reply-to
Content
> Does the technique you had in mind involve testing 1 ulp up or down to see whether its square is closer to the input?

Kinda sorta. Below is some code: it's essentially just pure integer operations, with a last minute conversion to float (implicit in the division in the case of the second branch). And it would need to be better tested, documented, and double-checked to be viable.


def isqrt_rto(n):
    """
    Square root of n, rounded to the nearest integer using round-to-odd.
    """
    a = math.isqrt(n)
    return a | (a*a != n)


def isqrt_frac_rto(n, m):
    """
    Square root of n/m, rounded to the nearest integer using round-to-odd.
    """
    quotient, remainder = divmod(isqrt_rto(4*n*m), 2*m)
    return quotient | bool(remainder)


def sqrt_frac(n, m):
    """
    Square root of n/m as a float, correctly rounded.
    """
    quantum = (n.bit_length() - m.bit_length() - 1) // 2 - 54
    if quantum >= 0:
        return float(isqrt_frac_rto(n, m << 2 * quantum) << quantum)
    else:
        return isqrt_frac_rto(n << -2 * quantum, m) / (1 << -quantum)
History
Date User Action Args
2021-11-23 17:47:27mark.dickinsonsetrecipients: + mark.dickinson, rhettinger, steven.daprano
2021-11-23 17:47:27mark.dickinsonsetmessageid: <1637689647.84.0.163294872354.issue45876@roundup.psfhosted.org>
2021-11-23 17:47:27mark.dickinsonlinkissue45876 messages
2021-11-23 17:47:27mark.dickinsoncreate