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