bugGNU Octave - Bugs: bug #45339, sind, cosd don't do range...

 
 

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

Submitter:  Charles Karney <karney>
Submitted:  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
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

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

(file #37084, file #37085)

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

(file #37079, file #37080, file #37081)

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: 

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

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:  

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.

Charles Karney <karney>

 

(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)

 

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)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 12 latest changes.

    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
        Attached File- Added cosdx.m, #37085
    2016-05-05 karney Attached File- Added sindx.m, #37079
        Attached File- Added atan2dx.m, #37080
        Attached File- Added cosdx.m, #37081
    2016-05-05 lachlan StatusPatch Submitted In Progress
    2016-05-03 lachlan StatusNone Patch Submitted
    2015-06-17 karney Attached File- Added sindx.m, #34242
        Attached File- Added cosdx.m, #34243
        Attached File- Added atan2dx.m, #34244

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code