On Wednesday, March 30, 2022 at 9:27:22 AM UTC+3, Terje Mathisen wrote:
> Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> >>
> >> Of course, it's 1ULP.
> >> But, according to my understanding, sqrt() is one of those very
> >> few primitives for which IEEE-754 not just recommends > correct
> >> rounding, but requires correct rounding.
> >
> > This is now https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101 .
> >
> Exact rounding rules are easy: All 5 core primitives
> (FADD/FSUB/FMUL/FDIV/FSQRT) for which it was known back in 1978 that it
> was relatively easy to always provide exact rounding, are in fact
> required to do so.
>
> Having a square root (for any supported fp size) which does not obey
> this is in fact a breaking bug.
>

Also various conversions to/from integer, right?
Probably, also binary scaling, including cases of normal-to-subnormal?

BTW, does it mean that "optional core" primitives, in particular, FMADD and FRSQRT,
if provided, are allowed to give inexact rounding?

> It is of course perfectly legal to provide an alternative "fastmath"
> library which disobeys these rules, but at that point you are not
> providing IEEE 754 math. I.e. Cray did this at some point, right?
>
> Similar issues with IBM hex math, but both of them was before 1978.
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

Thomas Koenig wrote:
> Terje Mathisen <terje.mathisen@tmsw.no> schrieb:
>> Thomas Koenig wrote:
>>> Michael S <already5chosen@yahoo.com> schrieb:
>>>>
>>>> Of course, it's 1ULP.
>>>> But, according to my understanding, sqrt() is one of those very
>>>> few primitives for which IEEE-754 not just recommends > correct
>>>> rounding, but requires correct rounding.
>>>
>>> This is now https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101 .
>>>
>>
>> Exact rounding rules are easy: All 5 core primitives
>> (FADD/FSUB/FMUL/FDIV/FSQRT) for which it was known back in 1978 that it
>> was relatively easy to always provide exact rounding, are in fact
>> required to do so.
>>
>> Having a square root (for any supported fp size) which does not obey
>> this is in fact a breaking bug.
>
> It's a wrong-code bug (and I have marked it as such). Compilers and
> libraries are known to have these, even thougn, in an ideal world,
> they should not exist.
>
> It is certainly too late to fix for gcc12 (regression-mode only),
> but gcc13 should be doable.
>
> Hmm... for those who are in the know about IEEE standards (which
> cost money :-(): Does sqrt need to follow rounding modes?
>
YES!

All those core ops have to abide by the current/actual rounding mode.

Trancendentals etc are different, there you will almost always try for
the best possible result, rounded according to nearest_or_even.

Terje

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

Michael S wrote:
> On Wednesday, March 30, 2022 at 9:27:22 AM UTC+3, Terje Mathisen wrote:
>> Thomas Koenig wrote:
>>> Michael S <already...@yahoo.com> schrieb:
>>>>
>>>> Of course, it's 1ULP.
>>>> But, according to my understanding, sqrt() is one of those very
>>>> few primitives for which IEEE-754 not just recommends > correct
>>>> rounding, but requires correct rounding.
>>>
>>> This is now https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101 .
>>>
>> Exact rounding rules are easy: All 5 core primitives
>> (FADD/FSUB/FMUL/FDIV/FSQRT) for which it was known back in 1978 that it
>> was relatively easy to always provide exact rounding, are in fact
>> required to do so.
>>
>> Having a square root (for any supported fp size) which does not obey
>> this is in fact a breaking bug.
>>
>
> Also various conversions to/from integer, right?

Most int-to-FP conversions are exact, unless you are converting to float
or fp16, but if you have to round, then you must do so correctly.

float-to-int will by default use the current rounding mode, so usually
nearest-or-even, but due to very heavy use of C int casting most
architectures have special operations to do truncation or arbitrary
directed rounding.

> Probably, also binary scaling, including cases of normal-to-subnormal?

Converting from a normal double which becomes a subnormal float is an
example of those grandfathered "round either before or after discovering
that it is subnormal" exceptions, but the correct choise imho is to
first align the mantissa correctly, then round.
>
> BTW, does it mean that "optional core" primitives, in particular, FMADD and FRSQRT,
> if provided, are allowed to give inexact rounding?

FMAC/FMADD pretty much have to use the same rules as FADD/FMUL, i.e.
exact rounding but I don't remember if it is required. You could
obviously claim it as an extension, not part of the 754 instruction set.

See below:
>
>> It is of course perfectly legal to provide an alternative "fastmath"
>> library which disobeys these rules, but at that point you are not
>> providing IEEE 754 math. I.e. Cray did this at some point, right?

Terje

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

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
> >> 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 for 6.63 vs. 6.62 - I rounded down log[2](1/99) not to
> overstate the accuracy :-)
<
-6.630927

On Wednesday, March 30, 2022 at 1:16:06 AM UTC-5, Terje Mathisen wrote:
> MitchAlsup wrote:
> >> Do you have a numeric example? Any error at all in FDIV is simply a bug,
> > <
> > https://www.cl.cam.ac.uk/~jrh13/slides/gelato-25may05/slides.pdf
> That pdf shows 0.50x accuracy for most of the trancendentals, but not a
> single mention of doing the same on FDIV, rather the opposite:
<
0.50x is not correctly rounded. 0.500... is correctly rounded.
<
I can't find the paper on FDIV and SQRT from Itanium by Markstein and somebody
that I read (Oh so many) years ago.
>
> > Exploiting non-atomic division
> > On Itanium, there’s no atomic division operation, but frcpa returns a
> > reciprocal approximation good to about 8 bits.
> > Techniques largely due to Markstein allow this to be refined to an
> > IEEE-correct quotient just using standard fma operations.
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

On Wednesday, March 30, 2022 at 2:57:08 PM UTC+3, MitchAlsup wrote:
> On Wednesday, March 30, 2022 at 1:16:06 AM UTC-5, Terje Mathisen wrote:
> > MitchAlsup wrote:
> > >> Do you have a numeric example? Any error at all in FDIV is simply a bug,
> > > <
> > > https://www.cl.cam.ac.uk/~jrh13/slides/gelato-25may05/slides.pdf
> > That pdf shows 0.50x accuracy for most of the trancendentals, but not a
> > single mention of doing the same on FDIV, rather the opposite:
> <
> 0.50x is not correctly rounded. 0.500... is correctly rounded.
> <
> I can't find the paper on FDIV and SQRT from Itanium by Markstein and somebody
> that I read (Oh so many) years ago.

Did the book ever existed in electronic form? I mean, outside of HP?

> >
> > > Exploiting non-atomic division
> > > On Itanium, there’s no atomic division operation, but frcpa returns a
> > > reciprocal approximation good to about 8 bits.
> > > Techniques largely due to Markstein allow this to be refined to an
> > > IEEE-correct quotient just using standard fma operations.
> > Terje
> >
> > --
> > - <Terje.Mathisen at tmsw.no>
> > "almost all programming can be viewed as an exercise in caching"

On Wednesday, March 30, 2022 at 4:18:39 PM UTC+3, Michael S wrote:
> On Wednesday, March 30, 2022 at 2:57:08 PM UTC+3, MitchAlsup wrote:
> > On Wednesday, March 30, 2022 at 1:16:06 AM UTC-5, Terje Mathisen wrote:
> > > MitchAlsup wrote:
> > > >> Do you have a numeric example? Any error at all in FDIV is simply a bug,
> > > > <
> > > > https://www.cl.cam.ac.uk/~jrh13/slides/gelato-25may05/slides.pdf
> > > That pdf shows 0.50x accuracy for most of the trancendentals, but not a
> > > single mention of doing the same on FDIV, rather the opposite:
> > <
> > 0.50x is not correctly rounded. 0.500... is correctly rounded.
> > <
> > I can't find the paper on FDIV and SQRT from Itanium by Markstein and somebody
> > that I read (Oh so many) years ago.
> Did the book ever existed in electronic form? I mean, outside of HP?

Or, may be, you mean that?
https://www.researchgate.net/publication/2413714_High_Precision_Division_and_Square_Root

> > >
> > > > Exploiting non-atomic division
> > > > On Itanium, there’s no atomic division operation, but frcpa returns a
> > > > reciprocal approximation good to about 8 bits.
> > > > Techniques largely due to Markstein allow this to be refined to an
> > > > IEEE-correct quotient just using standard fma operations.
> > > Terje
> > >
> > > --
> > > - <Terje.Mathisen at tmsw.no>
> > > "almost all programming can be viewed as an exercise in caching"

MitchAlsup wrote:
> On Wednesday, March 30, 2022 at 1:16:06 AM UTC-5, Terje Mathisen wrote:
>> MitchAlsup wrote:
>>>> Do you have a numeric example? Any error at all in FDIV is simply a bug,
>>> <
>>> https://www.cl.cam.ac.uk/~jrh13/slides/gelato-25may05/slides.pdf
>> That pdf shows 0.50x accuracy for most of the trancendentals, but not a
>> single mention of doing the same on FDIV, rather the opposite:
> <
> 0.50x is not correctly rounded. 0.500... is correctly rounded.

Of course, wasn't that clear from what I wrote?

Trancendentals are OK with "best effort", core ops have to be perfect.

In fact, even 0.500... isn't automatically correct, you have to also
obey the "or_even" part of the rounding rule, but you know that at least
as well as me. :-)

Terje

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

On Wednesday, March 30, 2022 at 11:47:34 AM UTC+3, Thomas Koenig wrote:
> Terje Mathisen <terje.m...@tmsw.no> schrieb:
> > Thomas Koenig wrote:
> >> Michael S <already...@yahoo.com> schrieb:
> >>>
> >>> Of course, it's 1ULP.
> >>> But, according to my understanding, sqrt() is one of those very
> >>> few primitives for which IEEE-754 not just recommends > correct
> >>> rounding, but requires correct rounding.
> >>
> >> This is now https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101 .
> >>
> >
> > Exact rounding rules are easy: All 5 core primitives
> > (FADD/FSUB/FMUL/FDIV/FSQRT) for which it was known back in 1978 that it
> > was relatively easy to always provide exact rounding, are in fact
> > required to do so.
> >
> > Having a square root (for any supported fp size) which does not obey
> > this is in fact a breaking bug.
> It's a wrong-code bug (and I have marked it as such). Compilers and
> libraries are known to have these, even thougn, in an ideal world,
> they should not exist.
>
> It is certainly too late to fix for gcc12 (regression-mode only),
> but gcc13 should be doable.
>
> Hmm... for those who are in the know about IEEE standards (which
> cost money :-(): Does sqrt need to follow rounding modes?

It seems, I found 100 times more serious quadmath bug than mere bad rounding in sqrtq().
On my system fmaq() call appears to change global rounding mode!

Demonstration:

// qsum.c
#include <quadmath.h>

__float128 qsum(__float128 a, __float128 b)
{ return a + b;
} // end of qsum.c

// qadd.c
#include <stdio.h>
#include <quadmath.h>

__float128 qsum(__float128 a, __float128 b);
int main(int argz, char** argv)
{ __float128 x[2];
int idx = 0;
for (int arg_i = 1; arg_i < argz; ++arg_i) {
char* arg = argv[arg_i];
char* endp;
x[idx % 2] = strtoflt128(arg, &endp);
if (arg != endp) {
++idx;
if ((idx % 2)==0) {
char buf[256];
__float128 dummy = 0;
__float128 sum0 = qsum(x[0], x[1]);
dummy = fmaq(x[0], x[0], x[0]);
__float128 sum1 = qsum(x[0], x[1]);
quadmath_snprintf(buf, sizeof(buf), "%-+46.36Qe", x[0]);
printf("%s x[0]\n", buf);
quadmath_snprintf(buf, sizeof(buf), "%-+46.36Qe", x[1]);
printf("%s x[1]\n", buf);
quadmath_snprintf(buf, sizeof(buf), "%-+46.36Qe", sum0);
printf("%s x[0]+x[1].\n", buf);
quadmath_snprintf(buf, sizeof(buf), "%-+46.36Qe", sum1);
printf("%s x[0]+x[1] after fmaq().\n", buf);
quadmath_snprintf(buf, sizeof(buf), "%-+46.36Qe", dummy);
printf("%s dummy\n", buf);
}
}
}
return 0;
} // end of qadd.c

Compile as
gcc -Wall -O2 qsum.c qadd.c -lquadmath

$ ./a 1 -1e-100
+1.000000000000000000000000000000000000e+00 x[0]
-1.000000000000000000000000000000000035e-100 x[1]
+1.000000000000000000000000000000000000e+00 x[0]+x[1].
+9.999999999999999999999999999999999037e-01 x[0]+x[1] after fmaq().
+2.000000000000000000000000000000000000e+00 dummy

MitchAlsup <MitchAlsup@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?

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.

> 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?

Michael S <already5chosen@yahoo.com> schrieb:

> On my system fmaq() call appears to change global rounding mode!

> Compile as
> gcc -Wall -O2 qsum.c qadd.c -lquadmath
>
> $ ./a 1 -1e-100
> +1.000000000000000000000000000000000000e+00 x[0]
> -1.000000000000000000000000000000000035e-100 x[1]
> +1.000000000000000000000000000000000000e+00 x[0]+x[1].
> +9.999999999999999999999999999999999037e-01 x[0]+x[1] after fmaq().
> +2.000000000000000000000000000000000000e+00 dummy

On Linux, I get

$ ./a.out 1 1e-100
+1.000000000000000000000000000000000000e+00 x[0]
+1.000000000000000000000000000000000035e-100 x[1]
+1.000000000000000000000000000000000000e+00 x[0]+x[1].
+1.000000000000000000000000000000000000e+00 x[0]+x[1] after fmaq().
+2.000000000000000000000000000000000000e+00 dummy

so the behavior appears to be platform-dependent. It might be
best to contact the maintainers of the platform you use on Windows
(Cygwin, MinGW or something else?)

On Wednesday, March 30, 2022 at 10:48:53 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On my system fmaq() call appears to change global rounding mode!
> > Compile as
> > gcc -Wall -O2 qsum.c qadd.c -lquadmath
> >
> > $ ./a 1 -1e-100
> > +1.000000000000000000000000000000000000e+00 x[0]
> > -1.000000000000000000000000000000000035e-100 x[1]
> > +1.000000000000000000000000000000000000e+00 x[0]+x[1].
> > +9.999999999999999999999999999999999037e-01 x[0]+x[1] after fmaq().
> > +2.000000000000000000000000000000000000e+00 dummy
> On Linux, I get
>
> $ ./a.out 1 1e-100
> +1.000000000000000000000000000000000000e+00 x[0]
> +1.000000000000000000000000000000000035e-100 x[1]
> +1.000000000000000000000000000000000000e+00 x[0]+x[1].
> +1.000000000000000000000000000000000000e+00 x[0]+x[1] after fmaq().
> +2.000000000000000000000000000000000000e+00 dummy
>

Wrong sign in command line parameter.
Should be './a.out 1 -1e-100'

> so the behavior appears to be platform-dependent. It might be
> best to contact the maintainers of the platform you use on Windows
> (Cygwin, MinGW or something else?)

Mingw64, a.k.a MSYS2.
I was never able to figure out who are maintainers.

On Wednesday, March 30, 2022 at 11:01:26 PM UTC+3, Michael S wrote:
> On Wednesday, March 30, 2022 at 10:48:53 PM UTC+3, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > On my system fmaq() call appears to change global rounding mode!
> > > Compile as
> > > gcc -Wall -O2 qsum.c qadd.c -lquadmath
> > >
> > > $ ./a 1 -1e-100
> > > +1.000000000000000000000000000000000000e+00 x[0]
> > > -1.000000000000000000000000000000000035e-100 x[1]
> > > +1.000000000000000000000000000000000000e+00 x[0]+x[1].
> > > +9.999999999999999999999999999999999037e-01 x[0]+x[1] after fmaq().
> > > +2.000000000000000000000000000000000000e+00 dummy
> > On Linux, I get
> >
> > $ ./a.out 1 1e-100
> > +1.000000000000000000000000000000000000e+00 x[0]
> > +1.000000000000000000000000000000000035e-100 x[1]
> > +1.000000000000000000000000000000000000e+00 x[0]+x[1].
> > +1.000000000000000000000000000000000000e+00 x[0]+x[1] after fmaq().
> > +2.000000000000000000000000000000000000e+00 dummy
> >
> Wrong sign in command line parameter.
> Should be './a.out 1 -1e-100'
> > so the behavior appears to be platform-dependent. It might be
> > best to contact the maintainers of the platform you use on Windows
> > (Cygwin, MinGW or something else?)
> Mingw64, a.k.a MSYS2.
> I was never able to figure out who are maintainers.

Today I spent way to many hours installing Linux with dev environment
on top of WSL on WS2019. Neither Debian nor Ubuntu were able to get
their development packages through corporate crapware (firewall, antvirii
and who knows what else). At the end I sort of succeeded with OpenSUSE.
Sort of, because in order to open basic bach I still need Windows Admin privilege.
But that's off topic. The bottom line is that I was to compile and both quadmath
and mpfr.

Indeed, under Linux I don't see a broken fmaq().
But I see quite horrible precision bugs in sqrtq().
E.g. 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024 delivers result that is off by 32 ULP!

"Michael S" wrote in message
news:7029a173-963d-402b-a184-642120b5e1b8n@googlegroups.com...

> The reason for this is that Chebychev poly of order N is defined as the
> best
> approximation for poly of order N+1. So, when term N+1 in infinite Tailor
> series is much bigger than sum of terms N+2 to Infinity then Cheby is
> good, but
> when term N+1 is not dominant it is less good.

True if best means infinity norm, but there are other choices if your
standard is L1 norm or L2 norm.

On Friday, April 1, 2022 at 9:51:49 AM UTC+3, James Van Buskirk wrote:
> "Michael S" wrote in message
> news:7029a173-963d-402b...@googlegroups.com...
>
> > The reason for this is that Chebychev poly of order N is defined as the
> > best
> > approximation for poly of order N+1. So, when term N+1 in infinite Tailor
> > series is much bigger than sum of terms N+2 to Infinity then Cheby is
> > good, but
> > when term N+1 is not dominant it is less good.
>
> True if best means infinity norm,

I am used to call it "the best approximation in the uniform sense" or minimax.

> but there are other choices if your
> standard is L1 norm or L2 norm.

In context of approximation of special functions to given fixed precision, i.e.
to 0.5 ULP for sqrt() or just a bit over 0.5ULP for trigs, the minimax norm is far
more relevant then the other norms, you mentioned.
For inner steps of rsqrt(), minimax is not the absolutely best norm, but quite
close to the best.

Michael S <already5chosen@yahoo.com> schrieb:

> Indeed, under Linux I don't see a broken fmaq().
> But I see quite horrible precision bugs in sqrtq().
> E.g. 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024 delivers result that is off by 32 ULP!

Can you share your test program?

I have

$ cat xx.c
#include <quadmath.h>
#include <stdio.h>
#define MPFR_WANT_FLOAT128
#include <mpfr.h>

#define NDIG 113

int
main ()
{ volatile __float128 af = 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024Q, bf;
mpfr_t a, b, tst, diff;
mpfr_set_default_prec (NDIG);

mpfr_init (a); mpfr_init(b); mpfr_init(tst);
mpfr_set_float128 (a, af, MPFR_RNDN);
bf = sqrtq (af);
mpfr_sqrt (b, a, MPFR_RNDN);
mpfr_set_float128 (tst, bf, MPFR_RNDN);
mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
printf ("\n");
mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
printf ("\n");
mpfr_out_str (stdout, 10, 0, tst, MPFR_RNDN);
printf ("\n");
return 0;
} $ gcc xx.c -lmpfr -lgmp -lquadmath
$ ./a.out
5.70906089989805317627271688535438725e-309
7.55583277997736869996788472321132647e-155
7.55583277997736869996788472321132647e-155

which looks all right.

On Sunday, April 3, 2022 at 5:36:34 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > Indeed, under Linux I don't see a broken fmaq().
> > But I see quite horrible precision bugs in sqrtq().
> > E.g. 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024 delivers result that is off by 32 ULP!
> Can you share your test program?
>

#include <stdio.h>
#include <stdlib.h>
#include <quadmath.h>
#include <gmp.h>
#include <mpfr.h>

static void to_mpfr(mpfr_t dst, __float128 src)
{ char buf[256];
quadmath_snprintf(buf, sizeof(buf), "%Qa", src);
mpfr_set_str(dst, buf, 0, GMP_RNDN);
}

static __float128 CalcULP(__float128 x, int dir)
{ if (dir < 0)
return x - nextafterq(x, 0.0q);
else
return nextafterq(x, FLT128_MAX) - x;
}

int main(int argz, char** argv)
{ const int MPFR_NBITS = 240;
mpfr_t xx, resx, ref, diff, ULPx, reldiff;
mpfr_init2(xx, MPFR_NBITS);
mpfr_init2(resx, MPFR_NBITS);
mpfr_init2(ref, MPFR_NBITS);
mpfr_init2(diff, MPFR_NBITS);
mpfr_init2(ULPx, MPFR_NBITS);
mpfr_init2(reldiff, MPFR_NBITS);

for (int arg_i = 1; arg_i < argz; ++arg_i) {
char* arg = argv[arg_i];
char* endp;
__float128 x = strtoflt128(arg, &endp);
if (arg != endp) {
char xebuf[256], yebuf[256];
char xhbuf[256], yhbuf[256];
__float128 y = sqrtq(x);
quadmath_snprintf(xebuf, sizeof(xebuf), "%-+50.38Qe", x);
quadmath_snprintf(yebuf, sizeof(yebuf), "%-+50.38Qe", y);

quadmath_snprintf(xhbuf, sizeof(xhbuf), "%-+40.28Qa", x);
quadmath_snprintf(yhbuf, sizeof(yhbuf), "%-+40.28Qa", y);

to_mpfr(xx, x);
to_mpfr(resx, y);
mpfr_sqrt(ref, xx, GMP_RNDN);

__float128 ULP = CalcULP(y, mpfr_cmp(ref, resx));
to_mpfr(ULPx ,ULP);

mpfr_sub(diff, resx, ref, GMP_RNDN);
mpfr_div(reldiff, diff, ULPx, GMP_RNDN);

mpfr_printf(
"arg %s\n"
"x %s %s\n"
"res %s %s\n"
"xx %-+50.38Re %-+50.38Ra\n"
"resx %-+50.38Re %-+50.38Ra\n"
"ref %-+50.38Re %-+50.38Ra\n"
"diff %-+50.38Re %-+50.38Ra\n"
"dULP %-+50.38Rf\n"
"\n"
, arg
, xebuf, xhbuf
, yebuf, yhbuf
, xx, xx
, resx, resx
, ref, ref
, diff, diff
, reldiff
);
}
}
return 0;
}

Compiled as:
gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp

On ws2019+WSL+OpenSUSE:
../a.out 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
res +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de4p-512
xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
resx +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de40000000000p-512
ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
diff +4.63601724942796444562165627161223615263e-187 +0x2.04653571c81b52db61f9a0b94f6857b5c00000p-620
dULP +32.27470917173350947337067829258154575147

On Windows+MSYS2 this particular test case works fine:
./a 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
res +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc4p-512
xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
resx +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc40000000000p-512
ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
diff +3.94598895362939939982164165850203162318e-189 +0x4.653571c81b52db61f9a0b94f6857b5c0000000p-628
dULP +0.27470917173350947337067829258154575147

It's hard to believe that the problem is in WSL.
According to my understanding, it does not affect user mode programs.

cat /etc/os-release
NAME="openSUSE Tumbleweed"
# VERSION="20211027"
ID="opensuse-tumbleweed"
ID_LIKE="opensuse suse"
VERSION_ID="20211027"
PRETTY_NAME="openSUSE Tumbleweed"
ANSI_COLOR="0;32"
CPE_NAME="cpe:/o:opensuse:tumbleweed:20211027"
BUG_REPORT_URL="https://bugs.opensuse.org"
HOME_URL="https://www.opensuse.org/"
DOCUMENTATION_URL="https://en.opensuse.org/Portal:Tumbleweed"
LOGO="distributor-logo-Tumbleweed"

uname -r
4.4.0-17763-Microsoft

> I have
>
> $ cat xx.c
> #include <quadmath.h>
> #include <stdio.h>
> #define MPFR_WANT_FLOAT128
> #include <mpfr.h>
>
> #define NDIG 113
>
> int
> main ()
> {
> volatile __float128 af = 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024Q, bf;
> mpfr_t a, b, tst, diff;
> mpfr_set_default_prec (NDIG);
>
> mpfr_init (a); mpfr_init(b); mpfr_init(tst);
> mpfr_set_float128 (a, af, MPFR_RNDN);
> bf = sqrtq (af);
> mpfr_sqrt (b, a, MPFR_RNDN);
> mpfr_set_float128 (tst, bf, MPFR_RNDN);
> mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
> printf ("\n");
> mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
> printf ("\n");
> mpfr_out_str (stdout, 10, 0, tst, MPFR_RNDN);
> printf ("\n");
> return 0;
> }
> $ gcc xx.c -lmpfr -lgmp -lquadmath
> $ ./a.out
> 5.70906089989805317627271688535438725e-309
> 7.55583277997736869996788472321132647e-155
> 7.55583277997736869996788472321132647e-155
>
> which looks all right.

Indeed, it's o.k.
But on my OpenSUSE I got:
5.70906089989805317627271688535438725e-309
7.55583277997736869996788472321132647e-155
7.55583277997736869996788472321137244e-155

Hardly o.k. :(

On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:
> On Sunday, April 3, 2022 at 5:36:34 PM UTC+3, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > Indeed, under Linux I don't see a broken fmaq().
> > > But I see quite horrible precision bugs in sqrtq().
> > > E.g. 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024 delivers result that is off by 32 ULP!
> > Can you share your test program?
> >
> #include <stdio.h>
> #include <stdlib.h>
> #include <quadmath.h>
> #include <gmp.h>
> #include <mpfr.h>
> static void to_mpfr(mpfr_t dst, __float128 src)
> {
> char buf[256];
> quadmath_snprintf(buf, sizeof(buf), "%Qa", src);
> mpfr_set_str(dst, buf, 0, GMP_RNDN);
> }
> static __float128 CalcULP(__float128 x, int dir)
> {
> if (dir < 0)
> return x - nextafterq(x, 0.0q);
> else
> return nextafterq(x, FLT128_MAX) - x;
> }
>
> int main(int argz, char** argv)
> {
> const int MPFR_NBITS = 240;
> mpfr_t xx, resx, ref, diff, ULPx, reldiff;
> mpfr_init2(xx, MPFR_NBITS);
> mpfr_init2(resx, MPFR_NBITS);
> mpfr_init2(ref, MPFR_NBITS);
> mpfr_init2(diff, MPFR_NBITS);
> mpfr_init2(ULPx, MPFR_NBITS);
> mpfr_init2(reldiff, MPFR_NBITS);
> for (int arg_i = 1; arg_i < argz; ++arg_i) {
> char* arg = argv[arg_i];
> char* endp;
> __float128 x = strtoflt128(arg, &endp);
> if (arg != endp) {
> char xebuf[256], yebuf[256];
> char xhbuf[256], yhbuf[256];
> __float128 y = sqrtq(x);
> quadmath_snprintf(xebuf, sizeof(xebuf), "%-+50.38Qe", x);
> quadmath_snprintf(yebuf, sizeof(yebuf), "%-+50.38Qe", y);
>
> quadmath_snprintf(xhbuf, sizeof(xhbuf), "%-+40.28Qa", x);
> quadmath_snprintf(yhbuf, sizeof(yhbuf), "%-+40.28Qa", y);
>
> to_mpfr(xx, x);
> to_mpfr(resx, y);
> mpfr_sqrt(ref, xx, GMP_RNDN);
>
> __float128 ULP = CalcULP(y, mpfr_cmp(ref, resx));
> to_mpfr(ULPx ,ULP);
>
> mpfr_sub(diff, resx, ref, GMP_RNDN);
> mpfr_div(reldiff, diff, ULPx, GMP_RNDN);
>
> mpfr_printf(
> "arg %s\n"
> "x %s %s\n"
> "res %s %s\n"
> "xx %-+50.38Re %-+50.38Ra\n"
> "resx %-+50.38Re %-+50.38Ra\n"
> "ref %-+50.38Re %-+50.38Ra\n"
> "diff %-+50.38Re %-+50.38Ra\n"
> "dULP %-+50.38Rf\n"
> "\n"
> , arg
> , xebuf, xhbuf
> , yebuf, yhbuf
> , xx, xx
> , resx, resx
> , ref, ref
> , diff, diff
> , reldiff
> );
> }
> }
> return 0;
> }
>
>
> Compiled as:
> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp
>
> On ws2019+WSL+OpenSUSE:
> ./a.out 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> res +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de4p-512
> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
> resx +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de40000000000p-512
> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
> diff +4.63601724942796444562165627161223615263e-187 +0x2.04653571c81b52db61f9a0b94f6857b5c00000p-620
> dULP +32.27470917173350947337067829258154575147
>
> On Windows+MSYS2 this particular test case works fine:
> ./a 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> res +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc4p-512
> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
> resx +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc40000000000p-512
> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
> diff +3.94598895362939939982164165850203162318e-189 +0x4.653571c81b52db61f9a0b94f6857b5c0000000p-628
> dULP +0.27470917173350947337067829258154575147
>
> It's hard to believe that the problem is in WSL.
> According to my understanding, it does not affect user mode programs.
>
> cat /etc/os-release
> NAME="openSUSE Tumbleweed"
> # VERSION="20211027"
> ID="opensuse-tumbleweed"
> ID_LIKE="opensuse suse"
> VERSION_ID="20211027"
> PRETTY_NAME="openSUSE Tumbleweed"
> ANSI_COLOR="0;32"
> CPE_NAME="cpe:/o:opensuse:tumbleweed:20211027"
> BUG_REPORT_URL="https://bugs.opensuse.org"
> HOME_URL="https://www.opensuse.org/"
> DOCUMENTATION_URL="https://en.opensuse.org/Portal:Tumbleweed"
> LOGO="distributor-logo-Tumbleweed"
>
> uname -r
> 4.4.0-17763-Microsoft
> > I have
> >
> > $ cat xx.c
> > #include <quadmath.h>
> > #include <stdio.h>
> > #define MPFR_WANT_FLOAT128
> > #include <mpfr.h>
> >
> > #define NDIG 113
> >
> > int
> > main ()
> > {
> > volatile __float128 af = 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024Q, bf;
> > mpfr_t a, b, tst, diff;
> > mpfr_set_default_prec (NDIG);
> >
> > mpfr_init (a); mpfr_init(b); mpfr_init(tst);
> > mpfr_set_float128 (a, af, MPFR_RNDN);
> > bf = sqrtq (af);
> > mpfr_sqrt (b, a, MPFR_RNDN);
> > mpfr_set_float128 (tst, bf, MPFR_RNDN);
> > mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
> > printf ("\n");
> > mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
> > printf ("\n");
> > mpfr_out_str (stdout, 10, 0, tst, MPFR_RNDN);
> > printf ("\n");
> > return 0;
> > }
> > $ gcc xx.c -lmpfr -lgmp -lquadmath
> > $ ./a.out
> > 5.70906089989805317627271688535438725e-309
> > 7.55583277997736869996788472321132647e-155
> > 7.55583277997736869996788472321132647e-155
> >
> > which looks all right.
> Indeed, it's o.k.
> But on my OpenSUSE I got:
> 5.70906089989805317627271688535438725e-309
> 7.55583277997736869996788472321132647e-155
> 7.55583277997736869996788472321137244e-155
>
> Hardly o.k. :(

I managed to install on the same server Fedora Remix for WSL (OT, unlike other distros, installing this one was painless).
It has exactly the same problem as OpenSUSE.
So, may be, the suggestion that bug is somehow related to WSL Linux kernel now appears to me less crazy than 30 minutes ago.

Michael S <already5chosen@yahoo.com> schrieb:
> On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:

>> Compiled as:
>> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp

Hmm.. what version of gcc are you using? Does "ldd a.out" point
to any strange places for libquadmath?

What you could try (clutching at straws here for a bit) is to link
the static version of the libquadmath library, i.e.

gcc -Wall -O2 bug_rep.c -lmpfr -lgmp /usr/lib/gcc/x86_64-linux-gnu/9/libquadmath.a

(or where you many find a static libquadmath).

>> On ws2019+WSL+OpenSUSE:
>> ./a.out 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> res +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de4p-512
>> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
>> resx +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de40000000000p-512
>> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
>> diff +4.63601724942796444562165627161223615263e-187 +0x2.04653571c81b52db61f9a0b94f6857b5c00000p-620
>> dULP +32.27470917173350947337067829258154575147
>>
>> On Windows+MSYS2 this particular test case works fine:
>> ./a 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
>> res +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc4p-512
>> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
>> resx +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc40000000000p-512
>> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
>> diff +3.94598895362939939982164165850203162318e-189 +0x4.653571c81b52db61f9a0b94f6857b5c0000000p-628
>> dULP +0.27470917173350947337067829258154575147
>>
>> It's hard to believe that the problem is in WSL.
>> According to my understanding, it does not affect user mode programs.

It should not.

>> cat /etc/os-release
>> NAME="openSUSE Tumbleweed"
>> # VERSION="20211027"
>> ID="opensuse-tumbleweed"
>> ID_LIKE="opensuse suse"
>> VERSION_ID="20211027"
>> PRETTY_NAME="openSUSE Tumbleweed"
>> ANSI_COLOR="0;32"
>> CPE_NAME="cpe:/o:opensuse:tumbleweed:20211027"
>> BUG_REPORT_URL="https://bugs.opensuse.org"
>> HOME_URL="https://www.opensuse.org/"
>> DOCUMENTATION_URL="https://en.opensuse.org/Portal:Tumbleweed"
>> LOGO="distributor-logo-Tumbleweed"

Tumbleweed is of course a bit experimental but it should certainly
not be this badly broken.

>> uname -r
>> 4.4.0-17763-Microsoft
>> > I have
>> >
>> > $ cat xx.c
>> > #include <quadmath.h>
>> > #include <stdio.h>
>> > #define MPFR_WANT_FLOAT128
>> > #include <mpfr.h>
>> >
>> > #define NDIG 113
>> >
>> > int
>> > main ()
>> > {
>> > volatile __float128 af = 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024Q, bf;
>> > mpfr_t a, b, tst, diff;
>> > mpfr_set_default_prec (NDIG);
>> >
>> > mpfr_init (a); mpfr_init(b); mpfr_init(tst);
>> > mpfr_set_float128 (a, af, MPFR_RNDN);
>> > bf = sqrtq (af);
>> > mpfr_sqrt (b, a, MPFR_RNDN);
>> > mpfr_set_float128 (tst, bf, MPFR_RNDN);
>> > mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
>> > printf ("\n");
>> > mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
>> > printf ("\n");
>> > mpfr_out_str (stdout, 10, 0, tst, MPFR_RNDN);
>> > printf ("\n");
>> > return 0;
>> > }
>> > $ gcc xx.c -lmpfr -lgmp -lquadmath
>> > $ ./a.out
>> > 5.70906089989805317627271688535438725e-309
>> > 7.55583277997736869996788472321132647e-155
>> > 7.55583277997736869996788472321132647e-155
>> >
>> > which looks all right.
>> Indeed, it's o.k.
>> But on my OpenSUSE I got:
>> 5.70906089989805317627271688535438725e-309
>> 7.55583277997736869996788472321132647e-155
>> 7.55583277997736869996788472321137244e-155
>>
>> Hardly o.k. :(

That is indeed really strange.

> I managed to install on the same server Fedora Remix for WSL
> (OT, unlike other distros, installing this one was painless).
> It has exactly the same problem as OpenSUSE.

> So, may be, the suggestion that bug is somehow related to WSL
> Linux kernel now appears to me less crazy than 30 minutes ago.

That would be _really_ weird.

Maybe you can submit a bug report against Fedora or OpenSUSE,
they are at least likely to not throw it out of the window
right away.

On Sunday, April 3, 2022 at 9:29:14 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:
>
> >> Compiled as:
> >> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp
> Hmm.. what version of gcc are you using? Does "ldd a.out" point
> to any strange places for libquadmath?
>
> What you could try (clutching at straws here for a bit) is to link
> the static version of the libquadmath library, i.e.
>
> gcc -Wall -O2 bug_rep.c -lmpfr -lgmp /usr/lib/gcc/x86_64-linux-gnu/9/libquadmath.a
>
> (or where you many find a static libquadmath).

-lmpfr -lgmp /usr/lib/gcc/x86_64-redhat-linux/11/libquadmath.a /usr/lib64/libm.a

Same shite.

> >> On ws2019+WSL+OpenSUSE:
> >> ./a.out 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> res +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de4p-512
> >> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
> >> resx +7.55583277997736869996788472321137243638e-155 +0x1.0358a8295b66c7fefb50c1298de40000000000p-512
> >> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
> >> diff +4.63601724942796444562165627161223615263e-187 +0x2.04653571c81b52db61f9a0b94f6857b5c00000p-620
> >> dULP +32.27470917173350947337067829258154575147
> >>
> >> On Windows+MSYS2 this particular test case works fine:
> >> ./a 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> arg 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> x +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0eap-1024
> >> res +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc4p-512
> >> xx +5.70906089989805317627271688535438724887e-309 +0x1.06bc82f7b9d71dfcbddf2358a0ea0000000000p-1024
> >> resx +7.55583277997736869996788472321132647080e-155 +0x1.0358a8295b66c7fefb50c1298dc40000000000p-512
> >> ref +7.55583277997736869996788472321132607621e-155 +0x1.0358a8295b66c7fefb50c1298dc3b9aca8e37ep-512
> >> diff +3.94598895362939939982164165850203162318e-189 +0x4.653571c81b52db61f9a0b94f6857b5c0000000p-628
> >> dULP +0.27470917173350947337067829258154575147
> >>
> >> It's hard to believe that the problem is in WSL.
> >> According to my understanding, it does not affect user mode programs.
> It should not.
> >> cat /etc/os-release
> >> NAME="openSUSE Tumbleweed"
> >> # VERSION="20211027"
> >> ID="opensuse-tumbleweed"
> >> ID_LIKE="opensuse suse"
> >> VERSION_ID="20211027"
> >> PRETTY_NAME="openSUSE Tumbleweed"
> >> ANSI_COLOR="0;32"
> >> CPE_NAME="cpe:/o:opensuse:tumbleweed:20211027"
> >> BUG_REPORT_URL="https://bugs.opensuse.org"
> >> HOME_URL="https://www.opensuse.org/"
> >> DOCUMENTATION_URL="https://en.opensuse.org/Portal:Tumbleweed"
> >> LOGO="distributor-logo-Tumbleweed"
> Tumbleweed is of course a bit experimental but it should certainly
> not be this badly broken.
> >> uname -r
> >> 4.4.0-17763-Microsoft
> >> > I have
> >> >
> >> > $ cat xx.c
> >> > #include <quadmath.h>
> >> > #include <stdio.h>
> >> > #define MPFR_WANT_FLOAT128
> >> > #include <mpfr.h>
> >> >
> >> > #define NDIG 113
> >> >
> >> > int
> >> > main ()
> >> > {
> >> > volatile __float128 af = 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024Q, bf;
> >> > mpfr_t a, b, tst, diff;
> >> > mpfr_set_default_prec (NDIG);
> >> >
> >> > mpfr_init (a); mpfr_init(b); mpfr_init(tst);
> >> > mpfr_set_float128 (a, af, MPFR_RNDN);
> >> > bf = sqrtq (af);
> >> > mpfr_sqrt (b, a, MPFR_RNDN);
> >> > mpfr_set_float128 (tst, bf, MPFR_RNDN);
> >> > mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
> >> > printf ("\n");
> >> > mpfr_out_str (stdout, 10, 0, b, MPFR_RNDN);
> >> > printf ("\n");
> >> > mpfr_out_str (stdout, 10, 0, tst, MPFR_RNDN);
> >> > printf ("\n");
> >> > return 0;
> >> > }
> >> > $ gcc xx.c -lmpfr -lgmp -lquadmath
> >> > $ ./a.out
> >> > 5.70906089989805317627271688535438725e-309
> >> > 7.55583277997736869996788472321132647e-155
> >> > 7.55583277997736869996788472321132647e-155
> >> >
> >> > which looks all right.
> >> Indeed, it's o.k.
> >> But on my OpenSUSE I got:
> >> 5.70906089989805317627271688535438725e-309
> >> 7.55583277997736869996788472321132647e-155
> >> 7.55583277997736869996788472321137244e-155
> >>
> >> Hardly o.k. :(
> That is indeed really strange.
> > I managed to install on the same server Fedora Remix for WSL
> > (OT, unlike other distros, installing this one was painless).
> > It has exactly the same problem as OpenSUSE.
>
> > So, may be, the suggestion that bug is somehow related to WSL
> > Linux kernel now appears to me less crazy than 30 minutes ago.
> That would be _really_ weird.
>
> Maybe you can submit a bug report against Fedora or OpenSUSE,
> they are at least likely to not throw it out of the window
> right away.

I think, throwing it out of the window is exactly what they are going to do.
Unlike WSL on later versions of Win10 and on Win11, WSL on WS2019 is
a ultimate red-headed step child. It is not supported by Microsoft themselves
so hard to expect that Red Hat (IBM) or SUSE care about it.

Is there any use for approximately equals? And how to go about determining
approximate equality. 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.

On Monday, April 4, 2022 at 12:10:40 AM UTC+3, Michael S wrote:
> On Sunday, April 3, 2022 at 9:29:14 PM UTC+3, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:
> >
> > >> Compiled as:
> > >> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp
> > Hmm.. what version of gcc are you using? Does "ldd a.out" point
> > to any strange places for libquadmath?
> >
> > What you could try (clutching at straws here for a bit) is to link
> > the static version of the libquadmath library, i.e.
> >
> > gcc -Wall -O2 bug_rep.c -lmpfr -lgmp /usr/lib/gcc/x86_64-linux-gnu/9/libquadmath.a
> >
> > (or where you many find a static libquadmath).
> -lmpfr -lgmp /usr/lib/gcc/x86_64-redhat-linux/11/libquadmath.a /usr/lib64/libm.a
>
> Same shite.

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.

On Monday, April 4, 2022 at 1:22:36 AM UTC+3, Michael S wrote:
> On Monday, April 4, 2022 at 12:10:40 AM UTC+3, Michael S wrote:
> > On Sunday, April 3, 2022 at 9:29:14 PM UTC+3, Thomas Koenig wrote:
> > > Michael S <already...@yahoo.com> schrieb:
> > > > On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:
> > >
> > > >> Compiled as:
> > > >> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp
> > > Hmm.. what version of gcc are you using? Does "ldd a.out" point
> > > to any strange places for libquadmath?
> > >
> > > What you could try (clutching at straws here for a bit) is to link
> > > the static version of the libquadmath library, i.e.
> > >
> > > gcc -Wall -O2 bug_rep.c -lmpfr -lgmp /usr/lib/gcc/x86_64-linux-gnu/9/libquadmath.a
> > >
> > > (or where you many find a static libquadmath).
> > -lmpfr -lgmp /usr/lib/gcc/x86_64-redhat-linux/11/libquadmath.a /usr/lib64/libm.a
> >
> > Same shite.
> 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.

I think, it would be fare to say that the whole compiler environment is broken.
Having 'long double' with 15-bit exponent and 53-bit mantissa is, by itself, legal, if weird.
But then LDBL_MANT_DIG in float.h should be set to 53. An in this case it is set to 64.

On Monday, April 4, 2022 at 2:09:06 AM UTC+3, Michael S wrote:
> On Monday, April 4, 2022 at 1:22:36 AM UTC+3, Michael S wrote:
> > On Monday, April 4, 2022 at 12:10:40 AM UTC+3, Michael S wrote:
> > > On Sunday, April 3, 2022 at 9:29:14 PM UTC+3, Thomas Koenig wrote:
> > > > Michael S <already...@yahoo.com> schrieb:
> > > > > On Sunday, April 3, 2022 at 8:09:02 PM UTC+3, Michael S wrote:
> > > >
> > > > >> Compiled as:
> > > > >> gcc -Wall -O2 bug_rep.c -lquadmath -lmpfr -lgmp
> > > > Hmm.. what version of gcc are you using? Does "ldd a.out" point
> > > > to any strange places for libquadmath?
> > > >
> > > > What you could try (clutching at straws here for a bit) is to link
> > > > the static version of the libquadmath library, i.e.
> > > >
> > > > gcc -Wall -O2 bug_rep.c -lmpfr -lgmp /usr/lib/gcc/x86_64-linux-gnu/9/libquadmath.a
> > > >
> > > > (or where you many find a static libquadmath).
> > > -lmpfr -lgmp /usr/lib/gcc/x86_64-redhat-linux/11/libquadmath.a /usr/lib64/libm.a
> > >
> > > Same shite.
> > 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.
> I think, it would be fare to say that the whole compiler environment is broken.
> Having 'long double' with 15-bit exponent and 53-bit mantissa is, by itself, legal, if weird.
> But then LDBL_MANT_DIG in float.h should be set to 53. An in this case it is set to 64.

Just to affirm what should be already obvious to attentive readers:
If I am adding following lines at the beginning of the program then Linux under WSL
under WS2019 behaves identically to x86-64 Linux tested by Thomas Koenig.

unsigned short cw;
_FPU_GETCW(cw);
cw = cw | _FPU_EXTENDED;
_FPU_SETCW(cw);