Thu 05 May 2016 06:04:10 PM UTC, comment #6:
Attached are new versions of the routines which comply with the
IEEE 759-2008 conventions for sinPi, cosPi, atan2Pi. No code
change was necessary for atan2d. However, I added a comment
with the special cases to check (and similarly for sind and cosd).
(file #37079, file #37080, file #37081)
|
Wed 04 May 2016 12:26:03 PM UTC, comment #4:
I've checked the specifications of sinPi, cosPi, atan2Pi in IEEE Std 754
2008. I recommend making sind, cosd, atan2d behave in compatible ways
(and I don't think that the patches I provided ensure this).
Specifically
[sind is odd]
sind(n*180) = -0 if n is a negative integer
sind(n*180) = +0 if n is a positive integer
sind(+/-0) = +/-0
[cosd is even]
cosd(n*180+90) = +0 for all integer n
[symmetries of atan2d]
atan2d(+/-0, -0) = +/-180
atan2d(+/-0, +0) = +/-0
atan2d(+/-0, x) = +/-180 for x < 0
atan2d(+/-0, x) = +/-0 for x > 0
and similar special cases for x or y = +/-inf
|
Fri 17 Jul 2015 05:29:32 PM UTC, comment #2:
Your observations about the challenges of argument reduction for the
sine and cosine routines are entirely correct.
However when the arguments are in degrees, as they are for sind and
cosd, argument reduction is exact. Indeed, this is one reason to keep
angles as degrees instead of converting them to radians. It's also
possible to reduce arguments to the range [-45deg, 45deg] exactly, and
this allows the sind and cosd routines to obey identities like sind(1) =
cosd(89).
Having these properties built into the default implementations of sind
and cosd is an obvious win. (And, as I noted, the MATLAB
implementations work like this.)
I should note one defect of my implementation and that is that sind(180)
= cosd(90) = -0 (instead of +0). This is easily fixed (just add 0 to
the result). Unfortunately sind(-0) returns +0 instead of -0; but that
is a consequence of rem(-0,360) returning +0 instead of -0 (a bug
surely?).
|
Mon 06 Jul 2015 10:00:50 AM UTC, comment #1:
Welcome to the pitfalls of floating point arithmetic! ;)
Although it is true that in the realm of real numbers you would expect the sum
sind(1e9) + sind(80)
to evaluate to zero, this is not necessarily the case for floating point numbers. Even worse, the kind of breakage experienced is likely to be FPU-dependent, see section "8.3.10 Transcendental Instruction Accuracy" in Intel's Manual.
Since the trigonometric functions in a processor cannot really be 2π-periodic, but will have a period depending on the internals of the processor, making the suggested changes would break compatibility with the other trigonometric functions.
Of course, all basic trigonometric functions could reduce their arguments to some value within an interval considered safe, but then again, x87 FPUs do that internally, anyway, and what should you do if you really wanted to evaluate the sine of huge values without the reduction, which turn out to be quite different from the expected results:
In this light, it might even be better to toss out the line
found in the current version, make the basic trigonometric functions throw a warning for arguments outside the range considered safe on the current FPU, instead, and leave any reduction of the argument to people who know how to reduce the arguments correctly in their respective case.
|
Wed 17 Jun 2015 10:13:08 AM UTC, original submission:
One of the reasons for working with degrees instead of radians
is that range reduction is exact. However, the octave trig
functions in terms of degrees don't make use of this. Thus
sind(1e9) + sind(80) returns 2e-10 instead of 0.
I attach replacements for these functions which remove this
deficiency and obey identities such as cosd(89) = sind(1).
Although the problems the atan2d are nearly as bad, I also
provide a replacement for this function which enforces the
quadrant symmetries.
The Matlab implementation of sind, cosd, atan2d behave
similarly to these attached versions.
|