[scilab-Users] Atan function

[Michaël Baudin]

Hi,

I do not see any bug in Scilab in this thread except, perhaps, in Matlab.

At first, we can inspect the plot of the atan(y,x) function:

scf();
N=20;
y=linspace(-4,4,N);
x=linspace(-4,4,N);
Z=(feval(y,x,myatan))';
surf (y ,x ,Z);
h = gcf ();
cmap = graycolormap (10);
h. color_map = cmap ;
xtitle("atan(y,x)","y","x");

There is a discontinuity on the line y=0 for x<0.

The following session:

-->atan(0,+0)
  ans  =
     0.
-->atan(0,-0)
  ans  =
     3.1415927

is caused by the signed zeros of the IEEE floating point arithmetic. The 
values returned by atan are consistent with the mathematical limit of 
the atan(0,x) function.

The signed zero exists because there are two ways of representing zero 
in a floating point format:

  * (-1)^0 * 0 * 2^-1022
  * (-1)^1 * 0 * 2^-1022

In both case, the significant is zero and the exponent is set to the 
minimum. But in the first case, the sign bit is set to 0, while in the 
second it is set to 1. This topic is covered in "Floating point numbers 
in Scilab", section 3.8 "Signed zeros":

http://forge.scilab.org/index.php/p/docscifloat/downloads/

These signed zeros can be consistently used to make branch cuts of 
inverse trigonometric function consistent, as emphasized by Kahan in 
"Branch cuts for complex elementary functions, or much ado about 
nothing's sign bit".

Exactly the same session can be reproduced in Octave, although Scilab 
atan(y,x) is called atan2(y,x) in Octave:

octave-3.2.4.exe:5> atan2([+0 -0 +0 -0],[+0 +0 -0 -0])
ans =

    0.00000  -0.00000   3.14159  -3.14159

This is the same result as in Scilab:

-->atan([+0 -0 +0 -0],[+0 +0 -0 -0])
  ans  =
     0.    0.    3.1415927  - 3.1415927

These values are consistent with the values given in the ISO C99 
standard (see the section F.9.1.4 "The atan2 functions"). These values 
are also consistent with the wikipedia page:

http://en.wikipedia.org/wiki/Atan2

See the "Variations" section and the note on the FPATAN function.

On the other hand, it seems to me that Matlab has a chosen to set zero 
on output whenever one of the inputs x or y is zero:

 >> atan2([+0 -0 +0 -0],[+0 +0 -0 -0])
ans =
      0     0     0     0

My point of view is that this is a bug, but it seems that this is on 
purpose:

http://www.mathworks.com/support/solutions/en/data/1-1CBTV/index.html?product=ML

Most functions consistently take into account for signed zeros in 
Scilab. The fact that

-->-0==+0
  ans  =
   T

is not really surprising: it is consistent with the mathematical 
definition of the zero.

The "sign" function in Scilab does not allow to distinguish these two zeros:

-->sign(+0)
  ans  =
     0.
-->sign(-0)
  ans  =
     0.

This is consistent with the mathematical definition of the "sign" 
function, which is zero when x is zero.

The missing function in Scilab is "signbit", which returns the sign bit 
of the given number. A implementation is available in the floating point 
module :

http://forge.scilab.org/index.php/p/floatingpoint/source/tree/HEAD/macros/flps_signbit.sci

This implementation is based on the fact that the most simple way to 
know the sign bit of a real number is to invert it, as shown in the 
following session:

-->ieee(2)
-->1/-0
  ans  =
   - Inf
-->1/+0
  ans  =
     Inf

Using this function is simple with ATOMS:

-->atomsInstall("floatingpoint");
-->atomsLoad("floatingpoint");
-->b = flps_signbit ( +0 )
  b  =
   F
-->b = flps_signbit ( -0 )
  b  =
   T

There are a couple of functions also missing in this category, 
including: copysign, hypot, isfinite, pow2.

Best regards,

Michaël Baudin

PS
Here are a couple of plots.

// 2D plots
scf();
subplot(2,2,1);
z2=atan(y,-4*ones(y));
plot(y,z2)
xtitle("atan(y,-4)","y");
//
subplot(2,2,2);
z2=atan(y,-2*ones(y));
plot(y,z2)
xtitle("atan(y,-2)","y");
//
subplot(2,2,3);
z2=atan(y,+2*ones(y));
plot(y,z2)
xtitle("atan(y,2)","y");
//
subplot(2,2,4);
z2=atan(y,+4*ones(y));
plot(y,z2)
xtitle("atan(y,4)","y");



Le 11/05/2011 15:49, Samuel GOUGEON a écrit :
>  ----- Message d'origine -----
> De : Serge Steer
> Date : 10/05/2011 22:57:
>> Le 10/05/2011 21:18, Samuel GOUGEON a écrit :
>>>  Look at that:
>>>
>>> -->atan(0, 0)
>>>  ans  =
>>>     0.00000000D+00
>>>
>>> -->atan(0, -0)
>>>  ans  =
>>>     3.14159265D+00
>>>
>>> --> -0 == 0
>>>  ans  =
>>>   T
>>>
>>> :-(( Funny, isn't it ?!
>> This later result is really quite  funny, but the atan result is just 
>> due to the branch cut definition for that function
> The "special" branching cut is used only when atan() is used with a 
> single argument.
> It is not the case here. The help page states: The range of atan(y,x) 
> is (-pi, pi].
> while : atan(-0,-0)== -%pi  => %T
>
> By the way, how is it possible to explain that the sign of -0 looks to 
> be recorded
> in floating point mode?? This "feature" makes Scilab somewhat 
> unconsistent:
> -->floor(-0)
>  ans  =
>     0.
> // -1 would instead be expected!
> -->sign(-0)
>  ans  =
>     0.
> -->sign(-number_properties("tiniest"))
>  ans  =
>   - 1.
> This is why atan(y,x) should return always the same result if x and y 
> are reals and null.
> -->help minus // does not document this special feature, no more than
> -->help sign
>
> Finally, to answer to Laurent:
>
> How can I guess that variable x is -0 and not 0?
>
> is_0minus = (x==0 & atan(0,x)==%pi);
>
> but we should hope that this puzzling feature will be fixed.
>
> Samuel
>


-- 
Michaël Baudin
Ingénieur de développement
michael.baudin at scilab.org
-------------------------
Consortium Scilab - Digiteo
Domaine de Voluceau - Rocquencourt
B.P. 105 - 78153 Le Chesnay Cedex
Tel. : 01 39 63 56 87 - Fax : 01 39 63 55 94