[Octave-bug-tracker] [bug #45339] sind, cosd don't do range reduction accurately

[Alexander Klein]

Follow-up Comment #1, bug #45339 (project octave):

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
<http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-software-developer-manual-325462.pdf>.

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.

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?45339>

_______________________________________________
  Nachricht gesendet von/durch Savannah
  http://savannah.gnu.org/