On Monday, March 28, 2022 at 1:14:52 AM UTC-5, James Van Buskirk wrote:
> "MitchAlsup" wrote in message
> news:600fb4d3-5e5f-490f...@googlegroups.com...
> > First 3 iterations are SP, next iteration is DP, next 2 are QP.
> Have you investigated the potential of using variable numbers of
> terms from the series
>
> x/(1-(1-x*D)) = x*sum([((1-x*D)**n,n=0,∞)])
> x/sqrt(1-(1-x**2*D)) =
> x*sum([(gamma(n+0.5)/(sqrt(pi)*gamma(n+1.0))*(1-x**2*D)**n,n=0,∞)])
<
Cursarilly: I develop HW algorithms, these are not constrained like SW algorithms
are. In HW, one often uses an algorithm that would be less optimal in SW imple-
mentation, but HW run faster--for example: Goldschmidt is "like" Newton-Raphson
but the multiplies are independent rather than dependent, so if you can keep the
multiplier tree occupied, it will be faster than N-R. SW is always bound by the
latency of the multiplier function unit, HW only by the multiplier tree itself. So,
Goldschmidt runs 1 iteration every 2 cycles, N-R runs 1 iteration every 2×multiplier
latency (generally 4 or 5). About 4× faster per iteration. HW has no trouble converting
from one precision to another, so 3 iterations in SP, 1 in DP, and 2 in QP is straight
forward (no conversion instructions.)
>
> to perhaps remove an iteration somewhere?
<
HW removes cycles directly, SW attempts to use tricks to remove 1 iteration..

On Friday, March 25, 2022 at 4:22:59 AM UTC-6, Elijah Stone wrote:
> The itanium had only an fp reciprocal; no division, that had to be done in
> software.

Yes, and that was a mistake when the Itanium did it, and it would be
a mistake for any other processor intended for general-purpose use.

Instructions take time to decode and fetch and so on. So if you
have no divide instruction, there's no way that you're going to be
able to do start one division going in the pipeline every single cycle.
Which is what a decently fast processor ought to be able to do, at
least as a peak rate.

John Savard

On Friday, March 25, 2022 at 4:32:10 PM UTC-6, MitchAlsup quoted, in part:
> That is basically why we don't do
> CORDIC anymore.

Hey, CORDIC was a _great_ algorithm. For computers that didn't
have a _multiply_ instruction in hardware. So of course it's not
useful on computers that can even do division and floating-point
arithmetic in hardware.

Pocket calculators, with four-bit processors that indeed don't
have hardware multiply, still do CORDIC.

John Savard

On Monday, March 28, 2022 at 9:02:37 PM UTC+3, Quadibloc wrote:
> On Friday, March 25, 2022 at 4:22:59 AM UTC-6, Elijah Stone wrote:
> > The itanium had only an fp reciprocal; no division, that had to be done in
> > software.
>
> Yes, and that was a mistake when the Itanium did it, and it would be
> a mistake for any other processor intended for general-purpose use.
>
> Instructions take time to decode and fetch and so on. So if you
> have no divide instruction, there's no way that you're going to be
> able to do start one division going in the pipeline every single cycle.
> Which is what a decently fast processor ought to be able to do, at
> least as a peak rate.
>
> John Savard

According to this definition, decently fast processors still do not exist.

On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> > Michael S <already...@yahoo.com> schrieb:
> > > On my system(s) I see another very strange timing effect:
> > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> > Source for the sqrt routine in libquadmath is here:
> >
> > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> >
> > Hmm.. they first do a test if it is within double precision range,
> > with two Newton iterations, then a test if it is within long double
> > range with a single Newton iteration, and if it is outside then
> > they pick apart the number and run two Newton iterations.
> >
> > Were the numbers outside the range of double inside the range
> > of long double?
> Yes.
> I didn't test "outside long double".
> If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> > That could explain things (and suggest
> > an improvement).
>
> I think, the whole routine needs rewrite, rather than improvements.
> The original author probably was not thinking very hard.

Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
For example,
sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
The correct answer is 1.03168097468685393293622695023074831927

I wonder if on Linux result is also wrong.

On Monday, March 28, 2022 at 1:02:37 PM UTC-5, Quadibloc wrote:
> On Friday, March 25, 2022 at 4:22:59 AM UTC-6, Elijah Stone wrote:
> > The itanium had only an fp reciprocal; no division, that had to be done in
> > software.
>
> Yes, and that was a mistake when the Itanium did it, and it would be
> a mistake for any other processor intended for general-purpose use.
>
> Instructions take time to decode and fetch and so on. So if you
> have no divide instruction, there's no way that you're going to be
> able to do start one division going in the pipeline every single cycle.
> Which is what a decently fast processor ought to be able to do, at
> least as a peak rate.
<
Err, no...........
<
Division is not used often enough to warrant 1-cycle throughput.
Division has 15-20 cycle latency (64-bit, int or fp)
<
The cost to make reciprocation fully pipelined is about the same cost
as 6-7 FAMC units. The cost to make division fully pipelined is about
the same cost as 10-11 FMAC units.
>
> John Savard

On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> > On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> > > Michael S <already...@yahoo.com> schrieb:
> > > > On my system(s) I see another very strange timing effect:
> > > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> > > Source for the sqrt routine in libquadmath is here:
> > >
> > > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> > >
> > > Hmm.. they first do a test if it is within double precision range,
> > > with two Newton iterations, then a test if it is within long double
> > > range with a single Newton iteration, and if it is outside then
> > > they pick apart the number and run two Newton iterations.
> > >
> > > Were the numbers outside the range of double inside the range
> > > of long double?
> > Yes.
> > I didn't test "outside long double".
> > If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> > > That could explain things (and suggest
> > > an improvement).
> >
> > I think, the whole routine needs rewrite, rather than improvements.
> > The original author probably was not thinking very hard.
> Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> For example,
> sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> The correct answer is 1.03168097468685393293622695023074831927
<
This is a 1ULP error, perhaps simply improperly rounded (or the 1:million cases that require the second N-R
iteration after you have 111 accurate fraction bits.
<
BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).
>
> I wonder if on Linux result is also wrong.

On Tuesday, March 29, 2022 at 12:44:55 AM UTC+3, MitchAlsup wrote:
> On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> > On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> > > On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> > > > Michael S <already...@yahoo.com> schrieb:
> > > > > On my system(s) I see another very strange timing effect:
> > > > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > > > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> > > > Source for the sqrt routine in libquadmath is here:
> > > >
> > > > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> > > >
> > > > Hmm.. they first do a test if it is within double precision range,
> > > > with two Newton iterations, then a test if it is within long double
> > > > range with a single Newton iteration, and if it is outside then
> > > > they pick apart the number and run two Newton iterations.
> > > >
> > > > Were the numbers outside the range of double inside the range
> > > > of long double?
> > > Yes.
> > > I didn't test "outside long double".
> > > If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> > > > That could explain things (and suggest
> > > > an improvement).
> > >
> > > I think, the whole routine needs rewrite, rather than improvements.
> > > The original author probably was not thinking very hard.
> > Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> > For example,
> > sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> > The correct answer is 1.03168097468685393293622695023074831927
> <
> This is a 1ULP error, perhaps simply improperly rounded

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.

> (or the 1:million cases that require the second N-R
> iteration after you have 111 accurate fraction bits.

Unlikely.

> <
> BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).

Didn't they provide compiler option for exact division?

> >
> > I wonder if on Linux result is also wrong.

On Tuesday, March 29, 2022 at 1:49:45 AM UTC+3, Michael S wrote:
> On Tuesday, March 29, 2022 at 12:44:55 AM UTC+3, MitchAlsup wrote:
> > On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> > > On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> > > > On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> > > > > Michael S <already...@yahoo.com> schrieb:
> > > > > > On my system(s) I see another very strange timing effect:
> > > > > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > > > > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> > > > > Source for the sqrt routine in libquadmath is here:
> > > > >
> > > > > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> > > > >
> > > > > Hmm.. they first do a test if it is within double precision range,
> > > > > with two Newton iterations, then a test if it is within long double
> > > > > range with a single Newton iteration, and if it is outside then
> > > > > they pick apart the number and run two Newton iterations.
> > > > >
> > > > > Were the numbers outside the range of double inside the range
> > > > > of long double?
> > > > Yes.
> > > > I didn't test "outside long double".
> > > > If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> > > > > That could explain things (and suggest
> > > > > an improvement).
> > > >
> > > > I think, the whole routine needs rewrite, rather than improvements.
> > > > The original author probably was not thinking very hard.
> > > Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> > > For example,
> > > sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> > > The correct answer is 1.03168097468685393293622695023074831927
> > <
> > This is a 1ULP error, perhaps simply improperly rounded
> 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.
> > (or the 1:million cases that require the second N-R
> > iteration after you have 111 accurate fraction bits.
> Unlikely.

I checked it. An error in this case = 0.7499966 ULP.

> > <
> > BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).
> Didn't they provide compiler option for exact division?
> > >
> > > I wonder if on Linux result is also wrong.

On Monday, March 28, 2022 at 3:39:51 PM UTC-6, MitchAlsup wrote:
> On Monday, March 28, 2022 at 1:02:37 PM UTC-5, Quadibloc wrote:
> > On Friday, March 25, 2022 at 4:22:59 AM UTC-6, Elijah Stone wrote:
> > > The itanium had only an fp reciprocal; no division, that had to be done in
> > > software.
> >
> > Yes, and that was a mistake when the Itanium did it, and it would be
> > a mistake for any other processor intended for general-purpose use.
> >
> > Instructions take time to decode and fetch and so on. So if you
> > have no divide instruction, there's no way that you're going to be
> > able to do start one division going in the pipeline every single cycle.
> > Which is what a decently fast processor ought to be able to do, at
> > least as a peak rate.
> <
> Err, no...........
> <
> Division is not used often enough to warrant 1-cycle throughput.
> Division has 15-20 cycle latency (64-bit, int or fp)
> <
> The cost to make reciprocation fully pipelined is about the same cost
> as 6-7 FAMC units. The cost to make division fully pipelined is about
> the same cost as 10-11 FMAC units.

There are definitely some hidden assumptions in my statement.

Obviously, letting a CPU start a division every cycle is going
to use die area. If die area is a crippling constraint, such that
each feature has to be weighed against, 'would I be better in
leaving out this feature, so I could put more cores on the die',
then indeed one division per cycle is an extravagance.

I, on the other hand, am of this opinion:

Once you can put *one* core on a die, that's all you need.

If you want 16 cores, you can put sixteen sockets on your
motherboard, or 16 dies on a substrate... or four sockets
on your motherboard, and four dies on each substrate.

You can't get away with having less than _one_
core, and adequate L1 cache, on a die. That's it.

L2 cache is nice also, however.

Having multiple cores more tightly coupled by putting
them on the same die and sharing a cache produces,
I had thought, so little performance benefit that it's
far outweighed by making bigger cores capable of
doing things faster.

And if Intel and AMD stand up to Microsoft, and say that
it needs to change its licensing policy, so you
can run Windows Home on motherboards with
four or eight CPU chips on them, not just one, I think that's
overdue, and a simple way to get greater performance
without requiring Moore's Law to override the laws
of physics.

John Savard

MitchAlsup wrote:
> On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
>> On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
>>> On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
>>>> Michael S <already...@yahoo.com> schrieb:
>>>>> On my system(s) I see another very strange timing effect:
>>>>> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
>>>>> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
>>>> Source for the sqrt routine in libquadmath is here:
>>>>
>>>> https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
>>>>
>>>> Hmm.. they first do a test if it is within double precision range,
>>>> with two Newton iterations, then a test if it is within long double
>>>> range with a single Newton iteration, and if it is outside then
>>>> they pick apart the number and run two Newton iterations.
>>>>
>>>> Were the numbers outside the range of double inside the range
>>>> of long double?
>>> Yes.
>>> I didn't test "outside long double".
>>> If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
>>>> That could explain things (and suggest
>>>> an improvement).
>>>
>>> I think, the whole routine needs rewrite, rather than improvements.
>>> The original author probably was not thinking very hard.
>> Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
>> For example,
>> sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
>> The correct answer is 1.03168097468685393293622695023074831927
> <
> This is a 1ULP error, perhaps simply improperly rounded (or the 1:million cases that require the second N-R
> iteration after you have 111 accurate fraction bits.
> <
> BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).

Do you have a numeric example? Any error at all in FDIV is simply a bug,
except for subnormal results where rounding errors got grandfathered in.
(I.e. it is explicitly allowed to do the rounding before you realize
that the result is subnormal, so that you have to denormalize it,
leading to double rounding (or maybe even truncation?)

According to the intentions of the full standard, there can be no doubt
that you should first align the result properly (usually normalize it,
but for subnormal you use the temporary 1 bit to stop the normalization
process at the correct amount of shifting), then you examine the
SIGN/LSB/GUARD/STICKY bits to determine the rounding.

Terje

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

On Tuesday, March 29, 2022 at 7:57:35 AM UTC-5, Terje Mathisen wrote:
> MitchAlsup wrote:
> > On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> >> On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> >>> On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> >>>> Michael S <already...@yahoo.com> schrieb:
> >>>>> On my system(s) I see another very strange timing effect:
> >>>>> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> >>>>> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> >>>> Source for the sqrt routine in libquadmath is here:
> >>>>
> >>>> https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> >>>>
> >>>> Hmm.. they first do a test if it is within double precision range,
> >>>> with two Newton iterations, then a test if it is within long double
> >>>> range with a single Newton iteration, and if it is outside then
> >>>> they pick apart the number and run two Newton iterations.
> >>>>
> >>>> Were the numbers outside the range of double inside the range
> >>>> of long double?
> >>> Yes.
> >>> I didn't test "outside long double".
> >>> If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> >>>> That could explain things (and suggest
> >>>> an improvement).
> >>>
> >>> I think, the whole routine needs rewrite, rather than improvements.
> >>> The original author probably was not thinking very hard.
> >> Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> >> For example,
> >> sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> >> The correct answer is 1.03168097468685393293622695023074831927
> > <
> > This is a 1ULP error, perhaps simply improperly rounded (or the 1:million cases that require the second N-R
> > iteration after you have 111 accurate fraction bits.
> > <
> > BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).
> 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
<
> except for subnormal results where rounding errors got grandfathered in.
> (I.e. it is explicitly allowed to do the rounding before you realize
> that the result is subnormal, so that you have to denormalize it,
> leading to double rounding (or maybe even truncation?)
>
> According to the intentions of the full standard, there can be no doubt
> that you should first align the result properly (usually normalize it,
> but for subnormal you use the temporary 1 bit to stop the normalization
> process at the correct amount of shifting), then you examine the
> SIGN/LSB/GUARD/STICKY bits to determine the rounding.
> Terje
>
> --
> - <Terje.Mathisen at tmsw.no>
> "almost all programming can be viewed as an exercise in caching"

On Tuesday, March 29, 2022 at 7:20:32 PM UTC+3, MitchAlsup wrote:
> On Tuesday, March 29, 2022 at 7:57:35 AM UTC-5, Terje Mathisen wrote:
> > MitchAlsup wrote:
> > > On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> > >> On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> > >>> On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> > >>>> Michael S <already...@yahoo.com> schrieb:
> > >>>>> On my system(s) I see another very strange timing effect:
> > >>>>> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> > >>>>> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> > >>>> Source for the sqrt routine in libquadmath is here:
> > >>>>
> > >>>> https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> > >>>>
> > >>>> Hmm.. they first do a test if it is within double precision range,
> > >>>> with two Newton iterations, then a test if it is within long double
> > >>>> range with a single Newton iteration, and if it is outside then
> > >>>> they pick apart the number and run two Newton iterations.
> > >>>>
> > >>>> Were the numbers outside the range of double inside the range
> > >>>> of long double?
> > >>> Yes.
> > >>> I didn't test "outside long double".
> > >>> If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> > >>>> That could explain things (and suggest
> > >>>> an improvement).
> > >>>
> > >>> I think, the whole routine needs rewrite, rather than improvements.
> > >>> The original author probably was not thinking very hard.
> > >> Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> > >> For example,
> > >> sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> > >> The correct answer is 1.03168097468685393293622695023074831927
> > > <
> > > This is a 1ULP error, perhaps simply improperly rounded (or the 1:million cases that require the second N-R
> > > iteration after you have 111 accurate fraction bits.
> > > <
> > > BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).
> > 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

On which page?
So far I don't see anything suggesting that compiler does non-perfectly rounded division.

> <
> > except for subnormal results where rounding errors got grandfathered in.
> > (I.e. it is explicitly allowed to do the rounding before you realize
> > that the result is subnormal, so that you have to denormalize it,
> > leading to double rounding (or maybe even truncation?)
> >
> > According to the intentions of the full standard, there can be no doubt
> > that you should first align the result properly (usually normalize it,
> > but for subnormal you use the temporary 1 bit to stop the normalization
> > process at the correct amount of shifting), then you examine the
> > SIGN/LSB/GUARD/STICKY bits to determine the rounding.
> > Terje
> >
> > --
> > - <Terje.Mathisen at tmsw.no>
> > "almost all programming can be viewed as an exercise in caching"

MitchAlsup <MitchAlsup@aol.com> schrieb:
> On Monday, March 28, 2022 at 3:19:46 AM UTC-5, Michael S wrote:
>> On Monday, March 28, 2022 at 9:54:31 AM UTC+3, Terje Mathisen wrote:
>> > Michael S wrote:
>> > > On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
>> > >> Marcus <m.de...@this.bitsnbites.eu> schrieb:
>> > >>> Today I found a pretty nice looking 2nd order polynomial for
>> > >>> approximating 1/x in the interval [1.0, 2.0):
>> > >>>
>> > >>> (2*x^2 - 9*x + 13) / 6
>> > >> Interesting. Could be useful for the start of a Newton iteration
>> > >> for a reciprocal, where the iteration formula would be
>> > >> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
>> > >>>
>> > >>> Maybe this is known stuff, and I'm not sure that it would be useful for
>> > >>> anything, but I thought I'd share it anyway.
>> > >>>
>> > >>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
>> > >>> LUT.
>> > >>>
>> > >>> The coefficients are "nice" from an implementation perspective, if you
>> > >>> were to hard-wire the polynomial. It should be possible to implement the
>> > >>> dividend with a single multiplier and a few integer adders. The 3 in the
>> > >>> divisor is muddying the waters, though. Is there a fast way to divide by
>> > >>> three?
>> > >> The standard trick of dividing by three in fixed point: For
>> > >> 32-bit numbers, multiply the number by the magic number 2863311531
>> > >> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
>> > >> "multiply high" and shift right one bit). For floating point
>> > >> numbers, this would have to be adjusted somewhat.
>> > >>
>> > >
>> > > And that causes the whole calculation to be of approximately the same
>> > > computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
>> > > potentially much better coefficients than those suggested.
>> > > The suggestion of Marcus is very naive.
>> > > Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
>> > > Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
>> > > down the road.
>> > > The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
>> > > Or when you don't care about speed, but that case is sort of off topic.
>> > If you are doing SIMD, but without a reciprocal/reciprocal square root
>> > lookup opcode, then I would use the infamous invsqrt() trick.
><
>> Magic constant and such?
>> I don't think it's a good idea if you seek top performance.
>> One thing to realize is that if you went to trouble of using rsqrt() NR steps
>> based on custom polynomials, with different coefficients on each step, then it
>> is better to keep error one-sided on all steps except possibly the last one.
><
> Notice that the error term for N-R iterations flip the sign every iteration.
><
>> Apart from of being slightly more precise, such strategy has advantage of being
>> simpler to do with unsigned integer math.
><
> You need to choose the approximat such that the error is uniformly positive after
> the last iteration. This allows you to choose from {r+0, r+1, and r+2} for the properly
> rounded result; which is easier than choosing {r-1, r+0, r+1} for the properly rounded
> result.
><
>> > With the
>> > latest Cheby-style tweaks to all the constants, this delivers 10+ bits
>> > after a single iteration.
>> >
>> > This is significantly better than an in-register 4-bit/16-term table lookup.
>> >
>> > Anyway, I am sure Mitch has much better polynomials to do it all
>> > directly. :-)
>> Mitch appears to like Chebychev polynomials.
>> Likely, because they are easy to find out :-)
><
> Once you figure out how to calculate them, they are, indeed, easy to determine.
>>
>> They are very close to optimal (in sense of minimizing two-sided error, but
>> probably can be easily offseted to minimize one-sided error too) when the
>> interval of input error is already small. But when the interval is large, as it
>> is a case on the first step without lookup table and, may be, on the second
>> step, when the 1st step was of low order, then Chebychev polynomials are
>> measurably suboptimal vs true equiripple polys. The later can be found by Remez
>> exchange algorithm.
><
> Cheby polynomials are within 1 ULP of optimal, where Remez is better, it is
> a lot harder to calculate, and if you are working to the last bit of precision,
> even calculating Remez requires more precision. Cheby does not have this
> property.

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.

On Tuesday, March 29, 2022 at 11:45:08 AM UTC-5, Thomas Koenig wrote:
> MitchAlsup <Mitch...@aol.com> schrieb:
> > On Monday, March 28, 2022 at 3:19:46 AM UTC-5, Michael S wrote:
> >> On Monday, March 28, 2022 at 9:54:31 AM UTC+3, Terje Mathisen wrote:
> >> > Michael S wrote:
> >> > > On Monday, March 28, 2022 at 12:26:48 AM UTC+3, Thomas Koenig wrote:
> >> > >> Marcus <m.de...@this.bitsnbites.eu> schrieb:
> >> > >>> Today I found a pretty nice looking 2nd order polynomial for
> >> > >>> approximating 1/x in the interval [1.0, 2.0):
> >> > >>>
> >> > >>> (2*x^2 - 9*x + 13) / 6
> >> > >> Interesting. Could be useful for the start of a Newton iteration
> >> > >> for a reciprocal, where the iteration formula would be
> >> > >> x_n+1 = 2*x_n - a * x_n^2 (division-free), or x_n*(2-a*x_n).
> >> > >>>
> >> > >>> Maybe this is known stuff, and I'm not sure that it would be useful for
> >> > >>> anything, but I thought I'd share it anyway.
> >> > >>>
> >> > >>> The accuracy that you get is roughly 6-7 bits, so it's not better than a
> >> > >>> LUT.
> >> > >>>
> >> > >>> The coefficients are "nice" from an implementation perspective, if you
> >> > >>> were to hard-wire the polynomial. It should be possible to implement the
> >> > >>> dividend with a single multiplier and a few integer adders. The 3 in the
> >> > >>> divisor is muddying the waters, though. Is there a fast way to divide by
> >> > >>> three?
> >> > >> The standard trick of dividing by three in fixed point: For
> >> > >> 32-bit numbers, multiply the number by the magic number 2863311531
> >> > >> (0xAAAAAAAB) and shift right by 33 bits (in practice, use a
> >> > >> "multiply high" and shift right one bit). For floating point
> >> > >> numbers, this would have to be adjusted somewhat.
> >> > >>
> >> > >
> >> > > And that causes the whole calculation to be of approximately the same
> >> > > computational complexity as 2nd-order polynomials with arbitrary 32-bit coefficients,
> >> > > potentially much better coefficients than those suggested.
> >> > > The suggestion of Marcus is very naive.
> >> > > Any suggestion to do sqrt or rsqrt approximation without lookup table is a loose proposition.
> >> > > Even very small table at the first step, like 32 or 64 entries, helps a lot and saves many steps
> >> > > down the road.
> >> > > The only exception to that could be when you do very wide SIMD that has no tools for even small lookups.
> >> > > Or when you don't care about speed, but that case is sort of off topic.
> >> > If you are doing SIMD, but without a reciprocal/reciprocal square root
> >> > lookup opcode, then I would use the infamous invsqrt() trick.
> ><
> >> Magic constant and such?
> >> I don't think it's a good idea if you seek top performance.
> >> One thing to realize is that if you went to trouble of using rsqrt() NR steps
> >> based on custom polynomials, with different coefficients on each step, then it
> >> is better to keep error one-sided on all steps except possibly the last one.
> ><
> > Notice that the error term for N-R iterations flip the sign every iteration.
> ><
> >> Apart from of being slightly more precise, such strategy has advantage of being
> >> simpler to do with unsigned integer math.
> ><
> > You need to choose the approximat such that the error is uniformly positive after
> > the last iteration. This allows you to choose from {r+0, r+1, and r+2} for the properly
> > rounded result; which is easier than choosing {r-1, r+0, r+1} for the properly rounded
> > result.
> ><
> >> > With the
> >> > latest Cheby-style tweaks to all the constants, this delivers 10+ bits
> >> > after a single iteration.
> >> >
> >> > This is significantly better than an in-register 4-bit/16-term table lookup.
> >> >
> >> > Anyway, I am sure Mitch has much better polynomials to do it all
> >> > directly. :-)
> >> Mitch appears to like Chebychev polynomials.
> >> Likely, because they are easy to find out :-)
> ><
> > Once you figure out how to calculate them, they are, indeed, easy to determine.
> >>
> >> They are very close to optimal (in sense of minimizing two-sided error, but
> >> probably can be easily offseted to minimize one-sided error too) when the
> >> interval of input error is already small. But when the interval is large, as it
> >> is a case on the first step without lookup table and, may be, on the second
> >> step, when the 1st step was of low order, then Chebychev polynomials are
> >> measurably suboptimal vs true equiripple polys. The later can be found by Remez
> >> exchange algorithm.
> ><
> > Cheby polynomials are within 1 ULP of optimal, where Remez is better, it is
> > a lot harder to calculate, and if you are working to the last bit of precision,
> > even calculating Remez requires more precision. Cheby does not have this
> > property.
> 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
<
with 6.63 bits of accuracy.
<-------------------------------------------------------------------------------------------------------------------
Why does your Remez get less precision than Chebychev ?

MitchAlsup <MitchAlsup@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
accuracy (a weight of x), and the maximum relative error is indeed
1/99 = 0.01010101.. A low bound for the relative error is better
as a starting point for Newton-Raphson, for example.

Optimizting for absolute error would have given a different result.

Michael S <already5chosen@yahoo.com> schrieb:
> On Tuesday, March 29, 2022 at 12:44:55 AM UTC+3, MitchAlsup wrote:
>> On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
>> > On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
>> > > On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
>> > > > Michael S <already...@yahoo.com> schrieb:
>> > > > > On my system(s) I see another very strange timing effect:
>> > > > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
>> > > > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
>> > > > Source for the sqrt routine in libquadmath is here:
>> > > >
>> > > > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
>> > > >
>> > > > Hmm.. they first do a test if it is within double precision range,
>> > > > with two Newton iterations, then a test if it is within long double
>> > > > range with a single Newton iteration, and if it is outside then
>> > > > they pick apart the number and run two Newton iterations.
>> > > >
>> > > > Were the numbers outside the range of double inside the range
>> > > > of long double?
>> > > Yes.
>> > > I didn't test "outside long double".
>> > > If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
>> > > > That could explain things (and suggest
>> > > > an improvement).
>> > >
>> > > I think, the whole routine needs rewrite, rather than improvements.
>> > > The original author probably was not thinking very hard.
>> > Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
>> > For example,
>> > sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
>> > The correct answer is 1.03168097468685393293622695023074831927
>> <
>> This is a 1ULP error, perhaps simply improperly rounded
>
> 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 .

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)
<
> accuracy (a weight of x), and the maximum relative error is indeed
> 1/99 = 0.01010101.. A low bound for the relative error is better
> as a starting point for Newton-Raphson, for example.
>
> Optimizting for absolute error would have given a different result.

On Tuesday, March 29, 2022 at 10:46:25 PM UTC+3, Thomas Koenig wrote:
> Michael S <already...@yahoo.com> schrieb:
> > On Tuesday, March 29, 2022 at 12:44:55 AM UTC+3, MitchAlsup wrote:
> >> On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
> >> > On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
> >> > > On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
> >> > > > Michael S <already...@yahoo.com> schrieb:
> >> > > > > On my system(s) I see another very strange timing effect:
> >> > > > > The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
> >> > > > > it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
> >> > > > Source for the sqrt routine in libquadmath is here:
> >> > > >
> >> > > > https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
> >> > > >
> >> > > > Hmm.. they first do a test if it is within double precision range,
> >> > > > with two Newton iterations, then a test if it is within long double
> >> > > > range with a single Newton iteration, and if it is outside then
> >> > > > they pick apart the number and run two Newton iterations.
> >> > > >
> >> > > > Were the numbers outside the range of double inside the range
> >> > > > of long double?
> >> > > Yes.
> >> > > I didn't test "outside long double".
> >> > > If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
> >> > > > That could explain things (and suggest
> >> > > > an improvement).
> >> > >
> >> > > I think, the whole routine needs rewrite, rather than improvements.
> >> > > The original author probably was not thinking very hard.
> >> > Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
> >> > For example,
> >> > sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
> >> > The correct answer is 1.03168097468685393293622695023074831927
> >> <
> >> This is a 1ULP error, perhaps simply improperly rounded
> >
> > 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 .

Funny, tests vs MPFR is exactly what I was doing for a last hour.
But I didn't know about existence mpfr_get_float128 ().
Was converting back and force via strings. Like that:

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);
}

My method is different, slower but more robust.
I start from pseudo-random float128 x, calculate res=sqrtq(x), then convert to MPFR
with 256-bit mantissa, calculate ref=mpfr_sqrt(), convert res to MPFR,
do a diff and compare a diff vs 0.5 ULP.

But at this point I am quite convinced of what I said couple of days ago: quadmath implementation
is hopeless and should be completely rewritten.

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

(And for 6.63 vs. 6.62 - I rounded down log[2](1/99) not to
overstate the accuracy :-)

MitchAlsup wrote:
> On Tuesday, March 29, 2022 at 7:57:35 AM UTC-5, Terje Mathisen wrote:
>> MitchAlsup wrote:
>>> On Monday, March 28, 2022 at 4:00:40 PM UTC-5, Michael S wrote:
>>>> On Sunday, March 27, 2022 at 2:32:34 AM UTC+3, Michael S wrote:
>>>>> On Sunday, March 27, 2022 at 1:25:29 AM UTC+3, Thomas Koenig wrote:
>>>>>> Michael S <already...@yahoo.com> schrieb:
>>>>>>> On my system(s) I see another very strange timing effect:
>>>>>>> The speed of sqrtq() strongly depends on the range of the input and does it in unexpected manner:
>>>>>>> it is much faster (1.5x to 2x) when input is *outside* the range of double precision!
>>>>>> Source for the sqrt routine in libquadmath is here:
>>>>>>
>>>>>> https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libquadmath/math/sqrtq.c;h=56ea5d3243c06df2276689849b448eebacae6367;hb=refs/heads/master
>>>>>>
>>>>>> Hmm.. they first do a test if it is within double precision range,
>>>>>> with two Newton iterations, then a test if it is within long double
>>>>>> range with a single Newton iteration, and if it is outside then
>>>>>> they pick apart the number and run two Newton iterations.
>>>>>>
>>>>>> Were the numbers outside the range of double inside the range
>>>>>> of long double?
>>>>> Yes.
>>>>> I didn't test "outside long double".
>>>>> If "long double" is 80-bit x87 format then "outside" would be mostly sub-normals, right?
>>>>>> That could explain things (and suggest
>>>>>> an improvement).
>>>>>
>>>>> I think, the whole routine needs rewrite, rather than improvements.
>>>>> The original author probably was not thinking very hard.
>>>> Hmm, on my system it's not just slow. Sometimes it rounds wrongly.
>>>> For example,
>>>> sqrtq(1.06436563353081694552652292498817875110) => 1.03168097468685393293622695023074812668
>>>> The correct answer is 1.03168097468685393293622695023074831927
>>> <
>>> This is a 1ULP error, perhaps simply improperly rounded (or the 1:million cases that require the second N-R
>>> iteration after you have 111 accurate fraction bits.
>>> <
>>> BTW: Itanium only gave 0.502 ULP FDIV accuracy (DP).
>> 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:

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

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

On Wednesday, March 30, 2022 at 12:27:22 AM UTC-6, Terje Mathisen wrote:

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

And the thing is, though, that for division (and, likely, square root) if you use
a fast algorithm, exact rounding is no longer easy, but requires significant extra
hardware, or a detectable amount of extra time - and, yes, one extra cycle is
detectable.

For addition, subtraction, and multiplication - only - is the extra effort for exact
rounding something that takes only a few gates _and_ a few gate delays that
will normally be swallowed up in the length of a cycle.

John Savard

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?

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?

Terje is (or was until recently?) a member of IEEE-754 Standard Committee,
so you can take his answer (i.e. Yes) in the above post as authoritative.