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