bugGNU Octave - Bugs: bug #37379, atan2 and other trig functions...

 
 

bug #37379: atan2 and other trig functions lack documentation

Submitter:  Michael Godfrey <godfrey>
Submitted:  Tue 18 Sep 2012 12:56:24 AM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Documentation
Status:  Fixed Assigned to:  None
Originator Name:  godfrey Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 26 Jan 2014 04:09:36 AM UTC, comment #21: 

Thanks for taking a look. Closing this report as fixed.

Mike Miller <mtmiller>
Group Member
Sun 26 Jan 2014 02:54:55 AM UTC, comment #20: 

I think that this was resolved a long time ago.
Good to close this report.

Michael Godfrey <godfrey>
Group Member
Fri 24 Jan 2014 10:48:36 PM UTC, comment #19: 

I have no idea anymore about what was going on with this report.  Sorry.

--Rik

Rik <rik5>
Group administrator
Fri 24 Jan 2014 10:10:42 PM UTC, comment #18: 

This bug history looks like it has gone through several iterations. It's labeled as documentation but the most recent comments seem to be talking about changing behavior of the atan2 function. Michael or Rik, can one of you refresh your memories on what was left to work on for this bug? Focus it on improving atan2 mapper function? Documentation issues are all resolved? If there were other related issues let's open separate reports to keep things straight.

Mike Miller <mtmiller>
Group Member
Fri 28 Sep 2012 11:53:11 PM UTC, comment #17: 

I looked at lo-mappers.cc.  It appears that
something like what is below would be a start,
but a serious numerical C++ person should take a look.

This should go somewhere after line 249 in lo-mappers.cc:

i.e., after:  // (complex, complex) -> complex mappers.

Complex
atan2 (const Complex& x, const Complex& y)
{
  return 2 * atan((sqrt (x)^2 + (y)^2 - (x) ) / (y) );
}
=====================================

This does not handle the exceptions and quadrants
as far as I can tell, but it may be a start.

Uniformity would definitely be enhanced if atan2 can
be integrated with the other trig functions in
lo-mappers.cc.

Michael Godfrey <godfrey>
Group Member
Fri 28 Sep 2012 08:49:38 PM UTC, comment #16: 

It should probably be in C++.  The atan2 function is defined in libinterp/interpfcn/data.cc beginning at line 213.  It doesn't look like it would be difficult to translate your m-file code to C++ in this instance.

There is another possibility as well which is that atan2 could be shifted to a mapper function.  The mapper functions are in liboctave/numeric/lo-mappers.[ch].  If you look in lo-mappers.h you will see that there is a prototype declaration for each function, like asin, that might accept a complex value.  lo-mappers.cc then has a tiny function to implement it.  I would ask jwe on the maintainers list about this strategy, though, before going forward.

Rik <rik5>
Group administrator
Fri 28 Sep 2012 08:34:50 PM UTC, comment #15: 

Well, the code did not apparently make it into
the comment below.  In Any case, this code only
handles the exceptions for scalar inputs.

function res = atan2(y, x)

% Test for special cases and quadrants.
  if(x == 0)
    if(y > 0)
      res = pi/2;
    elseif (y < 0)
      res = -pi/2;
    else
      res = NaN;
    endif
    return;
  endif
  res = atan(y/x);
  if(x < 0)
    if(y >= 0)
      res = res + pi;
    else
      res = res - pi;
    endif
  endif
  return;
endfunction

Michael Godfrey <godfrey>
Group Member
Fri 28 Sep 2012 07:07:42 PM UTC, comment #14: 

Below is a try at making atan2 more compatible
with the other trig functions, which allow complex
arguments.  This needs more checking, and maybe
it should really be in C?


function res = atan2(y, x)

% Test for special cases and quadrants.
  if(x == 0)
    if(y > 0)
      res = pi/2;
    elseif (y < 0)
      res = -pi/2;
    else
      res = NaN;
    endif
    return;
  endif
  res = atan(y/x);
  if(x < 0)
    if(y >= 0)
      res = res + pi;
    else
      res = res - pi;
    endif
  endif
  return;
endfunction

Michael Godfrey <godfrey>
Group Member
Fri 28 Sep 2012 02:07:48 AM UTC, comment #13: 

I have looked at this a bit more.  The single
argument trig and inverse trig functions appear
to follow Kahan for complex arguments.  This
differs from (extends) the ISO/IEC definitions
for reals (with +- 0).  However, atan2 does not
do this.  Thus, for complex arguments
atan2(y,x) != atan(y/x) despite the help doc which
says that it is.  atan2 does not accept complex
arguments.

Maybe atan2(y,x) should correctly preserve the
signs of y and x and use atan?

I will look more at how Maple does this.

Michael Godfrey <godfrey>
Group Member
Wed 26 Sep 2012 05:23:44 PM UTC, comment #12: 

That makes sense to me.  After a bit of thought,
I will add this behaviour to the documentation.

Michael Godfrey <godfrey>
Group Member
Wed 26 Sep 2012 05:02:49 PM UTC, comment #11: 

I don't think we need to adhere to the GNU spec.  In particular, I think it would be better to follow what mathematicians do and if they extend the range of the function with branch cuts or whatever then we can too.  To my simple mind, the acos() and cos() functions are inverses and that is the property I expect to hold at all times.  The domain of the function is a secondary concern to me.

In this case, the inverse property is maintained even though the domain has been extended.


x = acos (-2)
x =  3.1416 - 1.3170i
cos (x)
ans = -2.0000e+00 + 2.1211e-16i



Rik <rik5>
Group administrator
Wed 26 Sep 2012 04:00:15 PM UTC, comment #10: 

We are not consistent.  We use M_PI in several places, but it is defined only in liboctave/numeric/lo-specfun.cc:50.  So it would appear that some places, say quadcc.cc, may be using another definition.  Of course, maybe it is being defined in a system header file.  I was only grepping in the Octave tree.  See the attached file for a list of occurrences.

(file #26645)

Rik <rik5>
Group administrator
Wed 26 Sep 2012 03:35:00 PM UTC, comment #9: 

The GNU libc doc says:

The math library normally defines M_PI to a double approximation of pi. If strict ISO and/or POSIX compliance are requested this constant is not defined, but you can easily define it yourself:

     #define M_PI 3.14159265358979323846264338327
==========================

Do we do this?

Also, the GNU doc says:
            The arc cosine function is defined mathematically only over the domain -1 to 1. If x is outside the domain, acos signals a domain error.

And similar for other inverse functions.  Octave currently does not return an error.  Instead it returns a complex result.
This is also Matlab behaviour.  Should we adhere to the GNU
spec?


Michael Godfrey <godfrey>
Group Member
Wed 26 Sep 2012 03:24:49 PM UTC, comment #8: 

Rik,

Thanks!  I was just doing a build.
I have some more to add, including what
I can find out about the differences between
libm and the ISO/EDC spec.

Michael Godfrey <godfrey>
Group Member
Wed 26 Sep 2012 03:10:46 PM UTC, comment #7: 

I pushed a changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/46dd555edd33) for a malformed @tex construct that was stopping the documentation build.  You'll want to do an 'hg pull; hg update' before continuing with the patch.

Rik <rik5>
Group administrator
Wed 26 Sep 2012 02:17:34 PM UTC, comment #6: 

Well, OK. I had not even compiled it.
But, if you could push the attached minor
update,then I will do the rest as one
more final patch.

Thanks

(file #26644)

Michael Godfrey <godfrey>
Group Member
Wed 26 Sep 2012 12:44:03 PM UTC, comment #5: 

I have pushed your patch:

http://hg.savannah.gnu.org/hgweb/octave/rev/0b6c29cb53d0

It looks like useful information to me. Thanks. I am marking this bug "in progress" since it looks like you want to keep working on it.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Wed 26 Sep 2012 01:59:10 AM UTC, comment #4: 

Attached is a first draft of updates to the
trig functions doc.

I need to add indexing terms, and it would be
good to have more cross-references...

Comments?

(file #26638)

Michael Godfrey <godfrey>
Group Member
Tue 25 Sep 2012 09:46:08 PM UTC, comment #3: 

Oops, 2 typos:

— atan2(y, ±0) returns -pi/2 for y < 0. both -pi/2 *
— atan2(y, ±0) returns pi/2 for y > 0. both pi/2 *

should be:
— atan2(y, ±0) returns -pi/2 for y < 0. both -pi/2
— atan2(y, ±0) returns pi/2 for y > 0. both pi/2

It appears that Matlab does not recognize -0.
Someone might check this with newer Matlabs.

Michael Godfrey <godfrey>
Group Member
Tue 25 Sep 2012 09:37:08 PM UTC, comment #2: 

matlab r2009b differs from Octave and ISO/IEC as follows:
(*  in the right hand column means Matlab differs)

F.9.1.4 The atan2 functions     Matlab r2009B Differs from Octave
— atan2(±0, -0) returns ±pi.            both 0           *
— atan2(±0, +0) returns ±0.             both 0           *
— atan2(±0, x) returns ±pi for x < 0.    +- pi
— atan2(±0, x) returns ±0 for x > 0.    both 0           *
— atan2(y, ±0) returns -pi/2 for y < 0. both -pi/2       *
— atan2(y, ±0) returns pi/2 for y > 0.  both  pi/2       *
— atan2(±y, -Inf) returns ±pi for finite y > 0.  +- pi
— atan2(±y, +Inf) returns ±0 for finite y > 0.
                                        both 0           *
— atan2(±Inf, x) returns ±pi/2 for finite x. +- pi/2
— atan2(±Inf, -Inf) returns ±3pi/4.          +- 3pi/4
— atan2(±Inf, +Inf) returns ±pi/4.           +- pi/4

Octave adhers to ISO/IEC.  By far the most serious Matlab
issue is atan2(±0, -0)

I will work on documentation patch for this and other trig
functions.

Michael Godfrey <godfrey>
Group Member
Sun 23 Sep 2012 06:09:36 PM UTC, comment #1: 

There was a huge discussion about atan2 and how to handle float values on the mailing list some years back.  It seems that Apple has a different interpretation of the ISO C standards than everyone else. 

But for doubles, the situation is clearer.  Octave picks up the C atan2 function and the underlying C library should be C99 compliant.  In that case, Octave's behavior is defined by the ISO/IEC 9899 standard.  First, that means that atan2 is defined over the range [-pi,pi] inclusive.  Second, the behavior of special values is mandated in section Annex F.9.1.4.  I'll quote them for ease.  I don't know if it is worth documenting them all inside Octave or whether you want to write a patch that refers to the ISO standard.

F.9.1.4 The atan2 functions
— atan2(±0, -0) returns ±pi.
— atan2(±0, +0) returns ±0.
— atan2(±0, x) returns ±pi for x < 0.
— atan2(±0, x) returns ±0 for x > 0.
— atan2(y, ±0) returns -pi/2 for y < 0.
— atan2(y, ±0) returns pi/2 for y > 0.
— atan2(±y, -Inf) returns ±pi for finite y > 0.
— atan2(±y, +Inf) returns ±0 for finite y > 0.
— atan2(±Inf, x) returns ±pi/2 for finite x.
— atan2(±Inf, -Inf) returns ±3pi/4.
— atan2(±Inf, +Inf) returns ±pi/4.

Rik <rik5>
Group administrator
Tue 18 Sep 2012 12:56:24 AM UTC, original submission:  

I have been using atan2(y,x) and finally found a problem that I
might have guessed at sooner.  But on checking the documentation
I found that there is practically none.

It appears that the atan2 function that is used is compliant with
some standards.  But, a few of the cases that caused me
trouble were not included anywhere that I could find.  I was not
aware that: atan2(0, -0)  =>  pi  for example.

So, some checking and documentation is in order.

The attached file shows what I found.  If anyone knows of
better documentation that would be great.  Otherwise the Wikipedia
page plus some additions from the attached file may work

Comments?

Michael Godfrey <godfrey>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #26645:  M_PI.grep added by rik5 (5KiB - application/octet-stream)
file #26644:  trig_doc_02.diff added by godfrey (987B - text/x-patch)
file #26638:  trig_doc_01.diff added by godfrey (1KiB - text/x-patch)
file #26568:  atan2_def added by godfrey (3KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by godfrey (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 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-01-26 mtmiller StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2014-01-24 mtmiller StatusIn Progress Need Info
    2012-09-26 rik5 Attached File- Added M_PI.grep, #26645
    2012-09-26 godfrey Attached File- Added trig_doc_02.diff, #26644
    2012-09-26 jordigh StatusNone In Progress
    2012-09-26 godfrey Attached File- Added trig_doc_01.diff, #26638
    2012-09-18 godfrey Attached File- Added atan2_def, #26568

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code