# bug #45339: sind, cosd don't do range reduction accurately

 Submitted by: Charles Karney Submitted on: Wed 17 Jun 2015 10:13:08 AM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Inaccurate Result Status: Patch Submitted Assigned to: None Originator Name: karney Open/Closed: Open Release: dev Operating System: GNU/Linux

You are not logged in

## Discussion

Tue 01 May 2018 10:21:27 AM UTC, comment #8:

I recommend changing all occurrences of floor(x + 0.5) to round(x) in my patched files.  I don't think this changes the results but it's cleaner and x = 0.5 - 0.5^54 is rounded correctly.

Charles Karney <karney>
Thu 05 May 2016 08:38:55 PM UTC, comment #7:

Bletch, sind and cosd need to work with complex arguments.  The
latest attachments provide this.  (atan2d is restricted to real
arguments.)

Charles Karney <karney>
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).

Charles Karney <karney>
Thu 05 May 2016 07:24:03 AM UTC, comment #5:

Would you be able to extend your patches to comply to the standard?  That would be great!

I'm changing the status to "in progress".

Lachlan Andrew <lachlan>
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

Charles Karney <karney>
Tue 03 May 2016 10:47:22 PM UTC, comment #3:

Status: None -> Patch submitted

Lachlan Andrew <lachlan>
Fri 17 Jul 2015 05:29:32 PM UTC, comment #2:

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

Charles Karney <karney>
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:

octave:179> l=10.^linspace(0,250,26),sin(2*pi*l)
l =

Columns 1 through 11:

1.0000e+00    1.0000e+10    1.0000e+20    1.0000e+30    1.0000e+40    1.0000e+50    1.0000e+60    1.0000e+70    1.0000e+80    1.0000e+90   1.0000e+100

Columns 12 through 22:

1.0000e+110   1.0000e+120   1.0000e+130   1.0000e+140   1.0000e+150   1.0000e+160   1.0000e+170   1.0000e+180   1.0000e+190   1.0000e+200   1.0000e+210

Columns 23 through 26:

1.0000e+220   1.0000e+230   1.0000e+240   1.0000e+250

ans =

Columns 1 through 12:

-2.4493e-16  -4.4787e-06  -7.2437e-01  -8.5359e-01   5.9989e-02  -1.5867e-01   3.7819e-01   9.9193e-01   2.3609e-01   9.9670e-01   2.8200e-01  -9.5447e-01

Columns 13 through 24:

2.6493e-01  -6.2556e-01  -8.2755e-01  -9.9650e-01  -9.3242e-01  -8.2853e-01   9.3635e-01   6.8579e-01   7.2881e-01  -5.1089e-01  -6.8245e-01   8.8131e-01

Columns 25 and 26:

9.2369e-01  -9.5642e-01

In this light, it might even be better to toss out the line

y(I == fix (I) & finite (I)) = 0;

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.

Alexander Klein <matalex>
Wed 17 Jun 2015 10:13:08 AM UTC, original submission:

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

The Matlab implementation of sind, cosd, atan2d behave
similarly to these attached versions.

Charles Karney <karney>

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #37084:  sindx.m added by karney (765B - text/x-objcsrc - sind and cosd generalized to complex arguments)
file #37085:  cosdx.m added by karney (757B - text/x-objcsrc - sind and cosd generalized to complex arguments)
file #37079:  sindx.m added by karney (620B - text/x-objcsrc - New versions of sind, cosd, atan2d which treat +/-0 following IEEE 754-2008)
file #37080:  atan2dx.m added by karney (975B - text/x-objcsrc - New versions of sind, cosd, atan2d which treat +/-0 following IEEE 754-2008)
file #37081:  cosdx.m added by karney (615B - text/x-objcsrc - New versions of sind, cosd, atan2d which treat +/-0 following IEEE 754-2008)
file #34242:  sindx.m added by karney (319B - text/x-objcsrc - Replacements for sind, cosd, and atan2d)
file #34243:  cosdx.m added by karney (323B - text/x-objcsrc - Replacements for sind, cosd, and atan2d)
file #34244:  atan2dx.m added by karney (480B - text/x-objcsrc - Replacements for sind, cosd, and atan2d)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by mtmiller (Updated the item)
• -email is unavailable- added by lachlan (Posted a comment)
• -email is unavailable- added by matalex (Posted a comment)
• -email is unavailable- added by karney (Submitted the item)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.

## History

Date Changed by Updated Field Previous Value => Replaced by
2018-05-01 mtmiller Release3.8.2 dev
2018-05-01 mtmiller StatusIn Progress Patch Submitted
2016-05-05 karney Attached File- Added sindx.m, #37084