From 64b7e52987665dcdc317b63f54fe89162c8d80be Mon Sep 17 00:00:00 2001 From: Vadim Kuznetsov Date: Sun, 1 Sep 2024 13:40:16 +0300 Subject: [PATCH 1/4] Switch real.cpp to C++11 implmentation --- src/math/real.cpp | 66 +++-------------------------------------------- src/math/real.h | 1 - 2 files changed, 4 insertions(+), 63 deletions(-) diff --git a/src/math/real.cpp b/src/math/real.cpp index 521c58a1..a6abbeb0 100644 --- a/src/math/real.cpp +++ b/src/math/real.cpp @@ -130,14 +130,7 @@ nr_double_t tanh (const nr_double_t arg) { \return arc hyperbolic cosine of z */ nr_double_t acosh (const nr_double_t arg) { -#ifdef HAVE_STD_ACOSH - // c++11 return std::acosh (arg); -#elif HAVE_ACOSH - return ::acosh (arg); -#else - return log (arg + sqrt (arg * arg - 1.0)); -#endif } /*! \brief Compute arc hyperbolic sine @@ -146,14 +139,7 @@ nr_double_t acosh (const nr_double_t arg) { */ nr_double_t asinh (const nr_double_t arg) { -#ifdef HAVE_STD_ASINH - // c++11 return std::asinh (arg); -#elif HAVE_ASINH - return ::asinh (arg); -#else - return log (arg + sqrt (arg * arg + 1.0)); -#endif } /*! \brief Compute arc hyperbolic tangent @@ -162,14 +148,7 @@ nr_double_t asinh (const nr_double_t arg) */ nr_double_t atanh (const nr_double_t arg) { -#ifdef HAVE_STD_ATANH - // c++11 return std::atanh (arg); -#elif HAVE_ATANH - return ::atanh (arg); -#else - return 0.5 * log ( 2.0 / (1.0 - arg) - 1.0); -#endif } @@ -211,22 +190,7 @@ nr_double_t sqrt (const nr_double_t d) { \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$ */ nr_double_t xhypot (const nr_double_t a, const nr_double_t b) { -#ifdef HAVE_STD_HYPOT - return std::hypot(a,b) // c++11 -#else - nr_double_t c = fabs (a); - nr_double_t d = fabs (b); - if (c > d) { - nr_double_t e = d / c; - return c * sqrt (1 + e * e); - } - else if (d == 0) - return 0; - else { - nr_double_t e = c / d; - return d * sqrt (1 + e * e); - } -#endif + return std::hypot(a,b); // c++11 } // @@ -234,11 +198,7 @@ nr_double_t xhypot (const nr_double_t a, const nr_double_t b) { // nr_double_t erf( nr_double_t arg) { -#ifdef HAVE_STD_ERF return std::erf (arg); // c++11 -#elif HAVE_ERF - return ::erf (arg); -#endif } @@ -253,31 +213,13 @@ nr_double_t floor( nr_double_t arg) { return std::floor(arg); } -nr_double_t fmod( nr_double_t arg) { -#ifdef HAVE_STD_TRUNC - return std::fmod(arg); -#else - return fmod(arg); -#endif -} nr_double_t trunc( nr_double_t arg) { -#ifdef HAVE_STD_TRUNC return std::trunc(arg); -#elif HAVE_TRUNC - return ::trunc (arg); -#else - return arg > 0 ? floor (arg) : floor (arg + 1); -#endif } + nr_double_t round( nr_double_t arg) { -#ifdef HAVE_STD_ROUND return std::round(arg); -#elif HAVE_ROUND - return ::round (arg); -#else - return (arg > 0) ? floor (arg + 0.5) : ceil (arg - 0.5); -#endif } @@ -514,7 +456,7 @@ nr_double_t conj (const nr_double_t r) { * \return input in degree (x)*180/pi */ nr_double_t rad2deg (const nr_double_t x) { - return (180.0 * (x) / pi); + return (180.0 * (x) / pi); } /*! @@ -523,7 +465,7 @@ nr_double_t rad2deg (const nr_double_t x) { * \return input in radian (x)*pi/180 */ nr_double_t deg2rad (const nr_double_t x) { - return (pi * (x) / 180.0); + return (pi * (x) / 180.0); } diff --git a/src/math/real.h b/src/math/real.h index f64bb765..57ea4ddb 100644 --- a/src/math/real.h +++ b/src/math/real.h @@ -88,7 +88,6 @@ nr_double_t erf(const nr_double_t ); // nr_double_t ceil(const nr_double_t ); nr_double_t floor(const nr_double_t ); -nr_double_t fmod(const nr_double_t ); //FIXME nr_double_t trunc(const nr_double_t ); // c++11 nr_double_t round(const nr_double_t ); // c++11 From 154678f3921e603f06097f8354cf86342d97c829 Mon Sep 17 00:00:00 2001 From: Vadim Kuznetsov Date: Sun, 1 Sep 2024 18:57:09 +0300 Subject: [PATCH 2/4] Switch complex.cpp to C++11 implmentation --- src/math/complex.cpp | 45 -------------------------------------------- 1 file changed, 45 deletions(-) diff --git a/src/math/complex.cpp b/src/math/complex.cpp index 1749824b..4088e3a4 100644 --- a/src/math/complex.cpp +++ b/src/math/complex.cpp @@ -82,16 +82,7 @@ nr_complex_t tan (const nr_complex_t z) { \return arc cosine of z */ nr_complex_t acos (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ACOS return std::acos (z); -#else - // missing on OSX 10.6 - nr_double_t r = real (z); - nr_double_t i = imag (z); - nr_complex_t y = sqrt (z * z - 1.0); - if (r * i < 0.0) y = -y; - return nr_complex_t (0, -1.0) * log (z + y); -#endif } @@ -100,14 +91,7 @@ nr_complex_t acos (const nr_complex_t z) { \return arc sine of z */ nr_complex_t asin (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ASIN return std::asin (z); -#else - // missing on OSX 10.6 - nr_double_t r = real (z); - nr_double_t i = imag (z); - return nr_complex_t (0.0, -1.0) * log (nr_complex_t (-i, r) + sqrt (1.0 - z * z)); -#endif } /*! \brief Compute complex arc tangent @@ -115,12 +99,7 @@ nr_complex_t asin (const nr_complex_t z) { \return arc tangent of z */ nr_complex_t atan (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ATAN return std::atan (z); -#else - // missing on OSX 10.6 - return nr_complex_t (0.0, -0.5) * log (nr_complex_t (0, 2) / (z + nr_complex_t (0, 1)) - 1.0); -#endif } @@ -160,11 +139,7 @@ nr_complex_t tanh (const nr_complex_t z) { \return arc hyperbolic cosine of z */ nr_complex_t acosh (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ACOSH return std::acosh (z); -#else - return log (z + sqrt (z * z - 1.0)); -#endif } @@ -173,11 +148,7 @@ nr_complex_t acosh (const nr_complex_t z) { \return arc hyperbolic sine of z */ nr_complex_t asinh (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ACOSH return std::asinh (z); -#else - return log (z + sqrt (z * z + 1.0)); -#endif } @@ -186,11 +157,7 @@ nr_complex_t asinh (const nr_complex_t z) { \return arc hyperbolic tangent of z */ nr_complex_t atanh (const nr_complex_t z) { -#ifdef HAVE_CXX_COMPLEX_ATANH return std::atanh (z); -#else - return 0.5 * log ( 2.0 / (1.0 - z) - 1.0); -#endif } @@ -392,12 +359,8 @@ nr_complex_t atan2 (const nr_complex_t y, const nr_complex_t x) */ nr_complex_t log2 (const nr_complex_t z) { -#ifndef HAVE_CXX_COMPLEX_LOG2 nr_double_t phi = std::arg (z); return nr_complex_t (std::log (std::abs (z)) * log2e, phi * log2e); -#else - return std::log2 (z); -#endif } /*!\brief complex signum function @@ -550,11 +513,7 @@ nr_complex_t limexp (const nr_complex_t z) */ nr_complex_t polar (const nr_double_t mag, const nr_double_t ang ) { -#ifdef HAVE_CXX_COMPLEX_POLAR return std::polar (mag, ang); -#else - return nr_complex_t (mag * cos (ang), mag * sin (ang)); -#endif } /*!\brief Extension of polar construction to complex @@ -565,11 +524,7 @@ nr_complex_t polar (const nr_double_t mag, const nr_double_t ang ) */ nr_complex_t polar (const nr_complex_t a, const nr_complex_t p) { -#ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX - return std::polar (a, p); -#else return a * exp(nr_complex_t(0.0, 1.0) * p); -#endif } From 2dc1a13ba8e2b3048a86ddc208bc5d984dcedef2 Mon Sep 17 00:00:00 2001 From: Vadim Kuznetsov Date: Tue, 3 Sep 2024 13:15:38 +0300 Subject: [PATCH 3/4] Make real.cpp inline and move to real.h --- src/math/CMakeLists.txt | 2 +- src/math/real.h | 447 +++++++++++++++++++++++++++++++++++----- 2 files changed, 395 insertions(+), 54 deletions(-) diff --git a/src/math/CMakeLists.txt b/src/math/CMakeLists.txt index 1520b359..5cfbebb8 100644 --- a/src/math/CMakeLists.txt +++ b/src/math/CMakeLists.txt @@ -1,7 +1,7 @@ include_directories( ${CMAKE_CURRENT_SOURCE_DIR}) set(MATH_SRC # cbesselj.cpp - complex.cpp fspecial.cpp matrix.cpp real.cpp) + complex.cpp fspecial.cpp matrix.cpp ) set(HEADERS complex.h matrix.h precision.h real.h) diff --git a/src/math/real.h b/src/math/real.h index 57ea4ddb..88ef515c 100644 --- a/src/math/real.h +++ b/src/math/real.h @@ -27,118 +27,459 @@ #define __REAL_H__ #include +#include #include +#include "consts.h" /**It is preferred to add all used functions into the qucs namespace. * Doing so one is forced do think about compatibility instead of using std directly. * Inline is optional at this moment * \todo test if inline indeed performance improves (optimization flags should inline them anyway) */ - namespace qucs { // // trigonometric // -nr_double_t cos (const nr_double_t); -nr_double_t sin (const nr_double_t); -nr_double_t tan (const nr_double_t); -nr_double_t acos (const nr_double_t); -nr_double_t asin (const nr_double_t); -nr_double_t atan (const nr_double_t); -nr_double_t atan2 (const nr_double_t, const nr_double_t); //not used?, only for complex +/*! \brief Compute cosine of an angle + \param[in] z angle in radians + \return cosine of z +*/ +nr_double_t inline cos (const nr_double_t arg) { + return std::cos (arg); +} + +/*! \brief Compute sine of an angle + \param[in] z angle in radians + \return sine of z +*/ +nr_double_t inline sin (const nr_double_t arg) { + return std::sin (arg); +} + +/*! \brief Compute tangent of an angle + \param[in] z angle in radians + \return tangent of z +*/ +nr_double_t inline tan (const nr_double_t arg) { + return std::tan (arg); +} + +/*! \brief Compute arc cosine + \param[in] z arc + \return arc cosine of z +*/ +nr_double_t inline acos (const nr_double_t arg) { + return std::acos (arg); +} + +/*! \brief Compute arc sine + \param[in] z arc + \return arc sine of z +*/ +nr_double_t inline asin (const nr_double_t arg) { + return std::asin (arg); +} + +/*! \brief Compute arc tangent + \param[in] z arc + \return arc tangent of z +*/ +nr_double_t inline atan (const nr_double_t arg) { + return std::atan (arg); +} + +/*! \brief Compute arc tangent with two parameters (fortran like function) + \param[in] x proportion of x-coordinate + \param[in] y proportion of y-coordinate + \return principal value of the arc tangent of y/x, expressed in radians. +*/ +nr_double_t inline atan2 (const nr_double_t x, const nr_double_t y) { + return std::atan2 (x,y); +} // // hyperbolic // -nr_double_t cosh (const nr_double_t); -nr_double_t sinh (const nr_double_t); -nr_double_t tanh (const nr_double_t); -nr_double_t acosh (const nr_double_t); // c++11 -nr_double_t asinh (const nr_double_t); // c++11 -nr_double_t atanh (const nr_double_t); // c++11, not used?, only for complex + +/*! \brief Compute hyperbolic cosine + \param[in] z arc + \return hyperbolic cosine of z +*/ +nr_double_t inline cosh (const nr_double_t arg) { + return std::cosh (arg); +} + +/*! \brief Compute hyperbolic sine + \param[in] z arc + \return hyperbolic sine of z +*/ +nr_double_t inline sinh (const nr_double_t arg) { + return std::sinh (arg); +} + +/*! \brief Compute hyperbolic tangent + \param[in] z arc + \return hyperbolic tangent of z +*/ +nr_double_t inline tanh (const nr_double_t arg) { + return std::tanh (arg); +} + +/*! \brief Compute arc hyperbolic cosine + \param[in] z arc + \return arc hyperbolic cosine of z +*/ +nr_double_t inline acosh (const nr_double_t arg) { + return std::acosh (arg); +} + +/*! \brief Compute arc hyperbolic sine + \param[in] z arc + \return arc hyperbolic sine of z +*/ +nr_double_t inline asinh (const nr_double_t arg) +{ + return std::asinh (arg); +} + +/*! \brief Compute arc hyperbolic tangent + \param[in] z arc + \return arc hyperbolic tangent of z +*/ +nr_double_t inline atanh (const nr_double_t arg) +{ + return std::atanh (arg); +} // // exponential and logarithmic functions // -nr_double_t exp (const nr_double_t); -nr_double_t log (const nr_double_t); -nr_double_t log10 (const nr_double_t); - +nr_double_t inline exp (const nr_double_t arg) { + return std::exp (arg); +} +nr_double_t inline log (const nr_double_t arg) { + return std::log (arg); +} +nr_double_t inline log10 (const nr_double_t arg) { + return std::log10 (arg); +} // // power functions // -nr_double_t pow (const nr_double_t, const nr_double_t ); -nr_double_t sqrt (const nr_double_t ); -nr_double_t xhypot (const nr_double_t, const nr_double_t ); // same hypot in c++11? +nr_double_t inline pow (const nr_double_t a, const nr_double_t b) +{ + return std::pow (a,b); +} + +nr_double_t inline sqrt (const nr_double_t d) { + return std::sqrt (d); +} + +/*!\brief Euclidean distance function + + The xhypot() function returns \f$\sqrt{a^2+b^2}\f$. + This is the length of the hypotenuse of a right-angle triangle with sides + of length a and b, or the distance + of the point (a,b) from the origin. + + \param[in] a first length + \param[in] b second length + \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$ +*/ +nr_double_t inline xhypot (const nr_double_t a, const nr_double_t b) { + return std::hypot(a,b); // c++11 +} // // error functions // -nr_double_t erf(const nr_double_t ); + +nr_double_t inline erf( nr_double_t arg) { + return std::erf (arg); // c++11 +} // // rounding and remainder functions // -nr_double_t ceil(const nr_double_t ); -nr_double_t floor(const nr_double_t ); -nr_double_t trunc(const nr_double_t ); // c++11 -nr_double_t round(const nr_double_t ); // c++11 +nr_double_t inline ceil( nr_double_t arg) { + return std::ceil(arg); +} + +nr_double_t inline floor( nr_double_t arg) { + return std::floor(arg); +} + + +nr_double_t inline trunc( nr_double_t arg) { + return std::trunc(arg); +} + +nr_double_t inline round( nr_double_t arg) { + return std::round(arg); +} + // // Qucs extra trigonometric helper // -nr_double_t coth (const nr_double_t ); -nr_double_t sech (const nr_double_t ); -nr_double_t cosech (const nr_double_t ); +nr_double_t inline coth (const nr_double_t d) { + return 1.0 / std::tanh (d); +} + +nr_double_t inline sech (const nr_double_t d) { + return (1.0 / std::cosh (d)); +} + +nr_double_t inline cosech (const nr_double_t d) { + return (1.0 / std::sinh (d)); +} // // Qucs extra math functions // -nr_double_t sqr (const nr_double_t ); -unsigned int sqr (unsigned int); -nr_double_t quadr (const nr_double_t ); -nr_double_t rad2deg (const nr_double_t ); -nr_double_t deg2rad (const nr_double_t x ); +/*!\brief Square a value -/*!\brief Compute the third power of input */ -static inline nr_double_t cubic (const nr_double_t x) { return (x * x * x); } + \param[in] r Real number + \return \f$x^2\f$ +*/ +nr_double_t inline sqr (const nr_double_t r) { + return r * r; +} -/*!\brief Convert Celsius to Kelvin */ -static inline nr_double_t celsius2kelvin (const nr_double_t x) { return (x - K); } +unsigned int inline sqr (unsigned int r) { + return r * r; +} -/*!\brief Convert Kelvin to Celsius */ -static inline nr_double_t kelvin2celsius (const nr_double_t x) { return (x + K); } + + +/*!\brief Quartic function + + \param[in] r Real number + \return \f$x^4\f$ +*/ +nr_double_t inline quadr (const nr_double_t r) { + return r * r * r * r; +} // -// extra math functions +// extra math functions // -nr_double_t limexp (const nr_double_t); -nr_double_t signum (const nr_double_t); -nr_double_t sign (const nr_double_t); -nr_double_t sinc (const nr_double_t); -nr_double_t fix (const nr_double_t); -nr_double_t step (const nr_double_t); -unsigned int factorial (unsigned int); + +/*!\brief Compute limited exponential + + Compute limited exponential: + \f[ + \begin{cases} + \exp r & \text{if } r < \text{M\_LIMEXP} \\ + \exp (\text{M\_LIMEXP})\left[1.0 + (r - \text{M\_LIMEXP})\right] & + \text{else} + \end{cases} + \f] + + #limitexp is a constant + \param[in] r real number + \return limited exponential of r + \todo Change limexp(real) limexp(complex) file order + \todo Document #limitexp +*/ +nr_double_t inline limexp (const nr_double_t r) { + return r < limitexp ? exp (r) : exp (limitexp) * (1.0 + (r - limitexp)); +} + +/*!\brief real signum function + + compute \f[ + \mathrm{signum}\;d= + = \begin{cases} + O & \text{if } d=0 \\ + 1 & \text{if } d>0 \\ + -1 & \text{if } d<0 + \end{cases} + \f] + \param[in] d real number + \return signum of d + \todo Move near complex signum +*/ +nr_double_t inline signum (const nr_double_t d) { + if (d == 0) return 0; + return d < 0 ? -1 : 1; +} + +/*!\brief real sign function + + compute \f[ + \mathrm{sign}\;d= + = \begin{cases} + 1 & \text{if } d\ge 0 \\ + -1 & \text{if } d<0 + \end{cases} + \f] + \param[in] d real number + \return sign of d + \todo Move near complex sign +*/ +nr_double_t inline sign (const nr_double_t d) { + return d < 0 ? -1 : 1; +} + +/*!\brief Real cardinal sinus + + Compute \f$\mathrm{sinc}\;d=\frac{\sin d}{d}\f$ + \param[in] d real number + \return cardianal sinus of s + \todo Why not inline +*/ +nr_double_t inline sinc (const nr_double_t d) { + if (d == 0) return 1; + return sin (d) / d; +} + +/*!\brief Fix function + + Fix is nearest integral value in direction of 0, + \f[ + \operatorname{fix} d=\begin{cases} + \operatorname{floor} d & \text{if } d > 0 \\ + \operatorname{ceil} d & \text{else} + \end{cases} + \f] + + \param[in] d real number + \return fixed complex number + \todo Why not inline? +*/ +nr_double_t inline fix (const nr_double_t d) { + return (d > 0) ? floor (d) : ceil (d); +} + +/*!\brief Heaviside step function + + The Heaviside step function, H, also called unit step function, + is a discontinuous function whose value is zero for negative argument and + one for positive argument. For zero by convention, H(0)=0.5 + \param[in] d Heaviside argument + \return Heaviside step + \todo Create Heaviside alias +*/ +nr_double_t inline step (const nr_double_t d) { + nr_double_t x = d; + if (x < 0.0) + x = 0.0; + else if (x > 0.0) + x = 1.0; + else + x = 0.5; + return x; +} + +/*!\brief Compute factorial n ie \$n!\$ + +*/ +unsigned int inline +factorial (unsigned int n) { + unsigned int result = 1; + + /* 13! > 2^32 */ + assert (n < 13); + + if (n == 0) + return 1; + + for (; n > 1; n--) + result = result * n; + + return result; +} // // overload complex manipulations on reals // -nr_double_t real (const nr_double_t); -nr_double_t imag (const nr_double_t); -nr_double_t norm (const nr_double_t); -nr_double_t conj (const nr_double_t); -nr_double_t abs (const nr_double_t); + + +/*!\brief Real part of real number + + \param[in] r Real number + \return Real part of r ie r +*/ +nr_double_t inline real (const nr_double_t r) { + return r; +} + +/*!\brief Imaginary part of complex number + + \param[in] r Real number + \return Imaginary part of r +*/ +nr_double_t inline imag ([[maybe_unused]] const nr_double_t r) { + return 0.0; +} + +/*!\brief Compute euclidean norm of real number + + Compute \f$r^2\f$ + \param[in] r Real number + \return Euclidean norm of r +*/ +nr_double_t inline norm (const nr_double_t r) { + return r * r; +} + +/*!\brief Compute complex modulus of real number + + \param[in] r Real number + \return Modulus of r +*/ +nr_double_t inline abs (const nr_double_t r) { + return std::abs (r); +} + +/*!\brief Conjugate of real number + + \param[in] r Real number + \return Conjugate of real r ie r +*/ +nr_double_t inline conj (const nr_double_t r) { + return r; +} + +/*! + * \brief rad2deg Convert radian to degree + * \param x input + * \return input in degree (x)*180/pi + */ +nr_double_t inline rad2deg (const nr_double_t x) { + return (180.0 * (x) / pi); +} + +/*! + * \brief deg2rad Convert radian to degree + * \param x input + * \return input in radian (x)*pi/180 + */ +nr_double_t inline deg2rad (const nr_double_t x) { + return (pi * (x) / 180.0); +} + +/*!\brief Compute the third power of input */ +static inline nr_double_t cubic (const nr_double_t x) { return (x * x * x); } + +/*!\brief Convert Celsius to Kelvin */ +static inline nr_double_t celsius2kelvin (const nr_double_t x) { return (x - K); } + +/*!\brief Convert Kelvin to Celsius */ +static inline nr_double_t kelvin2celsius (const nr_double_t x) { return (x + K); } + } // namespace qucs From 0aaef11aac0cc2101b2f2299f11c425b01756182 Mon Sep 17 00:00:00 2001 From: Vadim Kuznetsov Date: Tue, 3 Sep 2024 13:56:12 +0300 Subject: [PATCH 4/4] Make complex.cpp inline and move to complex.h --- src/math/complex.cpp | 754 ----------------------------------------- src/math/complex.h | 791 ++++++++++++++++++++++++++++++++++++++----- 2 files changed, 708 insertions(+), 837 deletions(-) diff --git a/src/math/complex.cpp b/src/math/complex.cpp index 4088e3a4..7f76e65a 100644 --- a/src/math/complex.cpp +++ b/src/math/complex.cpp @@ -46,621 +46,6 @@ namespace qucs { -// -// trigonometric complex -// - -/*! \brief Compute complex cosine - \param[in] z complex angle - \return cosine of z -*/ -nr_complex_t cos (const nr_complex_t z) { - return std::cos (z); -} - - -/*! \brief Compute complex sine - \param[in] z complex angle - \return sine of z -*/ -nr_complex_t sin (const nr_complex_t z) { - return std::sin (z); -} - - -/*! \brief Compute complex tangent - \param[in] z complex angle - \return tangent of z -*/ -nr_complex_t tan (const nr_complex_t z) { - return std::tan (z); -} - - -/*! \brief Compute complex arc cosine - \param[in] z complex arc - \return arc cosine of z -*/ -nr_complex_t acos (const nr_complex_t z) { - return std::acos (z); -} - - -/*! \brief Compute complex arc sine - \param[in] z complex arc - \return arc sine of z -*/ -nr_complex_t asin (const nr_complex_t z) { - return std::asin (z); -} - -/*! \brief Compute complex arc tangent - \param[in] z complex arc - \return arc tangent of z -*/ -nr_complex_t atan (const nr_complex_t z) { - return std::atan (z); -} - - -// -// hyperbolic complex -// - -/*! \brief Compute complex hyperbolic cosine - \param[in] z complex arc - \return hyperbolic cosine of z -*/ -nr_complex_t cosh (const nr_complex_t z) { - return std::cosh (z); -} - - -/*! \brief Compute complex hyperbolic sine - \param[in] z complex arc - \return hyperbolic sine of z -*/ -nr_complex_t sinh (const nr_complex_t z) { - return std::sinh (z); -} - - -/*! \brief Compute complex hyperbolic tangent - \param[in] z complex arc - \return hyperbolic tangent of z -*/ -nr_complex_t tanh (const nr_complex_t z) { - return std::tanh (z); -} - - -/*! \brief Compute complex arc hyperbolic cosine - \param[in] z complex arc - \return arc hyperbolic cosine of z -*/ -nr_complex_t acosh (const nr_complex_t z) { - return std::acosh (z); -} - - -/*! \brief Compute complex arc hyperbolic sine - \param[in] z complex arc - \return arc hyperbolic sine of z -*/ -nr_complex_t asinh (const nr_complex_t z) { - return std::asinh (z); -} - - -/*! \brief Compute complex arc hyperbolic tangent - \param[in] z complex arc - \return arc hyperbolic tangent of z -*/ -nr_complex_t atanh (const nr_complex_t z) { - return std::atanh (z); -} - - -// -// transcendentals overloads -// - -/*! \brief Compute complex exponential - \param[in] z complex number - \return exponential of z -*/ -nr_complex_t exp (const nr_complex_t z) -{ - nr_double_t mag = exp (real (z)); - return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z))); -} - -/*! \brief Compute principal value of natural logarithm of z - \param[in] z complex number - \return principal value of natural logarithm of z -*/ -nr_complex_t log (const nr_complex_t z) -{ - nr_double_t phi = arg (z); - return nr_complex_t (log (abs (z)), phi); -} - -/*! \brief Compute principal value of decimal logarithm of z - \param[in] z complex number - \return principal value of decimal logarithm of z -*/ -nr_complex_t log10 (const nr_complex_t z) -{ - nr_double_t phi = arg (z); - return nr_complex_t (log10 (abs (z)), phi * log10e); -} - - -/*!\brief Compute power function with real exponent - - \param[in] z complex mantisse - \param[in] d real exponent - \return z power d (\f$z^d\f$) -*/ -nr_complex_t pow (const nr_complex_t z, const nr_double_t d) { - return std::pow (z, d); -} - -/*!\brief Compute power function with complex exponent but real mantisse - - \param[in] d real mantisse - \param[in] z complex exponent - \return d power z (\f$d^z\f$) -*/ -nr_complex_t pow (const nr_double_t d, const nr_complex_t z) { - return std::pow (d, z); -} - -/*!\brief Compute complex power function - - \param[in] z1 complex mantisse - \param[in] z2 complex exponent - \return d power z (\f$z_1^{z_2}\f$) -*/ -nr_complex_t pow (const nr_complex_t z1, const nr_complex_t z2) { - return std::pow (z1, z2); -} - - -/*!\brief Compute principal value of square root - - Compute the square root of a given complex number (except negative - real), and with a branch cut along the negative real axis. - - \param[in] z complex number - \return principal value of square root z -*/ -nr_complex_t sqrt (const nr_complex_t z) -{ - return std::sqrt (z); -} - - -/*!\brief Compute euclidean norm of complex number - - Compute \f$(\Re\mathrm{e}\;z )^2+ (\Im\mathrm{m}\;z)^2=|z|^2\f$ - \param[in] z Complex number - \return Euclidean norm of z -*/ -nr_double_t norm (const nr_complex_t z) -{ - return std::norm (z); -} - - -// -// Qucs extra math functions -// - -/*!\brief Compute complex cotangent - - \param[in] z complex angle - \return cotangent of z -*/ -nr_complex_t cot (const nr_complex_t z) -{ - nr_double_t r = 2.0 * std::real (z); - nr_double_t i = 2.0 * std::imag (z); - return nr_complex_t (0.0, 1.0) + nr_complex_t (0.0, 2.0) / (std::polar (std::exp (-i), r) - 1.0); -} - -/*!\brief Compute complex arc cotangent - - \param[in] z complex arc - \return arc cotangent of z -*/ -nr_complex_t acot (const nr_complex_t z) -{ - return nr_complex_t (0.0, -0.5) * std::log (nr_complex_t (0, 2) / (z - nr_complex_t (0, 1)) + 1.0); -} - -/*!\brief Compute complex hyperbolic cotangent - - \param[in] z complex angle - \return hyperbolic cotangent of z -*/ -nr_complex_t coth (const nr_complex_t z) -{ - nr_double_t r = 2.0 * std::real (z); - nr_double_t i = 2.0 * std::imag (z); - return 1.0 + 2.0 / (std::polar (std::exp (r), i) - 1.0); -} - -/*!\brief Compute complex argument hyperbolic cotangent - - \param[in] z complex arc - \return argument hyperbolic cotangent of z -*/ -nr_complex_t acoth (const nr_complex_t z) -{ - return 0.5 * std::log (2.0 / (z - 1.0) + 1.0); -} - - -/*!\brief Compute complex hyperbolic secant - - \param[in] z complex angle - \return hyperbolic secant of z -*/ -nr_complex_t sech (const nr_complex_t z) -{ - return (1.0 / std::cosh (z)); -} - -/*!\brief Compute complex argument hyperbolic secant - - \param[in] z complex arc - \return argument hyperbolic secant of z - \todo for symmetry reason implement sech -*/ -nr_complex_t asech (const nr_complex_t z) -{ - return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z); -} - -/*!\brief Compute complex argument hyperbolic cosec - - \param[in] z complex arc - \return argument hyperbolic cosec of z -*/ -nr_complex_t cosech (const nr_complex_t z) -{ - return (1.0 / std::sinh(z)); -} - -/*!\brief Compute complex arc tangent fortran like function - - atan2 is a two-argument function that computes the arc tangent of y / x - given y and x, but with a range of \f$(-\pi;\pi]\f$ - - \param[in] z complex angle - \return arc tangent of z -*/ -nr_complex_t atan2 (const nr_complex_t y, const nr_complex_t x) -{ - nr_complex_t a = qucs::atan (y / x); // std::atan missing on OSX 10.6 - return real (x) > 0.0 ? a : -a; -} - - -// -// extra math -// - -/*!\brief Compute principal value of binary logarithm of z - - \param[in] z complex number - \return principal value of binary logarithm of z -*/ -nr_complex_t log2 (const nr_complex_t z) -{ - nr_double_t phi = std::arg (z); - return nr_complex_t (std::log (std::abs (z)) * log2e, phi * log2e); -} - -/*!\brief complex signum function - - compute \f[ - \mathrm{signum}\;z= \mathrm{signum} (re^{i\theta}) - = \begin{cases} - 0 & \text{if } z=0 \\ - e^{i\theta} & \text{else} - \end{cases} - \f] - \param[in] z complex number - \return signum of z - \todo Better implementation z/abs(z) is not really stable -*/ -nr_complex_t signum (const nr_complex_t z) -{ - if (z == 0.0) return 0; - return z / abs (z); -} - -/*!\brief complex sign function - - compute \f[ - \mathrm{sign}\;z= \mathrm{sign} (re^{i\theta}) - = \begin{cases} - 1 & \text{if } z=0 \\ - e^{i\theta} & \text{else} - \end{cases} - \f] - \param[in] z complex number - \return sign of z - \todo Better implementation z/abs(z) is not really stable -*/ -nr_complex_t sign (const nr_complex_t z) -{ - if (z == 0.0) return nr_complex_t (1); - return z / abs (z); -} - - -/*!\brief Cardinal sine - - Compute \f$\mathrm{sinc}\;z=\frac{\sin z}{z}\f$ - \param[in] z complex number - \return cardianal sine of z -*/ -nr_complex_t sinc (const nr_complex_t z) -{ - if (real(z) == 0.0 && imag(z)) return 1; - return std::sin (z) / z; -} - -/*!\brief Euclidean distance function for complex argument - - The xhypot() function returns \f$\sqrt{a^2+b^2}\f$. - This is the length of the hypotenuse of a right-angle triangle with sides - of length a and b, or the distance - of the point (a,b) from the origin. - - \param[in] a first length - \param[in] b second length - \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$ -*/ -nr_double_t xhypot (const nr_complex_t a, const nr_complex_t b) -{ - nr_double_t c = norm (a); - nr_double_t d = norm (b); - if (c > d) - return abs (a) * std::sqrt (1.0 + d / c); - else if (d == 0.0) - return 0.0; - else - return abs (b) * std::sqrt (1.0 + c / d); -} - -/*!\brief Euclidean distance function for a double b complex */ -nr_double_t xhypot (nr_double_t a, nr_complex_t b) -{ - return xhypot (nr_complex_t (a), b); -} - -/*!\brief Euclidean distance function for b double a complex */ -nr_double_t xhypot (nr_complex_t a, nr_double_t b) -{ - return xhypot (a, nr_complex_t (b)); -} - - -/*!\brief Complex round - Round is the nearest integral value - Apply round to real and imaginary part - \param[in] z complex number - \return rounded complex number -*/ -nr_complex_t round (const nr_complex_t z) -{ - nr_double_t zreal = real (z); - nr_double_t zimag = imag (z); - // qucs::round resolved for double - zreal = qucs::round (zreal); - zimag = qucs::round (zimag); - return nr_complex_t (zreal, zimag); -} - - -/*!\brief Complex trunc - Apply round to integer, towards zero to real and imaginary part - \param[in] z complex number - \return rounded complex number -*/ -nr_complex_t trunc (const nr_complex_t z) -{ - nr_double_t zreal = real (z); - nr_double_t zimag = imag (z); - // qucs::round resolved for double - zreal = qucs::trunc (zreal); - zimag = qucs::trunc (zimag); - return nr_complex_t (zreal, zimag); -} - - -/*!\brief Magnitude in dB - Compute \f$10\log_{10} |z|^2=20\log_{10} |z|\f$ - \param[in] z complex number - \return Magnitude in dB -*/ -nr_double_t dB (const nr_complex_t z) -{ - return 10.0 * std::log10 (std::norm (z)); -} - - -/*!\brief Compute limited complex exponential - \param[in] z complex number - \return limited exponential of z - \todo Change limexp(real) limexp(complex) file order -*/ -nr_complex_t limexp (const nr_complex_t z) -{ - nr_double_t mag = qucs::limexp (real (z)); - return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z))); -} - - -/*!\brief Construct a complex number using polar notation - \param[in] mag Magnitude - \param[in] ang Angle - \return complex number in rectangular form -*/ -nr_complex_t polar (const nr_double_t mag, const nr_double_t ang ) -{ - return std::polar (mag, ang); -} - -/*!\brief Extension of polar construction to complex - \param[in] a Magnitude - \param[in] p Angle - \return complex number in rectangular form - \bug Do not seems holomorph form of real polar -*/ -nr_complex_t polar (const nr_complex_t a, const nr_complex_t p) -{ - return a * exp(nr_complex_t(0.0, 1.0) * p); -} - - -/*!\brief Converts impedance to reflexion coefficient - \param[in] z impedance - \param[in] zref normalisation impedance - \return reflexion coefficient -*/ -nr_complex_t ztor (const nr_complex_t z, nr_complex_t zref) { - return (z - zref) / (z + zref); -} - -/*!\brief Converts reflexion coefficient to impedance - \param[in] r reflexion coefficient - \param[in] zref normalisation impedance - \return impedance -*/ -nr_complex_t rtoz (const nr_complex_t r, nr_complex_t zref) { - return zref * (1.0 + r) / (1.0 - r); -} - -/*!\brief Converts admittance to reflexion coefficient - \param[in] y admitance - \param[in] zref normalisation impedance - \return reflexion coefficient -*/ -nr_complex_t ytor (const nr_complex_t y, nr_complex_t zref) { - return (1.0 - y * zref) / (1.0 + y * zref); -} - -/*!\brief Converts reflexion coefficient to admittance - \param[in] r reflexion coefficient - \param[in] zref normalisation impedance - \return admittance -*/ -nr_complex_t rtoy (const nr_complex_t r, nr_complex_t zref) { - return (1.0 - r) / (1.0 + r) / zref; -} - - - - - -/*!\brief Complex floor - - floor is the largest integral value not greater than argument - Apply floor to real and imaginary part - \param[in] z complex number - \return floored complex number -*/ -nr_complex_t floor (const nr_complex_t z) { - return nr_complex_t (std::floor (real (z)), std::floor (imag (z))); -} - - -/*!\brief Complex ceil - Ceil is the smallest integral value not less than argument - Apply ceil to real and imaginary part - \param[in] z complex number - \return ceilled complex number -*/ -nr_complex_t ceil (const nr_complex_t z) { - return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z))); -} - -/*!\brief Complex fix - - Apply fix to real and imaginary part - \param[in] z complex number - \return fixed complex number - \todo why not using real fix -*/ -nr_complex_t fix (const nr_complex_t z) { - nr_double_t x = real (z); - nr_double_t y = imag (z); - x = (x > 0) ? std::floor (x) : std::ceil (x); - y = (y > 0) ? std::floor (y) : std::ceil (y); - return nr_complex_t (x, y); -} - - - -/*!\brief Complex fmod - Apply fmod to the complex z - \param[in] x complex number (numerator) - \param[in] y complex number (denominator) - \return return \f$x - n * y\f$ where n is the quotient of \f$x / y\f$, - rounded towards zero to an integer. -*/ -nr_complex_t fmod (const nr_complex_t x, const nr_complex_t y) { - nr_complex_t n = qucs::floor (x / y); - return x - n * y; -} - - -/*!\brief Square of complex number - - \param[in] z complex number - \return squared complex number -*/ -nr_complex_t sqr (const nr_complex_t z) { - nr_double_t r = real (z); - nr_double_t i = imag (z); - return nr_complex_t (r * r - i * i, 2 * r * i); -} - - - - - -/*!\brief Heaviside step function for complex number - - Apply Heaviside to real and imaginary part - \param[in] z Heaviside argument - \return Heaviside step - \todo Create Heaviside alias - \todo Why not using real heaviside -*/ -nr_complex_t step (const nr_complex_t z) -{ - nr_double_t x = real (z); - nr_double_t y = imag (z); - if (x < 0.0) - x = 0.0; - else if (x > 0.0) - x = 1.0; - else - x = 0.5; - if (y < 0.0) - y = 0.0; - else if (y > 0.0) - y = 1.0; - else - y = 0.5; - return nr_complex_t (x, y); -} //using namespace fspecial; @@ -711,43 +96,6 @@ nr_complex_t i0 (const nr_complex_t z) return nr_complex_t (fspecial::i0 (std::real (z)), 0); } - -/*!\brief Error function - - \param[in] z argument - \return Error function - \bug Not implemented -*/ -nr_complex_t erf (const nr_complex_t z) -{ -#ifdef HAVE_STD_ERF - nr_double_t zerf = std::erf (std::real (z)); // c++11 -#elif HAVE_ERF - nr_double_t zerf = ::erf (std::real (z)); -#else - nr_double_t zerf = fspecial::erf (std::real (z)); -#endif - return nr_complex_t (zerf, 0); -} - -/*!\brief Complementart error function - - \param[in] z argument - \return Complementary error function - \bug Not implemented -*/ -nr_complex_t erfc (const nr_complex_t z) -{ -#ifdef HAVE_STD_ERF - nr_double_t zerfc = std::erfc (std::real (z)); // c++11 -#elif HAVE_ERFC - nr_double_t zerfc = ::erfc (std::real (z)); -#else - nr_double_t zerfc = fspecial::erfc (std::real (z)); -#endif - return nr_complex_t (zerfc, 0); -} - /*!\brief Inverse of error function \param[in] z argument @@ -770,107 +118,5 @@ nr_complex_t erfcinv (const nr_complex_t z) return nr_complex_t (fspecial::erfcinv (std::real (z)), 0); } - - -// ======================== - - -/*!\brief Modulo - \todo Why not inline -*/ -nr_complex_t operator%(const nr_complex_t z1, const nr_complex_t z2) -{ - return z1 - z2 * floor (z1 / z2); -} - -/*!\brief Modulo - \todo Why not inline -*/ -nr_complex_t operator%(const nr_complex_t z1, const nr_double_t r2) -{ - return z1 - r2 * floor (z1 / r2); -} - -/*!\brief Modulo - \todo Why not inline -*/ -nr_complex_t operator%(const nr_double_t r1, const nr_complex_t z2) -{ - return r1 - z2 * floor (r1 / z2); -} - - -/*!\brief Equality of two complex - \todo Why not inline - \note Like equality of double this test - is meaningless in finite precision - Use instead fabs(x-x0) < tol -*/ -bool operator==(const nr_complex_t z1, const nr_complex_t z2) -{ - return (std::real (z1) == std::real (z2)) && (std::imag (z1) == std::imag (z2)); -} - -/*!\brief Inequality of two complex - \todo Why not inline - \note Like inequality of double this test - is meaningless in finite precision - Use instead fabs(x-x0) > tol -*/ -bool operator!=(const nr_complex_t z1, const nr_complex_t z2) -{ - return (std::real (z1) != std::real (z2)) || (std::imag (z1) != std::imag (z2)); -} - -/*!\brief Superior of equal - \todo Why not inline -*/ -bool operator>=(const nr_complex_t z1, const nr_complex_t z2) -{ - return norm (z1) >= norm (z2); -} - -/*!\brief Inferior of equal - \todo Why not inline -*/ -bool operator<=(const nr_complex_t z1, const nr_complex_t z2) -{ - return norm (z1) <= norm (z2); -} - -/*!\brief Superior - \todo Why not inline -*/ -bool operator>(const nr_complex_t z1, const nr_complex_t z2) -{ - return norm (z1) > norm (z2); -} - -/*!\brief Inferior - \todo Why not inline -*/ -bool operator<(const nr_complex_t z1, const nr_complex_t z2) -{ - return norm (z1) < norm (z2); -} - -/*! - * \brief rad2deg Convert radian to degree - * \param x input - * \return real(x)*180/pi - */ -nr_double_t rad2deg (const nr_complex_t x) { - return rad2deg (real(x)); -} - -/*! - * \brief rad2deg Convert radian to degree - * \param x input - * \return real(x)*pi/180 - */ -nr_double_t deg2rad (const nr_complex_t x) { - return deg2rad (real(x)); -} - } // namespace qucs diff --git a/src/math/complex.h b/src/math/complex.h index b7d0b68d..77b02ac0 100644 --- a/src/math/complex.h +++ b/src/math/complex.h @@ -26,8 +26,11 @@ #define __COMPLEX_H__ #include + +#include #include "real.h" + typedef std::complex nr_complex_t; // undefine this macro if it is defined already @@ -42,91 +45,614 @@ namespace qucs { // // trigonometric complex // -nr_complex_t cos (const nr_complex_t); -nr_complex_t sin (const nr_complex_t); -nr_complex_t tan (const nr_complex_t); -nr_complex_t acos (const nr_complex_t); //c++11 -nr_complex_t asin (const nr_complex_t); //c++11 -nr_complex_t atan (const nr_complex_t); //c++11 + +/*! \brief Compute complex cosine + \param[in] z complex angle + \return cosine of z +*/ +nr_complex_t inline cos (const nr_complex_t z) { + return std::cos (z); +} + + +/*! \brief Compute complex sine + \param[in] z complex angle + \return sine of z +*/ +nr_complex_t inline sin (const nr_complex_t z) { + return std::sin (z); +} + + +/*! \brief Compute complex tangent + \param[in] z complex angle + \return tangent of z +*/ +nr_complex_t inline tan (const nr_complex_t z) { + return std::tan (z); +} + + +/*! \brief Compute complex arc cosine + \param[in] z complex arc + \return arc cosine of z +*/ +nr_complex_t inline acos (const nr_complex_t z) { + return std::acos (z); +} + + +/*! \brief Compute complex arc sine + \param[in] z complex arc + \return arc sine of z +*/ +nr_complex_t inline asin (const nr_complex_t z) { + return std::asin (z); +} + +/*! \brief Compute complex arc tangent + \param[in] z complex arc + \return arc tangent of z +*/ +nr_complex_t inline atan (const nr_complex_t z) { + return std::atan (z); +} // // hyperbolic complex // -nr_complex_t cosh (const nr_complex_t); -nr_complex_t sinh (const nr_complex_t); -nr_complex_t tanh (const nr_complex_t); -nr_complex_t acosh (const nr_complex_t); //c++11 -nr_complex_t asinh (const nr_complex_t); //c++11 -nr_complex_t atanh (const nr_complex_t); //c++11 + +/*! \brief Compute complex hyperbolic cosine + \param[in] z complex arc + \return hyperbolic cosine of z +*/ +nr_complex_t inline cosh (const nr_complex_t z) { + return std::cosh (z); +} + + +/*! \brief Compute complex hyperbolic sine + \param[in] z complex arc + \return hyperbolic sine of z +*/ +nr_complex_t inline sinh (const nr_complex_t z) { + return std::sinh (z); +} + + +/*! \brief Compute complex hyperbolic tangent + \param[in] z complex arc + \return hyperbolic tangent of z +*/ +nr_complex_t inline tanh (const nr_complex_t z) { + return std::tanh (z); +} + + +/*! \brief Compute complex arc hyperbolic cosine + \param[in] z complex arc + \return arc hyperbolic cosine of z +*/ +nr_complex_t inline acosh (const nr_complex_t z) { + return std::acosh (z); +} + + +/*! \brief Compute complex arc hyperbolic sine + \param[in] z complex arc + \return arc hyperbolic sine of z +*/ +nr_complex_t inline asinh (const nr_complex_t z) { + return std::asinh (z); +} + + +/*! \brief Compute complex arc hyperbolic tangent + \param[in] z complex arc + \return arc hyperbolic tangent of z +*/ +nr_complex_t inline atanh (const nr_complex_t z) { + return std::atanh (z); +} // // transcendentals overloads // -nr_complex_t exp (const nr_complex_t); -nr_complex_t log (const nr_complex_t); -nr_complex_t log10 (const nr_complex_t); -nr_complex_t pow (const nr_complex_t, const nr_double_t); -nr_complex_t pow (const nr_double_t, const nr_complex_t); -nr_complex_t pow (const nr_complex_t, const nr_complex_t); -nr_complex_t sqrt (const nr_complex_t); -nr_double_t norm (const nr_complex_t); +/*! \brief Compute complex exponential + \param[in] z complex number + \return exponential of z +*/ +nr_complex_t inline exp (const nr_complex_t z) +{ + nr_double_t mag = exp (real (z)); + return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z))); +} + +/*! \brief Compute principal value of natural logarithm of z + \param[in] z complex number + \return principal value of natural logarithm of z +*/ +nr_complex_t inline log (const nr_complex_t z) +{ + nr_double_t phi = arg (z); + return nr_complex_t (log (abs (z)), phi); +} + +/*! \brief Compute principal value of decimal logarithm of z + \param[in] z complex number + \return principal value of decimal logarithm of z +*/ +nr_complex_t inline log10 (const nr_complex_t z) +{ + nr_double_t phi = arg (z); + return nr_complex_t (log10 (abs (z)), phi * log10e); +} + + +/*!\brief Compute power function with real exponent + + \param[in] z complex mantisse + \param[in] d real exponent + \return z power d (\f$z^d\f$) +*/ +nr_complex_t inline pow (const nr_complex_t z, const nr_double_t d) { + return std::pow (z, d); +} + +/*!\brief Compute power function with complex exponent but real mantisse + + \param[in] d real mantisse + \param[in] z complex exponent + \return d power z (\f$d^z\f$) +*/ +nr_complex_t inline pow (const nr_double_t d, const nr_complex_t z) { + return std::pow (d, z); +} + +/*!\brief Compute complex power function + + \param[in] z1 complex mantisse + \param[in] z2 complex exponent + \return d power z (\f$z_1^{z_2}\f$) +*/ +nr_complex_t inline pow (const nr_complex_t z1, const nr_complex_t z2) { + return std::pow (z1, z2); +} + + +/*!\brief Compute principal value of square root + + Compute the square root of a given complex number (except negative + real), and with a branch cut along the negative real axis. + + \param[in] z complex number + \return principal value of square root z +*/ +nr_complex_t inline sqrt (const nr_complex_t z) +{ + return std::sqrt (z); +} + + +/*!\brief Compute euclidean norm of complex number + + Compute \f$(\Re\mathrm{e}\;z )^2+ (\Im\mathrm{m}\;z)^2=|z|^2\f$ + \param[in] z Complex number + \return Euclidean norm of z +*/ +nr_double_t inline norm (const nr_complex_t z) +{ + return std::norm (z); +} -// -// Qucs extra trancendental functions -// -nr_complex_t cot (const nr_complex_t); -nr_complex_t acot (const nr_complex_t); -nr_complex_t coth (const nr_complex_t); -nr_complex_t acoth (const nr_complex_t); -nr_complex_t sech (const nr_complex_t); -nr_complex_t asech (const nr_complex_t); -nr_complex_t cosech (const nr_complex_t); -nr_complex_t atan2 (const nr_complex_t, const nr_complex_t); // -// extra math +// Qucs extra trancendental functions // -nr_complex_t log2 (const nr_complex_t); -nr_complex_t signum (const nr_complex_t); -nr_complex_t sign (const nr_complex_t); -nr_complex_t sinc (const nr_complex_t); -nr_double_t xhypot (const nr_complex_t, const nr_complex_t); -nr_double_t xhypot (const nr_double_t, const nr_complex_t); -nr_double_t xhypot (const nr_complex_t, const nr_double_t); - -nr_complex_t round (const nr_complex_t); -nr_complex_t trunc (const nr_complex_t); - - -nr_double_t dB (const nr_complex_t); +/*!\brief Compute complex cotangent + + \param[in] z complex angle + \return cotangent of z +*/ +nr_complex_t inline cot (const nr_complex_t z) +{ + nr_double_t r = 2.0 * std::real (z); + nr_double_t i = 2.0 * std::imag (z); + return nr_complex_t (0.0, 1.0) + nr_complex_t (0.0, 2.0) / (std::polar (std::exp (-i), r) - 1.0); +} + +/*!\brief Compute complex arc cotangent + + \param[in] z complex arc + \return arc cotangent of z +*/ +nr_complex_t inline acot (const nr_complex_t z) +{ + return nr_complex_t (0.0, -0.5) * std::log (nr_complex_t (0, 2) / (z - nr_complex_t (0, 1)) + 1.0); +} + +/*!\brief Compute complex hyperbolic cotangent + + \param[in] z complex angle + \return hyperbolic cotangent of z +*/ +nr_complex_t inline coth (const nr_complex_t z) +{ + nr_double_t r = 2.0 * std::real (z); + nr_double_t i = 2.0 * std::imag (z); + return 1.0 + 2.0 / (std::polar (std::exp (r), i) - 1.0); +} + +/*!\brief Compute complex argument hyperbolic cotangent + + \param[in] z complex arc + \return argument hyperbolic cotangent of z +*/ +nr_complex_t inline acoth (const nr_complex_t z) +{ + return 0.5 * std::log (2.0 / (z - 1.0) + 1.0); +} + + +/*!\brief Compute complex hyperbolic secant + + \param[in] z complex angle + \return hyperbolic secant of z +*/ +nr_complex_t inline sech (const nr_complex_t z) +{ + return (1.0 / std::cosh (z)); +} + +/*!\brief Compute complex argument hyperbolic secant + + \param[in] z complex arc + \return argument hyperbolic secant of z + \todo for symmetry reason implement sech +*/ +nr_complex_t inline asech (const nr_complex_t z) +{ + return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z); +} + +/*!\brief Compute complex argument hyperbolic cosec + + \param[in] z complex arc + \return argument hyperbolic cosec of z +*/ +nr_complex_t inline cosech (const nr_complex_t z) +{ + return (1.0 / std::sinh(z)); +} + +/*!\brief Compute complex arc tangent fortran like function + + atan2 is a two-argument function that computes the arc tangent of y / x + given y and x, but with a range of \f$(-\pi;\pi]\f$ + + \param[in] z complex angle + \return arc tangent of z +*/ +nr_complex_t inline atan2 (const nr_complex_t y, const nr_complex_t x) +{ + nr_complex_t a = qucs::atan (y / x); // std::atan missing on OSX 10.6 + return real (x) > 0.0 ? a : -a; +} -nr_complex_t limexp (const nr_complex_t); -nr_complex_t polar (const nr_double_t mag, const nr_double_t theta = 0.0); -/// \bug are these needed/used? -//nr_complex_t polar (const nr_double_t a, const nr_complex_t p); -//nr_complex_t polar (const nr_complex_t a, const nr_double_t p = 0.0); -nr_complex_t polar (const nr_complex_t a, const nr_complex_t p = 0.0); -nr_complex_t ztor (const nr_complex_t, const nr_complex_t zref = 50.0); -nr_complex_t rtoz (const nr_complex_t, const nr_complex_t zref = 50.0); -nr_complex_t ytor (const nr_complex_t, const nr_complex_t zref = 50.0); -nr_complex_t rtoy (const nr_complex_t, const nr_complex_t zref = 50.0); - - -nr_complex_t floor (const nr_complex_t ); -nr_complex_t ceil (const nr_complex_t ); -nr_complex_t fix (const nr_complex_t ); - -/// \todo add fmod to manual -nr_complex_t fmod (const nr_complex_t x, const nr_complex_t y); -nr_complex_t sqr (const nr_complex_t z); +// +// extra math +// +/*!\brief Compute principal value of binary logarithm of z + + \param[in] z complex number + \return principal value of binary logarithm of z +*/ +nr_complex_t inline log2 (const nr_complex_t z) +{ + nr_double_t phi = std::arg (z); + return nr_complex_t (std::log (std::abs (z)) * log2e, phi * log2e); +} + +/*!\brief complex signum function + + compute \f[ + \mathrm{signum}\;z= \mathrm{signum} (re^{i\theta}) + = \begin{cases} + 0 & \text{if } z=0 \\ + e^{i\theta} & \text{else} + \end{cases} + \f] + \param[in] z complex number + \return signum of z + \todo Better implementation z/abs(z) is not really stable +*/ +nr_complex_t inline signum (const nr_complex_t z) +{ + if (z == 0.0) return 0; + return z / abs (z); +} + +/*!\brief complex sign function + + compute \f[ + \mathrm{sign}\;z= \mathrm{sign} (re^{i\theta}) + = \begin{cases} + 1 & \text{if } z=0 \\ + e^{i\theta} & \text{else} + \end{cases} + \f] + \param[in] z complex number + \return sign of z + \todo Better implementation z/abs(z) is not really stable +*/ +nr_complex_t inline sign (const nr_complex_t z) +{ + if (z == 0.0) return nr_complex_t (1); + return z / abs (z); +} + + +/*!\brief Cardinal sine + + Compute \f$\mathrm{sinc}\;z=\frac{\sin z}{z}\f$ + \param[in] z complex number + \return cardianal sine of z +*/ +nr_complex_t inline sinc (const nr_complex_t z) +{ + if (real(z) == 0.0 && imag(z)) return 1; + return std::sin (z) / z; +} + +/*!\brief Euclidean distance function for complex argument + + The xhypot() function returns \f$\sqrt{a^2+b^2}\f$. + This is the length of the hypotenuse of a right-angle triangle with sides + of length a and b, or the distance + of the point (a,b) from the origin. + + \param[in] a first length + \param[in] b second length + \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$ +*/ +nr_double_t inline xhypot (const nr_complex_t a, const nr_complex_t b) +{ + nr_double_t c = norm (a); + nr_double_t d = norm (b); + if (c > d) + return abs (a) * std::sqrt (1.0 + d / c); + else if (d == 0.0) + return 0.0; + else + return abs (b) * std::sqrt (1.0 + c / d); +} + +/*!\brief Euclidean distance function for a double b complex */ +nr_double_t inline xhypot (nr_double_t a, nr_complex_t b) +{ + return xhypot (nr_complex_t (a), b); +} + +/*!\brief Euclidean distance function for b double a complex */ +nr_double_t inline xhypot (nr_complex_t a, nr_double_t b) +{ + return xhypot (a, nr_complex_t (b)); +} + + +/*!\brief Complex round + Round is the nearest integral value + Apply round to real and imaginary part + \param[in] z complex number + \return rounded complex number +*/ +nr_complex_t inline round (const nr_complex_t z) +{ + nr_double_t zreal = real (z); + nr_double_t zimag = imag (z); + // qucs::round resolved for double + zreal = qucs::round (zreal); + zimag = qucs::round (zimag); + return nr_complex_t (zreal, zimag); +} + + +/*!\brief Complex trunc + Apply round to integer, towards zero to real and imaginary part + \param[in] z complex number + \return rounded complex number +*/ +nr_complex_t inline trunc (const nr_complex_t z) +{ + nr_double_t zreal = real (z); + nr_double_t zimag = imag (z); + // qucs::round resolved for double + zreal = qucs::trunc (zreal); + zimag = qucs::trunc (zimag); + return nr_complex_t (zreal, zimag); +} + + +/*!\brief Magnitude in dB + Compute \f$10\log_{10} |z|^2=20\log_{10} |z|\f$ + \param[in] z complex number + \return Magnitude in dB +*/ +nr_double_t inline dB (const nr_complex_t z) +{ + return 10.0 * std::log10 (std::norm (z)); +} + + +/*!\brief Compute limited complex exponential + \param[in] z complex number + \return limited exponential of z + \todo Change limexp(real) limexp(complex) file order +*/ +nr_complex_t inline limexp (const nr_complex_t z) +{ + nr_double_t mag = qucs::limexp (real (z)); + return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z))); +} + +/*!\brief Construct a complex number using polar notation + \param[in] mag Magnitude + \param[in] ang Angle + \return complex number in rectangular form +*/ +nr_complex_t inline polar (const nr_double_t mag, const nr_double_t ang = 0.0) +{ + return std::polar (mag, ang); +} + +/*!\brief Extension of polar construction to complex + \param[in] a Magnitude + \param[in] p Angle + \return complex number in rectangular form + \bug Do not seems holomorph form of real polar +*/ +nr_complex_t inline polar (const nr_complex_t a, const nr_complex_t p = 0.0) +{ + return a * exp(nr_complex_t(0.0, 1.0) * p); +} + + +/*!\brief Converts impedance to reflexion coefficient + \param[in] z impedance + \param[in] zref normalisation impedance + \return reflexion coefficient +*/ +nr_complex_t inline ztor (const nr_complex_t z, nr_complex_t zref = 50.0) { + return (z - zref) / (z + zref); +} + +/*!\brief Converts reflexion coefficient to impedance + \param[in] r reflexion coefficient + \param[in] zref normalisation impedance + \return impedance +*/ +nr_complex_t inline rtoz (const nr_complex_t r, nr_complex_t zref = 50.0) { + return zref * (1.0 + r) / (1.0 - r); +} + +/*!\brief Converts admittance to reflexion coefficient + \param[in] y admitance + \param[in] zref normalisation impedance + \return reflexion coefficient +*/ +nr_complex_t inline ytor (const nr_complex_t y, nr_complex_t zref = 50.0) { + return (1.0 - y * zref) / (1.0 + y * zref); +} + +/*!\brief Converts reflexion coefficient to admittance + \param[in] r reflexion coefficient + \param[in] zref normalisation impedance + \return admittance +*/ +nr_complex_t inline rtoy (const nr_complex_t r, nr_complex_t zref = 50.0) { + return (1.0 - r) / (1.0 + r) / zref; +} + +/*!\brief Complex floor + + floor is the largest integral value not greater than argument + Apply floor to real and imaginary part + \param[in] z complex number + \return floored complex number +*/ +nr_complex_t inline floor (const nr_complex_t z) { + return nr_complex_t (std::floor (real (z)), std::floor (imag (z))); +} + + +/*!\brief Complex ceil + Ceil is the smallest integral value not less than argument + Apply ceil to real and imaginary part + \param[in] z complex number + \return ceilled complex number +*/ +nr_complex_t inline ceil (const nr_complex_t z) { + return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z))); +} + +/*!\brief Complex fix + + Apply fix to real and imaginary part + \param[in] z complex number + \return fixed complex number + \todo why not using real fix +*/ +nr_complex_t inline fix (const nr_complex_t z) { + nr_double_t x = real (z); + nr_double_t y = imag (z); + x = (x > 0) ? std::floor (x) : std::ceil (x); + y = (y > 0) ? std::floor (y) : std::ceil (y); + return nr_complex_t (x, y); +} + + + +/*!\brief Complex fmod + Apply fmod to the complex z + \param[in] x complex number (numerator) + \param[in] y complex number (denominator) + \return return \f$x - n * y\f$ where n is the quotient of \f$x / y\f$, + rounded towards zero to an integer. +*/ +nr_complex_t inline fmod (const nr_complex_t x, const nr_complex_t y) { + nr_complex_t n = qucs::floor (x / y); + return x - n * y; +} + + +/*!\brief Square of complex number + + \param[in] z complex number + \return squared complex number +*/ +nr_complex_t inline sqr (const nr_complex_t z) { + nr_double_t r = real (z); + nr_double_t i = imag (z); + return nr_complex_t (r * r - i * i, 2 * r * i); +} + + + + + +/*!\brief Heaviside step function for complex number + + Apply Heaviside to real and imaginary part + \param[in] z Heaviside argument + \return Heaviside step + \todo Create Heaviside alias + \todo Why not using real heaviside +*/ +nr_complex_t inline step (const nr_complex_t z) +{ + nr_double_t x = real (z); + nr_double_t y = imag (z); + if (x < 0.0) + x = 0.0; + else if (x > 0.0) + x = 1.0; + else + x = 0.5; + if (y < 0.0) + y = 0.0; + else if (y > 0.0) + y = 1.0; + else + y = 0.5; + return nr_complex_t (x, y); +} -nr_complex_t step (const nr_complex_t); // bessel functions @@ -136,27 +662,126 @@ nr_complex_t i0 (const nr_complex_t); // error functions -nr_complex_t erf (const nr_complex_t); -nr_complex_t erfc (const nr_complex_t); nr_complex_t erfinv (const nr_complex_t); //see fspecial nr_complex_t erfcinv (const nr_complex_t); //see fspecial - -nr_double_t rad2deg (const nr_complex_t ); -nr_double_t deg2rad (const nr_complex_t ); - -// modulo operators -nr_complex_t operator % (const nr_complex_t, const nr_complex_t); -nr_complex_t operator % (const nr_complex_t, const nr_double_t); -nr_complex_t operator % (const nr_double_t, const nr_complex_t); - -// comparisons -bool operator == (const nr_complex_t, const nr_complex_t); -bool operator != (const nr_complex_t, const nr_complex_t); -bool operator >= (const nr_complex_t, const nr_complex_t); -bool operator <= (const nr_complex_t, const nr_complex_t); -bool operator > (const nr_complex_t, const nr_complex_t); -bool operator < (const nr_complex_t, const nr_complex_t); +/*!\brief Error function + + \param[in] z argument + \return Error function + \bug Not implemented +*/ +nr_complex_t inline erf (const nr_complex_t z) +{ + nr_double_t zerf = std::erf (std::real (z)); // c++11 + return nr_complex_t (zerf, 0); +} + +/*!\brief Complementart error function + + \param[in] z argument + \return Complementary error function + \bug Not implemented +*/ +nr_complex_t inline erfc (const nr_complex_t z) +{ + nr_double_t zerfc = std::erfc (std::real (z)); // c++11 + return nr_complex_t (zerfc, 0); +} + + +/*! + * \brief rad2deg Convert radian to degree + * \param x input + * \return real(x)*180/pi + */ +nr_double_t inline rad2deg (const nr_complex_t x) { + return rad2deg (real(x)); +} + +/*! + * \brief rad2deg Convert radian to degree + * \param x input + * \return real(x)*pi/180 + */ +nr_double_t inline deg2rad (const nr_complex_t x) { + return deg2rad (real(x)); +} + + + +// ======================== + + +/*!\brief Modulo +*/ +nr_complex_t inline operator%(const nr_complex_t z1, const nr_complex_t z2) +{ + return z1 - z2 * floor (z1 / z2); +} + +/*!\brief Modulo +*/ +nr_complex_t inline operator%(const nr_complex_t z1, const nr_double_t r2) +{ + return z1 - r2 * floor (z1 / r2); +} + +/*!\brief Modulo +*/ +nr_complex_t inline operator%(const nr_double_t r1, const nr_complex_t z2) +{ + return r1 - z2 * floor (r1 / z2); +} + + +/*!\brief Equality of two complex + \note Like equality of double this test + is meaningless in finite precision + Use instead fabs(x-x0) < tol +*/ +bool inline operator==(const nr_complex_t z1, const nr_complex_t z2) +{ + return (std::real (z1) == std::real (z2)) && (std::imag (z1) == std::imag (z2)); +} + +/*!\brief Inequality of two complex + \note Like inequality of double this test + is meaningless in finite precision + Use instead fabs(x-x0) > tol +*/ +bool inline operator!=(const nr_complex_t z1, const nr_complex_t z2) +{ + return (std::real (z1) != std::real (z2)) || (std::imag (z1) != std::imag (z2)); +} + +/*!\brief Superior of equal +*/ +bool inline operator>=(const nr_complex_t z1, const nr_complex_t z2) +{ + return norm (z1) >= norm (z2); +} + +/*!\brief Inferior of equal +*/ +bool inline operator<=(const nr_complex_t z1, const nr_complex_t z2) +{ + return norm (z1) <= norm (z2); +} + +/*!\brief Superior +*/ +bool inline operator>(const nr_complex_t z1, const nr_complex_t z2) +{ + return norm (z1) > norm (z2); +} + +/*!\brief Inferior +*/ +bool inline operator<(const nr_complex_t z1, const nr_complex_t z2) +{ + return norm (z1) < norm (z2); +} } // namespace qucs