Terje Mathisen wrote:

>

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

The 2 term polynomial begins to be accurate when |x| < ~2^-19 which is far
bigger than 2^-1023. (2^-19)^3 is 2^-57 which has few digits in the resulting
fraction 2^-19; certainly no 3rd term is required.

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

Many subnormal operands yield constant results or themselves::

cos(sn) = 1.0
sin(sn) = sn
tan(sn) = sn

> Terje

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

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

At which cost in tables sizes?

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

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.

The definition of 'exact' should be:
For any finite-precision function foo(x) lets designate the same
mathematical function calculated with infinite precision as Foo(x).
Let's designate an operation of rounding of infinite-precision number to
desired finite precision as Rnd(). Rounding is done in to-nearest mode.
Unlike in the case of basic operations, ties are allowed to be broken in
any direction.
The result of y=foo(x) for finite-precision number x considered
exact if *at least one* two conditions is true:
(1=Y-clause) Rnd(Foo(x)) == y
(2=X-clause) There exist an infinite precision number X for which
both Foo(X) == y and Rnd(X) == x.

As follows from the (2), it is possible and not uncommon that more
than one finite-precision number y is accepted exact result of foo(x).

If Committee omits the 2nd clause then the whole proposition will be not
just not useful, but harmful.

Steven G. Kargl wrote:

> sin(x) is not sinpi(x). The conversion factor that you're missing
> is M_PI as in sinpi(x) = sin(M_PI*x).

When I faced this kind of accuracy/precision problem in my HW transcendentals,

To get a sufficiently correct reduced argument I had to multiply the fraction
of x by a 2/pi number selected such that the HoB was aligned to 2 (quadrant)
and the multiplied result had 51+53 bits of precision so that up to 51 leading
bits of the reduced argument (after the quadrant bits) could be skipped if 0.
This required 3 uses of the multiplier tree one produced the leading bits::

Lead bits |/
Middle bits / /
trail bits /|

which were arranged |/ /| into a single 104 bit (minimum) product. My current
implementation uses 128 bits. And my polynomial evaluation uses 58-bit argu-
ments (min, 64-bit current) at each iteration. And I still only get 0.502
(58-bit:: 0.5002 64-bit) precision.

Payne and Hanek argument reduction--because it does not have access to the
intermediate bits, needs 4 multiplies instead of 2 and very careful arithmetic
to preserve accuracy. I can do this in 2 trips through the multiplier array
(patented)

So, I used 159-bits of 2/pi in 2 multiplies over 2 trips through the array
and get perfect argument (DP) reduction in 4 cycles. Payne ad Hanek use 256-
bits of DP FP operands and 30+ instruction to do the same thing.

My multiplier tree is cut into 2 sections (just like one would do for dual SP)
but here I feed the top 2/pi bits into the left hand side and the bottom 2/pi
bits into the right hand side so both get computed simultaneously; the subsequent
cycle multiplies by the middle bits of 2/pi. The 2/pi table is indexed by exponent.

Terje Mathisen wrote:

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

The drift rate of the oscillators in the satellites is such it will neve be.
The clock drift is reset every orbit.

Also note: hackers are now using ground based GPS transmitters to alter
where your GPS calculates where you think you are. This is most annoying
around airports when planes use GPS to auto guide the planes to runways.

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

And yet, driving by Edwards AFB sometimes my car's GPS shows my 50m off the
interstate quality road, and sometimes not.

> Terje

Michael S wrote:

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

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

> At which cost in tables sizes?

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

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

> The definition of 'exact' should be:
> For any finite-precision function foo(x) lets designate the same
> mathematical function calculated with infinite precision as Foo(x).
> Let's designate an operation of rounding of infinite-precision number to
> desired finite precision as Rnd(). Rounding is done in to-nearest mode.
> Unlike in the case of basic operations, ties are allowed to be broken in
> any direction.
> The result of y=foo(x) for finite-precision number x considered
> exact if *at least one* two conditions is true:
> (1=Y-clause) Rnd(Foo(x)) == y
> (2=X-clause) There exist an infinite precision number X for which
> both Foo(X) == y and Rnd(X) == x.

In the second clause:: are we guaranteed that RND(Foo(X))= Y ??

> As follows from the (2), it is possible and not uncommon that more
> than one finite-precision number y is accepted exact result of foo(x).

> If Committee omits the 2nd clause then the whole proposition will be not
> just not useful, but harmful.

An interesting side effect of greater intermediate precision is the lack
of need to round prior to the final result. Thus, my sin(x) over its
entire calculation suffers exactly 1 rounding. Payne & Hanek does not
have this prperty.

And thirdly::

Let us postulate that the reduced argument r is indeed calculated with 0.5ULP
error.

What makes you think you can calculate the polynomial without introducing
any more error ??

On Thu, 14 Mar 2024 17:34:57 +0000, MitchAlsup1 wrote:

> Michael S wrote:
>
>> (2=X-clause) There exist an infinite precision number X for which
>> both Foo(X) == y and Rnd(X) == x.
>
> In the second clause:: are we guaranteed that RND(Foo(X))= Y ??

No idea why that’s relevant. Michael S was talking about “Rnd”, not “RND”.

When you say “RND”, I think of a bad random-number generator found in
various implementations of BASIC.

On 3/14/2024 10:28 AM, MitchAlsup1 wrote:
> Terje Mathisen wrote:
>
>> 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.
>
> The drift rate of the oscillators in the satellites is such it will neve
> be.
> The clock drift is reset every orbit.
>
> Also note: hackers are now using ground based GPS transmitters to alter
> where your GPS calculates where you think you are. This is most annoying
> around airports when planes use GPS to auto guide the planes to runways.

GPS Spoofing?

>
>> 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.
>
> And yet, driving by Edwards AFB sometimes my car's GPS shows my 50m off the
> interstate quality road, and sometimes not.
>
>> Terje

On 3/14/2024 1:30 PM, Lawrence D'Oliveiro wrote:
> On Thu, 14 Mar 2024 17:34:57 +0000, MitchAlsup1 wrote:
>
>> Michael S wrote:
>>
>>> (2=X-clause) There exist an infinite precision number X for which
>>> both Foo(X) == y and Rnd(X) == x.
>>
>> In the second clause:: are we guaranteed that RND(Foo(X))= Y ??
>
> No idea why that’s relevant. Michael S was talking about “Rnd”, not “RND”.
>
> When you say “RND”, I think of a bad random-number generator found in
> various implementations of BASIC.

PRNG, an LGC? Fwiw, check this shit out:

https://groups.google.com/g/comp.lang.c++/c/7u_rLgQe86k/m/fYU9SnuAFQAJ

;^D

Chris M. Thomasson wrote:

> On 3/14/2024 10:28 AM, MitchAlsup1 wrote:
>>
>> The drift rate of the oscillators in the satellites is such it will neve
>> be.
>> The clock drift is reset every orbit.
>>
>> Also note: hackers are now using ground based GPS transmitters to alter
>> where your GPS calculates where you think you are. This is most annoying
>> around airports when planes use GPS to auto guide the planes to runways.

> GPS Spoofing?

More so in Europe than USA, but it is coming. Mentor Pilot had a video on
this last week.

Lawrence D'Oliveiro wrote:

> On Thu, 14 Mar 2024 17:34:57 +0000, MitchAlsup1 wrote:

>> Michael S wrote:
>>
>>> (2=X-clause) There exist an infinite precision number X for which
>>> both Foo(X) == y and Rnd(X) == x.
>>
>> In the second clause:: are we guaranteed that RND(Foo(X))= Y ??

> No idea why that’s relevant. Michael S was talking about “Rnd”, not “RND”.

You have just selected yourself as someone I will never reply to again.

> When you say “RND”, I think of a bad random-number generator found in
> various implementations of BASIC.

On 3/14/2024 3:19 PM, MitchAlsup1 wrote:
> Lawrence D'Oliveiro wrote:
>
>> On Thu, 14 Mar 2024 17:34:57 +0000, MitchAlsup1 wrote:
>
>>> Michael S wrote:
>>>
>>>> (2=X-clause) There exist an infinite precision number X for which
>>>> both Foo(X) == y and Rnd(X) == x.
>>>
>>> In the second clause:: are we guaranteed that RND(Foo(X))= Y ??
>
>> No idea why that’s relevant. Michael S was talking about “Rnd”, not
>> “RND”.
>
> You have just selected yourself as someone I will never reply to again.

It might be an AI? Just have a very strange feeling... Humm...

>
>> When you say “RND”, I think of a bad random-number generator found in
>> various implementations of BASIC.

On 3/14/2024 3:21 PM, Chris M. Thomasson wrote:
> On 3/14/2024 3:19 PM, MitchAlsup1 wrote:
[...]

An AI set on ass mode?

On 3/14/2024 3:18 PM, MitchAlsup1 wrote:
> Chris M. Thomasson wrote:
>
>> On 3/14/2024 10:28 AM, MitchAlsup1 wrote:
>>>
>>> The drift rate of the oscillators in the satellites is such it will
>>> neve be.
>>> The clock drift is reset every orbit.
>>>
>>> Also note: hackers are now using ground based GPS transmitters to alter
>>> where your GPS calculates where you think you are. This is most annoying
>>> around airports when planes use GPS to auto guide the planes to runways.
>
>> GPS Spoofing?
>
> More so in Europe than USA, but it is coming. Mentor Pilot had a video on
> this last week.

Thanks. The device thinks its all right... However, its connected to a
hyper nefarious dynamic data provider. Yikes!

Michael S wrote:

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

>>

> Both you and Mitch look at it from wrong perspective.

For the record, I am only addressing SIN() or SINPI() as a HW instruction.

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

At the HW level, when someone executes a
SIN R9,R17
instruction, they should get a result that has been computed to very
high precision, including argument reductions (which may require more
bits in the reduced argument than fit in an IEEE container).

> The real concern of API designer should be with avoidance of loss of
> precision in preparation of inputs and use of outputs.

With an instruction performing all the work, this argument is moot.

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

NOT when the argument reduction process does not contain ROUNDINGs !!!
OR when the argument reduction process contains more precision bits than
fit in an IEEE container !!! I determined for SIN() this needed 57-bits
of intermediate reduced argument, for TAN() and ATAN() this requires
58-bits:: compared to the 53 available in normalized IEEE 754.

> An implementation, no matter how good, can't recover what's
> already lost.

My 66000 ISA has these as instructions (not function calls) and has
the execution width to avoid your perils.

> sinpi() is slightly better, but only slightly, not enough better to
> justify less natural semantics.

SINPI() is also an instruction.
SIN() takes 19 cycles whereas SINPI() takes 16.....

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

In the end, you either have a precise result, or you do not.

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

See USPTO 10,761,806 and 10,983,755 for your answer.

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.

Actually using and depending upon the result is more stupid.

OTOH, it is and have always been, a core principle of ieee754 that basic
operations (FADD/FSUB/FMUL/FDIV/FSQRT) shall assume that the inputs are
exact (no fractional ulp uncertainty), and that we from that starting
point must deliver a correctly rounded version of the infinitely precise
exact result of the operation.

Given the latter, it is in fact very tempting to see if that basic
result rule could be applied to more of the non-core operations, but I
cannot foresee any situation where I would use it myself: If I find
myself in a situation where the final fractional ulp is important, then
I would far rather switch to doing the operation in fp128.

Terje

Michael S wrote:
> On Fri, 23 Feb 2024 11:10:00 +0100
> Terje Mathisen <terje.mathisen@tmsw.no> wrote:
>
>> 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.
>>
>
> At which cost in tables sizes?
>
>
>> There is a suggestion on the table to make that a (probably optional
>> imho) feature for an upcoming ieee754 revision.
>>
>> Terje
>>
>
> 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.
>
> The definition of 'exact' should be:
> For any finite-precision function foo(x) lets designate the same
> mathematical function calculated with infinite precision as Foo(x).
> Let's designate an operation of rounding of infinite-precision number to
> desired finite precision as Rnd(). Rounding is done in to-nearest mode.
> Unlike in the case of basic operations, ties are allowed to be broken in
> any direction.
> The result of y=foo(x) for finite-precision number x considered
> exact if *at least one* two conditions is true:
> (1=Y-clause) Rnd(Foo(x)) == y
> (2=X-clause) There exist an infinite precision number X for which
> both Foo(X) == y and Rnd(X) == x.
>
> As follows from the (2), it is possible and not uncommon that more
> than one finite-precision number y is accepted exact result of foo(x).
>
> If Committee omits the 2nd clause then the whole proposition will be not
> just not useful, but harmful.
>

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

MitchAlsup1 wrote:
> Terje Mathisen wrote:
>> 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.
>
> The drift rate of the oscillators in the satellites is such it will neve
> be.
> The clock drift is reset every orbit.

Sub-meter (typically 2 cm at 10-20 Hz) GPS requires a nearby (static)
base station which can measure the individual pseudo-range errors from
each sat and send that to the measuring device (the rover) over a
separate channel.

If you record all the pseudorange values, then you can do the same with
post-processing, this can be useful for DIY surveying.
>
> Also note: hackers are now using ground based GPS transmitters to alter
> where your GPS calculates where you think you are. This is most annoying
> around airports when planes use GPS to auto guide the planes to runways.

The potential for this attack is one of the reasons the military signal
is encrypted, so that it cannot be spoofed.
>
>> 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.
>
> And yet, driving by Edwards AFB sometimes my car's GPS shows my 50m off the
> interstate quality road, and sometimes not.

50 m is a bit high, but still fairly typical for what commercial
receivers can do in the vicinity of refelcting surfaces like cliffs or
buildings.

Modern multi-system receivers are actually getting much better at
detecting and mitigating such problems. Starting with
GPS+Glonass+Galileo having worldwide coverage, then adding in the
Chinese and Japanese sats means that as long as you don't worry to much
about pwer usage, you can do quite well. My personal target is sub 3m
when under a wet forest canopy, that is good enough for orienteering map
field survey work.

Terje

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

MitchAlsup1 wrote:
> And thirdly::
>
> Let us postulate that the reduced argument r is indeed calculated with
> 0.5ULP
> error.
>
> What makes you think you can calculate the polynomial without introducing
> any more error ??

It should be obvious that any argument reduction step must return a
value with significantly higher precision than the input(s), and that
this higher precision value is then used in any polynomial evaluation.

With careful setup, it is very often possible to reduce the amount of
extended-procision work needed to just one or two steps, i.e. for the
classic Taylor sin(x) series, with x fairly small, the x^3 and higher
terms can make do with double precision, so that the final step is to
add the two parts of the leading x term: First the trailing part and
then when adding the upper 53 bits of x you get a single rounding at
this stage.

This is easier and better when done with 64-bit fixed-point values,
augemented with a few 128-bit operations.

Terje

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

On Thu, 14 Mar 2024 17:34:57 +0000
mitchalsup@aol.com (MitchAlsup1) wrote:

> Michael S wrote:
>
> > (2=X-clause) There exist an infinite precision number X for which
> > both Foo(X) == y and Rnd(X) == x.
>
> In the second clause:: are we guaranteed that RND(Foo(X))= Y ??
>

No, we are not.

> > As follows from the (2), it is possible and not uncommon that more
> > than one finite-precision number y is accepted exact result of
> > foo(x).
>
> > If Committee omits the 2nd clause then the whole proposition will
> > be not just not useful, but harmful.
>
> An interesting side effect of greater intermediate precision is the
> lack of need to round prior to the final result. Thus, my sin(x) over
> its entire calculation suffers exactly 1 rounding. Payne & Hanek does
> not have this prperty.

Which does not help to recover precision lost during rounding of x
that happened before your wonderful instruction.

On Fri, 15 Mar 2024 11:23:45 +0100
Terje Mathisen <terje.mathisen@tmsw.no> 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.
>
> Actually using and depending upon the result is more stupid.
>
> OTOH, it is and have always been, a core principle of ieee754 that
> basic operations (FADD/FSUB/FMUL/FDIV/FSQRT) shall assume that the
> inputs are exact (no fractional ulp uncertainty), and that we from
> that starting point must deliver a correctly rounded version of the
> infinitely precise exact result of the operation.
>
> Given the latter, it is in fact very tempting to see if that basic
> result rule could be applied to more of the non-core operations, but
> I cannot foresee any situation where I would use it myself: If I find
> myself in a situation where the final fractional ulp is important,
> then I would far rather switch to doing the operation in fp128.
>
> Terje
>

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

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.

On Fri, 15 Mar 2024 12:00:45 +0100, Terje Mathisen
<terje.mathisen@tmsw.no> wrote:

>MitchAlsup1 wrote:
>> Terje Mathisen wrote:
>>> 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.
>>
>> The drift rate of the oscillators in the satellites is such it will neve
>> be.
>> The clock drift is reset every orbit.
>
>Sub-meter (typically 2 cm at 10-20 Hz) GPS requires a nearby (static)
>base station which can measure the individual pseudo-range errors from
>each sat and send that to the measuring device (the rover) over a
>separate channel.
>
>If you record all the pseudorange values, then you can do the same with
>post-processing, this can be useful for DIY surveying.
>>
>> Also note: hackers are now using ground based GPS transmitters to alter
>> where your GPS calculates where you think you are. This is most annoying
>> around airports when planes use GPS to auto guide the planes to runways.
>
>The potential for this attack is one of the reasons the military signal
>is encrypted, so that it cannot be spoofed.
>>
>>> 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.
>>
>> And yet, driving by Edwards AFB sometimes my car's GPS shows my 50m off the
>> interstate quality road, and sometimes not.
>
>50 m is a bit high, but still fairly typical for what commercial
>receivers can do in the vicinity of refelcting surfaces like cliffs or
>buildings.

Also the AFB operates differential GPS for its runways. An OTS unit
might be confused by close proximity to the ground transmitter.

>Modern multi-system receivers are actually getting much better at
>detecting and mitigating such problems. Starting with
>GPS+Glonass+Galileo having worldwide coverage, then adding in the
>Chinese and Japanese sats means that as long as you don't worry to much
>about pwer usage, you can do quite well. My personal target is sub 3m
>when under a wet forest canopy, that is good enough for orienteering map
>field survey work.
>
>Terje

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. :^)

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?

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. :^)
>
> Now, wrt the results, arbitrary precision for trig is useful, in say...
> Deep fractal zooms...

Say I want at least 30 digits of precision for a certain calculation of
complex numbers at a certain zoom level. Along the lines of trying to
match convergents of continued fractions, take pi:

3.1415926535897932384...

two digits of decimal precision:

22/7 = 3.14_28571428571428571428571428571

three digits of decimal precision:

333/106 = 3.1415_094339622641509433962264151

six digits of decimal precision:

355/113 = 3.141592_9203539823008849557522124

On and on...

Well...

>
> 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 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. :^)
>
> 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?

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? (Since I'm posting to
comp.lang.c, I'll mention that "^" is meant to be exponentation, not
bitwise xor.)

--
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 3/15/2024 2:16 PM, Keith Thompson wrote:
> "Chris M. Thomasson" <chris.m.thomasson.1@gmail.com> writes:
>> 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. :^)
>>
>> 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?
>
> 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?

Nope. :^) Fwiw, my fractals deal with zooming in and the numbers
comprising their initial state is tiny compared to 2^53. For instance
take these initial states for my multijulia fractal:

____________
c0 = {0.5,0.0}, c1 = {-5.5,0.0}

c0 = {0.0,1.0}, c1 = {0.0,-1.0}

c0 = {0.726514,0.240242}, c1 = {0.171039,0.235043}

c0 = {-1.444991,0.139179}, c1 = {-0.063294,-1.401774}

c0 = {1,0}, c1 = {-1,0}

c0 = {-.75, .06 }, c1 = {-.45, .6 }

c0 = {-1,0}, c1 = {1,0}, c2 = {0,-1}, c3 = {0,1}
____________

Notice how the numbers are small? However, if we zoom in enough, then
2^53 might become an issue..

> (Since I'm posting to
> comp.lang.c, I'll mention that "^" is meant to be exponentation, not
> bitwise xor.)
>

Indeed.