On Wednesday, April 20, 2022 at 11:12:39 AM UTC+3, Terje Mathisen wrote:
> Michael S wrote:
> > On Wednesday, April 20, 2022 at 1:33:50 AM UTC+3, MitchAlsup wrote:
> >> On Tuesday, April 19, 2022 at 5:05:28 PM UTC-5, Michael S wrote:
> >>> On Tuesday, April 19, 2022 at 10:42:08 PM UTC+3, MitchAlsup wrote:
> >>>> On Tuesday, April 19, 2022 at 3:10:26 AM UTC-5, Terje Mathisen wrote:
> >>>>> Marcus wrote:
> >>>>>> One more thing that I discovered about a week ago:
> >>>>>>
> >>>>>> In order to support the IEEE 754 semantics of 1 / 0 == Inf, you actually
> >>>>>> *need* a special Newtom-Raphson step instruction that handles the
> >>>>>> special case where x = 0 and y = Inf.
> >>>>> The obvious way to handle this, if you have a little bit of HW help in
> >>>>> the form of an input number(s) classifier, is to take out all the
> >>>>> special cases and merge these results with what you get from the normal
> >>>>> algorithm.
> >>>>>>
> >>>>>> Consider the approximation of y = 1 / x.
> >>>>>>
> >>>>>> A single NR step then becomes: y' = y * (2 - x * y)
> >>>> <
> >>>> Note: y ~= 1/x -- thus algebraically x×y ~= 1
> >>>> Here: 0×Inf should be 1 not NaN.
> >>>> In fact LIM[x->0]( x×~1/x ) -> ~1
> >>> 0*Inf = NaN for the same reason 0/0=Nan and Inf/Inf=Nan.
> >>> In all 3 cases the results are "unknown" and FP numbers have no better encoding for "unknown" than NaN.
> >> <
> >> They may be unknown "in general" but in this particular N-R iteration the LIM is known to be 1.000000000000000000000000000.
> >> <
> >
> > The rules are "in general". They are not created to satisfy needs of one particular not extremely important algorithm.
> > BTW, after thinking a bit about it I feel that existing inf rules are rather "too practical", scarifying mathematical rigor
> > in favor of non-proven practicality.
> > In particular, being today (tonight) in rigorous mood, I'd prefer (+inf - Finite_Positive_Number) to be NaN rather than +inf.
> > Because we don't know if result is bigger than MAX_DBL or, may be, smaller.
> > The same goes for (inf / Finite_Number_Above_One).
> > Etc...
> Sorry, but this is simply wrong:
>
> What you suggest here is to abandon knowledge we do have in favor of
> making NaN even more sticky than it currently is.

Yes, unproven practicality vs rigor.
With existing rules we preserve (still, not always) supposedly practically important
information that there was an overflow in a chain of calculation.
But we are no longer guaranteed that +inf represent a number that is bigger
than maximal representative finite number.

We could hah had both if we had 2 or 3 different infs. Like:
zero-inf - there was division by zero in the chain
huge-inf - bigger than MAX_DBL
ovf-inf - there was overflow in the chain, but we no longer sure that the current value is > MAX_DBL
IEEE NaN coding space even has room for all this goodies. But it's too late..

>
> If you have an algorithm which can spuriously overflow in the middle,
> but which with a different evaluation order, or higher intermediate
> precision/exponent range, would be OK, then giving Inf as the result
> still provides important information: This calculation overflowed
> somewhere in the middle.
>
> Always returning NaN instead means that you have lost this info, now the
> NaN could have been due to many different errors, including 0/0 etc.

Not always.
With my last-nights's-rules under many circumstances (Inf + Finite) or (Inf * Finite)
or (Inf/Finite) would still produce Inf. But yes, if the chain is long it's unlikely.
Even more unlikely (less likely ?) than DFP retaining info about # of significant digits.

>
> All the 754 ops have been defined to return the best possible answer,
> when taken in isolation.
>
> The original (1978) idea was to handle exponent overflow (or underflow)
> via a trap handler that would save the current exponent, bias it enough
> to bring the current operation into legal range, and then inspect and
> correct the final result based on any such bias operations performed
> during the processing.
>
> Afaik, nobody is seriously doing this for production code.
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Michael S wrote:
> On Wednesday, April 20, 2022 at 11:02:51 AM UTC+3, Terje Mathisen wrote:
>> MitchAlsup wrote:
>>> On Tuesday, April 19, 2022 at 5:05:28 PM UTC-5, Michael S wrote:
>>>> On Tuesday, April 19, 2022 at 10:42:08 PM UTC+3, MitchAlsup wrote:
>>>>> On Tuesday, April 19, 2022 at 3:10:26 AM UTC-5, Terje Mathisen wrote:
>>>>>> Marcus wrote:
>>>>>>> One more thing that I discovered about a week ago:
>>>>>>>
>>>>>>> In order to support the IEEE 754 semantics of 1 / 0 == Inf, you actually
>>>>>>> *need* a special Newtom-Raphson step instruction that handles the
>>>>>>> special case where x = 0 and y = Inf.
>>>>>> The obvious way to handle this, if you have a little bit of HW help in
>>>>>> the form of an input number(s) classifier, is to take out all the
>>>>>> special cases and merge these results with what you get from the normal
>>>>>> algorithm.
>>>>>>>
>>>>>>> Consider the approximation of y = 1 / x.
>>>>>>>
>>>>>>> A single NR step then becomes: y' = y * (2 - x * y)
>>>>> <
>>>>> Note: y ~= 1/x -- thus algebraically x×y ~= 1
>>>>> Here: 0×Inf should be 1 not NaN.
>>>>> In fact LIM[x->0]( x×~1/x ) -> ~1
>>>> 0*Inf = NaN for the same reason 0/0=Nan and Inf/Inf=Nan.
>>>> In all 3 cases the results are "unknown" and FP numbers have no better encoding for "unknown" than NaN.
>>> <
>>> They may be unknown "in general" but in this particular N-R iteration the LIM is known to be 1.000000000000000000000000000.
>> Which is why you cannot use normal (individual/standard) FP ops to do
>> this, you need something higher level which knows about the special
>> rules in effect right here. I.e. like that ARM reciprocal step opcode.
>>
>> Alternatively, and the solution I prefer, you handle all special inputs
>> in parallel with the normal algorithm and merge at the end.
>>
>> I.e. if the divisor is 0 then the 754 standard specifies exactly what
>> the output should be for all the various dividend/divisor combinations:
>>
>> non-zero/zero -> inf with sign according to normal (XOR) rules.
>>
>> zero/zero -> NaN
>> inf/zero -> NaN
>
> Are you sure about it?
> I didn't look in the IEEE docs, but both my intuition and gcc compiler disagree.
> Both say:
> inf/zero -> inf with sign according to normal (XOR) rules.

Oops, sorry!

Yeah, dividing anything except Zero or NaN by zero gives you Inf.

Mea Culpa.

Terje

--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"

On Wednesday, April 20, 2022 at 2:42:28 PM UTC+3, Terje Mathisen wrote:
> Michael S wrote:
> > On Wednesday, April 20, 2022 at 11:02:51 AM UTC+3, Terje Mathisen wrote:
> >> MitchAlsup wrote:
> >>> On Tuesday, April 19, 2022 at 5:05:28 PM UTC-5, Michael S wrote:
> >>>> On Tuesday, April 19, 2022 at 10:42:08 PM UTC+3, MitchAlsup wrote:
> >>>>> On Tuesday, April 19, 2022 at 3:10:26 AM UTC-5, Terje Mathisen wrote:
> >>>>>> Marcus wrote:
> >>>>>>> One more thing that I discovered about a week ago:
> >>>>>>>
> >>>>>>> In order to support the IEEE 754 semantics of 1 / 0 == Inf, you actually
> >>>>>>> *need* a special Newtom-Raphson step instruction that handles the
> >>>>>>> special case where x = 0 and y = Inf.
> >>>>>> The obvious way to handle this, if you have a little bit of HW help in
> >>>>>> the form of an input number(s) classifier, is to take out all the
> >>>>>> special cases and merge these results with what you get from the normal
> >>>>>> algorithm.
> >>>>>>>
> >>>>>>> Consider the approximation of y = 1 / x.
> >>>>>>>
> >>>>>>> A single NR step then becomes: y' = y * (2 - x * y)
> >>>>> <
> >>>>> Note: y ~= 1/x -- thus algebraically x×y ~= 1
> >>>>> Here: 0×Inf should be 1 not NaN.
> >>>>> In fact LIM[x->0]( x×~1/x ) -> ~1
> >>>> 0*Inf = NaN for the same reason 0/0=Nan and Inf/Inf=Nan.
> >>>> In all 3 cases the results are "unknown" and FP numbers have no better encoding for "unknown" than NaN.
> >>> <
> >>> They may be unknown "in general" but in this particular N-R iteration the LIM is known to be 1.000000000000000000000000000.
> >> Which is why you cannot use normal (individual/standard) FP ops to do
> >> this, you need something higher level which knows about the special
> >> rules in effect right here. I.e. like that ARM reciprocal step opcode.
> >>
> >> Alternatively, and the solution I prefer, you handle all special inputs
> >> in parallel with the normal algorithm and merge at the end.
> >>
> >> I.e. if the divisor is 0 then the 754 standard specifies exactly what
> >> the output should be for all the various dividend/divisor combinations:
> >>
> >> non-zero/zero -> inf with sign according to normal (XOR) rules.
> >>
> >> zero/zero -> NaN
> >> inf/zero -> NaN
> >
> > Are you sure about it?
> > I didn't look in the IEEE docs, but both my intuition and gcc compiler disagree.
> > Both say:
> > inf/zero -> inf with sign according to normal (XOR) rules.
> Oops, sorry!
>
> Yeah, dividing anything except Zero or NaN by zero gives you Inf.
>
> Mea Culpa.
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

BTW,
I posted implementations of fast quad-precision square root and division
as followups to Thomas Koenig's bug report.
If your are interested, look here:
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101