On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
Lawrence D'Oliveiro <ldo@nz.invalid> wrote:

> On Thu, 22 Feb 2024 20:16:25 -0000 (UTC), Steven G. Kargl wrote:
>
> > Apparently, you missed the part about argument reduction. For
> > sinpi(x), it is fairly easy to reduce x = n + r with n an integer
> > and r in [0,1). For the extended interval, x in [0,2^23], there are
> > roughly 2^23 values with r = 0.5; and thus, sinpi(x) = 1 exactly.
> > There are no such values for sin(x), and argument reduction for
> > sin(x) is much more involved.
>
> You are working with approximations anyway. That those approximations
> happen to exactly equal some random value seems irrelevant.
>
> > As to real-world use, how about any physical phenomena where one is
> > interest in resonance frequencies of the system. For a simple
> > example see https://www.feynmanlectures.caltech.edu/I_49.html where
> > one might write f(x) = sin(kx) = sin(pi * (2*x/L)) with L a length
> > of say a clamped string.
>
> I don’t see a formula anything like that anywhere on that page. On
> the other hand, I see lots of use of “ω” symbols, representing
> frequencies in radians per second.
>
> > There are also uses with computing other functions, e.g., the true
> > gamma function via the Euler reflection formula.
> >
> > gamma(x) * gamma(1 - x) = pi / sin(pi * x) = pi / sinpi(x)
>
> π radians = half a circle. Are there any other examples of the
> usefulness of half-circles as an angle unit? As opposed to the dozens
> or hundreds of examples of the usefulness of radians as an angle unit?

In digital signal processing circle-based units are pretty much always
more natural than radians.
For specific case of 1/2 circle, I can't see where it can be used
directly.
From algorithmic perspective, full circle looks like the most obvious
choice.
From [binary floating point] numerical properties perspective, 1/8th of
the circle (==pi/4 radians = 45 degrees) is probably the best option
for a library routine, because for sin() its derivative at 0 is closest
to 1 among all powers of two which means that loss of precision
near 0 is very similar for input and for output. But this advantage
does not sound like particularly big deal.

Michael S wrote:

> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:

>> On Thu, 22 Feb 2024 20:16:25 -0000 (UTC), Steven G. Kargl wrote:
>>
>> > Apparently, you missed the part about argument reduction. For
>> > sinpi(x), it is fairly easy to reduce x = n + r with n an integer
>> > and r in [0,1).

You are better off with r in [-½,+½]

>> For the extended interval, x in [0,2^23], there are
>> > roughly 2^23 values with r = 0.5; and thus, sinpi(x) = 1 exactly.
>> > There are no such values for sin(x), and argument reduction for
>> > sin(x) is much more involved.

I have a patent on doing argument reduction for sin(x) in 4 cycles...that
is as good as Payne-Haneok argument reduction over -infinity..+infinity
>>
>> You are working with approximations anyway. That those approximations
>> happen to exactly equal some random value seems irrelevant.
>>
>> > As to real-world use, how about any physical phenomena where one is
>> > interest in resonance frequencies of the system. For a simple
>> > example see https://www.feynmanlectures.caltech.edu/I_49.html where
>> > one might write f(x) = sin(kx) = sin(pi * (2*x/L)) with L a length
>> > of say a clamped string.
>>
>> I don’t see a formula anything like that anywhere on that page. On
>> the other hand, I see lots of use of “ω” symbols, representing
>> frequencies in radians per second.
>>
>> > There are also uses with computing other functions, e.g., the true
>> > gamma function via the Euler reflection formula.
>> >
>> > gamma(x) * gamma(1 - x) = pi / sin(pi * x) = pi / sinpi(x)

I would not have expected for the magnitude of gamma(x) * gamma(1 - x)
to be always greater than than pi....but there it is.

>>
>> π radians = half a circle. Are there any other examples of the
>> usefulness of half-circles as an angle unit? As opposed to the dozens
>> or hundreds of examples of the usefulness of radians as an angle unit?

> In digital signal processing circle-based units are pretty much always
> more natural than radians.
> For specific case of 1/2 circle, I can't see where it can be used
> directly.
> From algorithmic perspective, full circle looks like the most obvious
> choice.
> From [binary floating point] numerical properties perspective, 1/8th of
> the circle (==pi/4 radians = 45 degrees) is probably the best option
> for a library routine, because for sin() its derivative at 0 is closest
> to 1 among all powers of two which means that loss of precision
> near 0 is very similar for input and for output. But this advantage
> does not sound like particularly big deal.

I should remind that my patents include methods for calculating sin()
and cos() in 18 cycles, sinpi() and cospi() in 14 cycles while achieving
an accuracy n the 0.5002-bits of error. All of this is in HW, and all
carried out in precision wider than IEEE 754 with a bunch of HW tricks
not available to SW.

At this performance level, let the algorithms wanting pi based trigonometry
have it and those wanting unit based trigonometry have it too. Let program-
mers use what is easiest for them.

On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:

> Michael S wrote:
>
>> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
>> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
>
>>> On Thu, 22 Feb 2024 20:16:25 -0000 (UTC), Steven G. Kargl wrote:
>>>
>>> > Apparently, you missed the part about argument reduction. For
>>> > sinpi(x), it is fairly easy to reduce x = n + r with n an integer
>>> > and r in [0,1).
>
> You are better off with r in [-½,+½]

Partial implementation details for single precision C code (trying
to keep this on topic for c.l.c) are at
https://cgit.freebsd.org/src/tree/lib/msun/src/s_sinpif.c

One uses sinpi(-|x|) = -sinpi(|x|), which natually gets one to
r = x - floor(x) in [0,1). One then uses kernels for sin() and
cos() to do the actual computation.

r in [0,0.25) is ksin(r).
r in [0.25,0.5) is kcos(0.5-r).
r in [0.5,0.75) is kcos(r - 0.5).
r in [0.75,1) is ksin(1 - r).

>>> For the extended interval, x in [0,2^23], there are
>>> > roughly 2^23 values with r = 0.5; and thus, sinpi(x) = 1 exactly.
>>> > There are no such values for sin(x), and argument reduction for
>>> > sin(x) is much more involved.
>
> I have a patent on doing argument reduction for sin(x) in 4 cycles...that
> is as good as Payne-Haneok argument reduction over -infinity..+infinity

By a strange coincidence, I have the Payn-Hanek paper on the
top of stack of papers on my desk. Still need to internalize
the details.

(much trimmed for brevity)

>> In digital signal processing circle-based units are pretty much always
>> more natural than radians.
>> For specific case of 1/2 circle, I can't see where it can be used
>> directly.
>> From algorithmic perspective, full circle looks like the most obvious
>> choice.
>> From [binary floating point] numerical properties perspective, 1/8th of
>> the circle (==pi/4 radians = 45 degrees) is probably the best option
>> for a library routine, because for sin() its derivative at 0 is closest
>> to 1 among all powers of two which means that loss of precision
>> near 0 is very similar for input and for output. But this advantage
>> does not sound like particularly big deal.
>
> I should remind that my patents include methods for calculating sin()
> and cos() in 18 cycles, sinpi() and cospi() in 14 cycles while achieving
> an accuracy n the 0.5002-bits of error. All of this is in HW, and all
> carried out in precision wider than IEEE 754 with a bunch of HW tricks
> not available to SW.
>
> At this performance level, let the algorithms wanting pi based trigonometry
> have it and those wanting unit based trigonometry have it too. Let program-
> mers use what is easiest for them.

Agreed a programmer should use what is required by the problem
that they are solving. I'll note that SW implementations have
their sets of tricks (e.g., use of double-double arithmetic to
achieve double precision).

--
steve

On 2/22/2024 4:13 PM, Steven G. Kargl wrote:
> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>
>> Michael S wrote:
>>
>>> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
>>> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
>>
>>>> On Thu, 22 Feb 2024 20:16:25 -0000 (UTC), Steven G. Kargl wrote:
>>>>
>>>>> Apparently, you missed the part about argument reduction. For
>>>>> sinpi(x), it is fairly easy to reduce x = n + r with n an integer
>>>>> and r in [0,1).
>>
>> You are better off with r in [-½,+½]
>
> Partial implementation details for single precision C code (trying
> to keep this on topic for c.l.c) are at
> https://cgit.freebsd.org/src/tree/lib/msun/src/s_sinpif.c
>
> One uses sinpi(-|x|) = -sinpi(|x|), which natually gets one to
> r = x - floor(x) in [0,1). One then uses kernels for sin() and
> cos() to do the actual computation.
>
> r in [0,0.25) is ksin(r).
> r in [0.25,0.5) is kcos(0.5-r).
> r in [0.5,0.75) is kcos(r - 0.5).
> r in [0.75,1) is ksin(1 - r).
>
>>>> For the extended interval, x in [0,2^23], there are
>>>>> roughly 2^23 values with r = 0.5; and thus, sinpi(x) = 1 exactly.
>>>>> There are no such values for sin(x), and argument reduction for
>>>>> sin(x) is much more involved.
>>
>> I have a patent on doing argument reduction for sin(x) in 4 cycles...that
>> is as good as Payne-Haneok argument reduction over -infinity..+infinity
>
> By a strange coincidence, I have the Payn-Hanek paper on the
> top of stack of papers on my desk. Still need to internalize
> the details.
>
> (much trimmed for brevity)
>
>
>>> In digital signal processing circle-based units are pretty much always
>>> more natural than radians.
>>> For specific case of 1/2 circle, I can't see where it can be used
>>> directly.
>>> From algorithmic perspective, full circle looks like the most obvious
>>> choice.
>>> From [binary floating point] numerical properties perspective, 1/8th of
>>> the circle (==pi/4 radians = 45 degrees) is probably the best option
>>> for a library routine, because for sin() its derivative at 0 is closest
>>> to 1 among all powers of two which means that loss of precision
>>> near 0 is very similar for input and for output. But this advantage
>>> does not sound like particularly big deal.
>>
>> I should remind that my patents include methods for calculating sin()
>> and cos() in 18 cycles, sinpi() and cospi() in 14 cycles while achieving
>> an accuracy n the 0.5002-bits of error. All of this is in HW, and all
>> carried out in precision wider than IEEE 754 with a bunch of HW tricks
>> not available to SW.
>>
>> At this performance level, let the algorithms wanting pi based trigonometry
>> have it and those wanting unit based trigonometry have it too. Let program-
>> mers use what is easiest for them.
>
> Agreed a programmer should use what is required by the problem
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Indeed! :^)

> that they are solving. I'll note that SW implementations have
> their sets of tricks (e.g., use of double-double arithmetic to
> achieve double precision).
>

Use some form of arbitrary precision if you want to get to:

A 10^4141 zoom mset zoom!? It must took many hours to render... ;^)

https://youtu.be/Xjy_HSUujaw

On 2/22/2024 4:49 PM, Chris M. Thomasson wrote:
[...]> Use some form of arbitrary precision if you want to get to:
>
> A 10^4141 zoom mset zoom!? It must took many hours to render... ;^)
>
> https://youtu.be/Xjy_HSUujaw

A 10^275 zoom:

https://youtu.be/0jGaio87u3A

Steven G. Kargl wrote:

> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:

>> Michael S wrote:
>>
>>> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
>>> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
>>
>>>> On Thu, 22 Feb 2024 20:16:25 -0000 (UTC), Steven G. Kargl wrote:
>>>>
>>>> > Apparently, you missed the part about argument reduction. For
>>>> > sinpi(x), it is fairly easy to reduce x = n + r with n an integer
>>>> > and r in [0,1).
>>
>> You are better off with r in [-½,+½]

> Partial implementation details for single precision C code (trying
> to keep this on topic for c.l.c) are at
> https://cgit.freebsd.org/src/tree/lib/msun/src/s_sinpif.c

> One uses sinpi(-|x|) = -sinpi(|x|), which natually gets one to
> r = x - floor(x) in [0,1). One then uses kernels for sin() and
> cos() to do the actual computation.

> r in [0,0.25) is ksin(r).
> r in [0.25,0.5) is kcos(0.5-r).
> r in [0.5,0.75) is kcos(r - 0.5).
> r in [0.75,1) is ksin(1 - r).

By "better off" I mean you have ½-bit of greater reduced argument precision
when the reduced argument is centered on zero.

>>>> For the extended interval, x in [0,2^23], there are
>>>> > roughly 2^23 values with r = 0.5; and thus, sinpi(x) = 1 exactly.
>>>> > There are no such values for sin(x), and argument reduction for
>>>> > sin(x) is much more involved.
>>
>> I have a patent on doing argument reduction for sin(x) in 4 cycles...that
>> is as good as Payne-Haneok argument reduction over -infinity..+infinity

> By a strange coincidence, I have the Payn-Hanek paper on the
> top of stack of papers on my desk. Still need to internalize
> the details.

Here is the version found in the SUN transcendentals::

static void ReduceFull(double *xp, int *a, double x)
{ Double X = { x };
int ec = X.s.exponent - (1023+33);
int k = (ec + 26) * (607*4) >> 16;
int m = 27*k - ec;
int offset = m >> 3;
x *= 0x1p-400;
double xDekker = x * (0x1p27 + 1);
double x0 = xDekker - (xDekker - x);
double x1 = x - x0;
const double *p0 = &TwoOverPiWithOffset[offset][k]; // 180 DP FP numbers
const double fp0 = p0[0];
const double fp1 = p0[1];
const double fp2 = p0[2];
const double fp3 = p0[3];
const double f0 = x1 * fp0 + fp1 * x0;
double f = x1 * fp1 + fp2 * x0;
const double fi = f0 + f;
static const double IntegerBias = 0x1.8p52;
Double Fi = { fi + IntegerBias };
*a = Fi.s.significand2;
double fint = Fi.d - IntegerBias;
const double fp4 = p0[4];
const double fp5 = p0[5];
const double fp6 = p0[6];
f = f0 - fint + f;
f += x1 * fp2 + fp3 * x0;
f += x1 * fp3 + fp4 * x0;
f += x1 * fp4 + fp5 * x0;
f += x1 * fp5 + fp6 * x0;
*xp = f * 0x3.243F6A8885A3p-1;
}

Which I can do in 4 cycles in HW !! {{I happen to have this on hand to explain how
I can do this in HW in 4 cycles......}}

I should also note that the conversion of f back into *xp looses ¼-bit of precision.

> (much trimmed for brevity)

>>> In digital signal processing circle-based units are pretty much always
>>> more natural than radians.
>>> For specific case of 1/2 circle, I can't see where it can be used
>>> directly.
>>> From algorithmic perspective, full circle looks like the most obvious
>>> choice.
>>> From [binary floating point] numerical properties perspective, 1/8th of
>>> the circle (==pi/4 radians = 45 degrees) is probably the best option
>>> for a library routine, because for sin() its derivative at 0 is closest
>>> to 1 among all powers of two which means that loss of precision
>>> near 0 is very similar for input and for output. But this advantage
>>> does not sound like particularly big deal.
>>
>> I should remind that my patents include methods for calculating sin()
>> and cos() in 18 cycles, sinpi() and cospi() in 14 cycles while achieving
>> an accuracy n the 0.5002-bits of error. All of this is in HW, and all
>> carried out in precision wider than IEEE 754 with a bunch of HW tricks
>> not available to SW.
>>
>> At this performance level, let the algorithms wanting pi based trigonometry
>> have it and those wanting unit based trigonometry have it too. Let program-
>> mers use what is easiest for them.

> Agreed a programmer should use what is required by the problem
> that they are solving. I'll note that SW implementations have
> their sets of tricks (e.g., use of double-double arithmetic to
> achieve double precision).

To get near IEEE desired precision, one HAS TO use more than 754 precision.

Chris M. Thomasson wrote:

> On 2/22/2024 4:49 PM, Chris M. Thomasson wrote:
> [...]> Use some form of arbitrary precision if you want to get to:
>>
>> A 10^4141 zoom mset zoom!? It must took many hours to render... ;^)
>>
>> https://youtu.be/Xjy_HSUujaw

> A 10^275 zoom:

> https://youtu.be/0jGaio87u3A

Argument reduction zoom::

COS( 6381956970095103×2^797) = -4.68716592425462761112×10-19

You can argue all you want about whether this has any significance, but
it happens to be the value closest to that radian argument when performed
with infinite precision.

Michael S wrote:
> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
>>
>> π radians = half a circle. Are there any other examples of the
>> usefulness of half-circles as an angle unit? As opposed to the dozens
>> or hundreds of examples of the usefulness of radians as an angle unit?
>
> In digital signal processing circle-based units are pretty much always
> more natural than radians.
> For specific case of 1/2 circle, I can't see where it can be used
> directly.
> From algorithmic perspective, full circle looks like the most obvious
> choice.
> From [binary floating point] numerical properties perspective, 1/8th of
> the circle (==pi/4 radians = 45 degrees) is probably the best option
> for a library routine, because for sin() its derivative at 0 is closest
> to 1 among all powers of two which means that loss of precision
> near 0 is very similar for input and for output. But this advantage
> does not sound like particularly big deal.

ieee754 defines sinpi() and siblings, but imho it really doesn't matter
if you use circles, half-circles (i.e. sinpi) or some other binary
fraction of a circle: Argument reduction for huge inputs are just as
easy, you might just have to multiply by the corresponding power of two
(i.e. adjust the exponent) before extracting the fractional term.

For sinpi(x) I could do it like this:

if (abs(x) >= two_to_52nd_power) error("Zero significant bits.");
ix = int(x);
x_reduced = x - (double) (ix & ~1);
if (x_reduced < 0.0) x_reduced += 2.0;

but it is probably better to return a value in the [-1.0 .. 1.0> range?

Terje

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

MitchAlsup1 wrote:
> Steven G. Kargl wrote:
>> Agreed a programmer should use what is required by the problem
>> that they are solving.  I'll note that SW implementations have
>> their sets of tricks (e.g., use of double-double arithmetic to
>> achieve double precision).
>
> To get near IEEE desired precision, one HAS TO use more than 754 precision.

There are groups who have shown that exactly rounded trancendental
functions are in fact achievable with maybe 3X reduced performance.

There is a suggestion on the table to make that a (probably optional
imho) feature for an upcoming ieee754 revision.

Terje

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

On Fri, 23 Feb 2024 11:01:02 +0100
Terje Mathisen <terje.mathisen@tmsw.no> wrote:

> Michael S wrote:
> > On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
> > Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
> >>
> >> π radians = half a circle. Are there any other examples of the
> >> usefulness of half-circles as an angle unit? As opposed to the
> >> dozens or hundreds of examples of the usefulness of radians as an
> >> angle unit?
> >
> > In digital signal processing circle-based units are pretty much
> > always more natural than radians.
> > For specific case of 1/2 circle, I can't see where it can be used
> > directly.
> > From algorithmic perspective, full circle looks like the most
> > obvious choice.
> > From [binary floating point] numerical properties perspective,
> > 1/8th of the circle (==pi/4 radians = 45 degrees) is probably the
> > best option for a library routine, because for sin() its derivative
> > at 0 is closest to 1 among all powers of two which means that loss
> > of precision near 0 is very similar for input and for output. But
> > this advantage does not sound like particularly big deal.
>
> ieee754 defines sinpi() and siblings, but imho it really doesn't
> matter if you use circles, half-circles (i.e. sinpi) or some other
> binary fraction of a circle: Argument reduction for huge inputs are
> just as easy, you might just have to multiply by the corresponding
> power of two (i.e. adjust the exponent) before extracting the
> fractional term.
>
> For sinpi(x) I could do it like this:
>
> if (abs(x) >= two_to_52nd_power) error("Zero significant bits.");
> ix = int(x);
> x_reduced = x - (double) (ix & ~1);
> if (x_reduced < 0.0) x_reduced += 2.0;
>
> but it is probably better to return a value in the [-1.0 .. 1.0>
> range?
>
> Terje
>

Both you and Mitch look at it from wrong perspective.
When we define a library API, an ease of implementation of the library
function should be pretty low on our priority scale. As long as
reasonably precise calculation is theoretically possible, we should
give credit to intelligence of implementor, that's all.
The real concern of API designer should be with avoidance of loss of
precision in preparation of inputs and use of outputs.
In specific case of y=sin2pi(x), it is x that is more problematic,
because near 0 it starts to lose precision 3 octaves before y. In
subnormal range we lose ~2.5 bits of precision in preparation of the
argument. An implementation, no matter how good, can't recover what's
already lost.
sinpi() is slightly better, but only slightly, not enough better to
justify less natural semantics.

My yesterday suggestion is a step in right direction, but today I think
that it is not sufficiently radical step. In specific case of sin/cos
there is no good reason to match loss of precision on input with loss of
precision.
<Thinking out load>
There are advantages in matched loss of precision for input and for
output when both input and output occupy full range of real numbers
(ignoring sign). Output of sin/cos does not occupy o full range. But
for tan() it does. So, may be, there is a one good reason for matching
input with output for sin/cos - consistency with tan.
</Thinking out load>
So, ignoring tan(), what is really an optimal input scaling for sin/cos
inputs? Today I think that it is a scaling in which full circle
corresponds to 2**64. With such scaling you never lose any input
precision to subnoramals before precision of the result is lost
completely.
Now, one could ask "Why 2**64, why not 2**56 that has the same
property?". My answer is "Because 2**64 is a scaling that is most
convenient for preparation of trig arguments in fix-point [on 64-bit
computer]." I.e. apart from being a good scaling for avoidance of loss
of precision in tiny range, it happens to be the best scaling for
interoperability with fixed point.
That is my answer today. Will it hold tomorrow? Tomorrow will tell.

On Fri, 23 Feb 2024 00:13:02 -0000 (UTC)
"Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> wrote:

> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>
> >
> > I have a patent on doing argument reduction for sin(x) in 4
> > cycles...that is as good as Payne-Haneok argument reduction over
> > -infinity..+infinity
>
> By a strange coincidence, I have the Payn-Hanek paper on the
> top of stack of papers on my desk. Still need to internalize
> the details.
>

Payne-Hanek is a right answer to a wrong question.
In science/engineering floating-point number x represents a range
[x-ULP:x+ULP] with hopefully much higher probability of being in range
[x-ULP/2:x+ULP/2]. However in this reduced range probability is flat,
an "exact" x is no worse and no better than any other point within this
range.
Which means that when we look for y=sin(x) then for
scientific/engineering purposes any answer in range*
[sin(x-ULP/2):sin(x+ULP/2)] is as "right" as any other answer. The
right answer to the question "What is a sin() of IEEE-754 binary64
number 1e17?" is "Anything in range [-1:1]". The wise and useful answer
to the same question is "You're holding it wrong".

* in paragraph Sin(x) designates result calculated with infinite
precision.

Michael S wrote:

> On Fri, 23 Feb 2024 00:13:02 -0000 (UTC)
> "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> wrote:

>> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>>
>> >
>> > I have a patent on doing argument reduction for sin(x) in 4
>> > cycles...that is as good as Payne-Haneok argument reduction over
>> > -infinity..+infinity
>>
>> By a strange coincidence, I have the Payn-Hanek paper on the
>> top of stack of papers on my desk. Still need to internalize
>> the details.
>>

> Payne-Hanek is a right answer to a wrong question.

If someone looking for fast sin(x) you are correct, for a person looking
for the best possible (correctly rounded) result, you are not. In any
event I have driven that monstrosity down from 50-odd cycles to 4--
making it at least palatable.

> In science/engineering floating-point number x represents a range
> [x-ULP:x+ULP] with hopefully much higher probability of being in range
> [x-ULP/2:x+ULP/2]. However in this reduced range probability is flat,
> an "exact" x is no worse and no better than any other point within this
> range.
> Which means that when we look for y=sin(x) then for
> scientific/engineering purposes any answer in range*
> [sin(x-ULP/2):sin(x+ULP/2)] is as "right" as any other answer. The
> right answer to the question "What is a sin() of IEEE-754 binary64
> number 1e17?" is "Anything in range [-1:1]". The wise and useful answer
> to the same question is "You're holding it wrong".

There are verification tests that know the exact value of both x and sin(x).
For these uses you argument is not valid--however I am willing to grant
that outside of verification, and in actual use, your argument holds as
well ad the noise in the data provides.

> * in paragraph Sin(x) designates result calculated with infinite
> precision.

On 2/22/2024 6:42 PM, MitchAlsup1 wrote:
> Chris M. Thomasson wrote:
>
>> On 2/22/2024 4:49 PM, Chris M. Thomasson wrote:
>> [...]> Use some form of arbitrary precision if you want to get to:
>>>
>>> A 10^4141 zoom mset zoom!? It must took many hours to render... ;^)
>>>
>>> https://youtu.be/Xjy_HSUujaw
>
>> A 10^275 zoom:
>
>> https://youtu.be/0jGaio87u3A
>
> Argument reduction zoom::
>
> COS( 6381956970095103×2^797) = -4.68716592425462761112×10-19
>
> You can argue all you want about whether this has any significance, but
> it happens to be the value closest to that radian argument when performed
> with infinite precision.

Well, I use a lot of trig wrt my work. One example is getting at the
roots of complex numbers, polar form to rectangular form, ect... So
zooming in "deep" on say, something like my MultiJulia IFS work is going
to require arbitrary precision trig...

https://paulbourke.net/fractals/multijulia

Btw, my friend Paul took the time to create this nice write up and all
of the renders of my MultiJulia IFS experiment.

On Fri, 23 Feb 2024 20:02:00 +0000, MitchAlsup1 wrote:

> Michael S wrote:
>
>> On Fri, 23 Feb 2024 00:13:02 -0000 (UTC)
>> "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> wrote:
>
>>> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>>>
>>>
>>> > I have a patent on doing argument reduction for sin(x) in 4
>>> > cycles...that is as good as Payne-Haneok argument reduction over
>>> > -infinity..+infinity
>>>
>>> By a strange coincidence, I have the Payn-Hanek paper on the top of
>>> stack of papers on my desk. Still need to internalize the details.
>>>
>>>
>> Payne-Hanek is a right answer to a wrong question.
>
> If someone looking for fast sin(x) you are correct, for a person looking
> for the best possible (correctly rounded) result, you are not. In any
> event I have driven that monstrosity down from 50-odd cycles to 4--
> making it at least palatable.
>
>> In science/engineering floating-point number x represents a range
>> [x-ULP:x+ULP] with hopefully much higher probability of being in range
>> [x-ULP/2:x+ULP/2]. However in this reduced range probability is flat,
>> an "exact" x is no worse and no better than any other point within this
>> range.
>> Which means that when we look for y=sin(x) then for
>> scientific/engineering purposes any answer in range*
>> [sin(x-ULP/2):sin(x+ULP/2)] is as "right" as any other answer. The
>> right answer to the question "What is a sin() of IEEE-754 binary64
>> number 1e17?" is "Anything in range [-1:1]". The wise and useful answer
>> to the same question is "You're holding it wrong".
>
> There are verification tests that know the exact value of both x and
> sin(x).
> For these uses you argument is not valid--however I am willing to grant
> that outside of verification, and in actual use, your argument holds as
> well ad the noise in the data provides.

Fortunately, for 32-bit single precision floating point, one can
exhaustively test single-argument functions. For the SW implementation
on FreeBSD, exhaustive testing shows

% tlibm sinpi -fPED -x 0 -X 0x1p23 -s 0
Interval tested for sinpif: [0,8.38861e+06]
1000000 calls, 0.015453 secs, 0.01545 usecs/call
ulp <= 0.5: 99.842% 1253631604 | 99.842% 1253631604
0.5 < ulp < 0.6: 0.158% 1989420 | 100.000% 1255621024
Max ulp: 0.583666 at 6.07950642e-05

% tlibm sin -fPED -x 0 -X 0x1p23 -s 0
Interval tested for sinf: [0,8.38861e+06]
1000000 calls, 0.019520 secs, 0.01952 usecs/call
ulp <= 0.5: 99.995% 1249842491 | 99.995% 1249842491
0.5 < ulp < 0.6: 0.005% 60102 | 100.000% 1249902593
Max ulp: 0.500900 at 3.68277281e+05

The speed test is for 1M values evenly distributed over
the interval. The difference in speed for sinpi vs sin
is due to the argument reduction.

Note 1: libm built with clang/llvm with only -O2 on a
AMD Phenom II X2 560 Processor (3300.14-MHz K8-class CPU).

Note 2: subnormal values for x are test not included in
the count as subnormal are tested with a different approach.

--
steve

Steven G. Kargl wrote:

> On Fri, 23 Feb 2024 20:02:00 +0000, MitchAlsup1 wrote:

>> Michael S wrote:
>>
>>> On Fri, 23 Feb 2024 00:13:02 -0000 (UTC)
>>> "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> wrote:
>>
>>>> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>>>>
>>>>
>>>> > I have a patent on doing argument reduction for sin(x) in 4
>>>> > cycles...that is as good as Payne-Haneok argument reduction over
>>>> > -infinity..+infinity
>>>>
>>>> By a strange coincidence, I have the Payn-Hanek paper on the top of
>>>> stack of papers on my desk. Still need to internalize the details.
>>>>
>>>>
>>> Payne-Hanek is a right answer to a wrong question.
>>
>> If someone looking for fast sin(x) you are correct, for a person looking
>> for the best possible (correctly rounded) result, you are not. In any
>> event I have driven that monstrosity down from 50-odd cycles to 4--
>> making it at least palatable.
>>
>>> In science/engineering floating-point number x represents a range
>>> [x-ULP:x+ULP] with hopefully much higher probability of being in range
>>> [x-ULP/2:x+ULP/2]. However in this reduced range probability is flat,
>>> an "exact" x is no worse and no better than any other point within this
>>> range.
>>> Which means that when we look for y=sin(x) then for
>>> scientific/engineering purposes any answer in range*
>>> [sin(x-ULP/2):sin(x+ULP/2)] is as "right" as any other answer. The
>>> right answer to the question "What is a sin() of IEEE-754 binary64
>>> number 1e17?" is "Anything in range [-1:1]". The wise and useful answer
>>> to the same question is "You're holding it wrong".
>>
>> There are verification tests that know the exact value of both x and
>> sin(x).
>> For these uses you argument is not valid--however I am willing to grant
>> that outside of verification, and in actual use, your argument holds as
>> well ad the noise in the data provides.

> Fortunately, for 32-bit single precision floating point, one can

All of the numerics I have spoken about are DP (64-bit) values.

> exhaustively test single-argument functions. For the SW implementation
> on FreeBSD, exhaustive testing shows

> % tlibm sinpi -fPED -x 0 -X 0x1p23 -s 0
> Interval tested for sinpif: [0,8.38861e+06]
> 1000000 calls, 0.015453 secs, 0.01545 usecs/call
> ulp <= 0.5: 99.842% 1253631604 | 99.842% 1253631604
> 0.5 < ulp < 0.6: 0.158% 1989420 | 100.000% 1255621024
> Max ulp: 0.583666 at 6.07950642e-05

> % tlibm sin -fPED -x 0 -X 0x1p23 -s 0
> Interval tested for sinf: [0,8.38861e+06]
> 1000000 calls, 0.019520 secs, 0.01952 usecs/call
> ulp <= 0.5: 99.995% 1249842491 | 99.995% 1249842491
> 0.5 < ulp < 0.6: 0.005% 60102 | 100.000% 1249902593
> Max ulp: 0.500900 at 3.68277281e+05

I find it interesting that sinpi() has worse numeric error than sin().

I also find it interesting that the highest error is in the easy part
of polynomial evaluation without argument reduction whereas sin() has
its worst result in a region where significant argument reduction has
transpired.

{{Seems like another case of faster poor answers top slower correct
answers.}}

> The speed test is for 1M values evenly distributed over
> the interval. The difference in speed for sinpi vs sin
> is due to the argument reduction.

But result precision suffers nonetheless.

> Note 1: libm built with clang/llvm with only -O2 on a
> AMD Phenom II X2 560 Processor (3300.14-MHz K8-class CPU).

> Note 2: subnormal values for x are test not included in
> the count as subnormal are tested with a different approach.

On Fri, 23 Feb 2024 22:29:17 +0000, MitchAlsup1 wrote:

> I find it interesting that sinpi() has worse numeric error than sin().

Another argument for sticking with the simpler solution? One set of
radian-specific functions plus conversion factors, instead of a separate
set of unit-specific functions for every unit?

On Fri, 23 Feb 2024 22:29:17 +0000, MitchAlsup1 wrote:

> Steven G. Kargl wrote:
>
>> On Fri, 23 Feb 2024 20:02:00 +0000, MitchAlsup1 wrote:
>
>>> Michael S wrote:
>>>
>>>> On Fri, 23 Feb 2024 00:13:02 -0000 (UTC)
>>>> "Steven G. Kargl" <sgk@REMOVEtroutmask.apl.washington.edu> wrote:
>>>
>>>>> On Thu, 22 Feb 2024 22:57:20 +0000, MitchAlsup1 wrote:
>>>>>
>>>>>
>>>>> > I have a patent on doing argument reduction for sin(x) in 4
>>>>> > cycles...that is as good as Payne-Haneok argument reduction over
>>>>> > -infinity..+infinity
>>>>>
>>>>> By a strange coincidence, I have the Payn-Hanek paper on the top of
>>>>> stack of papers on my desk. Still need to internalize the details.
>>>>>
>>>>>
>>>> Payne-Hanek is a right answer to a wrong question.
>>>
>>> If someone looking for fast sin(x) you are correct, for a person
>>> looking for the best possible (correctly rounded) result, you are not.
>>> In any event I have driven that monstrosity down from 50-odd cycles to
>>> 4-- making it at least palatable.
>>>
>>>> In science/engineering floating-point number x represents a range
>>>> [x-ULP:x+ULP] with hopefully much higher probability of being in
>>>> range [x-ULP/2:x+ULP/2]. However in this reduced range probability is
>>>> flat, an "exact" x is no worse and no better than any other point
>>>> within this range.
>>>> Which means that when we look for y=sin(x) then for
>>>> scientific/engineering purposes any answer in range*
>>>> [sin(x-ULP/2):sin(x+ULP/2)] is as "right" as any other answer. The
>>>> right answer to the question "What is a sin() of IEEE-754 binary64
>>>> number 1e17?" is "Anything in range [-1:1]". The wise and useful
>>>> answer to the same question is "You're holding it wrong".
>>>
>>> There are verification tests that know the exact value of both x and
>>> sin(x).
>>> For these uses you argument is not valid--however I am willing to
>>> grant that outside of verification, and in actual use, your argument
>>> holds as well ad the noise in the data provides.
>
>> Fortunately, for 32-bit single precision floating point, one can
>
> All of the numerics I have spoken about are DP (64-bit) values.

Sure, but it isn't possible to exhaustively DP on current hardware
(available to me). Simple testing can give one a warm fuzzy feeling.
With 10M values in the range [0,0x1p53], the SW sinpi seems to be better
than the SW sin.

% tlibm sin -d -x 0 -X 0x1p53 -N 10 -s 0
Interval tested for sin: [0,9.0072e+15]
10000000 calls, 2.594275 secs, 0.25943 usecs/call
count: 10000000
xm = 3.4119949728897955e+15, /* 0x43283e61, 0xf8ab4587 */
libm = 7.2135591711063873e-01, /* 0x3fe71559, 0x0118855e */
mpfr = 7.2135591711063862e-01, /* 0x3fe71559, 0x0118855d */
ULP = 0.77814

% tlibm sinpi -d -x 0 -X 0x1p53 -N 10 -s 0
Interval tested for sinpi: [0,9.0072e+15]
10000000 calls, 0.184541 secs, 0.01845 usecs/call
count: 10000000
xm = 1.5384297865527400e+12, /* 0x42766318, 0xf99b8bd7 */
libm = 7.2898962872051931e-01, /* 0x3fe753e2, 0x0ecf4a3e */
mpfr = 7.2898962872051942e-01, /* 0x3fe753e2, 0x0ecf4a3f */
ULP = 0.69476

Note, the timing difference is again dominated by argument reduction.

Also note, that each of the above tests took ~180 cpu seconds. At
the moment, tlibm will not use multiple cores or cpus on other nodes.

>> exhaustively test single-argument functions. For the SW implementation
>> on FreeBSD, exhaustive testing shows
>
>> % tlibm sinpi -fPED -x 0 -X 0x1p23 -s 0 Interval tested for sinpif:
>> [0,8.38861e+06]
>> 1000000 calls, 0.015453 secs, 0.01545 usecs/call
>> ulp <= 0.5: 99.842% 1253631604 | 99.842% 1253631604
>> 0.5 < ulp < 0.6: 0.158% 1989420 | 100.000% 1255621024 Max ulp:
>> 0.583666 at 6.07950642e-05
>
>> % tlibm sin -fPED -x 0 -X 0x1p23 -s 0 Interval tested for sinf:
>> [0,8.38861e+06]
>> 1000000 calls, 0.019520 secs, 0.01952 usecs/call
>> ulp <= 0.5: 99.995% 1249842491 | 99.995% 1249842491
>> 0.5 < ulp < 0.6: 0.005% 60102 | 100.000% 1249902593 Max ulp:
>> 0.500900 at 3.68277281e+05
>
> I find it interesting that sinpi() has worse numeric error than sin().
>
> I also find it interesting that the highest error is in the easy part of
> polynomial evaluation without argument reduction whereas sin() has its
> worst result in a region where significant argument reduction has
> transpired.
>
> {{Seems like another case of faster poor answers top slower correct
> answers.}}

For my purposes, a max ulp 0.538 is good enough. The time needed
to reduce this to 0.5009 is better spent on other unimplemented
libm functions in Freebsd's libm or implemented functions with
much worse ULP

>> The speed test is for 1M values evenly distributed over the interval.
>> The difference in speed for sinpi vs sin is due to the argument
>> reduction.
>
> But result precision suffers nonetheless.

Argument reduction itself isn't the culprit. Go read the FreeBSD source
code. What is done once one has the reduced argument is the issue. I
suspect that any patch you come up with would be gladly accepted.

--
steve

On Fri, 23 Feb 2024 22:39:19 +0000, Lawrence D'Oliveiro wrote:

> On Fri, 23 Feb 2024 22:29:17 +0000, MitchAlsup1 wrote:
>
>> I find it interesting that sinpi() has worse numeric error than sin().
>
> Another argument for sticking with the simpler solution? One set of
> radian-specific functions plus conversion factors, instead of a separate
> set of unit-specific functions for every unit?

Here's a counter argument for your simple conversion
factors. Changing FreeBSD's libm to return sinf(M_PI*x)
for sinpif(x), one observes

%./tlibm sinpi -f -x 0 -X 0x1p22 -PED
Interval tested for sinpif: [0,4.1943e+06]
ulp <= 0.5: 83.174% 1037372501 | 83.174% 1037372501
0.5 < ulp < 0.6: 0.937% 11690904 | 84.111% 1049063405
0.6 < ulp < 0.7: 0.689% 8587516 | 84.800% 1057650921
0.7 < ulp < 0.8: 0.497% 6200902 | 85.297% 1063851823
0.8 < ulp < 0.9: 0.323% 4023726 | 85.620% 1067875549
0.9 < ulp < 1.0: 0.157% 1952256 | 85.776% 1069827805
1.0 < ulp < 1.5: 0.347% 4332879 | 86.124% 1074160684
1.5 < ulp < 2.0: 0.247% 3083647 | 86.371% 1077244331
2.0 < ulp < 3.0: 0.360% 4491936 | 86.731% 1081736267
3.0 < ulp < 0.0: 13.269% 165496149 | 100.000% 1247232416
Max ulp: 8388278.500000 at 5.65300049e+03

% ./tlibm sinpi -f -a 5.65300049e+03
x = 5.65300049e+03f, /* 0x45b0a801 */
libm = -2.51050433e-03f, /* 0xbb248746 */
mpfr = -1.53398013e-03f, /* 0xbac90fd5 */
ULP = 8388278.50000

I can assure one of the libm and mpfr values, when x = 5653.00049,
is quite wrong.

--
steve

On Sat, 24 Feb 2024 04:03:08 -0000 (UTC), Steven G. Kargl wrote:

> I can assure one of the libm and mpfr values, when x = 5653.00049,
> is quite wrong.

That’s 9 figures, but single-precision IEEE754 floats only allow about 7.

On Sat, 24 Feb 2024 04:47:29 +0000, Lawrence D'Oliveiro wrote:

> On Sat, 24 Feb 2024 04:03:08 -0000 (UTC), Steven G. Kargl wrote:
>
>> I can assure one of the libm and mpfr values, when x = 5653.00049,
>> is quite wrong.
>
> That’s 9 figures, but single-precision IEEE754 floats only allow about
> 7.

ROTFL. I suggest (re-)reading the entire section

IEEE Std 754TM-2008

5.12.2 External decimal character sequences representing
finite numbers

However, I'll give you the salient part on p.32

For the purposes of discussing the limits on correctly
rounded conversion, define the following quantities:

for binary32, Pmin (binary32) = 9
....

Conversions from a supported binary format bf to an external
character sequence and back again results in a copy of the
original number so long as there are at least Pmin (bf)
significant digits specified and the rounding-direction
attributes in effect during the two conversions are round to
nearest rounding-direction attributes.

You also seem to miss that I gave you the bit pattern as
an unsigned 32-bit int. There are no extra binary digits.

At this point, I think it might be prudent for you to
read Goldberg's paper [1] and IEEE 754.

[1] https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

I'm done responding to your trolls and apologies other c.l.c
posters.

--
steve

On Sat, 24 Feb 2024 04:03:08 -0000 (UTC), Steven G. Kargl wrote:

> % ./tlibm sinpi -f -a 5.65300049e+03
> x = 5.65300049e+03f, /* 0x45b0a801 */
> libm = -2.51050433e-03f, /* 0xbb248746 */
> mpfr = -1.53398013e-03f, /* 0xbac90fd5 */
> ULP = 8388278.50000
>
> I can assure one of the libm and mpfr values, when x = 5653.00049,
> is quite wrong.

They’re all wrong. Python says (in double-precision), that
math.sin(5653.00049) is -0.9566595256716378.

Even GCC’s sinf function says the value is -0.9565167, which is only a
slight discrepancy in the fourth decimal place, so even in single-
precision it is doing a darn sight better than you have managed above.

So where is there a version of the code you’re using, using different
angle units, that gets closer to the right answer?

On 2024-02-24, Steven G. Kargl <sgk@REMOVEtroutmask.apl.washington.edu> wrote:
> On Sat, 24 Feb 2024 04:47:29 +0000, Lawrence D'Oliveiro wrote:
>
>> On Sat, 24 Feb 2024 04:03:08 -0000 (UTC), Steven G. Kargl wrote:
>>
>>> I can assure one of the libm and mpfr values, when x = 5653.00049,
>>> is quite wrong.
>>
>> That’s 9 figures, but single-precision IEEE754 floats only allow about
>> 7.
>
> ROTFL. I suggest (re-)reading the entire section

[ snip ]

> for binary32, Pmin (binary32) = 9

Also, from the comp.lang.c point of view, in the C language we have
FLT_DIG and FLT_DECIMAL_DIG in <float.h>.

On IEEE 754 systems these are 6 and 9.

FLT_DIG gives you, roughly speaking, how many decimal digits of
precision a float can preserve in the decimal -> float -> decimal
conversion sequence

FLT_DECIMAL_DIG indicates how many digits are required to exactly
preserve a float value in decimal form: float -> decimal -> float.

--
TXR Programming Language: http://nongnu.org/txr
Cygnal: Cygwin Native Application Library: http://kylheku.com/cygnal
Mastodon: @Kazinator@mstdn.ca

On Sat, 24 Feb 2024 05:38:53 +0000, Lawrence D'Oliveiro wrote:

> On Sat, 24 Feb 2024 04:03:08 -0000 (UTC), Steven G. Kargl wrote:
>
>> % ./tlibm sinpi -f -a 5.65300049e+03
>> x = 5.65300049e+03f, /* 0x45b0a801 */
>> libm = -2.51050433e-03f, /* 0xbb248746 */
>> mpfr = -1.53398013e-03f, /* 0xbac90fd5 */
>> ULP = 8388278.50000
>>
>> I can assure one of the libm and mpfr values, when x = 5653.00049,
>> is quite wrong.
>
> They’re all wrong. Python says (in double-precision), that
> math.sin(5653.00049) is -0.9566595256716378.
>
> Even GCC’s sinf function says the value is -0.9565167, which is only a
> slight discrepancy in the fourth decimal place, so even in single-
> precision it is doing a darn sight better than you have managed above.
>
> So where is there a version of the code you’re using, using different
> angle units, that gets closer to the right answer?

sin(x) is not sinpi(x). The conversion factor that you're missing
is M_PI as in sinpi(x) = sin(M_PI*x). The whole point of the
half-cycle trig functions is that one does not need to do the
multiplication by pi, or more importantly the argument reduction
modulo pi.

% ./tlibm sin -f -a 5.65300049e+3
x = 5.65300049e+03f, /* 0x45b0a801 */
libm = -9.56659019e-01f, /* 0xbf74e79b */
mpfr = -9.56659019e-01f, /* 0xbf74e79b */
ULP = 0.10338

--
steve

Michael S wrote:
> On Fri, 23 Feb 2024 11:01:02 +0100
> Terje Mathisen <terje.mathisen@tmsw.no> wrote:
>
>> Michael S wrote:
>>> On Thu, 22 Feb 2024 21:04:52 -0000 (UTC)
>>> Lawrence D'Oliveiro <ldo@nz.invalid> wrote:
>>>>
>>>> π radians = half a circle. Are there any other examples of the
>>>> usefulness of half-circles as an angle unit? As opposed to the
>>>> dozens or hundreds of examples of the usefulness of radians as an
>>>> angle unit?
>>>
>>> In digital signal processing circle-based units are pretty much
>>> always more natural than radians.
>>> For specific case of 1/2 circle, I can't see where it can be used
>>> directly.
>>> From algorithmic perspective, full circle looks like the most
>>> obvious choice.
>>> From [binary floating point] numerical properties perspective,
>>> 1/8th of the circle (==pi/4 radians = 45 degrees) is probably the
>>> best option for a library routine, because for sin() its derivative
>>> at 0 is closest to 1 among all powers of two which means that loss
>>> of precision near 0 is very similar for input and for output. But
>>> this advantage does not sound like particularly big deal.
>>
>> ieee754 defines sinpi() and siblings, but imho it really doesn't
>> matter if you use circles, half-circles (i.e. sinpi) or some other
>> binary fraction of a circle: Argument reduction for huge inputs are
>> just as easy, you might just have to multiply by the corresponding
>> power of two (i.e. adjust the exponent) before extracting the
>> fractional term.
>>
>> For sinpi(x) I could do it like this:
>>
>> if (abs(x) >= two_to_52nd_power) error("Zero significant bits.");
>> ix = int(x);
>> x_reduced = x - (double) (ix & ~1);
>> if (x_reduced < 0.0) x_reduced += 2.0;
>>
>> but it is probably better to return a value in the [-1.0 .. 1.0>
>> range?
>>
>> Terje
>>
>
> Both you and Mitch look at it from wrong perspective.
> When we define a library API, an ease of implementation of the library
> function should be pretty low on our priority scale. As long as
> reasonably precise calculation is theoretically possible, we should
> give credit to intelligence of implementor, that's all.
> The real concern of API designer should be with avoidance of loss of
> precision in preparation of inputs and use of outputs.
> In specific case of y=sin2pi(x), it is x that is more problematic,
> because near 0 it starts to lose precision 3 octaves before y. In
> subnormal range we lose ~2.5 bits of precision in preparation of the
> argument. An implementation, no matter how good, can't recover what's
> already lost.
> sinpi() is slightly better, but only slightly, not enough better to
> justify less natural semantics.
>
> My yesterday suggestion is a step in right direction, but today I think
> that it is not sufficiently radical step. In specific case of sin/cos
> there is no good reason to match loss of precision on input with loss of
> precision.
> <Thinking out load>
> There are advantages in matched loss of precision for input and for
> output when both input and output occupy full range of real numbers
> (ignoring sign). Output of sin/cos does not occupy o full range. But
> for tan() it does. So, may be, there is a one good reason for matching
> input with output for sin/cos - consistency with tan.
> </Thinking out load>
> So, ignoring tan(), what is really an optimal input scaling for sin/cos
> inputs? Today I think that it is a scaling in which full circle
> corresponds to 2**64. With such scaling you never lose any input
> precision to subnoramals before precision of the result is lost
> completely.
> Now, one could ask "Why 2**64, why not 2**56 that has the same
> property?". My answer is "Because 2**64 is a scaling that is most
> convenient for preparation of trig arguments in fix-point [on 64-bit
> computer]." I.e. apart from being a good scaling for avoidance of loss
> of precision in tiny range, it happens to be the best scaling for
> interoperability with fixed point.
> That is my answer today. Will it hold tomorrow? Tomorrow will tell.

This is, except for today being a 64-bit world as opposed to 30 years
earlier, exactly the same reasoning Garmin's programmers used to decide
that all their lat/long calculations would use 2^32 as the full circle.

With a signed 32-bit int you get a resolution of 40e6m / 2^32 = 0.0093m
or 9.3 mm, which they considered was more than good enough back in the
days of SA and its ~100 RMS noise, and even after Clinton got rid of
that (May 2 2000 or 2001?), sub-cm GPS is very rarely available.

Doing the same with 64-bit means that you get a resolution of 2.17e-12 m
which is 2.17 picometers or 0.0217 Å, so significantly smaller than a
single H atom which is about 1 Å in size.

Terje

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

Steven G. Kargl wrote:
> Fortunately, for 32-bit single precision floating point, one can
> exhaustively test single-argument functions. For the SW implementation
> on FreeBSD, exhaustive testing shows
>
> % tlibm sinpi -fPED -x 0 -X 0x1p23 -s 0
> Interval tested for sinpif: [0,8.38861e+06]
> 1000000 calls, 0.015453 secs, 0.01545 usecs/call
> ulp <= 0.5: 99.842% 1253631604 | 99.842% 1253631604
> 0.5 < ulp < 0.6: 0.158% 1989420 | 100.000% 1255621024
> Max ulp: 0.583666 at 6.07950642e-05

That seems worse than it should be, I would expect better results than
the 0.5009 quoted for regular sin() below, and still with the same or
better performance. The only obvious explanation would be that sin(x)
over +/- 45 degrees would have a first Taylor term of 1.0*x, and you
could probably force the first Cheby term to be the same, so by adding
that as the last stage you get improved precision for free?

>
> % tlibm sin -fPED -x 0 -X 0x1p23 -s 0
> Interval tested for sinf: [0,8.38861e+06]
> 1000000 calls, 0.019520 secs, 0.01952 usecs/call
> ulp <= 0.5: 99.995% 1249842491 | 99.995% 1249842491
> 0.5 < ulp < 0.6: 0.005% 60102 | 100.000% 1249902593
> Max ulp: 0.500900 at 3.68277281e+05
>
> The speed test is for 1M values evenly distributed over
> the interval. The difference in speed for sinpi vs sin
> is due to the argument reduction.
>
> Note 1: libm built with clang/llvm with only -O2 on a
> AMD Phenom II X2 560 Processor (3300.14-MHz K8-class CPU).
>
> Note 2: subnormal values for x are test not included in
> the count as subnormal are tested with a different approach.
>

sin(x) for Subnormal x seems like it should be trivial: If you just
return x the maximum error of a Taylor series with two terms
(x - x^3/3!)
would be that second term which is by definition zero since x is already
subnormal, right?

If you want to perform well on a platform where subnormal ops traps to a
much slower hardware or software path, then it might make sense to start
by checking for |x| < 2^-340 or so?

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