Re: [PATCH] Fix old thinko in formula to compute sweight in numeric_sqrt().

Hi,

On Tue, Jan 31, 2023, at 14:40, Dean Rasheed wrote:
> That's still not right. If you want a proper mathematically justified
> formula, it's fairly easy to derive.
...
> or equivalently, in code with truncated integer division:
>
>     if (arg.weight >= 0)
>         sweight = arg.weight * DEC_DIGITS / 2 + 1;
>     else
>         sweight = 1 - (1 - arg.weight * DEC_DIGITS) / 2;

Beautiful! Thank you for magnificent analysis and extraordinary good explanation, now I finally get it.

> When DEC_DIGITS = 4, this formula also reduces to sweight = 2 *
> arg.weight + 1, but neither gcc nor clang is smart enough to spot that
> (clang doesn't simplify your formula either, BTW).

Oh, that's a shame. :-(

> So even though I
> believe that the above is mathematically correct, and won't change any
> results for DEC_DIGITS = 4, I'm still hesitant to use it, because it
> will have a (small) performance impact, and I don't believe it does
> anything to improve code readability (and certainly not without an
> explanatory comment).

I also think the performance impact no matter how small isn't worth it,
but a comment based on your comments would be very valuable IMO.

Below is an attempt at summarising your text, and to avoid the performance impact,
maybe an #if so we get the general correct formula for DEC_DIGITS 1 or 2,
and the reduced hand-optimised form for DEC_DIGITS 4?
That could also improve readabilty, since readers perhaps more easily would see
the relation between sweight and arg.weight, for the only DEC_DIGITS case we care about.

Suggestion:

    /*
     * Here we approximate the decimal weight of the square root (sweight),
     * given the NBASE-weight (arg.weight) of the input argument.
     *
     * The lower bound of the decimal weight of the input argument is used to
     * calculate the decimal weight of the square root, with integer division
     * being truncated.
     * 
     * The general formula is:
     *
     *     sweight = floor((n+1) / 2)
     * 
     * In our case, since the base is NBASE, not 10, and since we only
     * require an approximation, we don't take the trouble to compute n
     * exactly, we just use the fact that it lies in the range
     * 
     *     arg.weight * DEC_DIGITS + 1 <= n <= (arg.weight + 1) * DEC_DIGITS
     *
     * Since we want to ensure at least a certain number of significant
     * digits in the result, we're only interested in the lower bound.
     * Plugging that into the formula above gives:
     * 
     *     sweight >= floor(arg.weight * DEC_DIGITS / 2 + 1)
     *
     * Which leads us to the formula below with truncated integer division.
     */
#if DEC_DIGITS == 1 || DEC_DIGITS == 2

    if (arg.weight >= 0)
        sweight = arg.weight * DEC_DIGITS / 2 + 1;
    else
        sweight = 1 - (1 - arg.weight * DEC_DIGITS) / 2;

#elif DEC_DIGITS == 4

    /*
     * Neither gcc nor clang is smart enough to spot that
     * the formula above neatly reduces to the below
     * when DEC_DIGITS == 4.
     */
    sweight = 2 * arg.weight + 1;

#else
#error unsupported NBASE
#endif

/Joel