On 3/15/2024 2:26 PM, Chris M. Thomasson wrote:
[...]
> Notice how the numbers are small? However, if we zoom in enough, then
> 2^53 might become an issue..
[...]

Argh! I mean if we zoom in enough, then only 2^53 of precision might
become an issue... We can go in really deep:

A deep zoom might be 10^275. Somebody did one to the 10^4141! I feel
like including a youtube link here, but I will refrain out of respect
for you.

On 3/15/2024 1:59 PM, Chris M. Thomasson wrote:
> On 3/15/2024 3:23 AM, Terje Mathisen wrote:
>> Michael, I for the main part agree with you here, i.e. calculating
>> sin(x) with x larger than 2^53 or so, is almost certainly stupid.
> [...]
>
> ;^D tooooooo big. :^)

Actually, I can think of one of my fractals that tries to scale back
numbers that get too big.

https://www.shadertoy.com/view/XtscDl

When a complex number breaches an escape time barrier, my code simply
scales it back, say z = z * .242, and tries the fractal again. It has
the effect of exploding the fractal out beyond its original boundaries...

>
> Now, wrt the results, arbitrary precision for trig is useful, in say...
> Deep fractal zooms...
>
> Zooming in really deep in say something like this, well the precision of
> trig can become an issue:
>
> https://paulbourke.net/fractals/multijulia/
>
> Trig would be used, say, in rectangular to-from polar forms wrt getting
> the n-ary roots of a complex number?
>

"Chris M. Thomasson" <chris.m.thomasson.1@gmail.com> writes:
> On 3/15/2024 2:26 PM, Chris M. Thomasson wrote:
> [...]
>> Notice how the numbers are small? However, if we zoom in enough,
>> then 2^53 might become an issue..
> [...]
>
> Argh! I mean if we zoom in enough, then only 2^53 of precision might
> become an issue... We can go in really deep:
>
> A deep zoom might be 10^275. Somebody did one to the 10^4141! I feel
> like including a youtube link here, but I will refrain out of respect
> for you.

I don't object to YouTube links as long as they're reasonably relevant
and there's enough information for readers to decide whether it's worth
their time to open them.

--
Keith Thompson (The_Other_Keith) Keith.S.Thompson+u@gmail.com
Working, but not speaking, for Medtronic
void Void(void) { Void(); } /* The recursive call of the void */

Keith Thompson wrote:

> I can see how computing sin(x) with high precision for "reasonable"
> values of x would be useful, but does any of that benefit from being
> able to compute sin(2^53) accurately?

Because accurate argument reduction reduces the burden on the programmer
to remain within his sandbox. Since early 1990s, when Payne and Hanek
argument reduction became widely available, there is no reason NOT to
perform argument reduction with great accuracy. NOTE: you don't HAVE to
use Payne and Hanek when simpler strategies word adequately; but when
they do not, you do not have to "Throw up your hands" and run away.

Michael S wrote:

> To make it less tempting, you could try to push for inclusion of
> rsqrt() into basic set. Long overdue, IMHO.

Many Newton-Raphson iterations converge more rapidly with RSQRT than
with SQRT, for reasons I never looked that deeply into.

> Right now, I can't think of any other transcendental that I really want
> to elevate to higher status. It seems to me that elevation of log2(x)
> and of 2**x will do no harm, but I am not sure about usefulness.

Ln2 and EXP2 are the basis of My 66000 log and EXP families, performed in
HW function unit.

Ln2 has the property that a binade [1..2) maps directly to another binade
[2..4); exp2 has a similar property. Ln2 is easier to get good accuracy
from than Ln or Log as a starting point. But in any event:: you are always
within a high precision multiply of the correct result = round(53×106);

MitchAlsup1 wrote:
> Michael S wrote:
>
>> To make it less tempting, you could try to push for inclusion of
>> rsqrt() into basic set. Long overdue, IMHO.
>
> Many Newton-Raphson iterations converge more rapidly with RSQRT than
> with SQRT, for reasons I never looked that deeply into.

The classic NR iteration for rsqrt() uses only fmul & fadd/fsub, while
the naive sqrt() version needs a divide in every iteration.

Personally I prefer to use the y = rsqrt(x) form and then just multiply
by x in the end:

y = sqrt(x) = x/sqrt(x) = rsqrt(x)*x

In order to get the rounding correct from this form you can probably
just square the final result and compare this with the original input?
>
>> Right now, I can't think of any other transcendental that I really want
>> to elevate to higher status. It seems to me that elevation of log2(x)
>> and of 2**x will do no harm, but I am not sure about usefulness.
>
> Ln2 and EXP2 are the basis of My 66000 log and EXP families, performed in
> HW function unit.
>
> Ln2 has the property that a binade [1..2) maps directly to another binade
> [2..4); exp2 has a similar property. Ln2 is easier to get good accuracy
> from than Ln or Log as a starting point. But in any event:: you are
> always within a high precision multiply of the correct result =
> round(53×106);

Exactly the same reasoning I used to compute rsqrt() over the [1.0 -
4.0> range, or (simplifying the exp logic) over [0.5 - 2.0> so that the
last exp bit maps to the first vs second half of the range. Back before
we got approximate rsqrt() lookup opcodes, I wrote code that took the
last exp bit and the first n mantissa bits and looked up an optimal
starting point for that value + 0.5 ulp. Getting 11-12 bits this way
means approximatly correct float after a single iteration and useful
double after two.

Terje

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

On Sat, 16 Mar 2024 01:16:25 +0000
mitchalsup@aol.com (MitchAlsup1) wrote:

> Keith Thompson wrote:
>
> > I can see how computing sin(x) with high precision for "reasonable"
> > values of x would be useful, but does any of that benefit from being
> > able to compute sin(2^53) accurately?
>
> Because accurate argument reduction reduces the burden on the
> programmer to remain within his sandbox.

Not really.
sin(2**53) is not what programmer wanted, because 2**53 is not the x he
really wanted. It has nothing to do with calculation of sin().

> Since early 1990s, when
> Payne and Hanek argument reduction became widely available, there is
> no reason NOT to perform argument reduction with great accuracy.
> NOTE: you don't HAVE to use Payne and Hanek when simpler strategies
> word adequately; but when they do not, you do not have to "Throw up
> your hands" and run away.

Michael S wrote:

> On Sat, 16 Mar 2024 01:16:25 +0000
> mitchalsup@aol.com (MitchAlsup1) wrote:

>> Keith Thompson wrote:
>>
>> > I can see how computing sin(x) with high precision for "reasonable"
>> > values of x would be useful, but does any of that benefit from being
>> > able to compute sin(2^53) accurately?
>>
>> Because accurate argument reduction reduces the burden on the
>> programmer to remain within his sandbox.

> Not really.

Say you are a programmer and you receive a value like 2^53 from an
Input read and you wan the most accurate possible SIN( of that ).

With high precision argument reduction you have to do (wait for it)
NoThInG !!.

Without high precision argument reduction yo have to do a lot of
programming--programming that I insist should have been done by
the programmer responsible for sin(x).

> sin(2**53) is not what programmer wanted, because 2**53 is not the x he
> really wanted. It has nothing to do with calculation of sin().

And EvErYtHiNg to do with high precision argument reduction; and
maybe the programmer has no requirement on what he is allowed to
read in from a file and pass onto sin(x).

>> Since early 1990s, when
>> Payne and Hanek argument reduction became widely available, there is
>> no reason NOT to perform argument reduction with great accuracy.
>> NOTE: you don't HAVE to use Payne and Hanek when simpler strategies
>> word adequately; but when they do not, you do not have to "Throw up
>> your hands" and run away.

mitchalsup@aol.com (MitchAlsup1) writes:
>Michael S wrote:
>
>> On Sat, 16 Mar 2024 01:16:25 +0000
>> mitchalsup@aol.com (MitchAlsup1) wrote:
>
>>> Keith Thompson wrote:
>>>
>>> > I can see how computing sin(x) with high precision for "reasonable"
>>> > values of x would be useful, but does any of that benefit from being
>>> > able to compute sin(2^53) accurately?
>>>
>>> Because accurate argument reduction reduces the burden on the
>>> programmer to remain within his sandbox.
>
>> Not really.
>
>Say you are a programmer and you receive a value like 2^53 from an
>Input read and you wan the most accurate possible SIN( of that ).

That generally seems to be an highly unlikely scenario.

On Sat, 16 Mar 2024 17:22:50 +0000
mitchalsup@aol.com (MitchAlsup1) wrote:

> Michael S wrote:
>
> > On Sat, 16 Mar 2024 01:16:25 +0000
> > mitchalsup@aol.com (MitchAlsup1) wrote:
>
> >> Keith Thompson wrote:
> >>
> >> > I can see how computing sin(x) with high precision for
> >> > "reasonable" values of x would be useful, but does any of that
> >> > benefit from being able to compute sin(2^53) accurately?
> >>
> >> Because accurate argument reduction reduces the burden on the
> >> programmer to remain within his sandbox.
>
> > Not really.
>
> Say you are a programmer and you receive a value like 2^53 from an
> Input read and you wan the most accurate possible SIN( of that ).
>
> With high precision argument reduction you have to do (wait for it)
> NoThInG !!.
>
> Without high precision argument reduction yo have to do a lot of
> programming--programming that I insist should have been done by
> the programmer responsible for sin(x).
>

I never met programmers that wants do anything of this sort.
Programmers that I meet in real world could sometimes want something
else. They could want, for example, to draw 1000 points of wave of blue
LED with frequency 663.42 THz (663.42e12 Hz) with step of 10 nm (1e-8 m)
starting at distance of 1 million km. Well, not really, but at least I
can imagine them wanting to do it*.
Will double precision sin() routine that does perfect argument
reduction be any more helpful for their task then sin() implementation
that makes no effort at all beyound simple modulo 2*pi? No, they both
would be the same.
May be, sinpi() that does a perfect reduction more or less by
definition, be more helpful?
For naive folks - not at all.
For advanced folks - yes, but just a little. It still wouldn't be easy.
And among programmers that I know naives outnumber advanced by more
than factor of 10.

-----
* - For more real case, they could want to draw the same wave at
distance of 2m using single-precision arithmetic.

On 3/15/2024 2:13 PM, Chris M. Thomasson wrote:
[...]
> three digits of decimal precision:

OOPS! I mean four:

>
> 333/106 = 3.1415_094339622641509433962264151
[...]

mitchalsup@aol.com (MitchAlsup1) writes:
> Michael S wrote:
>> On Sat, 16 Mar 2024 01:16:25 +0000
>> mitchalsup@aol.com (MitchAlsup1) wrote:
>>> Keith Thompson wrote:
>>> > I can see how computing sin(x) with high precision for
>>> > "reasonable" values of x would be useful, but does any of that
>>> > benefit from being able to compute sin(2^53) accurately?
>>> Because accurate argument reduction reduces the burden on the
>>> programmer to remain within his sandbox.
>
>> Not really.
>
> Say you are a programmer and you receive a value like 2^53 from an
> Input read and you wan the most accurate possible SIN( of that ).

I can't think of a scenario where that would be useful (other than just
doing it for the sake of doing it).

If 2^53 represents a physical quantity, how likely is the actual value
to be known within ±π (+/i pi for those who prefer ASCII)?

If you can get better precision without too much extra cost, that's
great. I don't know enough to have an opinion about what the best
tradeoff is, but I presume it's going to be different depending on the
application.

Here's a C program that shows how precise sin(2^53) can be for types
float, double, and long double (I used gcc and glibc). The nextafter
functions are used to compute the nearest representable number. For
long double, the value of sin() changes by about 1 part in 1600, which
seems decent, but it's not nearly as precise as for values around 1.0.
For float and double, the imprecision of the argument is enough to make
the result practically meaningless.

#include <math.h>
#include <stdio.h>
#include <limits.h>
#include <float.h>
int main(void) {
{
printf("float (%zu bits, %d mantissa bits)\n", CHAR_BIT * sizeof (float), FLT_MANT_DIG);
const float x = (float)(1LL<<53);
const float y = nextafterf(x, x*2);
printf("%.8f %.8f\n", x, sinf(x));
printf("%.8f %.8f\n", y, sinf(y));
}
putchar('\n');
{
printf("double (%zu bits, %d mantissa bits)\n", CHAR_BIT * sizeof (double), DBL_MANT_DIG);
const double x = (double)(1LL<<53);
const double y = nextafter(x, x*2);
printf("%.8f %.8f\n", x, sin(x));
printf("%.8f %.8f\n", y, sin(y));
}
putchar('\n');
{
printf("long double (%zu bits, %d mantissa bits)\n", CHAR_BIT * sizeof (long double), LDBL_MANT_DIG);
const long double x = (long double)(1LL<<53);
const long double y = nextafterl(x, x*2);
printf("%.8Lf %.8Lf\n", x, sinl(x));
printf("%.8Lf %.8Lf\n", y, sinl(y));
}
}

Output:

float (32 bits, 24 mantissa bits)
9007199254740992.00000000 -0.84892595
9007200328482816.00000000 -0.34159181

double (64 bits, 53 mantissa bits)
9007199254740992.00000000 -0.84892596
9007199254740994.00000000 -0.12729655

long double (128 bits, 64 mantissa bits)
9007199254740992.00000000 -0.84892596
9007199254740992.00097656 -0.84944168

--
Keith Thompson (The_Other_Keith) Keith.S.Thompson+u@gmail.com
Working, but not speaking, for Medtronic
void Void(void) { Void(); } /* The recursive call of the void */

On 16/03/2024 23:19, Keith Thompson wrote:
> mitchalsup@aol.com (MitchAlsup1) writes:

>> Say you are a programmer and you receive a value like 2^53 from an
>> Input read and you wan the most accurate possible SIN( of that ).
>
> I can't think of a scenario where that would be useful (other than just
> doing it for the sake of doing it).
>
> If 2^53 represents a physical quantity, how likely is the actual value
> to be known within ±π (+/i pi for those who prefer ASCII)?
>
> If you can get better precision without too much extra cost, that's
> great. I don't know enough to have an opinion about what the best
> tradeoff is, but I presume it's going to be different depending on the
> application.
>
> Here's a C program that shows how precise sin(2^53) can be for types
> float, double, and long double (I used gcc and glibc). The nextafter
> functions are used to compute the nearest representable number. For
> long double, the value of sin() changes by about 1 part in 1600, which
> seems decent, but it's not nearly as precise as for values around 1.0.
> For float and double, the imprecision of the argument is enough to make
> the result practically meaningless.
>
> ...
> Output:
>
> float (32 bits, 24 mantissa bits)
> 9007199254740992.00000000 -0.84892595
> 9007200328482816.00000000 -0.34159181
>
> double (64 bits, 53 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740994.00000000 -0.12729655
>
> long double (128 bits, 64 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740992.00097656 -0.84944168

Is this output supposed to be different between gcc -O0 and gcc -O3?

bart <bc@freeuk.com> writes:
> On 16/03/2024 23:19, Keith Thompson wrote:
>> mitchalsup@aol.com (MitchAlsup1) writes:
>
>>> Say you are a programmer and you receive a value like 2^53 from an
>>> Input read and you wan the most accurate possible SIN( of that ).
>> I can't think of a scenario where that would be useful (other than
>> just
>> doing it for the sake of doing it).
>> If 2^53 represents a physical quantity, how likely is the actual
>> value
>> to be known within ±π (+/i pi for those who prefer ASCII)?
>> If you can get better precision without too much extra cost, that's
>> great. I don't know enough to have an opinion about what the best
>> tradeoff is, but I presume it's going to be different depending on the
>> application.
>> Here's a C program that shows how precise sin(2^53) can be for types
>> float, double, and long double (I used gcc and glibc). The nextafter
>> functions are used to compute the nearest representable number. For
>> long double, the value of sin() changes by about 1 part in 1600, which
>> seems decent, but it's not nearly as precise as for values around 1.0.
>> For float and double, the imprecision of the argument is enough to make
>> the result practically meaningless.
>>
>> ...
>> Output:
>> float (32 bits, 24 mantissa bits)
>> 9007199254740992.00000000 -0.84892595
>> 9007200328482816.00000000 -0.34159181
>> double (64 bits, 53 mantissa bits)
>> 9007199254740992.00000000 -0.84892596
>> 9007199254740994.00000000 -0.12729655
>> long double (128 bits, 64 mantissa bits)
>> 9007199254740992.00000000 -0.84892596
>> 9007199254740992.00097656 -0.84944168

>
> Is this output supposed to be different between gcc -O0 and gcc -O3?

No, why do you ask?

On my system, the output is identical with gcc and clang, and tcc at all
optimization levels, with both glibc and musl. I do get slightly
different results for long double on Cygwin vs. Ubuntu (Cygwin uses
newlib, not glibc).

--
Keith Thompson (The_Other_Keith) Keith.S.Thompson+u@gmail.com
Working, but not speaking, for Medtronic
void Void(void) { Void(); } /* The recursive call of the void */

On 17/03/2024 01:38, Keith Thompson wrote:
> bart <bc@freeuk.com> writes:
>> On 16/03/2024 23:19, Keith Thompson wrote:
>>> mitchalsup@aol.com (MitchAlsup1) writes:
>>
>>>> Say you are a programmer and you receive a value like 2^53 from an
>>>> Input read and you wan the most accurate possible SIN( of that ).
>>> I can't think of a scenario where that would be useful (other than
>>> just
>>> doing it for the sake of doing it).
>>> If 2^53 represents a physical quantity, how likely is the actual
>>> value
>>> to be known within ±π (+/i pi for those who prefer ASCII)?
>>> If you can get better precision without too much extra cost, that's
>>> great. I don't know enough to have an opinion about what the best
>>> tradeoff is, but I presume it's going to be different depending on the
>>> application.
>>> Here's a C program that shows how precise sin(2^53) can be for types
>>> float, double, and long double (I used gcc and glibc). The nextafter
>>> functions are used to compute the nearest representable number. For
>>> long double, the value of sin() changes by about 1 part in 1600, which
>>> seems decent, but it's not nearly as precise as for values around 1.0.
>>> For float and double, the imprecision of the argument is enough to make
>>> the result practically meaningless.
>>>
>>> ...
>>> Output:
>>> float (32 bits, 24 mantissa bits)
>>> 9007199254740992.00000000 -0.84892595
>>> 9007200328482816.00000000 -0.34159181
>>> double (64 bits, 53 mantissa bits)
>>> 9007199254740992.00000000 -0.84892596
>>> 9007199254740994.00000000 -0.12729655
>>> long double (128 bits, 64 mantissa bits)
>>> 9007199254740992.00000000 -0.84892596
>>> 9007199254740992.00097656 -0.84944168
>
>>
>> Is this output supposed to be different between gcc -O0 and gcc -O3?
>
> No, why do you ask?
>
> On my system, the output is identical with gcc and clang, and tcc at all
> optimization levels, with both glibc and musl. I do get slightly
> different results for long double on Cygwin vs. Ubuntu (Cygwin uses
> newlib, not glibc).

Because I get the results shown below. Even if the library is different
from yours, I would have assumed matching results.

------------------------------

C:\c>gcc c.c -O3

C:\c>a
float (32 bits, 24 mantissa bits)
9007199254740992.00000000 -0.84892595
9007200328482816.00000000 -0.34159181

double (64 bits, 53 mantissa bits)
9007199254740992.00000000 -0.84892596
9007199254740994.00000000 -0.12729655

long double (128 bits, 64 mantissa bits)
9007199254740992.00000000 -0.84892596
9007199254740992.00097656 -0.84944168

C:\c>gcc c.c -O0

C:\c>a
float (32 bits, 24 mantissa bits)
9007199254740992.00000000 -0.84893209
9007200328482816.00000000 -0.34160271

double (64 bits, 53 mantissa bits)
9007199254740992.00000000 -0.84893209
9007199254740994.00000000 -0.12728505

long double (128 bits, 64 mantissa bits)
9007199254740992.00000000 -0.84893209
9007199254740992.00097656 -0.84944780

C:\c>gcc --version
gcc (tdm64-1) 10.3.0

bart <bc@freeuk.com> writes:
> On 17/03/2024 01:38, Keith Thompson wrote:
>> bart <bc@freeuk.com> writes:
>>> On 16/03/2024 23:19, Keith Thompson wrote:
>>>> mitchalsup@aol.com (MitchAlsup1) writes:
>>>
>>>>> Say you are a programmer and you receive a value like 2^53 from an
>>>>> Input read and you wan the most accurate possible SIN( of that ).
>>>> I can't think of a scenario where that would be useful (other than
>>>> just
>>>> doing it for the sake of doing it).
>>>> If 2^53 represents a physical quantity, how likely is the actual
>>>> value
>>>> to be known within ±π (+/i pi for those who prefer ASCII)?
>>>> If you can get better precision without too much extra cost, that's
>>>> great. I don't know enough to have an opinion about what the best
>>>> tradeoff is, but I presume it's going to be different depending on the
>>>> application.
>>>> Here's a C program that shows how precise sin(2^53) can be for types
>>>> float, double, and long double (I used gcc and glibc). The nextafter
>>>> functions are used to compute the nearest representable number. For
>>>> long double, the value of sin() changes by about 1 part in 1600, which
>>>> seems decent, but it's not nearly as precise as for values around 1.0.
>>>> For float and double, the imprecision of the argument is enough to make
>>>> the result practically meaningless.
>>>>
>>>> ...
>>>> Output:
>>>> float (32 bits, 24 mantissa bits)
>>>> 9007199254740992.00000000 -0.84892595
>>>> 9007200328482816.00000000 -0.34159181
>>>> double (64 bits, 53 mantissa bits)
>>>> 9007199254740992.00000000 -0.84892596
>>>> 9007199254740994.00000000 -0.12729655
>>>> long double (128 bits, 64 mantissa bits)
>>>> 9007199254740992.00000000 -0.84892596
>>>> 9007199254740992.00097656 -0.84944168
>>
>>>
>>> Is this output supposed to be different between gcc -O0 and gcc -O3?
>> No, why do you ask?
>> On my system, the output is identical with gcc and clang, and tcc at
>> all
>> optimization levels, with both glibc and musl. I do get slightly
>> different results for long double on Cygwin vs. Ubuntu (Cygwin uses
>> newlib, not glibc).
>
>
> Because I get the results shown below. Even if the library is
> different from yours, I would have assumed matching results.
>
>
> ------------------------------
>
> C:\c>gcc c.c -O3
>
> C:\c>a
> float (32 bits, 24 mantissa bits)
> 9007199254740992.00000000 -0.84892595
> 9007200328482816.00000000 -0.34159181
>
> double (64 bits, 53 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740994.00000000 -0.12729655
>
> long double (128 bits, 64 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740992.00097656 -0.84944168
>
> C:\c>gcc c.c -O0
>
> C:\c>a
> float (32 bits, 24 mantissa bits)
> 9007199254740992.00000000 -0.84893209
> 9007200328482816.00000000 -0.34160271
>
> double (64 bits, 53 mantissa bits)
> 9007199254740992.00000000 -0.84893209
> 9007199254740994.00000000 -0.12728505
>
> long double (128 bits, 64 mantissa bits)
> 9007199254740992.00000000 -0.84893209
> 9007199254740992.00097656 -0.84944780
>
> C:\c>gcc --version
> gcc (tdm64-1) 10.3.0

I see similar results (possibly identical, I haven't checked) with
i686-w64-mingw32-gcc on Cygwin. I think it uses the same, or a similar,
runtime library as the tdm64 implementation you're using.

Comparing the generated assembly, at -O0 it generates a call to "sin",
while at -O1 and above it replaces the expression with a compile-time
constant. Apparently that compile-time constant matches the run-time
computed value for glibc, but not for whatever tdm64 uses.

I don't *think* it's a conformance issue as long as the precision is
within the standard's (fairly vague) requirements. When I modify the
program to start with sin(1.0) rather than sin(2**53), the output is
identical at all optimization levels.

Apparently glibc and the library used by tgm64 behave differently when
computing the sin of very large values. The result for sin(2**53) in
long double differs by about one part in 2**17 between -O0 and -O1, or
between gcc/glibc and tgm64.

--
Keith Thompson (The_Other_Keith) Keith.S.Thompson+u@gmail.com
Working, but not speaking, for Medtronic
void Void(void) { Void(); } /* The recursive call of the void */

On Sat, 16 Mar 2024 16:19:11 -0700
Keith Thompson <Keith.S.Thompson+u@gmail.com> wrote:

> mitchalsup@aol.com (MitchAlsup1) writes:
> > Michael S wrote:
> >> On Sat, 16 Mar 2024 01:16:25 +0000
> >> mitchalsup@aol.com (MitchAlsup1) wrote:
> >>> Keith Thompson wrote:
> >>> > I can see how computing sin(x) with high precision for
> >>> > "reasonable" values of x would be useful, but does any of that
> >>> > benefit from being able to compute sin(2^53) accurately?
> >>> Because accurate argument reduction reduces the burden on the
> >>> programmer to remain within his sandbox.
> >
> >> Not really.
> >
> > Say you are a programmer and you receive a value like 2^53 from an
> > Input read and you wan the most accurate possible SIN( of that ).
>
> I can't think of a scenario where that would be useful (other than
> just doing it for the sake of doing it).
>
> If 2^53 represents a physical quantity, how likely is the actual value
> to be known within ±π (+/i pi for those who prefer ASCII)?
>
> If you can get better precision without too much extra cost, that's
> great. I don't know enough to have an opinion about what the best
> tradeoff is, but I presume it's going to be different depending on the
> application.
>
> Here's a C program that shows how precise sin(2^53) can be for types
> float, double, and long double (I used gcc and glibc). The nextafter
> functions are used to compute the nearest representable number. For
> long double, the value of sin() changes by about 1 part in 1600, which
> seems decent, but it's not nearly as precise as for values around 1.0.
> For float and double, the imprecision of the argument is enough to
> make the result practically meaningless.
>
> #include <math.h>
> #include <stdio.h>
> #include <limits.h>
> #include <float.h>
> int main(void) {
> {
> printf("float (%zu bits, %d mantissa bits)\n", CHAR_BIT *
> sizeof (float), FLT_MANT_DIG); const float x = (float)(1LL<<53);
> const float y = nextafterf(x, x*2);
> printf("%.8f %.8f\n", x, sinf(x));
> printf("%.8f %.8f\n", y, sinf(y));
> }
> putchar('\n');
> {
> printf("double (%zu bits, %d mantissa bits)\n", CHAR_BIT *
> sizeof (double), DBL_MANT_DIG); const double x = (double)(1LL<<53);
> const double y = nextafter(x, x*2);
> printf("%.8f %.8f\n", x, sin(x));
> printf("%.8f %.8f\n", y, sin(y));
> }
> putchar('\n');
> {
> printf("long double (%zu bits, %d mantissa bits)\n", CHAR_BIT
> * sizeof (long double), LDBL_MANT_DIG); const long double x = (long
> double)(1LL<<53); const long double y = nextafterl(x, x*2);
> printf("%.8Lf %.8Lf\n", x, sinl(x));
> printf("%.8Lf %.8Lf\n", y, sinl(y));
> }
> }
>
> Output:
>
> float (32 bits, 24 mantissa bits)
> 9007199254740992.00000000 -0.84892595
> 9007200328482816.00000000 -0.34159181
>
> double (64 bits, 53 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740994.00000000 -0.12729655
>
> long double (128 bits, 64 mantissa bits)
> 9007199254740992.00000000 -0.84892596
> 9007199254740992.00097656 -0.84944168
>

As written, your example does not emphasize that the problem has
nothing to do with implementation of sinX() library routine.
It's best illustrated by followup conversation with bart, IMHO 100%
O.T.
To make the point more clear I'd rather change it to following form:

#include <math.h>
#include <stdio.h>
#include <limits.h>
#include <float.h>

void foo(long double x1, long double x2)
{ const double y1 = (double)sinl(x1);
const double y2 = (double)sinl(x2);
printf("%.20Le %.17f\n", x1, y1);
printf("%.20Le %.17f\n", x2, y2);
}

int main(void) {
const float x0 = (float)(1LL<<53);
{
printf("float (%zu bits, %d mantissa bits)\n", CHAR_BIT * sizeof
(float), FLT_MANT_DIG); const float x1 = x0;
const float x2 = nextafterf(x1, FLT_MAX);
foo(x1, x2);
}
putchar('\n');
{
printf("double (%zu bits, %d mantissa bits)\n", CHAR_BIT * sizeof
(double), DBL_MANT_DIG); const double x1 = x0;
const double x2 = nextafter(x1, FLT_MAX);
foo(x1, x2);
}
putchar('\n');
{
printf("long double (%zu bits, %d mantissa bits)\n", CHAR_BIT *
sizeof (long double), LDBL_MANT_DIG); const long double x1 = x0;
const long double x2 = nextafterl(x1, FLT_MAX);
foo(x1, x2);
}
}

On Sun, 17 Mar 2024 11:06:21 +0200
Michael S <already5chosen@yahoo.com> wrote:

> On Sat, 16 Mar 2024 16:19:11 -0700
> Keith Thompson <Keith.S.Thompson+u@gmail.com> wrote:
>
> > mitchalsup@aol.com (MitchAlsup1) writes:
> > > Michael S wrote:
> > >> On Sat, 16 Mar 2024 01:16:25 +0000
> > >> mitchalsup@aol.com (MitchAlsup1) wrote:
> > >>> Keith Thompson wrote:
> > >>> > I can see how computing sin(x) with high precision for
> > >>> > "reasonable" values of x would be useful, but does any of that
> > >>> > benefit from being able to compute sin(2^53) accurately?
> > >>> Because accurate argument reduction reduces the burden on the
> > >>> programmer to remain within his sandbox.
> > >
> > >> Not really.
> > >
> > > Say you are a programmer and you receive a value like 2^53 from an
> > > Input read and you wan the most accurate possible SIN( of that ).
> > >
> >
> > I can't think of a scenario where that would be useful (other than
> > just doing it for the sake of doing it).
> >
> > If 2^53 represents a physical quantity, how likely is the actual
> > value to be known within ±π (+/i pi for those who prefer ASCII)?
> >
> > If you can get better precision without too much extra cost, that's
> > great. I don't know enough to have an opinion about what the best
> > tradeoff is, but I presume it's going to be different depending on
> > the application.
> >
> > Here's a C program that shows how precise sin(2^53) can be for types
> > float, double, and long double (I used gcc and glibc). The
> > nextafter functions are used to compute the nearest representable
> > number. For long double, the value of sin() changes by about 1
> > part in 1600, which seems decent, but it's not nearly as precise as
> > for values around 1.0. For float and double, the imprecision of the
> > argument is enough to make the result practically meaningless.
> >

<snip>

> >
>
> As written, your example does not emphasize that the problem has
> nothing to do with implementation of sinX() library routine.
> It's best illustrated by followup conversation with bart, IMHO 100%
> O.T.
> To make the point more clear I'd rather change it to following form:
>
<snip>
>

BTW, personally I'd prefer to illustrate futility of Payne-Hanek or
equivalent reduction with following example:

#include <math.h>
#include <stdio.h>
#include <limits.h>
#include <float.h>

void foo(long double x)
{ const double y = (double)sinl(x);
printf("%.20Le %25.16La %.17f\n", x, x, y);
}

int main(void) {
const float a0 = 0x1p56;
const float b0 = 11;
{
printf("%-12s (%3zu bits, %d mantissa bits) ",
"float", CHAR_BIT * sizeof (float), FLT_MANT_DIG);
const float a = a0;
const float b = b0;
const float x = a / b;
foo(x);
}
{
printf("%-12s (%3zu bits, %d mantissa bits) ",
"double", CHAR_BIT * sizeof (double), DBL_MANT_DIG);
const double a = a0;
const double b = b0;
const double x = a / b;
foo(x);
}
{
printf("%-12s (%3zu bits, %d mantissa bits) ",
"long double", CHAR_BIT * sizeof (long double), LDBL_MANT_DIG);
const long double a = a0;
const long double b = b0;
const long double x = a / b;
foo(x);
}
}

On 17/03/2024 03:57, Keith Thompson wrote:
> bart <bc@freeuk.com> writes:
>> On 17/03/2024 01:38, Keith Thompson wrote:
>>> bart <bc@freeuk.com> writes:
>>>> On 16/03/2024 23:19, Keith Thompson wrote:
>>>>> mitchalsup@aol.com (MitchAlsup1) writes:
>>>>
>>>>>> Say you are a programmer and you receive a value like 2^53 from an
>>>>>> Input read and you wan the most accurate possible SIN( of that ).
>>>>> I can't think of a scenario where that would be useful (other than
>>>>> just
>>>>> doing it for the sake of doing it).
>>>>> If 2^53 represents a physical quantity, how likely is the actual
>>>>> value
>>>>> to be known within ±π (+/i pi for those who prefer ASCII)?
>>>>> If you can get better precision without too much extra cost, that's
>>>>> great. I don't know enough to have an opinion about what the best
>>>>> tradeoff is, but I presume it's going to be different depending on the
>>>>> application.
>>>>> Here's a C program that shows how precise sin(2^53) can be for types
>>>>> float, double, and long double (I used gcc and glibc). The nextafter
>>>>> functions are used to compute the nearest representable number. For
>>>>> long double, the value of sin() changes by about 1 part in 1600, which
>>>>> seems decent, but it's not nearly as precise as for values around 1.0.
>>>>> For float and double, the imprecision of the argument is enough to make
>>>>> the result practically meaningless.
>>>>>
>>>>> ...
>>>>> Output:
>>>>> float (32 bits, 24 mantissa bits)
>>>>> 9007199254740992.00000000 -0.84892595
>>>>> 9007200328482816.00000000 -0.34159181
>>>>> double (64 bits, 53 mantissa bits)
>>>>> 9007199254740992.00000000 -0.84892596
>>>>> 9007199254740994.00000000 -0.12729655
>>>>> long double (128 bits, 64 mantissa bits)
>>>>> 9007199254740992.00000000 -0.84892596
>>>>> 9007199254740992.00097656 -0.84944168
>>>
>>>>
>>>> Is this output supposed to be different between gcc -O0 and gcc -O3?
>>> No, why do you ask?
>>> On my system, the output is identical with gcc and clang, and tcc at
>>> all
>>> optimization levels, with both glibc and musl. I do get slightly
>>> different results for long double on Cygwin vs. Ubuntu (Cygwin uses
>>> newlib, not glibc).
>>
>>
>> Because I get the results shown below. Even if the library is
>> different from yours, I would have assumed matching results.
>>
>>
>> ------------------------------
>>
>> C:\c>gcc c.c -O3
>>
>> C:\c>a
>> float (32 bits, 24 mantissa bits)
>> 9007199254740992.00000000 -0.84892595
>> 9007200328482816.00000000 -0.34159181
>>
>> double (64 bits, 53 mantissa bits)
>> 9007199254740992.00000000 -0.84892596
>> 9007199254740994.00000000 -0.12729655
>>
>> long double (128 bits, 64 mantissa bits)
>> 9007199254740992.00000000 -0.84892596
>> 9007199254740992.00097656 -0.84944168
>>
>> C:\c>gcc c.c -O0
>>
>> C:\c>a
>> float (32 bits, 24 mantissa bits)
>> 9007199254740992.00000000 -0.84893209
>> 9007200328482816.00000000 -0.34160271
>>
>> double (64 bits, 53 mantissa bits)
>> 9007199254740992.00000000 -0.84893209
>> 9007199254740994.00000000 -0.12728505
>>
>> long double (128 bits, 64 mantissa bits)
>> 9007199254740992.00000000 -0.84893209
>> 9007199254740992.00097656 -0.84944780
>>
>> C:\c>gcc --version
>> gcc (tdm64-1) 10.3.0
>
> I see similar results (possibly identical, I haven't checked) with
> i686-w64-mingw32-gcc on Cygwin. I think it uses the same, or a similar,
> runtime library as the tdm64 implementation you're using.
>
> Comparing the generated assembly, at -O0 it generates a call to "sin",
> while at -O1 and above it replaces the expression with a compile-time
> constant. Apparently that compile-time constant matches the run-time
> computed value for glibc, but not for whatever tdm64 uses.
>
> I don't *think* it's a conformance issue as long as the precision is
> within the standard's (fairly vague) requirements. When I modify the
> program to start with sin(1.0) rather than sin(2**53), the output is
> identical at all optimization levels.
>
> Apparently glibc and the library used by tgm64 behave differently when
> computing the sin of very large values. The result for sin(2**53) in
> long double differs by about one part in 2**17 between -O0 and -O1, or
> between gcc/glibc and tgm64.
>

That is interesting. It might not be a conformance issue, but it is a
divergence from the intention of the compiler. gcc makes use of a very
comprehensive library that emulates a whole range of floating point
hardware implementations to try to get bit-perfect compile-time versions
of run-time calculations, precisely so that it can do compile-time
calculations that fully match what you would get at run-time, regardless
of compilation flags and optimisation.

So if the program is producing different results for run-time and
compile-time versions of "sin", then there is a mismatch between the
configuration of the compiler and the run-time library it is using. I
don't know enough about the configuration options for gcc (there are a
/lot/ of them), but I know it can be configured to assume you are using
glibc or a glibc-compatible library. Perhaps the mingw build here is
configured with this option, but is actually using a run-time library
with a slightly implementation of these functions (such as an MS DLL C
library). If that is the case, then this is a mistake in the
configuration of gcc here, and I would think it could be reported back
to the mingw people as a bug.

On 3/15/2024 3:48 PM, Keith Thompson wrote:
> "Chris M. Thomasson" <chris.m.thomasson.1@gmail.com> writes:
>> On 3/15/2024 2:26 PM, Chris M. Thomasson wrote:
>> [...]
>>> Notice how the numbers are small? However, if we zoom in enough,
>>> then 2^53 might become an issue..
>> [...]
>>
>> Argh! I mean if we zoom in enough, then only 2^53 of precision might
>> become an issue... We can go in really deep:
>>
>> A deep zoom might be 10^275. Somebody did one to the 10^4141! I feel
>> like including a youtube link here, but I will refrain out of respect
>> for you.
>
> I don't object to YouTube links as long as they're reasonably relevant
> and there's enough information for readers to decide whether it's worth
> their time to open them.
>

Fair enough. :^)

On 3/17/2024 1:41 PM, Chris M. Thomasson wrote:
> On 3/15/2024 3:48 PM, Keith Thompson wrote:
>> "Chris M. Thomasson" <chris.m.thomasson.1@gmail.com> writes:
>>> On 3/15/2024 2:26 PM, Chris M. Thomasson wrote:
>>> [...]
>>>> Notice how the numbers are small? However, if we zoom in enough,
>>>> then 2^53 might become an issue..
>>> [...]
>>>
>>> Argh! I mean if we zoom in enough, then only 2^53 of precision might
>>> become an issue... We can go in really deep:
>>>
>>> A deep zoom might be 10^275. Somebody did one to the 10^4141! I feel
>>> like including a youtube link here, but I will refrain out of respect
>>> for you.
>>
>> I don't object to YouTube links as long as they're reasonably relevant
>> and there's enough information for readers to decide whether it's worth
>> their time to open them.
>>
>
> Fair enough. :^)

Fwiw, another algorithm of mine that would benefit from high
precision... The fun part is that a point can be placed anywhere in the
field. So, it can be a big one. Depending on the users placement, so to
speak:

https://gallery.bridgesmathart.org/exhibitions/2024-joint-mathematics-meetings/chris-m-thomasson

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

> At which cost in tables sizes?

My impression was that it wasn't costly in that respect, but since my
recollection seems to be off on the performance, maybe it's off here
as well.

> The critical point here is definition of what considered exact. If
> 'exact' is measured only on y side of y=foo(x), disregarding
> possible imprecision on the x side then you are very likely to end up
> with results that are slower to calculate, but not at all more useful
> from point of view of engineer or physicist. Exactly like Payne-Hanek
> or Mitch's equivalent of Payne-Hanek.

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.

Stefan

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×

>> At which cost in tables sizes?

Making transcendental faster is always a tradeoff between table size
and speed. SIN() COS() can use 10-term polynomials when the reduced
argument is -¼pi..+¼pi, on the other hand, when the reduced argument
is ±0.008 a 3 term polynomial is just as accurate but you need 128×3
tables.

> My impression was that it wasn't costly in that respect, but since my
> recollection seems to be off on the performance, maybe it's off here
> as well.

ITANIC did rather well, here, because it had 2 FMAC units and could use
both for the polynomials.

>> The critical point here is definition of what considered exact. If
>> 'exact' is measured only on y side of y=foo(x), disregarding
>> possible imprecision on the x side then you are very likely to end up
>> with results that are slower to calculate, but not at all more useful
>> from point of view of engineer or physicist. Exactly like Payne-Hanek
>> or Mitch's equivalent of Payne-Hanek.

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

> Stefan

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