Mercurial > octave
view liboctave/numeric/lo-mappers.h @ 31205:b0e90ca8e679 stable
quad2d: fix unintended complex conjugate return (bug #62972)
quad2d: use .' instead of ' to avoid complex conjugate in q and qerr.
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Sun, 28 Aug 2022 12:21:17 -0400 |
parents | 796f54d4ddbf |
children | e88a07dec498 |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1996-2022 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// #if ! defined (octave_lo_mappers_h) #define octave_lo_mappers_h 1 #include "octave-config.h" #include <cmath> #include <limits> #include "lo-ieee.h" #include "oct-cmplx.h" #include "oct-inttypes-fwd.h" namespace octave { namespace math { extern OCTAVE_API bool isna (double x); extern OCTAVE_API bool isna (float x); extern OCTAVE_API bool isna (const Complex& x); extern OCTAVE_API bool isna (const FloatComplex& x); extern OCTAVE_API bool is_NaN_or_NA (const Complex& x); extern OCTAVE_API bool is_NaN_or_NA (const FloatComplex& x); inline double copysign (double x, double y) { return std::copysign (x, y); } inline float copysign (float x, float y) { return std::copysignf (x, y); } inline double signbit (double x) { return std::signbit (x); } inline float signbit (float x) { return std::signbit (x); } // Test for negative sign. extern OCTAVE_API bool negative_sign (double x); extern OCTAVE_API bool negative_sign (float x); // Test for positive sign. inline bool positive_sign (double x) { return ! negative_sign (x); } inline bool positive_sign (float x) { return ! negative_sign (x); } extern OCTAVE_API Complex acos (const Complex& x); extern OCTAVE_API FloatComplex acos (const FloatComplex& x); extern OCTAVE_API Complex asin (const Complex& x); extern OCTAVE_API FloatComplex asin (const FloatComplex& x); inline Complex atan (const Complex& x) { return std::atan (x); } inline FloatComplex atan (const FloatComplex& x) { return std::atan (x); } // The C++ standard would normally return a std::complex value for conj // even when the input is fully real. Octave overrides this. inline double conj (double x) { return x; } inline float conj (float x) { return x; } template <typename T> std::complex<T> conj (const std::complex<T>& x) { return std::conj (x); } inline double log2 (double x) { return std::log2 (x); } inline float log2 (float x) { return std::log2f (x); } extern OCTAVE_API Complex log2 (const Complex& x); extern OCTAVE_API FloatComplex log2 (const FloatComplex& x); extern OCTAVE_API double log2 (double x, int& exp); extern OCTAVE_API float log2 (float x, int& exp); extern OCTAVE_API Complex log2 (const Complex& x, int& exp); extern OCTAVE_API FloatComplex log2 (const FloatComplex& x, int& exp); inline double exp2 (double x) { return std::exp2 (x); } inline float exp2 (float x) { return std::exp2f (x); } template <typename T> std::complex<T> ceil (const std::complex<T>& x) { return std::complex<T> (std::ceil (std::real (x)), std::ceil (std::imag (x))); } template <typename T> std::complex<T> trunc (const std::complex<T>& x) { return std::complex<T> (std::trunc (std::real (x)), std::trunc (std::imag (x))); } // Provide alias for trunc under the more familiar name of fix. inline double fix (double x) { return std::trunc (x); } inline float fix (float x) { return std::trunc (x); } template <typename T> std::complex<T> fix (const std::complex<T>& x) { return trunc (x); } template <typename T> std::complex<T> floor (const std::complex<T>& x) { return std::complex<T> (std::floor (std::real (x)), std::floor (std::imag (x))); } inline double round (double x) { return std::round (x); } inline float round (float x) { return std::roundf (x); } template <typename T> std::complex<T> round (const std::complex<T>& x) { return std::complex<T> (round (std::real (x)), round (std::imag (x))); } inline double roundb (double x) { double t = round (x); if (fabs (x - t) == 0.5) t = 2 * std::trunc (0.5 * t); return t; } inline float roundb (float x) { float t = round (x); if (fabsf (x - t) == 0.5f) t = 2 * std::trunc (0.5f * t); return t; } template <typename T> std::complex<T> roundb (const std::complex<T>& x) { return std::complex<T> (roundb (std::real (x)), roundb (std::imag (x))); } extern OCTAVE_API double frexp (double x, int *expptr); extern OCTAVE_API float frexp (float x, int *expptr); inline bool isnan (bool) { return false; } inline bool isnan (char) { return false; } inline bool isnan (double x) { return std::isnan (x); } inline bool isnan (float x) { return std::isnan (x); } // FIXME: Do we need the isnan overload for complex? template <typename T> bool isnan (const std::complex<T>& x) { return (isnan (std::real (x)) || isnan (std::imag (x))); } inline bool isfinite (double x) { return std::isfinite (x); } inline bool isfinite (float x) { return std::isfinite (x); } // FIXME: Do we need isfinite overload for complex? template <typename T> bool isfinite (const std::complex<T>& x) { return (isfinite (std::real (x)) && isfinite (std::imag (x))); } inline bool isinf (double x) { return std::isinf (x); } inline bool isinf (float x) { return std::isinf (x); } template <typename T> bool isinf (const octave_int<T>&) { return false; } // FIXME: Do we need isinf overload for complex? template <typename T> bool isinf (const std::complex<T>& x) { return (isinf (std::real (x)) || isinf (std::imag (x))); } // Some useful tests, that are commonly repeated. // Test for a finite integer. // FIXME: Benchmark whether trunc might be faster than round here. inline bool isinteger (double x) { return isfinite (x) && x == round (x); } inline bool isinteger (float x) { return isfinite (x) && x == round (x); } inline double signum (double x) { double tmp = 0.0; if (x < 0.0) tmp = -1.0; else if (x > 0.0) tmp = 1.0; return isnan (x) ? numeric_limits<double>::NaN () : tmp; } inline float signum (float x) { float tmp = 0.0f; if (x < 0.0f) tmp = -1.0f; else if (x > 0.0f) tmp = 1.0f; return isnan (x) ? numeric_limits<float>::NaN () : tmp; } template <typename T> std::complex<T> signum (const std::complex<T>& x) { T tmp = abs (x); return tmp == 0 ? 0.0 : x / tmp; } // Convert X to the nearest integer value. Should not pass NaN to // this function. // For integer types? Hmm. Need to be sure T is an integer type... template <typename T> T x_nint (T x) { return x; } template <> inline double x_nint (double x) { return (isfinite (x) ? std::floor (x + 0.5) : x); } template <> inline float x_nint (float x) { return (isfinite (x) ? std::floor (x + 0.5f) : x); } extern OCTAVE_API octave_idx_type nint_big (double x); extern OCTAVE_API octave_idx_type nint_big (float x); extern OCTAVE_API int nint (double x); extern OCTAVE_API int nint (float x); template <typename T> T mod (T x, T y) { T retval; if (y == 0) retval = x; else { T q = x / y; if (x_nint (y) != y && (std::abs ((q - x_nint (q)) / x_nint (q)) < std::numeric_limits<T>::epsilon ())) retval = 0; else { T n = std::floor (q); // Prevent use of extra precision. volatile T tmp = y * n; retval = x - tmp; } } if (x != y && y != 0) retval = copysign (retval, y); return retval; } template <typename T> T rem (T x, T y) { T retval; if (y == 0) retval = numeric_limits<T>::NaN (); else { T q = x / y; if (x_nint (y) != y && (std::abs ((q - x_nint (q)) / x_nint (q)) < std::numeric_limits<T>::epsilon ())) retval = 0; else { T n = std::trunc (q); // Prevent use of extra precision. volatile T tmp = y * n; retval = x - tmp; } } if (x != y && y != 0) retval = copysign (retval, x); return retval; } // Generic min, max definitions template <typename T> T min (T x, T y) { return x <= y ? x : y; } template <typename T> T max (T x, T y) { return x >= y ? x : y; } // This form is favorable. GCC will translate (x <= y ? x : y) without a // jump, hence the only conditional jump involved will be the first // (isnan), infrequent and hence friendly to branch prediction. inline double min (double x, double y) { return isnan (y) ? x : (x <= y ? x : y); } inline double max (double x, double y) { return isnan (y) ? x : (x >= y ? x : y); } inline float min (float x, float y) { return isnan (y) ? x : (x <= y ? x : y); } inline float max (float x, float y) { return isnan (y) ? x : (x >= y ? x : y); } inline std::complex<double> min (const std::complex<double>& x, const std::complex<double>& y) { return abs (x) <= abs (y) ? x : (isnan (x) ? x : y); } inline std::complex<float> min (const std::complex<float>& x, const std::complex<float>& y) { return abs (x) <= abs (y) ? x : (isnan (x) ? x : y); } inline std::complex<double> max (const std::complex<double>& x, const std::complex<double>& y) { return abs (x) >= abs (y) ? x : (isnan (x) ? x : y); } inline std::complex<float> max (const std::complex<float>& x, const std::complex<float>& y) { return abs (x) >= abs (y) ? x : (isnan (x) ? x : y); } template <typename T> inline octave_int<T> min (const octave_int<T>& x, const octave_int<T>& y) { return xmin (x, y); } template <typename T> inline octave_int<T> max (const octave_int<T>& x, const octave_int<T>& y) { return xmax (x, y); } // These map reals to Complex. extern OCTAVE_API Complex rc_acos (double); extern OCTAVE_API FloatComplex rc_acos (float); extern OCTAVE_API Complex rc_acosh (double); extern OCTAVE_API FloatComplex rc_acosh (float); extern OCTAVE_API Complex rc_asin (double); extern OCTAVE_API FloatComplex rc_asin (float); extern OCTAVE_API Complex rc_atanh (double); extern OCTAVE_API FloatComplex rc_atanh (float); extern OCTAVE_API Complex rc_log (double); extern OCTAVE_API FloatComplex rc_log (float); extern OCTAVE_API Complex rc_log2 (double); extern OCTAVE_API FloatComplex rc_log2 (float); extern OCTAVE_API Complex rc_log10 (double); extern OCTAVE_API FloatComplex rc_log10 (float); extern OCTAVE_API Complex rc_sqrt (double); extern OCTAVE_API FloatComplex rc_sqrt (float); } } #endif