Mercurial > octave-nkf

changeset 19276:ba7e42dea4b2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix returned phase of asin, acos outside principal branch (bug #43349). * NEWS: Announce change in phase for inputs outside the principal branch [-1,1]. * mappers.cc (Fasin, Facos): Add BIST tests to check values outside principal branch. * lo-mappers.cc (acos (const Complex& x), asin (const Complex& x), acos (const FloatComplex& x),asin (const FloatComplex& x): Update double and float versions of asin, acos to check if input is real (imag == 0). If so, avoid intermediate calculation that produces -0i which affects phase (but not magnitude) of result.
author Rik <rik@octave.org>
date Thu, 09 Oct 2014 16:14:46 -0700
parents 64dc954cf06c
children 91a6f06c5052
files NEWS libinterp/corefcn/mappers.cc liboctave/numeric/lo-mappers.cc
diffstat 3 files changed, 99 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS	Tue Oct 07 22:37:39 2014 +0530
+++ b/NEWS	Thu Oct 09 16:14:46 2014 -0700
@@ -31,6 +31,13 @@
     triangulation is always returned.  The delaunay3 function which
     handles 3-D inputs has been deprecated in favor of delaunay.
 
+ ** The trigonometric functions asin and acos return different phase values
+    from previous versions of Octave when the input is outside the principal
+    branch ([-1, 1]).  If the real portion of the input is greater than 1 then
+    the limit from below is taken.  If the real portion is less than 1 then the
+    limit from above is taken.  This criteria is consistent with several other
+    numerical analysis software packages.
+
  ** Integer formats used in the printf family of functions now work for
     64-bit integers and are more compatible with Matlab when printing
     non-integer values.  Now instead of truncating, Octave will switch
--- a/libinterp/corefcn/mappers.cc	Tue Oct 07 22:37:39 2014 +0530
+++ b/libinterp/corefcn/mappers.cc	Thu Oct 09 16:14:46 2014 -0700
@@ -114,6 +114,20 @@
 %! v = single ([0, pi/6, pi/4, pi/3, pi/2, 2*pi/3, 3*pi/4, 5*pi/6, pi]);
 %! assert (acos (x), v, sqrt (eps ("single")));
 
+## Test values on either side of branch cut
+%!test
+%! rval = 0;
+%! ival = 1.31695789692481635;
+%! obs = acos ([2, 2-i*eps, 2+i*eps]);
+%! exp = [rval + ival*i, rval + ival*i, rval - ival*i];
+%! assert (obs, exp, 2*eps);
+%! rval = pi;
+%! obs = acos ([-2, -2-i*eps, -2+i*eps]);
+%! exp = [rval - ival*i, rval + ival*i, rval - ival*i];
+%! assert (obs, exp, 2*eps);
+%! assert (acos ([2 0]),  [ival*i, pi/2], 2*eps);
+%! assert (acos ([2 0i]), [ival*i, pi/2], 2*eps);
+
 %!error acos ()
 %!error acos (1, 2)
 */
@@ -236,12 +250,32 @@
 }
 
 /*
-%!test
+%!shared rt2, rt3
 %! rt2 = sqrt (2);
 %! rt3 = sqrt (3);
+
+%!test
 %! x = [0, 1/2, rt2/2, rt3/2, 1, rt3/2, rt2/2, 1/2, 0];
 %! v = [0, pi/6, pi/4, pi/3, pi/2, pi/3, pi/4, pi/6, 0];
-%! assert (all (abs (asin (x) - v) < sqrt (eps)));
+%! assert (asin (x), v, sqrt (eps));
+
+%!test
+%! x = single ([0, 1/2, rt2/2, rt3/2, 1, rt3/2, rt2/2, 1/2, 0]);
+%! v = single ([0, pi/6, pi/4, pi/3, pi/2, pi/3, pi/4, pi/6, 0]);
+%! assert (asin (x), v, sqrt (eps ("single")));
+
+## Test values on either side of branch cut
+%!test
+%! rval = pi/2;
+%! ival = 1.31695789692481635;
+%! obs = asin ([2, 2-i*eps, 2+i*eps]);
+%! exp = [rval - ival*i, rval - ival*i, rval + ival*i];
+%! assert (obs, exp, 2*eps);
+%! obs = asin ([-2, -2-i*eps, -2+i*eps]);
+%! exp = [-rval + ival*i, -rval - ival*i, -rval + ival*i];
+%! assert (obs, exp, 2*eps);
+%! assert (asin ([2 0]),  [rval - ival*i, 0], 2*eps);
+%! assert (asin ([2 0i]), [rval - ival*i, 0], 2*eps);
 
 %!error asin ()
 %!error asin (1, 2)
--- a/liboctave/numeric/lo-mappers.cc	Tue Oct 07 22:37:39 2014 +0530
+++ b/liboctave/numeric/lo-mappers.cc	Thu Oct 09 16:14:46 2014 -0700
@@ -184,7 +184,20 @@
 {
   static Complex i (0, 1);
 
-  return -i * (log (x + i * (sqrt (1.0 - x*x))));
+  Complex tmp;
+
+  if (imag (x) == 0.0)
+    {
+      // If the imaginary part of X is 0, then avoid generating an
+      // imaginary part of -0 for the expression 1-x*x.
+      // This effectively chooses the same phase of the branch cut as Matlab.
+      double xr = real (x);
+      tmp = Complex (1.0 - xr*xr);
+    }
+  else
+    tmp = 1.0 - x*x;
+
+  return -i * log (x + i * sqrt (tmp));
 }
 
 Complex
@@ -198,7 +211,20 @@
 {
   static Complex i (0, 1);
 
-  return -i * log (i*x + sqrt (1.0 - x*x));
+  Complex tmp;
+
+  if (imag (x) == 0.0)
+    {
+      // If the imaginary part of X is 0, then avoid generating an
+      // imaginary part of -0 for the expression 1-x*x.
+      // This effectively chooses the same phase of the branch cut as Matlab.
+      double xr = real (x);
+      tmp = Complex (1.0 - xr*xr);
+    }
+  else
+    tmp = 1.0 - x*x;
+
+  return -i * log (i*x + sqrt (tmp));
 }
 
 Complex
@@ -401,7 +427,20 @@
 {
   static FloatComplex i (0, 1);
 
-  return -i * (log (x + i * (sqrt (static_cast<float>(1.0) - x*x))));
+  FloatComplex tmp;
+
+  if (imag (x) == 0.0f)
+    {
+      // If the imaginary part of X is 0, then avoid generating an
+      // imaginary part of -0 for the expression 1-x*x.
+      // This effectively chooses the same phase of the branch cut as Matlab.
+      float xr = real (x);
+      tmp = FloatComplex (1.0f - xr*xr);
+    }
+  else
+    tmp = 1.0f - x*x;
+
+  return -i * log (x + i * sqrt (tmp));
 }
 
 FloatComplex
@@ -415,7 +454,20 @@
 {
   static FloatComplex i (0, 1);
 
-  return -i * log (i*x + sqrt (static_cast<float>(1.0) - x*x));
+  FloatComplex tmp;
+
+  if (imag (x) == 0.0f)
+    {
+      // If the imaginary part of X is 0, then avoid generating an
+      // imaginary part of -0 for the expression 1-x*x.
+      // This effectively chooses the same phase of the branch cut as Matlab.
+      float xr = real (x);
+      tmp = FloatComplex (1.0f - xr*xr);
+    }
+  else
+    tmp = 1.0f - x*x;
+
+  return -i * log (i*x + sqrt (tmp));
 }
 
 FloatComplex