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


Reply via email to