Mercurial > hg > octave-nkf

diff liboctave/numeric/lo-mappers.cc @ 19434: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 161ebb78ac1b
children 5b263e517c95
line wrap: on
line diff
--- a/liboctave/numeric/lo-mappers.cc
+++ b/liboctave/numeric/lo-mappers.cc
@@ -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