robf...@gmail.com <robfi680@gmail.com> schrieb:
>
> Is there any use for approximately equals?

Absolutely. In numerical code,

if (abs(a-b) < eps) call hooray("Heureka!")

is ubiquitous.

> And how to go about determining
> approximate equality.

That is rather difficult and depends a lot on your application.
Assume you are looking for a solution to f(x) = 0, which is
well-behaved (i.e. has a large derivative around the root), then
a few times machine precision would serve.

If you are looking for a minimum which has the "usual" behavior
like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
precise - hoping for better accurady than sqrt(epsilon), where
1+epsilon is the smallest number which does not equal 1, will
not be fulfilled.

>I was thinking of an equality operator that takes in a
> number of significant bits that must be equal. So, if there was a double
> precision comparison, it could say equal to within 23 significant bits.

Even an abdiff instruction like r = abs(a-b) could save an
instruction.

Or, depending on how many bits there are for specifying
a register and if there are implicit flags, an instruction
that sets a flag depending on the sign of abs(a-b) < eps
could be helpful.

Also quite useful could be a hypot function, which would calculate
sqrt(a**2 + b** 2) to machine accuracy. People are still writing
papers about how to do it in the absence of such an instruction
(Article 9 in ACM Transactions on Mathematical Software 47, 1
(2021))

Thomas Koenig wrote:
> robf...@gmail.com <robfi680@gmail.com> schrieb:
>>
>> Is there any use for approximately equals?
>
> Absolutely. In numerical code,
>
> if (abs(a-b) < eps) call hooray("Heureka!")
>
> is ubiquitous.
>
>> And how to go about determining
>> approximate equality.
>
> That is rather difficult and depends a lot on your application.
> Assume you are looking for a solution to f(x) = 0, which is
> well-behaved (i.e. has a large derivative around the root), then
> a few times machine precision would serve.
>
> If you are looking for a minimum which has the "usual" behavior
> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
> precise - hoping for better accurady than sqrt(epsilon), where
> 1+epsilon is the smallest number which does not equal 1, will
> not be fulfilled.
>
>> I was thinking of an equality operator that takes in a
>> number of significant bits that must be equal. So, if there was a double
>> precision comparison, it could say equal to within 23 significant bits.

Actually, having a compare operation which was capable of returning the
number of identical bits, even when near a flipping boundary (0.99999..
vs 1.0000...) would be quite useful, and probably faster than the
traditional relative offset measure:

rel_err = (estimate-exact)/exact.

>
> Even an abdiff instruction like r = abs(a-b) could save an
> instruction.
>
> Or, depending on how many bits there are for specifying
> a register and if there are implicit flags, an instruction
> that sets a flag depending on the sign of abs(a-b) < eps
> could be helpful.
>
> Also quite useful could be a hypot function, which would calculate
> sqrt(a**2 + b** 2) to machine accuracy. People are still writing
> papers about how to do it in the absence of such an instruction
> (Article 9 in ACM Transactions on Mathematical Software 47, 1
> (2021))

The new Augmented(Addition/Multiplication) fp operations would make this
a lot easier, i.e. everything inside the sqrt() call would be exact, you
just need to create a not-quite-quad sqrt function to get exact rounding.

Terje

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

On 4/4/2022 10:19 AM, Terje Mathisen wrote:
> Thomas Koenig wrote:
>> robf...@gmail.com <robfi680@gmail.com> schrieb:
>>>
>>> Is there any use for approximately equals?
>>
>> Absolutely.  In numerical code,
>>
>>    if (abs(a-b) < eps) call hooray("Heureka!")
>>
>> is ubiquitous.
>>
>>> And how to go about determining
>>> approximate equality.
>>
>> That is rather difficult and depends a lot on your application.
>> Assume you are looking for a solution to f(x) = 0, which is
>> well-behaved (i.e. has a large derivative around the root), then
>> a few times machine precision would serve.
>>
>> If you are looking for a minimum which has the "usual" behavior
>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
>> precise - hoping for better accurady than sqrt(epsilon), where
>> 1+epsilon is the smallest number which does not equal 1, will
>> not be fulfilled.
>>
>>> I was thinking of an equality operator that takes in a
>>> number of significant bits that must be equal. So, if there was a double
>>> precision comparison, it could say equal to within 23 significant bits.
>
> Actually, having a compare operation which was capable of returning the
> number of identical bits, even when near a flipping boundary (0.99999..
> vs 1.0000...) would be quite useful, and probably faster than the
> traditional relative offset measure:
>
>  rel_err = (estimate-exact)/exact.

I noticed that there is space for that in the results register of the
FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
in hardware, and it can be done in parallel with determining the other
bits in the result, perhaps you can convince Mitch to add it.

--
- Stephen Fuld
(e-mail address disguised to prevent spam)

On Wednesday, March 30, 2022 at 1:19:52 PM UTC-5, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > On Wednesday, March 30, 2022 at 12:46:22 AM UTC-5, Thomas Koenig wrote:
> >> MitchAlsup <Mitch...@aol.com> schrieb:
> >> > On Tuesday, March 29, 2022 at 1:55:03 PM UTC-5, Thomas Koenig wrote:
> >> >> MitchAlsup <Mitch...@aol.com> schrieb:
> >> >> > On Tuesday, March 29, 2022 at 11:45:08 AM UTC-5, Thomas Koenig wrote:
> >> >>
> >> >> >> I've looked at the optimum Remez polynomial for 1/x in the range
> >> >> >> of 1..2. It is 70/33 - 16/11 * x + 32/99 * x**2, with a weight of x
> >> >> >> (so optimizing for the relative error).
> >> >> >>
> >> >> >> The maximum relative error is 1/99, reached at four points in the
> >> >> >> interval [1..2], so around 6.62 bits of minimum accuracy.
> >> >> ><
> >> >> > I stated way above::---------------------------------------------------------------------------------------
> >> >> > Equivalent to::
> >> >> ><
> >> >> > .33333×x^2 - 1.5×x + 2.1666
> >> >> ><
> >> >> > The Chebychev 2nd order Coefficients on the interval [1..2) are:
> >> >> ><
> >> >> > 0.32323232×x^2 -0.48484848×x + 0.66666667
> This is the formula we are discussing.
> >> >> There is somethig wrong with your formula. f(1) would be
> >> >> 0.32323232 - 0.48484848 + 0.66666667 = 0.50505051, which does
> >> >> not even closely approximate 1/1 = 1. Is there some rescaling
> >> >> somewhere?
> >> >> ><
> >> >> > with 6.63 bits of accuracy.
> >> >> ><-------------------------------------------------------------------------------------------------------------------
> >> >> > Why does your Remez get less precision than Chebychev ?
> >> >> Absolute or relative precision? As I wrote above, I used relative
> >> ><
> >> > Mine was relative (too)
> >> Easy enough to check. What was the actual formula you
> >> arrived at? Not the one you wrote above, obviously.
> ><
> >=IF(AA29=0,-99,LOG(ABS(AA29),2))
> ... and this is a formula for calculating an error.
>
> We seem to have a disconnect here somewhere.
>
> We were discussing a second-degree polynomial approximation to 1/x.
> The formula you gave doesn't have anything close to 6 bits of
> accuracy, it has _zero_ bits of accuracy (as I showed for the case
> of x=1 above).
>
> So I'm not sure what you are in fact approximating. I would, however,
> like to compare your Chebyshev formula for 1/x with the Remez formula
> I gave above.
>
> Could you give that Chebyshev formula for 1/x?
<
p(x) = 0.32323232×x^2 -0.48484848×x + 0.66666667
<
on the interval 1.0..2.0
<
Which is compared to 1000 randomized points against 1/x
<
The highest error encountered has 6.63 bits of precision.
<---------------------------------
But the x I put into the polynomial was the distance from the mid-point of the interval = 1.5
<---------------------------------
Making the formula::
<
p(x) = 0.32323232×(x-1.5)^2 -0.48484848×(x-1.5) + 0.66666667
<
x in the interval {1.0..2.0} polynomial argument in the range {-0.5..+0.5}
<
<
<
It is so easy to get lost in eXcel spreadsheet equations.

On Wednesday, March 30, 2022 at 2:48:44 PM UTC-5, Michael S wrote:
> On Wednesday, March 30, 2022 at 9:19:52 PM UTC+3, Thomas Koenig wrote:
> > MitchAlsup <Mitch...@aol.com> schrieb:
> > > On Wednesday, March 30, 2022 at 12:46:22 AM UTC-5, Thomas Koenig wrote:
> > >> MitchAlsup <Mitch...@aol.com> schrieb:
> > >> > On Tuesday, March 29, 2022 at 1:55:03 PM UTC-5, Thomas Koenig wrote:
> > >> >> MitchAlsup <Mitch...@aol.com> schrieb:
> > >> >> > On Tuesday, March 29, 2022 at 11:45:08 AM UTC-5, Thomas Koenig wrote:
> > >> >>
> > >> >> >> I've looked at the optimum Remez polynomial for 1/x in the range
> > >> >> >> of 1..2. It is 70/33 - 16/11 * x + 32/99 * x**2, with a weight of x
> > >> >> >> (so optimizing for the relative error).
> > >> >> >>
> > >> >> >> The maximum relative error is 1/99, reached at four points in the
> > >> >> >> interval [1..2], so around 6.62 bits of minimum accuracy.
> > >> >> ><
> > >> >> > I stated way above::---------------------------------------------------------------------------------------
> > >> >> > Equivalent to::
> > >> >> ><
> > >> >> > .33333×x^2 - 1.5×x + 2.1666
> > >> >> ><
> > >> >> > The Chebychev 2nd order Coefficients on the interval [1..2) are:
> > >> >> ><
> > >> >> > 0.32323232×x^2 -0.48484848×x + 0.66666667
> > This is the formula we are discussing.
> > >> >> There is somethig wrong with your formula. f(1) would be
> > >> >> 0.32323232 - 0.48484848 + 0.66666667 = 0.50505051, which does
> > >> >> not even closely approximate 1/1 = 1. Is there some rescaling
> > >> >> somewhere?
> > >> >> ><
> > >> >> > with 6.63 bits of accuracy.
> > >> >> ><-------------------------------------------------------------------------------------------------------------------
> > >> >> > Why does your Remez get less precision than Chebychev ?
> > >> >> Absolute or relative precision? As I wrote above, I used relative
> > >> ><
> > >> > Mine was relative (too)
> > >> Easy enough to check. What was the actual formula you
> > >> arrived at? Not the one you wrote above, obviously.
> > ><
> > >=IF(AA29=0,-99,LOG(ABS(AA29),2))
> > ... and this is a formula for calculating an error.
> >
> > We seem to have a disconnect here somewhere.
> >
> > We were discussing a second-degree polynomial approximation to 1/x.
> > The formula you gave doesn't have anything close to 6 bits of
> > accuracy, it has _zero_ bits of accuracy (as I showed for the case
> > of x=1 above).
> >
> He gave poly on interval [2:1] instead of [1:2].
> After affine transformation (in Matlab/Octave polyaffine(pp, [3 -1]) ) it's
> the same poly as yours, except that Mitch rounded coefficients to 8 decimal places.
<
I choose only to "print" the first 8 decimal digits, eXcel choose to round.
<
> > So I'm not sure what you are in fact approximating. I would, however,
> > like to compare your Chebyshev formula for 1/x with the Remez formula
> > I gave above.
> >
> > Could you give that Chebyshev formula for 1/x?

On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
> > Thomas Koenig wrote:
> >> robf...@gmail.com <robf...@gmail.com> schrieb:
> >>>
> >>> Is there any use for approximately equals?
> >>
> >> Absolutely. In numerical code,
> >>
> >> if (abs(a-b) < eps) call hooray("Heureka!")
> >>
> >> is ubiquitous.
> >>
> >>> And how to go about determining
> >>> approximate equality.
> >>
> >> That is rather difficult and depends a lot on your application.
> >> Assume you are looking for a solution to f(x) = 0, which is
> >> well-behaved (i.e. has a large derivative around the root), then
> >> a few times machine precision would serve.
> >>
> >> If you are looking for a minimum which has the "usual" behavior
> >> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
> >> precise - hoping for better accurady than sqrt(epsilon), where
> >> 1+epsilon is the smallest number which does not equal 1, will
> >> not be fulfilled.
> >>
> >>> I was thinking of an equality operator that takes in a
> >>> number of significant bits that must be equal. So, if there was a double
> >>> precision comparison, it could say equal to within 23 significant bits.
> >
> > Actually, having a compare operation which was capable of returning the
> > number of identical bits, even when near a flipping boundary (0.99999..
> > vs 1.0000...) would be quite useful, and probably faster than the
> > traditional relative offset measure:
> >
> > rel_err = (estimate-exact)/exact.
> I noticed that there is space for that in the results register of the
> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
> in hardware, and it can be done in parallel with determining the other
> bits in the result, perhaps you can convince Mitch to add it.
>
Yes, space if there, in fact only 20-ish of the 64-bits are defined.
<
But I am not sure at all how to do a divide in the FCMP instruction--
That is I don't see what you want put there.
>
> --
> - Stephen Fuld
> (e-mail address disguised to prevent spam)

On 4/4/2022 3:32 PM, MitchAlsup wrote:
> On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
>> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
>>> Thomas Koenig wrote:
>>>> robf...@gmail.com <robf...@gmail.com> schrieb:
>>>>>
>>>>> Is there any use for approximately equals?
>>>>
>>>> Absolutely. In numerical code,
>>>>
>>>> if (abs(a-b) < eps) call hooray("Heureka!")
>>>>
>>>> is ubiquitous.
>>>>
>>>>> And how to go about determining
>>>>> approximate equality.
>>>>
>>>> That is rather difficult and depends a lot on your application.
>>>> Assume you are looking for a solution to f(x) = 0, which is
>>>> well-behaved (i.e. has a large derivative around the root), then
>>>> a few times machine precision would serve.
>>>>
>>>> If you are looking for a minimum which has the "usual" behavior
>>>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
>>>> precise - hoping for better accurady than sqrt(epsilon), where
>>>> 1+epsilon is the smallest number which does not equal 1, will
>>>> not be fulfilled.
>>>>
>>>>> I was thinking of an equality operator that takes in a
>>>>> number of significant bits that must be equal. So, if there was a double
>>>>> precision comparison, it could say equal to within 23 significant bits.
>>>
>>> Actually, having a compare operation which was capable of returning the
>>> number of identical bits, even when near a flipping boundary (0.99999..
>>> vs 1.0000...) would be quite useful, and probably faster than the
>>> traditional relative offset measure:
>>>
>>> rel_err = (estimate-exact)/exact.
>> I noticed that there is space for that in the results register of the
>> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
>> in hardware, and it can be done in parallel with determining the other
>> bits in the result, perhaps you can convince Mitch to add it.
>>
> Yes, space if there, in fact only 20-ish of the 64-bits are defined.
> <
> But I am not sure at all how to do a divide in the FCMP instruction--
> That is I don't see what you want put there.

I may have misunderstood Terje. I thought he wanted the number of
mantissa bits that are identical between the two operands. This would
just be an integer value up to 52.

--
- Stephen Fuld
(e-mail address disguised to prevent spam)

On Monday, April 4, 2022 at 5:41:12 PM UTC-5, Stephen Fuld wrote:
> On 4/4/2022 3:32 PM, MitchAlsup wrote:
> > On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
> >> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
> >>> Thomas Koenig wrote:
> >>>> robf...@gmail.com <robf...@gmail.com> schrieb:
> >>>>>
> >>>>> Is there any use for approximately equals?
> >>>>
> >>>> Absolutely. In numerical code,
> >>>>
> >>>> if (abs(a-b) < eps) call hooray("Heureka!")
> >>>>
> >>>> is ubiquitous.
> >>>>
> >>>>> And how to go about determining
> >>>>> approximate equality.
> >>>>
> >>>> That is rather difficult and depends a lot on your application.
> >>>> Assume you are looking for a solution to f(x) = 0, which is
> >>>> well-behaved (i.e. has a large derivative around the root), then
> >>>> a few times machine precision would serve.
> >>>>
> >>>> If you are looking for a minimum which has the "usual" behavior
> >>>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
> >>>> precise - hoping for better accurady than sqrt(epsilon), where
> >>>> 1+epsilon is the smallest number which does not equal 1, will
> >>>> not be fulfilled.
> >>>>
> >>>>> I was thinking of an equality operator that takes in a
> >>>>> number of significant bits that must be equal. So, if there was a double
> >>>>> precision comparison, it could say equal to within 23 significant bits.
> >>>
> >>> Actually, having a compare operation which was capable of returning the
> >>> number of identical bits, even when near a flipping boundary (0.99999..
> >>> vs 1.0000...) would be quite useful, and probably faster than the
> >>> traditional relative offset measure:
> >>>
> >>> rel_err = (estimate-exact)/exact.
> >> I noticed that there is space for that in the results register of the
> >> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
> >> in hardware, and it can be done in parallel with determining the other
> >> bits in the result, perhaps you can convince Mitch to add it.
> >>
> > Yes, space if there, in fact only 20-ish of the 64-bits are defined.
> > <
> > But I am not sure at all how to do a divide in the FCMP instruction--
> > That is I don't see what you want put there.
> I may have misunderstood Terje. I thought he wanted the number of
> mantissa bits that are identical between the two operands. This would
> just be an integer value up to 52.
<
64 XOR gate 1-delay
64-bit FF1.... 4-delays
64-unary to 6-bit binary converter ... 4-delays
So, it seems to fit
<
But some lower end models are dependent on the comparison not taking more
than 8 gates because CMP-BB is CoIssued.
<
Still, its close to fitting.
> --
> - Stephen Fuld
> (e-mail address disguised to prevent spam)

On Tuesday, April 5, 2022 at 1:41:12 AM UTC+3, Stephen Fuld wrote:
> On 4/4/2022 3:32 PM, MitchAlsup wrote:
> > On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
> >> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
> >>> Thomas Koenig wrote:
> >>>> robf...@gmail.com <robf...@gmail.com> schrieb:
> >>>>>
> >>>>> Is there any use for approximately equals?
> >>>>
> >>>> Absolutely. In numerical code,
> >>>>
> >>>> if (abs(a-b) < eps) call hooray("Heureka!")
> >>>>
> >>>> is ubiquitous.
> >>>>
> >>>>> And how to go about determining
> >>>>> approximate equality.
> >>>>
> >>>> That is rather difficult and depends a lot on your application.
> >>>> Assume you are looking for a solution to f(x) = 0, which is
> >>>> well-behaved (i.e. has a large derivative around the root), then
> >>>> a few times machine precision would serve.
> >>>>
> >>>> If you are looking for a minimum which has the "usual" behavior
> >>>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
> >>>> precise - hoping for better accurady than sqrt(epsilon), where
> >>>> 1+epsilon is the smallest number which does not equal 1, will
> >>>> not be fulfilled.
> >>>>
> >>>>> I was thinking of an equality operator that takes in a
> >>>>> number of significant bits that must be equal. So, if there was a double
> >>>>> precision comparison, it could say equal to within 23 significant bits.
> >>>
> >>> Actually, having a compare operation which was capable of returning the
> >>> number of identical bits, even when near a flipping boundary (0.99999..
> >>> vs 1.0000...) would be quite useful, and probably faster than the
> >>> traditional relative offset measure:
> >>>
> >>> rel_err = (estimate-exact)/exact.
> >> I noticed that there is space for that in the results register of the
> >> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
> >> in hardware, and it can be done in parallel with determining the other
> >> bits in the result, perhaps you can convince Mitch to add it.
> >>
> > Yes, space if there, in fact only 20-ish of the 64-bits are defined.
> > <
> > But I am not sure at all how to do a divide in the FCMP instruction--
> > That is I don't see what you want put there.
> I may have misunderstood Terje. I thought he wanted the number of
> mantissa bits that are identical between the two operands. This would
> just be an integer value up to 52.
> --

I.e. __clz(bits_of(a) xor bits_of(b)) ?
I'm not sure how useful it is. Such criterion misses case of very close FP numbers that belong to different octaves.

__clz(abs(bits_of(a) - bits_of(b))) is more useful but even here we have sore cases of tiny positives vs tiny negatives that
are going to report big difference when in reality the difference is small.

So, what we really want is a bit more complicated. During subtraction we want to treat bits of FP numbers as integers in sign-magnitude
format. Then abs() and finally CLZ.

Thinking about it, nearly all required logic should be preset anyway, as part of FP adder, but it can be located in such way that
if we reuse it then our proximity estimator will have a "FP-like" latency of 2-3 clock cycles, instead of desirable latency of 1.

> - Stephen Fuld
> (e-mail address disguised to prevent spam)

On 4/4/2022 4:22 PM, Michael S wrote:
> On Tuesday, April 5, 2022 at 1:41:12 AM UTC+3, Stephen Fuld wrote:
>> On 4/4/2022 3:32 PM, MitchAlsup wrote:
>>> On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
>>>> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
>>>>> Thomas Koenig wrote:
>>>>>> robf...@gmail.com <robf...@gmail.com> schrieb:
>>>>>>>
>>>>>>> Is there any use for approximately equals?
>>>>>>
>>>>>> Absolutely. In numerical code,
>>>>>>
>>>>>> if (abs(a-b) < eps) call hooray("Heureka!")
>>>>>>
>>>>>> is ubiquitous.
>>>>>>
>>>>>>> And how to go about determining
>>>>>>> approximate equality.
>>>>>>
>>>>>> That is rather difficult and depends a lot on your application.
>>>>>> Assume you are looking for a solution to f(x) = 0, which is
>>>>>> well-behaved (i.e. has a large derivative around the root), then
>>>>>> a few times machine precision would serve.
>>>>>>
>>>>>> If you are looking for a minimum which has the "usual" behavior
>>>>>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
>>>>>> precise - hoping for better accurady than sqrt(epsilon), where
>>>>>> 1+epsilon is the smallest number which does not equal 1, will
>>>>>> not be fulfilled.
>>>>>>
>>>>>>> I was thinking of an equality operator that takes in a
>>>>>>> number of significant bits that must be equal. So, if there was a double
>>>>>>> precision comparison, it could say equal to within 23 significant bits.
>>>>>
>>>>> Actually, having a compare operation which was capable of returning the
>>>>> number of identical bits, even when near a flipping boundary (0.99999..
>>>>> vs 1.0000...) would be quite useful, and probably faster than the
>>>>> traditional relative offset measure:
>>>>>
>>>>> rel_err = (estimate-exact)/exact.
>>>> I noticed that there is space for that in the results register of the
>>>> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
>>>> in hardware, and it can be done in parallel with determining the other
>>>> bits in the result, perhaps you can convince Mitch to add it.
>>>>
>>> Yes, space if there, in fact only 20-ish of the 64-bits are defined.
>>> <
>>> But I am not sure at all how to do a divide in the FCMP instruction--
>>> That is I don't see what you want put there.
>> I may have misunderstood Terje. I thought he wanted the number of
>> mantissa bits that are identical between the two operands. This would
>> just be an integer value up to 52.
>> --
>
> I.e. __clz(bits_of(a) xor bits_of(b)) ?
> I'm not sure how useful it is. Such criterion misses case of very close FP numbers that belong to different octaves.
>
> __clz(abs(bits_of(a) - bits_of(b))) is more useful but even here we have sore cases of tiny positives vs tiny negatives that
> are going to report big difference when in reality the difference is small.
>
> So, what we really want is a bit more complicated. During subtraction we want to treat bits of FP numbers as integers in sign-magnitude
> format. Then abs() and finally CLZ.
>
> Thinking about it, nearly all required logic should be preset anyway, as part of FP adder, but it can be located in such way that
> if we reuse it then our proximity estimator will have a "FP-like" latency of 2-3 clock cycles, instead of desirable latency of 1.

I am not qualified to talk about the numerical analysis parts of this.
I was only reacting to Terje's and I assumed the utility of what he
asked for. Whether it is the right solution or not - I can only say
"Let's you and him fight!". :-)

--
- Stephen Fuld
(e-mail address disguised to prevent spam)

Stephen Fuld wrote:
> On 4/4/2022 3:32 PM, MitchAlsup wrote:
>> On Monday, April 4, 2022 at 2:00:18 PM UTC-5, Stephen Fuld wrote:
>>> On 4/4/2022 10:19 AM, Terje Mathisen wrote:
>>>> Thomas Koenig wrote:
>>>>> robf...@gmail.com <robf...@gmail.com> schrieb:
>>>>>>
>>>>>> Is there any use for approximately equals?
>>>>>
>>>>> Absolutely.  In numerical code,
>>>>>
>>>>>     if (abs(a-b) < eps) call hooray("Heureka!")
>>>>>
>>>>> is ubiquitous.
>>>>>
>>>>>> And how to go about determining
>>>>>> approximate equality.
>>>>>
>>>>> That is rather difficult and depends a lot on your application.
>>>>> Assume you are looking for a solution to f(x) = 0, which is
>>>>> well-behaved (i.e. has a large derivative around the root), then
>>>>> a few times machine precision would serve.
>>>>>
>>>>> If you are looking for a minimum which has the "usual" behavior
>>>>> like f(x) = a*(x-x0)**2 + b, then the optimum will be far less
>>>>> precise - hoping for better accurady than sqrt(epsilon), where
>>>>> 1+epsilon is the smallest number which does not equal 1, will
>>>>> not be fulfilled.
>>>>>
>>>>>> I was thinking of an equality operator that takes in a
>>>>>> number of significant bits that must be equal. So, if there was a
>>>>>> double
>>>>>> precision comparison, it could say equal to within 23 significant
>>>>>> bits.
>>>>
>>>> Actually, having a compare operation which was capable of returning the
>>>> number of identical bits, even when near a flipping boundary (0.99999..
>>>> vs 1.0000...) would be quite useful, and probably faster than the
>>>> traditional relative offset measure:
>>>>
>>>>   rel_err = (estimate-exact)/exact.
>>> I noticed that there is space for that in the results register of the
>>> FMCP instruction in Mitch's My 66000. Since it seem pretty simple to do
>>> in hardware, and it can be done in parallel with determining the other
>>> bits in the result, perhaps you can convince Mitch to add it.
>>>
>> Yes, space if there, in fact only 20-ish of the 64-bits are defined.
>> <
>> But I am not sure at all how to do a divide in the FCMP instruction--
>> That is I don't see what you want put there.
>
> I may have misunderstood Terje.  I thought he wanted the number of
> mantissa bits that are identical between the two operands. This would
> just be an integer value up to 52.

A bit harder: I want an approximate binary log of the difference between
two fp numbers, scaled by the first. I.e. for comparing 1.000 vs 0.999 I
want something around 10 (exact value=9.965784285) indicating the number
of significant bits.

So
double similarity(double a, double b)
{ // Return approximately log2(a/(a-b)) = log2(a)-log2(a-b)
return log2(a)-log2(abs(a-b))
}

This function would return 0 when b is twice as large as a, i.e. zero
equal mantissa bits, then negative numbers for being further away in
magnitude, with similarity(1,100) returning -6.629, or close to the
negative of similarity(1,1.01) which is 6.644. This indicates that the
exponents differ by 6 or 7 bits which is correct.

Terje

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

Michael S wrote:
> I.e. __clz(bits_of(a) xor bits_of(b)) ? I'm not sure how useful it
> is. Such criterion misses case of very close FP numbers that belong
> to different octaves.
>
> __clz(abs(bits_of(a) - bits_of(b))) is more useful but even here we
> have sore cases of tiny positives vs tiny negatives that are going to
> report big difference when in reality the difference is small.
>
> So, what we really want is a bit more complicated. During subtraction
> we want to treat bits of FP numbers as integers in sign-magnitude
> format. Then abs() and finally CLZ.
>
> Thinking about it, nearly all required logic should be preset anyway,
> as part of FP adder, but it can be located in such way that if we
> reuse it then our proximity estimator will have a "FP-like" latency
> of 2-3 clock cycles, instead of desirable latency of 1.

This is exactly right, and will at least be far faster than the real
thing which requires first a real subtraction, then one or two binary
logs. :-)

Terje

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

On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
<already5chosen@yahoo.com> wrote:

>I'm starting to understand the problem.
>It is related to implementation of long double, more specifically, to
>settings of x87 control world. Somehow, under WSL, x87 control world
>is set to 53-bit precision (default Windows settings). Of course,
>quadmath, being Linux-originated, expects x87 control word to be set
>to 64-bit precision. Who is at fault, kernel (Microsoft) or userland
>(Fedora/SUSE) ? I'd guess, one is going to blame another and vice
>versa.

Blame should fall on the library - if code needs control flags set in
some particular way, it should make sure they are set correctly.
Relying on system defaults because they happen to line up with
expectations is just lazy.

Also the Windows default x87 setting IS 64-bit (full width) precision.
There are differences in how Windows and Linux handle save/restore of
x87 state, but if it's all working correctly, I don't see how that
would cause any problem.

I know nothing of quadmath, so this may be a stupid question ... but
documentation seems to indicate the library dates from ~2010. x86
compilers were using SSE (SIMD) for floating point long before that -
at least since Pentium 4. So why does quadmath use the x87?

George

> So
> double similarity(double a, double b)
> { // Return approximately log2(a/(a-b)) = log2(a)-log2(a-b)
> return log2(a)-log2(abs(a-b))
> }

Matches my intuition, indeed.

> This function would return 0 when b is twice as large as a, i.e. zero equal
> mantissa bits, then negative numbers for being further away in magnitude,
> with similarity(1,100) returning -6.629, or close to the negative of
> similarity(1,1.01) which is 6.644. This indicates that the exponents differ
> by 6 or 7 bits which is correct.

I suspect it would be good enough to focus on positive values (and
return 0 rather than a negative value). But admittedly, negative values
are probably a small matter of subtracting the two exponents, which
should be easy to do in a single cycle.

Regarding stuffing it into My66k's generic FMCP instruction, I'm not
sure it's such a great idea since you then have to extract the
corresponding bits to make use of them. AFAIK the current FMCP returns
a vector of booleans (and these can be extracted "for free" as part of
the PRED instruction, for example) so this "small int" would feel out
of place.

Stefan

George Neuner <gneuner2@comcast.net> writes:
>On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
><already5chosen@yahoo.com> wrote:
>
>>I'm starting to understand the problem.
>>It is related to implementation of long double, more specifically, to
>>settings of x87 control world. Somehow, under WSL, x87 control world
>>is set to 53-bit precision (default Windows settings). Of course,
>>quadmath, being Linux-originated, expects x87 control word to be set
>>to 64-bit precision. Who is at fault, kernel (Microsoft) or userland
>>(Fedora/SUSE) ? I'd guess, one is going to blame another and vice
>>versa.
>
>Blame should fall on the library - if code needs control flags set in
>some particular way, it should make sure they are set correctly.
>Relying on system defaults because they happen to line up with
>expectations is just lazy.

The ABI
<https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
clearly specifies in Table 3.3 that at process initialization the PC
bits in the x87 Floating-Point control word are set to 11 ("Double
extended precision", i.e., 64-bit rounding of the mantissa).

If WSL does not satisfy this requirement, it's obviously WSL that's to
blame. If a user-level program relies on the specification instead of
unnecessarily setting the bits itself, one may consider it lazy, but I
consider it smart not to do unnecessary busywork.

Apparently Microsoft was not willing to do the grunt work of
implementing every bit of the Linux kernel/user interface, that's why
they gave up on the WSL approach and switched to the VM approach of
WSL2. Are Microsoft lazy, or are they smart?

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

MitchAlsup <MitchAlsup@aol.com> writes:

> On Wednesday, March 30, 2022 at 1:19:52 PM UTC-5, Thomas Koenig wrote:

[...]

>> Could you give that Chebyshev formula for 1/x?
>
> p(x) = 0.32323232 x^2 -0.48484848 x + 0.66666667
> on the interval 1.0..2.0
> Which is compared to 1000 randomized points against 1/x
> The highest error encountered has 6.63 bits of precision.

If p(x) =
0.32323232 * (x-1.5) * (x-1.5)
- 0.48484848 * (x-1.5)
+ 0.66666667

observe that, at x = 1, so 1/x = 1, the expression

- log( (1 - p(1)) / 1 ) / log( 2 )

is less than 6.6294, and that,
at x = 2, so 1/x = 0.5, the expression

- log( (p(2) - 0.5) / 0.5) / log( 2 )

also is less than 6.6294

On Tuesday, April 5, 2022 at 7:32:10 PM UTC+3, Anton Ertl wrote:
> George Neuner <gneu...@comcast.net> writes:
> >On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
> ><already...@yahoo.com> wrote:
> >
> >>I'm starting to understand the problem.
> >>It is related to implementation of long double, more specifically, to
> >>settings of x87 control world. Somehow, under WSL, x87 control world
> >>is set to 53-bit precision (default Windows settings). Of course,
> >>quadmath, being Linux-originated, expects x87 control word to be set
> >>to 64-bit precision. Who is at fault, kernel (Microsoft) or userland
> >>(Fedora/SUSE) ? I'd guess, one is going to blame another and vice
> >>versa.
> >
> >Blame should fall on the library - if code needs control flags set in
> >some particular way, it should make sure they are set correctly.
> >Relying on system defaults because they happen to line up with
> >expectations is just lazy.
> The ABI
> <https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
> clearly specifies in Table 3.3 that at process initialization the PC
> bits in the x87 Floating-Point control word are set to 11 ("Double
> extended precision", i.e., 64-bit rounding of the mantissa).
>
> If WSL does not satisfy this requirement, it's obviously WSL that's to
> blame. If a user-level program relies on the specification instead of
> unnecessarily setting the bits itself, one may consider it lazy, but I
> consider it smart not to do unnecessary busywork.
>

Except that WSL didn't promise ABI-level compatibility with iAMD64 Linux.
All it promised is a decent source-level compatibility, say, to the same or slightly
better level as between iAMD64 Linux and aarch64 Linux.
And, in this case, they didn't deliver on the promise.
The correct and the simplest thing to do for them would be, probably, to force
"long double==double" at the level of compiler and standard headers.
Then, at least for this library, everything will work correctly, same way it is
working correctly on aarch64 or POWER-LE.

> Apparently Microsoft was not willing to do the grunt work of
> implementing every bit of the Linux kernel/user interface, that's why
> they gave up on the WSL approach and switched to the VM approach of
> WSL2. Are Microsoft lazy, or are they smart?

Both, IMHO, but more the former.

>
> - anton
> --
> 'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
> Mitch Alsup, <c17fcd89-f024-40e7...@googlegroups.com>

On Tuesday, April 5, 2022 at 4:26:24 PM UTC+3, George Neuner wrote:
> On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
> <already...@yahoo.com> wrote:
>
> >I'm starting to understand the problem.
> >It is related to implementation of long double, more specifically, to
> >settings of x87 control world. Somehow, under WSL, x87 control world
> >is set to 53-bit precision (default Windows settings). Of course,
> >quadmath, being Linux-originated, expects x87 control word to be set
> >to 64-bit precision. Who is at fault, kernel (Microsoft) or userland
> >(Fedora/SUSE) ? I'd guess, one is going to blame another and vice
> >versa.
>
> Blame should fall on the library - if code needs control flags set in
> some particular way, it should make sure they are set correctly.
> Relying on system defaults because they happen to line up with
> expectations is just lazy.

The author of the library relied on autoconf infrastructure that
is supposed to provides HAVE_SQRTL.
I was not able to find out what HAVE_SQRTL should mean, but
obviously it does not mean what the author thought it means.
But even if the author was more pedantic and did things in a Standard C
way, i.e. doing something like '#if LDBL_MANT_DIG >= 64' it would
still produce a wrong result, because LDBL_MANT_DIG is defined as 64!

>
> Also the Windows default x87 setting IS 64-bit (full width) precision.
> There are differences in how Windows and Linux handle save/restore of
> x87 state, but if it's all working correctly, I don't see how that
> would cause any problem.
>
>
> I know nothing of quadmath, so this may be a stupid question ... but
> documentation seems to indicate the library dates from ~2010.

But very likely is built on earlier work.

> x86 compilers were using SSE (SIMD) for floating point long before that -
> at least since Pentium 4. So why does quadmath use the x87?
>

That, may be, an interesting and important question, but it is
related to performance rather than correctness. I.e. it is of lower importance.

Also, it's not fare to say that quadmath *use x87*.
It tries to use whatever compiler provides for a 'long double' type on given platform.
Which, on x86-64 Linux and Windows/MSYS2 happens to be x87 and on "native
x86-64 Windows" happens to be SSE2.
But on WSL it happenned to be 'castrated x87'.

> George

Michael S <already5chosen@yahoo.com> writes:
>On Tuesday, April 5, 2022 at 7:32:10 PM UTC+3, Anton Ertl wrote:
>> The ABI
>> <https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
>> clearly specifies in Table 3.3 that at process initialization the PC
>> bits in the x87 Floating-Point control word are set to 11 ("Double
>> extended precision", i.e., 64-bit rounding of the mantissa).
>>
>> If WSL does not satisfy this requirement, it's obviously WSL that's to
>> blame. If a user-level program relies on the specification instead of
>> unnecessarily setting the bits itself, one may consider it lazy, but I
>> consider it smart not to do unnecessary busywork.
>>
>
>Except that WSL didn't promise ABI-level compatibility with iAMD64 Linux.

They didn't? If that is the case, then anyone who runs a Linux binary
on WSL is to blame.

However, I certainly had the impression that binary compatibility was
the point and the goal of WSL. Admittedly, my interest in Windows 10
features is close to zero, so I could easily have misunderstood
something.

>All it promised is a decent source-level compatibility, say, to the same or slightly
>better level as between iAMD64 Linux and aarch64 Linux.

[Citation needed] Assuming this claim is true, they made a whole
toolchain for WSL, and every user of WSL would have to rebuild every
piece of code with that toolchain? I never heard of that. It also
looks to me like something that very few would use; those who wanted
to build Unix programs for Windows already had Cygwin, and that
produces applications that run across many more Windows versions, and
without the user having to install WSL. I don't see a reason why
Microsoft would go for source-level compatibility, especially given
that they provided little or no support to Cygwin AFAIK.

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

On Wednesday, April 6, 2022 at 10:15:04 AM UTC+3, Anton Ertl wrote:
> Michael S <already...@yahoo.com> writes:
> >On Tuesday, April 5, 2022 at 7:32:10 PM UTC+3, Anton Ertl wrote:
> >> The ABI
> >> <https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
> >> clearly specifies in Table 3.3 that at process initialization the PC
> >> bits in the x87 Floating-Point control word are set to 11 ("Double
> >> extended precision", i.e., 64-bit rounding of the mantissa).
> >>
> >> If WSL does not satisfy this requirement, it's obviously WSL that's to
> >> blame. If a user-level program relies on the specification instead of
> >> unnecessarily setting the bits itself, one may consider it lazy, but I
> >> consider it smart not to do unnecessary busywork.
> >>
> >
> >Except that WSL didn't promise ABI-level compatibility with iAMD64 Linux..
> They didn't? If that is the case, then anyone who runs a Linux binary
> on WSL is to blame.
>
> However, I certainly had the impression that binary compatibility was
> the point and the goal of WSL. Admittedly, my interest in Windows 10
> features is close to zero, so I could easily have misunderstood
> something.
> >All it promised is a decent source-level compatibility, say, to the same or slightly
> >better level as between iAMD64 Linux and aarch64 Linux.
> [Citation needed]

It's hard to find any official statements about anything WSL.
Even docs are often wrong. It looks like the whole enterprise was
never thought out to the ultimate end.
I wouldn't be surprised that [if WSL2 works well then] in a year or two
Microsoft will try to pretend that WSL never existed.

On one hand, they did promise to run "Linux binaries".
On the other hand, they admitted that syscall-level compatibility is not 100%.
Which means that you can't reliably run Linux binaries, doesn't it?

> Assuming this claim is true, they made a whole
> toolchain for WSL, and every user of WSL would have to rebuild every
> piece of code with that toolchain?

It seems, you're right.
$ gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/10/lto-wrapper
OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa:hsa
OFFLOAD_TARGET_DEFAULT=1
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Debian 10.2.1-6' --with-bugurl=file:///usr/share/doc/gcc-10/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-10 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-10-Km9U7s/gcc-10-10.2.1/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-10-Km9U7s/gcc-10-10.2.1/debian/tmp-gcn/usr,hsa --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu --with-build-config=bootstrap-lto-lean --enable-link-mutex
Thread model: posix
Supported LTO compression algorithms: zlib zstd
gcc version 10.2.1 20210110 (Debian 10.2.1-6)

So, the tool is a regular x86_64 Linux tool.

> I never heard of that. It also
> looks to me like something that very few would use; those who wanted
> to build Unix programs for Windows already had Cygwin, and that
> produces applications that run across many more Windows versions, and
> without the user having to install WSL.

Cygwin is compatible, but very slow.
Its user interface (cat&past, drag&drop etc) feels old.
And its package manager sucks.

MSYS2 is quite fast and its user interface is not bad.
But it is very significantly incompatible.
And its package manager sucks.
Or, may be, by now, after many breakages and fixes, package manager
by itself does not suck, but relatively to popular distros too few
packages available in pre-build form.

WSL was supposed to be as fast as MSYS2, at least as long as files are accessed
on its native filesystem and even access to files on Windows side was supposed
to be not much slower than MSYS2.
Its user interface is about the same as MSYS2.
And WSL was supposed to be at least as compatible at source level as Cygwin..
But most importantly, it was supposed to give you package management of major distros.

> I don't see a reason why
> Microsoft would go for source-level compatibility, especially given
> that they provided little or no support to Cygwin AFAIK.

It looks like you are correct with regard to Microsoft's intentions.
At least since ~2017.

> - anton
> --
> 'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
> Mitch Alsup, <c17fcd89-f024-40e7...@googlegroups.com>

On Tue, 05 Apr 2022 16:16:51 GMT, anton@mips.complang.tuwien.ac.at
(Anton Ertl) wrote:

>George Neuner <gneuner2@comcast.net> writes:
>
>>Blame should fall on the library - if code needs control flags set in
>>some particular way, it should make sure they are set correctly.
>>Relying on system defaults because they happen to line up with
>>expectations is just lazy.
>
>The ABI
><https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
>clearly specifies in Table 3.3 that at process initialization the PC
>bits in the x87 Floating-Point control word are set to 11 ("Double
>extended precision", i.e., 64-bit rounding of the mantissa).
>
>If WSL does not satisfy this requirement, it's obviously WSL that's to
>blame. If a user-level program relies on the specification instead of
>unnecessarily setting the bits itself, one may consider it lazy, but I
>consider it smart not to do unnecessary busywork.

It's hardly "busywork". We're not talking about a program - we're
talking about a /library/, which can't know whether the user program
has changed from the default.

A library MUST explicitly set the environment it needs, and it must
save/restore the program environment in case that is different.

George

In article <vofo4hh9npgd0vaefo51khtndt80g440if@4ax.com>,
gneuner2@comcast.net (George Neuner) wrote:
> On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
> <already5chosen@yahoo.com> wrote:
> >I'm starting to understand the problem.
> >It is related to implementation of long double, more specifically,
> >to settings of x87 control world. Somehow, under WSL, x87 control
> >world is set to 53-bit precision (default Windows settings). Of
> >course, quadmath, being Linux-originated, expects x87 control word
> >to be set to 64-bit precision. Who is at fault, kernel (Microsoft)
> >or userland (Fedora/SUSE)? I'd guess, one is going to blame
> >another and vice versa.
>
> Blame should fall on the library - if code needs control flags set
> in some particular way, it should make sure they are set correctly.
> Relying on system defaults because they happen to line up with
> expectations is just lazy.

Changing those flags takes time - on the microsecond scale,because it
often requires the pipeline to empty before the change - so you don't
want to be doing it at every entry to a quad-precision library, where
you'd hope that some operations would be timed on the nanosecond scale.
It's probably better to document what's needed, so that the application
can control things.

> Also the Windows default x87 setting IS 64-bit (full width)
> precision.

Default in what circumstances? In the Microsoft C/C++ run-time
environment, the default has been ordinary double precision, with a
53-bit mantissa, since 1996 to my certain knowledge. The hardware default
for the x87 registers is long double, with a 64-bit mantissa.

I've had painful experiences with the precision being set wrong in both
directions:

A customer for the modeller I work on had bought in a Visual Basic
interpreter from a third party supplier, which they used as a macro
language in their application. (This was before Microsoft offered an
embeddable interpreter.) This third-party interpreter changed the default
precision to long double "because it was more accurate" although they
didn't seem to have any idea about its effects. Because this was new to
us, it took us a week or so to figure out what the problem was, during
which the customer was chasing us twice a day with a lot of strong
language. When we could prove what the problem was, they rolled all that
anger back up, called the head of the Basic company, and let it all off
at him. The interpreter was changed, swiftly, and gave no further trouble.

A couple of years later, I was updating the x86 Linux build standard.
Linux had shown lots of fiddling small differences in test results, which
were a nuisance. The C/C++ on Linux was inheriting the hardware default
of long double evaluation, but variables saves into memory were saved as
doubles. This adds enough jitter to be noticeable: setting it to use
double evaluation made it far more consistent with other platforms. We
did that in the test harness and documented it for customers.

A few years after that, someone was writing example code for Microsoft's
DirectX viewing library. But everything you viewed in it was distorted,
and rapidly got smushed into really weird shapes. I grabbed the processor
manual, broke into the debugger, checked the control flags, found they
were set to single-precision evaluation, and got accused of witchcraft.
DirectX on 32-bit x86 sets the precision to single on entry, and by
default doesn't change it back. This naturally messes up any other
libraries that were expecting double. Microsoft could not explain why
they were doing this, but I suspect there was some assembler in there,
written by someone who had left. They did point out an option that would
make DirectX save and restore the floating-point controls, but warned us
that it would make things slow. They hadn't understood the timescales:
humans interact with DirectX on timescales of tenths of seconds, so a
couple of microseconds is irrelevant.

John

On Wednesday, April 6, 2022 at 4:01:27 PM UTC-5, John Dallman wrote:
> In article <vofo4hh9npgd0vaef...@4ax.com>,
> gneu...@comcast.net (George Neuner) wrote:
> > On Sun, 3 Apr 2022 15:22:34 -0700 (PDT), Michael S
> > <already...@yahoo.com> wrote:
> > >I'm starting to understand the problem.
> > >It is related to implementation of long double, more specifically,
> > >to settings of x87 control world. Somehow, under WSL, x87 control
> > >world is set to 53-bit precision (default Windows settings). Of
> > >course, quadmath, being Linux-originated, expects x87 control word
> > >to be set to 64-bit precision. Who is at fault, kernel (Microsoft)
> > >or userland (Fedora/SUSE)? I'd guess, one is going to blame
> > >another and vice versa.
> >
> > Blame should fall on the library - if code needs control flags set
> > in some particular way, it should make sure they are set correctly.
> > Relying on system defaults because they happen to line up with
> > expectations is just lazy.
<
> Changing those flags takes time - on the microsecond scale,because it
> often requires the pipeline to empty before the change - so you don't
> want to be doing it at every entry to a quad-precision library, where
> you'd hope that some operations would be timed on the nanosecond scale.
> It's probably better to document what's needed, so that the application
> can control things.
<
Yes, the x87 flag design is a true debacle..........
<
> > Also the Windows default x87 setting IS 64-bit (full width)
> > precision.
> Default in what circumstances? In the Microsoft C/C++ run-time
> environment, the default has been ordinary double precision, with a
> 53-bit mantissa, since 1996 to my certain knowledge. The hardware default
> for the x87 registers is long double, with a 64-bit mantissa.
>
> I've had painful experiences with the precision being set wrong in both
> directions:
<
>--------------------------------------Computer architecture Lesson-------------------------------------
>
> A customer for the modeller I work on had bought in a Visual Basic
> interpreter from a third party supplier, which they used as a macro
> language in their application. (This was before Microsoft offered an
> embeddable interpreter.) This third-party interpreter changed the default
> precision to long double "because it was more accurate" although they
> didn't seem to have any idea about its effects. Because this was new to
> us, it took us a week or so to figure out what the problem was, during
> which the customer was chasing us twice a day with a lot of strong
> language. When we could prove what the problem was, they rolled all that
> anger back up, called the head of the Basic company, and let it all off
> at him. The interpreter was changed, swiftly, and gave no further trouble.
>
>
> A couple of years later, I was updating the x86 Linux build standard.
> Linux had shown lots of fiddling small differences in test results, which
> were a nuisance. The C/C++ on Linux was inheriting the hardware default
> of long double evaluation, but variables saves into memory were saved as
> doubles. This adds enough jitter to be noticeable: setting it to use
> double evaluation made it far more consistent with other platforms. We
> did that in the test harness and documented it for customers.
>
> A few years after that, someone was writing example code for Microsoft's
> DirectX viewing library. But everything you viewed in it was distorted,
> and rapidly got smushed into really weird shapes. I grabbed the processor
> manual, broke into the debugger, checked the control flags, found they
> were set to single-precision evaluation, and got accused of witchcraft.
> DirectX on 32-bit x86 sets the precision to single on entry, and by
> default doesn't change it back. This naturally messes up any other
> libraries that were expecting double. Microsoft could not explain why
> they were doing this, but I suspect there was some assembler in there,
> written by someone who had left. They did point out an option that would
> make DirectX save and restore the floating-point controls, but warned us
> that it would make things slow. They hadn't understood the timescales:
> humans interact with DirectX on timescales of tenths of seconds, so a
> couple of microseconds is irrelevant.
<
>-------------------------------------------End of Lesson---------------------------------------------------------
<
3 (more) anecdotal reasons one does not want to change FP widths with
mode bits. You want to change FP widths by using a different OpCode !!
<
>--------------------------------------------------------------------------------------------------------------------------
<
But notice that the FPU can accept a new set of mode bits every cycle !!
This is the only way for HyperThreading to work (reasonably) !!!
<
But nobody ever got around to fixing how slow updating the mode bits is.
>
> John

George Neuner <gneuner2@comcast.net> writes:
>On Tue, 05 Apr 2022 16:16:51 GMT, anton@mips.complang.tuwien.ac.at
>(Anton Ertl) wrote:
>>The ABI
>><https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
>>clearly specifies in Table 3.3 that at process initialization the PC
>>bits in the x87 Floating-Point control word are set to 11 ("Double
>>extended precision", i.e., 64-bit rounding of the mantissa).
>>
>>If WSL does not satisfy this requirement, it's obviously WSL that's to
>>blame. If a user-level program relies on the specification instead of
>>unnecessarily setting the bits itself, one may consider it lazy, but I
>>consider it smart not to do unnecessary busywork.
>
>It's hardly "busywork". We're not talking about a program - we're
>talking about a /library/, which can't know whether the user program
>has changed from the default.

We were actually talking about the program that worked nicely on Linux,
and failed on WSL.

But yes, the library may not be doing what it should do in a program
that actively sets the precision to some other values than long
double.

I can imagine some measures that would deal with the problem in the
usual case: E.g., on first invocation the function for setting the
precision control changes the vector of the quad-math functions to
versions that set the precision control for their own needs (and save
and restore the x87 control word). However, that would not work for
programs that don't use these functions, but set the control word
directly (in assembly language), so unless the library is documented
as requiring the function calls for changing the precision control,
the library would still be deficient.

- anton
--
'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
Mitch Alsup, <c17fcd89-f024-40e7-a594-88a85ac10d20o@googlegroups.com>

On Thursday, April 7, 2022 at 12:03:32 PM UTC+3, Anton Ertl wrote:
> George Neuner <gneu...@comcast.net> writes:
> >On Tue, 05 Apr 2022 16:16:51 GMT, an...@mips.complang.tuwien.ac.at
> >(Anton Ertl) wrote:
> >>The ABI
> >><https://raw.githubusercontent.com/wiki/hjl-tools/x86-psABI/x86-64-psABI-1.0.pdf>
> >>clearly specifies in Table 3.3 that at process initialization the PC
> >>bits in the x87 Floating-Point control word are set to 11 ("Double
> >>extended precision", i.e., 64-bit rounding of the mantissa).
> >>
> >>If WSL does not satisfy this requirement, it's obviously WSL that's to
> >>blame. If a user-level program relies on the specification instead of
> >>unnecessarily setting the bits itself, one may consider it lazy, but I
> >>consider it smart not to do unnecessary busywork.
> >
> >It's hardly "busywork". We're not talking about a program - we're
> >talking about a /library/, which can't know whether the user program
> >has changed from the default.
> We were actually talking about the program that worked nicely on Linux,
> and failed on WSL.
>

More I think about it, less I see a justification for WSL control word defaults.

> But yes, the library may not be doing what it should do in a program
> that actively sets the precision to some other values than long
> double.
>
> I can imagine some measures that would deal with the problem in the
> usual case: E.g., on first invocation the function for setting the
> precision control changes the vector of the quad-math functions to
> versions that set the precision control for their own needs (and save
> and restore the x87 control word). However, that would not work for
> programs that don't use these functions, but set the control word
> directly (in assembly language), so unless the library is documented
> as requiring the function calls for changing the precision control,
> the library would still be deficient.
>

In this particular case, the only winning strategy is to refuse to play.
I.e. quadmath library should be and could be coded without any use of 80-bit
FP and with very minimalist use of 64-bit FP. On modern 64-bit x86 (and,
I suppose, on modern ARM and POWER) it's not only the most robust way,
precision-wise, but also the fastest.

But I can imagine other cases where judicious use of 80-bit FP is really
beneficiary.

> - anton
> --
> 'Anyone trying for "industrial quality" ISA should avoid undefined behavior.'
> Mitch Alsup, <c17fcd89-f024-40e7...@googlegroups.com>