>>>> There are groups who have shown that exactly rounded trancendental
>>>> functions are in fact achievable with maybe 3X reduced performance.
>> That much? I had the impression it was significantly cheaper.
> The J. M. Muller book indicates about 2× to 2.5×

The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project claims
to get much better performance (basically, in the same ballpark as
not-correctly-rounded implementations).

[ Their key insight is the idea that to get correct rounding, you
shouldn't try to compute the best approximation of the exact result
and then round, but you should instead try to compute any
approximation whose rounding gives the correct result. ]

My impression was that their performance was good enough that the case
for not-correctly-rounded implementations becomes very weak.

>> I don't know what are/were the motivations for the people working on
>> exact transcendentals, but they have applications unrelated to the fact
>> that they're "better": the main benefit (from this here PL guy) is that
>> it gives them a reliable, reproducible semantics.
>> Bit-for-bit reproducibility makes several things much easier.
>
> Consider moving an application which uses libm from machine to machine.
> When libm is correctly rounded, there is no issue at all; not so other-
> wise.

Exactly!
[ Or should I say "Correctly rounded!"? ]

Stefan

On Wed, 20 Mar 2024 09:54:36 -0400
Stefan Monnier <monnier@iro.umontreal.ca> wrote:

> >>>> There are groups who have shown that exactly rounded
> >>>> trancendental functions are in fact achievable with maybe 3X
> >>>> reduced performance.
> >> That much? I had the impression it was significantly cheaper.
> > The J. M. Muller book indicates about 2× to 2.5×
>
> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project
> claims to get much better performance (basically, in the same
> ballpark as not-correctly-rounded implementations).
>

I had only read the 1st page.
It sounds like they are not particularly interested in IEEE binary64
which appears to be the major point of interest of math libs of
conventional languages.

> [ Their key insight is the idea that to get correct rounding, you
> shouldn't try to compute the best approximation of the exact result
> and then round, but you should instead try to compute any
> approximation whose rounding gives the correct result. ]
>
> My impression was that their performance was good enough that the case
> for not-correctly-rounded implementations becomes very weak.
>

It all depend of what you compare against.
For scalar call for majority of transcendental functions on IEEE-754
list, it's probably very easy to get correctly rounded binary32 results
in approximately the same time as results calculated with max. err of,
say, 0.75 ULP. Especially so if target machine has fast binary64
arithmetic.
But in practice when we use lower (than binary64) precision we often
care about vector performance rather than scalar.
I.e. we care little about speed of sinf(), but want ippsTone_32f() as
fast as possible. In case you wonder, this function is part Intel
Performance Primitives and it is FAST. Writing correctly rounded
function that approaches the speed of this *almost* correctly
rounded routine (I think, for sane input ranges it's better than
0.55 ULP) would not be easy at all!

> >> I don't know what are/were the motivations for the people working
> >> on exact transcendentals, but they have applications unrelated to
> >> the fact that they're "better": the main benefit (from this here
> >> PL guy) is that it gives them a reliable, reproducible semantics.
> >> Bit-for-bit reproducibility makes several things much easier.
> >
> > Consider moving an application which uses libm from machine to
> > machine. When libm is correctly rounded, there is no issue at all;
> > not so other- wise.
>
> Exactly!
> [ Or should I say "Correctly rounded!"? ]
>
>
> Stefan

You like this proposal because you are implementer of the language/lib.
It makes your regression tests easier. And it's good challenge.
I don't like it because I am primarily user of the language/lib. My
floating-point tests have zero chance of repeatability of this sort for
a thousand of other reasons.
I don't want to pay for correct rounding of transcendental functions.
Neither in speed and especially nor in tables footprint. Not even a
little. Because for me there are no advantages.
Now, there are things that I am ready to pay for. E.g. preservation of
mathematical properties of original exact function. I.e. if original is
monotonic on certain interval then I do want at least weak monotonicity
of approximation. If original is even (or odd) I want the same for
approximation. If original never exceed 1, I want the same
for approximation. Etc... But correct rounding is not on the list.

>> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project
>> claims to get much better performance (basically, in the same
>> ballpark as not-correctly-rounded implementations).
> I had only read the 1st page.
> It sounds like they are not particularly interested in IEEE binary64
> which appears to be the major point of interest of math libs of
> conventional languages.

Indeed, I hoped by now they had attacked the world of 64bit floats, but
it looks like it's not the case yet. In theory the same idea is
applicable there and should give similar results, but of course it's
harder for two reasons:

- The search space to consider is much larger, making it much more
costly to do exhaustive testing (and to collect the info they need to
choose/compute the polynomials).
- You can't cheaply use higher-precision floats for internal computations.

> It all depend of what you compare against.
> For scalar call for majority of transcendental functions on IEEE-754
> list, it's probably very easy to get correctly rounded binary32 results
> in approximately the same time as results calculated with max. err of,
> say, 0.75 ULP. Especially so if target machine has fast binary64
> arithmetic.

IIUC that was not the case before their work: it was "easy" to get the
correct result in 99% of the cases, but covering all 100% of the cases
used to be costly because those few cases needed a lot more
internal precision.

> I don't want to pay for correct rounding of transcendental functions.
> Neither in speed and especially nor in tables footprint. Not even a
> little. Because for me there are no advantages.

For that reason, correctly rounded results were nowhere near the menu,
indeed. But the proposition of Rlibm is that maybe we can have our cake
and eat it too, because (if all goes well) you don't have to pay for it
in speed nor in table size.

I'm hoping that pans out :-)

> Now, there are things that I am ready to pay for. E.g. preservation of
> mathematical properties of original exact function. I.e. if original is
> monotonic on certain interval then I do want at least weak monotonicity
> of approximation. If original is even (or odd) I want the same for
> approximation. If original never exceed 1, I want the same
> for approximation. Etc... But correct rounding is not on the list.

But correct rounding would give you those properties.

Stefan

On Wed, 20 Mar 2024 18:21:47 +0200, Michael S wrote:

> On Wed, 20 Mar 2024 09:54:36 -0400
> Stefan Monnier <monnier@iro.umontreal.ca> wrote:
>
>> >>>> There are groups who have shown that exactly rounded
>> >>>> trancendental functions are in fact achievable with maybe 3X
>> >>>> reduced performance.
>> >> That much? I had the impression it was significantly cheaper.
>> > The J. M. Muller book indicates about 2× to 2.5×
>>
>> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project
>> claims to get much better performance (basically, in the same
>> ballpark as not-correctly-rounded implementations).
>>
>
> I had only read the 1st page.
> It sounds like they are not particularly interested in IEEE binary64
> which appears to be the major point of interest of math libs of
> conventional languages.
>

I skimmed their logf(x) implementation. Their technique will
fall a part for binary64 (and other higher precisions). With
logf(x), they combine an argument step with table look-up and
a 5th-order polynomial approximation. The polynomial is computed
in double precision, and provides the correct rounding. For
binary64, the polynomial would require many more terms and
double-double or binary128 evaluation.

--
steve

Stefan Monnier wrote:

>>>>> There are groups who have shown that exactly rounded trancendental
>>>>> functions are in fact achievable with maybe 3X reduced performance.
>>> That much? I had the impression it was significantly cheaper.
>> The J. M. Muller book indicates about 2× to 2.5×

> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project claims
> to get much better performance (basically, in the same ballpark as
> not-correctly-rounded implementations).

> [ Their key insight is the idea that to get correct rounding, you
> shouldn't try to compute the best approximation of the exact result
> and then round, but you should instead try to compute any
> approximation whose rounding gives the correct result. ]

This sounds like the "Minefield Method"

> My impression was that their performance was good enough that the case
> for not-correctly-rounded implementations becomes very weak.

>>> I don't know what are/were the motivations for the people working on
>>> exact transcendentals, but they have applications unrelated to the fact
>>> that they're "better": the main benefit (from this here PL guy) is that
>>> it gives them a reliable, reproducible semantics.
>>> Bit-for-bit reproducibility makes several things much easier.
>>
>> Consider moving an application which uses libm from machine to machine.
>> When libm is correctly rounded, there is no issue at all; not so other-
>> wise.

> Exactly!
> [ Or should I say "Correctly rounded!"? ]

> Stefan

Michael S wrote:

> On Wed, 20 Mar 2024 09:54:36 -0400
> Stefan Monnier <monnier@iro.umontreal.ca> wrote:

>> [ Their key insight is the idea that to get correct rounding, you
>> shouldn't try to compute the best approximation of the exact result
>> and then round, but you should instead try to compute any
>> approximation whose rounding gives the correct result. ]
>>
>> My impression was that their performance was good enough that the case
>> for not-correctly-rounded implementations becomes very weak.
>>

> It all depend of what you compare against.
> For scalar call for majority of transcendental functions on IEEE-754
> list, it's probably very easy to get correctly rounded binary32 results
> in approximately the same time as results calculated with max. err of,
> say, 0.75 ULP. Especially so if target machine has fast binary64
> arithmetic.

> But in practice when we use lower (than binary64) precision we often
> care about vector performance rather than scalar.
> I.e. we care little about speed of sinf(), but want ippsTone_32f() as
> fast as possible. In case you wonder, this function is part Intel
> Performance Primitives and it is FAST. Writing correctly rounded
> function that approaches the speed of this *almost* correctly
> rounded routine (I think, for sane input ranges it's better than
> 0.55 ULP) would not be easy at all!

I challenge ANY software version of SIN() correctly rounded or not
to compete with my <patented> HW implementations for speed (or even
power).

>> >> I don't know what are/were the motivations for the people working
>> >> on exact transcendentals, but they have applications unrelated to
>> >> the fact that they're "better": the main benefit (from this here
>> >> PL guy) is that it gives them a reliable, reproducible semantics.
>> >> Bit-for-bit reproducibility makes several things much easier.
>> >
>> > Consider moving an application which uses libm from machine to
>> > machine. When libm is correctly rounded, there is no issue at all;
>> > not so other- wise.
>>
>> Exactly!
>> [ Or should I say "Correctly rounded!"? ]
>>
>>
>> Stefan

> You like this proposal because you are implementer of the language/lib.
> It makes your regression tests easier. And it's good challenge.
> I don't like it because I am primarily user of the language/lib. My
> floating-point tests have zero chance of repeatability of this sort for
> a thousand of other reasons.

> I don't want to pay for correct rounding of transcendental functions.

Even when the HW is 7× faster than SW algorithms ??

> Neither in speed and especially nor in tables footprint. Not even a
> little. Because for me there are no advantages.

My HW algorithms use no memory (cache, DRAM, or even LDs.)

> Now, there are things that I am ready to pay for. E.g. preservation of
> mathematical properties of original exact function.

You know, of course, that incorrectly rounded SIN() and COS() do not
maintain the property of SIN()^2+COS()^2 == 1.0

> I.e. if original is
> monotonic on certain interval then I do want at least weak monotonicity
> of approximation.

My HW algorithms have a numeric proof of this maintenance.

> If original is even (or odd) I want the same for
> approximation. If original never exceed 1, I want the same
> for approximation. Etc... But correct rounding is not on the list.

SIN(x) can be incorrectly rounded to be greater than 1.0.

Still want incorrect rounding--or just a polynomial that does not
have SI(X) > 1.0 ??

>> [ Their key insight is the idea that to get correct rounding, you
>> shouldn't try to compute the best approximation of the exact result
>> and then round, but you should instead try to compute any
>> approximation whose rounding gives the correct result. ]
>
> This sounds like the "Minefield Method"

Indeed it is.

Stefan

Stefan Monnier wrote:

>>> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project
>>> claims to get much better performance (basically, in the same
>>> ballpark as not-correctly-rounded implementations).
>> I had only read the 1st page.
>> It sounds like they are not particularly interested in IEEE binary64
>> which appears to be the major point of interest of math libs of
>> conventional languages.

> Indeed, I hoped by now they had attacked the world of 64bit floats, but
> it looks like it's not the case yet. In theory the same idea is
> applicable there and should give similar results, but of course it's
> harder for two reasons:

> - The search space to consider is much larger, making it much more
> costly to do exhaustive testing (and to collect the info they need to
> choose/compute the polynomials).
> - You can't cheaply use higher-precision floats for internal computations.

You can in HW, just not in SW. That is a strong algorithms to migrate
these functions from SW procedures to HW instructions. You do not need
FP128, just FP with a 64-bit fraction.

>> It all depend of what you compare against.

OK:: SIN(x) takes 19-cycles and is 0.5002 accurate with Payne & Hanek argument
reduction--Against any SW implementation.

>> For scalar call for majority of transcendental functions on IEEE-754
>> list, it's probably very easy to get correctly rounded binary32 results
>> in approximately the same time as results calculated with max. err of,
>> say, 0.75 ULP.

Yes, at the cost of 4× larger tables.

> Especially so if target machine has fast binary64
>> arithmetic.

Obviously.

> IIUC that was not the case before their work: it was "easy" to get the
> correct result in 99% of the cases, but covering all 100% of the cases
> used to be costly because those few cases needed a lot more
> internal precision.

Muller indicates one typically need 2×n+6 to 2×n+12 bits to get correct
roundings 100% of the time. FP128 only has 2×n+3 and is insufficient by
itself.

>> I don't want to pay for correct rounding of transcendental functions.
>> Neither in speed and especially nor in tables footprint. Not even a
>> little. Because for me there are no advantages.

> For that reason, correctly rounded results were nowhere near the menu,
> indeed. But the proposition of Rlibm is that maybe we can have our cake
> and eat it too, because (if all goes well) you don't have to pay for it
> in speed nor in table size.

> I'm hoping that pans out :-)

>> Now, there are things that I am ready to pay for. E.g. preservation of
>> mathematical properties of original exact function. I.e. if original is
>> monotonic on certain interval then I do want at least weak monotonicity
>> of approximation. If original is even (or odd) I want the same for
>> approximation. If original never exceed 1, I want the same
>> for approximation. Etc... But correct rounding is not on the list.

> But correct rounding would give you those properties.

> Stefan

Steven G. Kargl wrote:

> On Wed, 20 Mar 2024 18:21:47 +0200, Michael S wrote:

>> On Wed, 20 Mar 2024 09:54:36 -0400
>> Stefan Monnier <monnier@iro.umontreal.ca> wrote:
>>
>>> >>>> There are groups who have shown that exactly rounded
>>> >>>> trancendental functions are in fact achievable with maybe 3X
>>> >>>> reduced performance.
>>> >> That much? I had the impression it was significantly cheaper.
>>> > The J. M. Muller book indicates about 2× to 2.5×
>>>
>>> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project
>>> claims to get much better performance (basically, in the same
>>> ballpark as not-correctly-rounded implementations).
>>>
>>
>> I had only read the 1st page.
>> It sounds like they are not particularly interested in IEEE binary64
>> which appears to be the major point of interest of math libs of
>> conventional languages.
>>

> I skimmed their logf(x) implementation. Their technique will
> fall a part for binary64 (and other higher precisions). With
> logf(x), they combine an argument step with table look-up and
> a 5th-order polynomial approximation. The polynomial is computed
> in double precision, and provides the correct rounding. For
> binary64, the polynomial would require many more terms and
> double-double or binary128 evaluation.

The first 7 terms of the DP polynomial do not require FP128, whereas
the last 3 do/will.

Whereas: HW with 64-bits of fraction does not need FP128 at all.
64-bits of fraction with 64-bit coefficients added into a 128-bit
accumulator is more than sufficient.

Note FP64 has 53-bits of fraction.....

On Wed, 20 Mar 2024 20:33:44 +0000
mitchalsup@aol.com (MitchAlsup1) wrote:

> Michael S wrote:
>
> > On Wed, 20 Mar 2024 09:54:36 -0400
> > Stefan Monnier <monnier@iro.umontreal.ca> wrote:
>
> >> [ Their key insight is the idea that to get correct rounding, you
> >> shouldn't try to compute the best approximation of the exact
> >> result and then round, but you should instead try to compute any
> >> approximation whose rounding gives the correct result. ]
> >>
> >> My impression was that their performance was good enough that the
> >> case for not-correctly-rounded implementations becomes very weak.
> >>
>
> > It all depend of what you compare against.
> > For scalar call for majority of transcendental functions on IEEE-754
> > list, it's probably very easy to get correctly rounded binary32
> > results in approximately the same time as results calculated with
> > max. err of, say, 0.75 ULP. Especially so if target machine has
> > fast binary64 arithmetic.
>
> > But in practice when we use lower (than binary64) precision we often
> > care about vector performance rather than scalar.
> > I.e. we care little about speed of sinf(), but want ippsTone_32f()
> > as fast as possible. In case you wonder, this function is part Intel
> > Performance Primitives and it is FAST. Writing correctly rounded
> > function that approaches the speed of this *almost* correctly
> > rounded routine (I think, for sane input ranges it's better than
> > 0.55 ULP) would not be easy at all!
>
> I challenge ANY software version of SIN() correctly rounded or not
> to compete with my <patented> HW implementations for speed (or even
> power).
>

Before you post this response, you could as well look at what
ippsTone_32f() is doing. Hint - it's not generic scalar sin().
IMHO, for long enough vector and on modern enough Intel or AMD CPU it
will very easily beat any scalar-oriented binary64-oriented HW
implementation of sin() or cos().
This function is not about latency. It's about throughput.

AFAIR, youu were quite surprised by speed (throughput) of another IPP
primitive, ippsSin_64f_A53() when I posted results of timing
measurement here less than 2 yeears ago. So, before you issue a
challenge, just take into account that ippsTone_32f() is both more
specialized than ippsSin_64f_A53() and has much lower precision. So,
while I didn't test, I expect that it is much much faster.

Stefan Monnier wrote:
>>>>> There are groups who have shown that exactly rounded trancendental
>>>>> functions are in fact achievable with maybe 3X reduced performance.
>>> That much? I had the impression it was significantly cheaper.
>> The J. M. Muller book indicates about 2× to 2.5×
>
> The [Rlibm](https://people.cs.rutgers.edu/~sn349/rlibm/) project claims
> to get much better performance (basically, in the same ballpark as
> not-correctly-rounded implementations).
>
> [ Their key insight is the idea that to get correct rounding, you
> shouldn't try to compute the best approximation of the exact result
> and then round, but you should instead try to compute any
> approximation whose rounding gives the correct result. ]

That is indeed interesting. However, it is also very interesting that
they only do this for 32-bit or less. I.e. the domains for which it is
almost trivially easy to verify the results by checking all possible
inputs. :-)
>
> My impression was that their performance was good enough that the case
> for not-correctly-rounded implementations becomes very weak.

I agree in priniciple: If you can use polys that are not much more
complicated than the min-max/cheby case, and which always round to the
desired values, then that's an obvious good.
...
Reading the full log2f() code, it seems that they can use double for all
evaluations, and with a single (!) mantissa exception, this produces the
correct results for all inputs and all rounding modes. :-)

I.e. with 53 bits to work with, giving up about one ulp for the range
reduction, the 52 remaining bits corresponds to 2n+6 bits available for
the 5-term poly to evaluate a final float result.

When asking for perfectly rounded trancendentals, with approximately the
same runtime, the float case is a form of cheating, simply because
current FPUs tend to run float and double at the same speed.

OTOH, I'm still persuaded that for a float library, using this approach
might in fact be included in the 754 standard.

Doing the same for double by "simply" doing all operations with fp128
variables would definitely take significantly longer, and verification
would be an "interesting" problem. (Interesting in the "multiple PhDs"
domain.)

Terje

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

MitchAlsup1 wrote:
> Stefan Monnier wrote:
>> IIUC that was not the case before their work: it was "easy" to get the
>> correct result in 99% of the cases, but covering all 100% of the cases
>> used to be costly because those few cases needed a lot more
>> internal precision.
>
> Muller indicates one typically need 2×n+6 to 2×n+12 bits to get correct
> roundings 100% of the time. FP128 only has 2×n+3 and is insufficient by
> itself.

I agree with everything else you've written about this subject, but
afair, fp128 is using 1:15:112 while double is of course 1:10:53.

On the one hand, 53*2+6 -> 112, on the other hand (if we include the
hidden bits) we get 54*2+5 -> 113.

So significantly more than 2n+3 but not enough on its own to guarantee
correct rounding.

As Michael S have mentioned, we want these algorithms to work for
vector/SIMD calculations, and at that point both lookup tables and
widening the size of the temporaries are very costly.

Terje

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

On Thu, 21 Mar 2024 08:52:18 +0100
Terje Mathisen <terje.mathisen@tmsw.no> wrote:

> MitchAlsup1 wrote:
> > Stefan Monnier wrote:
> >> IIUC that was not the case before their work: it was "easy" to get
> >> the correct result in 99% of the cases, but covering all 100% of
> >> the cases used to be costly because those few cases needed a lot
> >> more internal precision.
> >
> > Muller indicates one typically need 2×n+6 to 2×n+12 bits to get
> > correct roundings 100% of the time. FP128 only has 2×n+3 and is
> > insufficient by itself.
>
> I agree with everything else you've written about this subject, but
> afair, fp128 is using 1:15:112 while double is of course 1:10:53.
>

IEEE-754 binary64 is 1:11:52 :-)

But anyway I am skeptical about Miller's rules of thumb.
I'd expect that different transcendental functions would exercise
non-trivially different behaviors, mostly because they have different
relationships between input and output ranges. Some of them compress
wider inputs into narrower output and some do the opposite.
Yet another factor is luck.

Besides, I see nothing special about binary128 as a helper format.
It is not supported on wast majority of HW, And even when it is
supported, like on IBM POWER, for majority of operations it is slower
than emulated 128-bit fixed-point. Fix-point is more work for coder, but
sounds like more sure path to success.

> On the one hand, 53*2+6 -> 112, on the other hand (if we include the
> hidden bits) we get 54*2+5 -> 113.
>
> So significantly more than 2n+3 but not enough on its own to
> guarantee correct rounding.
>
> As Michael S have mentioned, we want these algorithms to work for
> vector/SIMD calculations, and at that point both lookup tables and
> widening the size of the temporaries are very costly.
>
> Terje
>

Michael S wrote:

> On Thu, 21 Mar 2024 08:52:18 +0100
> Terje Mathisen <terje.mathisen@tmsw.no> wrote:

>> MitchAlsup1 wrote:
>> > Stefan Monnier wrote:
>> >> IIUC that was not the case before their work: it was "easy" to get
>> >> the correct result in 99% of the cases, but covering all 100% of
>> >> the cases used to be costly because those few cases needed a lot
>> >> more internal precision.
>> >
>> > Muller indicates one typically need 2×n+6 to 2×n+12 bits to get
>> > correct roundings 100% of the time. FP128 only has 2×n+3 and is
>> > insufficient by itself.
>>
>> I agree with everything else you've written about this subject, but
>> afair, fp128 is using 1:15:112 while double is of course 1:10:53.
>>

> IEEE-754 binary64 is 1:11:52 :-)

> But anyway I am skeptical about Miller's rules of thumb.

See Chapter 10 "Elementary Functions: Algorithms and Implementation"
Jean-Michael Muller {you can find a free copy on the net} and in
particular tables 10.5-10.14 -- quoting from the book::

"In his PhD dissertation [206], Lef`evre gives algorithms for finding
the worst cases of the table maker’s dilemma. These algorithms use linear
approximations and massive parallelism. A recent presentation of these
algorithms, with some improvements, can be found in [207]. We have run
Lef`evre’s algorithms to find worst cases in double-precision for the
most common functions and domains. The results obtained so far, first
published in [208] are given in Tables 10.5 to 10.14

They are NOT rules of thumb !! But go read it for yourself.

> I'd expect that different transcendental functions would exercise
> non-trivially different behaviors, mostly because they have different
> relationships between input and output ranges. Some of them compress
> wider inputs into narrower output and some do the opposite.
> Yet another factor is luck.

> Besides, I see nothing special about binary128 as a helper format.
> It is not supported on wast majority of HW, And even when it is
> supported, like on IBM POWER, for majority of operations it is slower
> than emulated 128-bit fixed-point. Fix-point is more work for coder, but
> sounds like more sure path to success.

>> On the one hand, 53*2+6 -> 112, on the other hand (if we include the
>> hidden bits) we get 54*2+5 -> 113.
>>
>> So significantly more than 2n+3 but not enough on its own to
>> guarantee correct rounding.
>>
>> As Michael S have mentioned, we want these algorithms to work for
>> vector/SIMD calculations, and at that point both lookup tables and
>> widening the size of the temporaries are very costly.
>>
>> Terje
>>

On 2/24/2024 7:26 PM, Chris M. Thomasson wrote:
> On 2/24/2024 2:52 PM, Lawrence D'Oliveiro wrote:
>> On Sat, 24 Feb 2024 12:40:52 -0800, Chris M. Thomasson wrote:
>>
[...]
> think of the following interesting aspect:
>
> https://forums.parallax.com/discussion/147522/dog-leg-hypotenuse-approximation
>
> Vs "needing" to get much more accurate results, for say medical imaging
> wrt volumetric renders...

Another "high cycle" fractal of mine that benefits from higher
precision, under iteration:

https://paulbourke.net/fractals/cubicjulia

Michael S wrote:
> On Thu, 21 Mar 2024 08:52:18 +0100
> Terje Mathisen <terje.mathisen@tmsw.no> wrote:
>
>> MitchAlsup1 wrote:
>>> Stefan Monnier wrote:
>>>> IIUC that was not the case before their work: it was "easy" to get
>>>> the correct result in 99% of the cases, but covering all 100% of
>>>> the cases used to be costly because those few cases needed a lot
>>>> more internal precision.
>>>
>>> Muller indicates one typically need 2×n+6 to 2×n+12 bits to get
>>> correct roundings 100% of the time. FP128 only has 2×n+3 and is
>>> insufficient by itself.
>>
>> I agree with everything else you've written about this subject, but
>> afair, fp128 is using 1:15:112 while double is of course 1:10:53.
>>
>
> IEEE-754 binary64 is 1:11:52 :-)

Oops! Mea Culpa!
I _do_ know that double has a 10-bit exponent bias (1023), so it has to
be 11 bits wide. :-(

>
> But anyway I am skeptical about Miller's rules of thumb.
> I'd expect that different transcendental functions would exercise
> non-trivially different behaviors, mostly because they have different
> relationships between input and output ranges. Some of them compress
> wider inputs into narrower output and some do the opposite.
> Yet another factor is luck.

I agree, this is a per-function problem, with some being substantially
harder than others.
>
> Besides, I see nothing special about binary128 as a helper format.
> It is not supported on wast majority of HW, And even when it is
> supported, like on IBM POWER, for majority of operations it is slower
> than emulated 128-bit fixed-point. Fix-point is more work for coder, but
> sounds like more sure path to success.

In my own code (since I don't have Mitch's ability to use much wider
internal fp formats) I also prefer 64-bit u64 as the working chunk size.

Almost 30 years ago, during the FDIV workaround, I needed a q&d way to
verify that our fpatan2 replacement was correct, so what I did was to
write a 1:31:96 format library over a weekend.

Yeah, it was much more exponent than needed, but with only 32-bit
registers available it was much easier to get the asm correct.

For the fpatan2 I used a dead simple approach with little range
reduction, just a longish Taylor series (i.e. no Cheby optimizations).

I had previously written 3-4 different iomplementations of arbitrary
precision atan2() when I wanted to calculatie as many digits of pi as
possible, so I just reused one of those algorithms.

Terje

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