Mercurial > hg > octave-nkf
changeset 19434:ba7e42dea4b2
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 +++ b/NEWS @@ -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 +++ b/libinterp/corefcn/mappers.cc @@ -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 +++ 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