Mercurial > octave
changeset 21945:e9765b62d4e8
Use C++11 standard complex trig functions when available (bug #44310, bug #45507)
* acinclude.m4 (OCTAVE_CHECK_FUNC_COMPLEX): New macro.
* configure.ac: Use it to test for complex variants of trig and hyperbolic
functions.
* lo-mappers.cc (octave::math::acos, octave::math::asin, octave::math::atan):
Use standard library functions when available.
* lo-specfun.cc (octave::math::acosh, octave::math::asinh, octave::math::atanh):
Likewise.
* mappers.cc: Add BIST tests for affected functions.
| author | Mike Miller <mtmiller@octave.org> |
|---|---|
| date | Fri, 17 Jun 2016 14:24:27 -0700 |
| parents | 4b2eab5d2a6a |
| children | 835d070ede9c |
| files | configure.ac libinterp/corefcn/mappers.cc liboctave/numeric/lo-mappers.cc liboctave/numeric/lo-specfun.cc m4/acinclude.m4 |
| diffstat | 5 files changed, 138 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/configure.ac Fri Jun 17 16:51:09 2016 -0400 +++ b/configure.ac Fri Jun 17 14:24:27 2016 -0700 @@ -2578,6 +2578,15 @@ ;; esac +## Look in <complex> for C++ variants of math functions that we need. + +OCTAVE_CHECK_FUNC_COMPLEX(acos) +OCTAVE_CHECK_FUNC_COMPLEX(acosh) +OCTAVE_CHECK_FUNC_COMPLEX(asin) +OCTAVE_CHECK_FUNC_COMPLEX(asinh) +OCTAVE_CHECK_FUNC_COMPLEX(atan) +OCTAVE_CHECK_FUNC_COMPLEX(atanh) + ## Check for nonstandard, but common math functions, that we need. dnl Use multiple AC_CHECKs to avoid line continuations '\' in list
--- a/libinterp/corefcn/mappers.cc Fri Jun 17 16:51:09 2016 -0400 +++ b/libinterp/corefcn/mappers.cc Fri Jun 17 14:24:27 2016 -0700 @@ -125,6 +125,11 @@ %! assert (acos ([2 0]), [ival*i, pi/2], 2*eps); %! assert (acos ([2 0i]), [ival*i, pi/2], 2*eps); +## Test large magnitude arguments (bug #45507) +%! x = [1, -1, i, -i] .* 1e150; +%! v = [0, pi, pi/2, pi/2]; +%! assert (real (acos (x)), v); + %!error acos () %!error acos (1, 2) */ @@ -164,6 +169,12 @@ %! assert (acosh (single (10i)), re + i*im, 5*eps ("single")); %! assert (acosh (single (-10i)), re - i*im, 5*eps ("single")); +## Test large magnitude arguments (bug #45507) +%!test +%! x = [1, -1, i, -i] .* 1e150; +%! v = [0, pi, pi/2, -pi/2]; +%! assert (imag (acosh (x)), v); + %!error acosh () %!error acosh (1, 2) */ @@ -277,6 +288,11 @@ %! assert (asin ([2 0]), [rval - ival*i, 0], 2*eps); %! assert (asin ([2 0i]), [rval - ival*i, 0], 2*eps); +## Test large magnitude arguments (bug #45507) +%! x = [1, -1, i, -i] .* 1e150; +%! v = [pi/2, -pi/2, 0, -0]; +%! assert (real (asin (x)), v); + %!error asin () %!error asin (1, 2) */ @@ -305,6 +321,11 @@ %! x = single ([0, i, 0, -i]); %! assert (asinh (x), v, sqrt (eps ("single"))); +## Test large magnitude arguments (bug #45507) +%! x = [1, -1, i, -i] .* 1e150; +%! v = [0, 0, pi/2, -pi/2]; +%! assert (imag (asinh (x)), v); + %!error asinh () %!error asinh (1, 2) */ @@ -337,6 +358,12 @@ %! x = single ([0, rt3/3, 1, rt3, -rt3, -1, -rt3/3, 0]); %! assert (atan (x), v, sqrt (eps ("single"))); +## Test large magnitude arguments (bug #44310, bug #45507) +%! x = [1, -1, i, -i] .* 1e150; +%! v = [pi/2, -pi/2, pi/2, -pi/2]; +%! assert (real (atan (x)), v); +%! assert (imag (atan (x)), [0, 0, 0, 0], eps); + %!error atan () %!error atan (1, 2) */ @@ -365,6 +392,12 @@ %! x = single ([0, 0]); %! assert (atanh (x), v, sqrt (eps ("single"))); +## Test large magnitude arguments (bug #44310, bug #45507) +%! x = [1, -1, i, -i] .* 1e150; +%! v = [pi/2, pi/2, pi/2, -pi/2]; +%! assert (imag (atanh (x)), v); +%! assert (real (atanh (x)), [0, 0, 0, 0], eps); + %!error atanh () %!error atanh (1, 2) */
--- a/liboctave/numeric/lo-mappers.cc Fri Jun 17 16:51:09 2016 -0400 +++ b/liboctave/numeric/lo-mappers.cc Fri Jun 17 14:24:27 2016 -0700 @@ -88,6 +88,14 @@ Complex acos (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ACOS) + Complex y = std::acos (x); + + if (imag (x) == 0.0 && real (x) > 1.0) + return conj (y); + else + return y; +#else static Complex i (0, 1); Complex tmp; @@ -104,11 +112,20 @@ tmp = 1.0 - x*x; return -i * log (x + i * sqrt (tmp)); +#endif } FloatComplex acos (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ACOS) + FloatComplex y = std::acos (x); + + if (imag (x) == 0.0f && real (x) > 1.0f) + return conj (y); + else + return y; +#else static FloatComplex i (0, 1); FloatComplex tmp; @@ -125,11 +142,20 @@ tmp = 1.0f - x*x; return -i * log (x + i * sqrt (tmp)); +#endif } Complex asin (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ASIN) + Complex y = std::asin (x); + + if (imag (x) == 0.0 && real (x) > 1.0) + return conj (y); + else + return y; +#else static Complex i (0, 1); Complex tmp; @@ -146,11 +172,20 @@ tmp = 1.0 - x*x; return -i * log (i*x + sqrt (tmp)); +#endif } FloatComplex asin (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ASIN) + FloatComplex y = std::asin (x); + + if (imag (x) == 0.0f && real (x) > 1.0f) + return conj (y); + else + return y; +#else static FloatComplex i (0, 1); FloatComplex tmp; @@ -167,22 +202,31 @@ tmp = 1.0f - x*x; return -i * log (i*x + sqrt (tmp)); +#endif } Complex atan (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ATAN) + return std::atan (x); +#else static Complex i (0, 1); return i * log ((i + x) / (i - x)) / 2.0; +#endif } FloatComplex atan (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ATAN) + return std::atan (x); +#else static FloatComplex i (0, 1); return i * log ((i + x) / (i - x)) / 2.0f; +#endif } double log2 (double x) { return std::log2 (x); }
--- a/liboctave/numeric/lo-specfun.cc Fri Jun 17 16:51:09 2016 -0400 +++ b/liboctave/numeric/lo-specfun.cc Fri Jun 17 14:24:27 2016 -0700 @@ -229,13 +229,21 @@ Complex acosh (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ACOSH) + return std::acosh (x); +#else return log (x + sqrt (x + 1.0) * sqrt (x - 1.0)); +#endif } FloatComplex acosh (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ACOSH) + return std::acosh (x); +#else return log (x + sqrt (x + 1.0f) * sqrt (x - 1.0f)); +#endif } double @@ -265,13 +273,21 @@ Complex asinh (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ASINH) + return std::asinh (x); +#else return log (x + sqrt (x*x + 1.0)); +#endif } FloatComplex asinh (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ASINH) + return std::asinh (x); +#else return log (x + sqrt (x*x + 1.0f)); +#endif } double @@ -301,13 +317,21 @@ Complex atanh (const Complex& x) { +#if defined (HAVE_COMPLEX_STD_ATANH) + return std::atanh (x); +#else return log ((1.0 + x) / (1.0 - x)) / 2.0; +#endif } FloatComplex atanh (const FloatComplex& x) { +#if defined (HAVE_COMPLEX_STD_ATANH) + return std::atanh (x); +#else return log ((1.0f + x) / (1.0f - x)) / 2.0f; +#endif } double
--- a/m4/acinclude.m4 Fri Jun 17 16:51:09 2016 -0400 +++ b/m4/acinclude.m4 Fri Jun 17 14:24:27 2016 -0700 @@ -320,6 +320,34 @@ fi ]) dnl +dnl Check whether a complex-valued function is available in <complex>. +dnl Will define HAVE_COMPLEX_STD_FUNC if the function is available in the +dnl std namespace and is callable on both std::complex<double> and +dnl std::complex<float>. The return type of the function is expected to +dnl be of the same std::complex<T> type. +dnl +AC_DEFUN([OCTAVE_CHECK_FUNC_COMPLEX], [ + ac_safe=`echo "$1" | $SED 'y% ./+-:=%___p___%'` + + AC_CACHE_CHECK([for std::$1 in <complex>], + [octave_cv_func_complex_std_$ac_safe], + [AC_LANG_PUSH(C++) + AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ + #include <complex> + ]], [[ + std::complex<double> z = std::$1 (std::complex<double> (1.0, 1.0)); + std::complex<float> c = std::$1 (std::complex<float> (1.0, 1.0)); + ]])], + [eval "octave_cv_func_complex_std_$ac_safe=yes"], + [eval "octave_cv_func_complex_std_$ac_safe=no"]) + AC_LANG_POP(C++) + ]) + if eval "test \"`echo '$octave_cv_func_complex_std_'$ac_safe`\" = yes"; then + AC_DEFINE(AS_TR_CPP([[HAVE_COMPLEX_STD_][$1]]), 1, + [Define to 1 if <complex> provides std::$1(std::complex<T>).]) + fi +]) +dnl dnl Check whether Qscintilla has version 2.6.0 or later dnl FIXME: This test uses a version number. It potentially could dnl be re-written to actually call the function, but is it worth it?