Qucs-core  0.0.19
complex.cpp
Go to the documentation of this file.
00001 /*
00002  * complex.cpp - complex number class implementation
00003  *
00004  * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org>
00005  * Copyright (C) 2006 Gunther Kraut <gn.kraut@t-online.de>
00006  *
00007  * This is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2, or (at your option)
00010  * any later version.
00011  *
00012  * This software is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this package; see the file COPYING.  If not, write to
00019  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00020  * Boston, MA 02110-1301, USA.
00021  *
00022  * $Id$
00023  *
00024  */
00025 
00030 #if HAVE_CONFIG_H
00031 # include <config.h>
00032 #endif
00033 
00034 #include <cmath>
00035 #include <assert.h>
00036 #include <errno.h>
00037 #include <stdio.h>
00038 #include <stdlib.h>
00039 
00040 #include "constants.h"
00041 #include "precision.h"
00042 #include "complex.h"
00043 #include "consts.h"
00044 #include "fspecial.h"
00045 
00046 
00047 namespace qucs {
00048 
00049 //
00050 // trigonometric complex
00051 //
00052 
00057 nr_complex_t cos (const nr_complex_t z) {
00058     return std::cos (z);
00059 }
00060 
00061 
00066 nr_complex_t sin (const nr_complex_t z) {
00067     return std::sin (z);
00068 }
00069 
00070 
00075 nr_complex_t tan (const nr_complex_t z) {
00076     return std::tan (z);
00077 }
00078 
00079 
00084 nr_complex_t acos (const nr_complex_t z) {
00085 #ifdef HAVE_CXX_COMPLEX_ACOS
00086     return std::acos (z);
00087 #else
00088     // missing on OSX 10.6
00089     nr_double_t r = real (z);
00090     nr_double_t i = imag (z);
00091     nr_complex_t y = sqrt (z * z - 1.0);
00092     if (r * i < 0.0) y = -y;
00093     return nr_complex_t (0, -1.0) * log (z + y);
00094 #endif
00095 }
00096 
00097 
00102 nr_complex_t asin (const nr_complex_t z) {
00103 #ifdef HAVE_CXX_COMPLEX_ASIN
00104     return std::asin (z);
00105 #else
00106     // missing on OSX 10.6
00107     nr_double_t r = real (z);
00108     nr_double_t i = imag (z);
00109     return nr_complex_t (0.0, -1.0) * log (nr_complex_t (-i, r) + sqrt (1.0 - z * z));
00110 #endif
00111 }
00112 
00117 nr_complex_t atan (const nr_complex_t z) {
00118 #ifdef HAVE_CXX_COMPLEX_ATAN
00119     return std::atan (z);
00120 #else
00121     // missing on OSX 10.6
00122     return nr_complex_t (0.0, -0.5) * log (nr_complex_t (0, 2) / (z + nr_complex_t (0, 1)) - 1.0);
00123 #endif
00124 }
00125 
00126 
00127 //
00128 // hyperbolic complex
00129 //
00130 
00135 nr_complex_t cosh (const nr_complex_t z) {
00136     return std::cosh (z);
00137 }
00138 
00139 
00144 nr_complex_t sinh (const nr_complex_t z) {
00145     return std::sinh (z);
00146 }
00147 
00148 
00153 nr_complex_t tanh (const nr_complex_t z) {
00154     return std::tanh (z);
00155 }
00156 
00157 
00162 nr_complex_t acosh (const nr_complex_t z) {
00163 #ifdef HAVE_CXX_COMPLEX_ACOSH
00164     return std::acosh (z);
00165 #else
00166     return log (z + sqrt (z * z - 1.0));
00167 #endif
00168 }
00169 
00170 
00175 nr_complex_t asinh (const nr_complex_t z) {
00176 #ifdef HAVE_CXX_COMPLEX_ACOSH
00177     return std::asinh (z);
00178 #else
00179     return log (z + sqrt (z * z + 1.0));
00180 #endif
00181 }
00182 
00183 
00188 nr_complex_t atanh (const nr_complex_t z) {
00189 #ifdef HAVE_CXX_COMPLEX_ATANH
00190     return std::atanh (z);
00191 #else
00192     return 0.5 * log ( 2.0 / (1.0 - z) - 1.0);
00193 #endif
00194 }
00195 
00196 
00197 //
00198 // transcendentals overloads
00199 //
00200 
00205 nr_complex_t exp (const nr_complex_t z)
00206 {
00207     nr_double_t mag = exp (real (z));
00208     return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
00209 }
00210 
00215 nr_complex_t log (const nr_complex_t z)
00216 {
00217     nr_double_t phi = arg (z);
00218     return nr_complex_t (log (abs (z)), phi);
00219 }
00220 
00225 nr_complex_t log10 (const nr_complex_t z)
00226 {
00227     nr_double_t phi = arg (z);
00228     return nr_complex_t (log10 (abs (z)), phi * log10e);
00229 }
00230 
00231 
00238 nr_complex_t pow (const nr_complex_t z, const nr_double_t d) {
00239     return std::pow (z, d);
00240 }
00241 
00248 nr_complex_t pow (const nr_double_t d, const nr_complex_t z) {
00249     return std::pow (d, z);
00250 }
00251 
00258 nr_complex_t pow (const nr_complex_t z1, const nr_complex_t z2) {
00259     return std::pow (z1, z2);
00260 }
00261 
00262 
00271 nr_complex_t sqrt (const nr_complex_t z)
00272 {
00273     return std::sqrt (z);
00274 }
00275 
00276 
00283 nr_double_t norm (const nr_complex_t z)
00284 {
00285     return std::norm (z);
00286 }
00287 
00288 
00289 //
00290 //  Qucs extra math functions
00291 //
00292 
00298 nr_complex_t cot (const nr_complex_t z)
00299 {
00300     nr_double_t r = 2.0 * std::real (z);
00301     nr_double_t i = 2.0 * std::imag (z);
00302     return nr_complex_t (0.0, 1.0) + nr_complex_t (0.0, 2.0) / (std::polar (std::exp (-i), r) - 1.0);
00303 }
00304 
00310 nr_complex_t acot (const nr_complex_t z)
00311 {
00312     return nr_complex_t (0.0, -0.5) * std::log (nr_complex_t (0, 2) / (z - nr_complex_t (0, 1)) + 1.0);
00313 }
00314 
00320 nr_complex_t coth (const nr_complex_t z)
00321 {
00322     nr_double_t r = 2.0 * std::real (z);
00323     nr_double_t i = 2.0 * std::imag (z);
00324     return 1.0 + 2.0 / (std::polar (std::exp (r), i) - 1.0);
00325 }
00326 
00332 nr_complex_t acoth (const nr_complex_t z)
00333 {
00334     return 0.5 * std::log (2.0 / (z - 1.0) + 1.0);
00335 }
00336 
00337 
00343 nr_complex_t sech (const nr_complex_t z)
00344 {
00345     return (1.0 / std::cosh (z));
00346 }
00347 
00354 nr_complex_t asech (const nr_complex_t z)
00355 {
00356     return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z);
00357 }
00358 
00364 nr_complex_t     cosech (const nr_complex_t z)
00365 {
00366     return (1.0 / std::sinh(z));
00367 }
00368 
00377 nr_complex_t atan2 (const nr_complex_t y, const nr_complex_t x)
00378 {
00379     nr_complex_t a = qucs::atan (y / x);  // std::atan missing on OSX 10.6
00380     return real (x) > 0.0 ? a : -a;
00381 }
00382 
00383 
00384 //
00385 // extra math
00386 //
00387 
00393 nr_complex_t log2 (const nr_complex_t z)
00394 {
00395 #ifndef HAVE_CXX_COMPLEX_LOG2
00396     nr_double_t phi = std::arg (z);
00397     return nr_complex_t (std::log (std::abs (z)) * log2e, phi * log2e);
00398 #else
00399     return std::log2 (z);
00400 #endif
00401 }
00402 
00416 nr_complex_t signum (const nr_complex_t z)
00417 {
00418     if (z == 0.0) return 0;
00419     return z / abs (z);
00420 }
00421 
00435 nr_complex_t sign (const nr_complex_t z)
00436 {
00437     if (z == 0.0) return nr_complex_t (1);
00438     return z / abs (z);
00439 }
00440 
00441 
00448 nr_complex_t    sinc (const nr_complex_t z)
00449 {
00450     if (real(z) == 0.0 && imag(z)) return 1;
00451     return std::sin (z) / z;
00452 }
00453 
00465 nr_double_t xhypot (const nr_complex_t a, const nr_complex_t b)
00466 {
00467     nr_double_t c = norm (a);
00468     nr_double_t d = norm (b);
00469     if (c > d)
00470         return abs (a) * std::sqrt (1.0 + d / c);
00471     else if (d == 0.0)
00472         return 0.0;
00473     else
00474         return abs (b) * std::sqrt (1.0 + c / d);
00475 }
00476 
00478 nr_double_t xhypot (nr_double_t a, nr_complex_t b)
00479 {
00480     return xhypot (nr_complex_t (a), b);
00481 }
00482 
00484 nr_double_t xhypot (nr_complex_t a, nr_double_t b)
00485 {
00486     return xhypot (a, nr_complex_t (b));
00487 }
00488 
00489 
00496 nr_complex_t round (const nr_complex_t z)
00497 {
00498     nr_double_t zreal = real (z);
00499     nr_double_t zimag = imag (z);
00500     // qucs::round resolved for double
00501     zreal = qucs::round (zreal);
00502     zimag = qucs::round (zimag);
00503     return nr_complex_t (zreal, zimag);
00504 }
00505 
00506 
00512 nr_complex_t trunc (const nr_complex_t z)
00513 {
00514   nr_double_t zreal = real (z);
00515   nr_double_t zimag = imag (z);
00516   // qucs::round resolved for double
00517   zreal = qucs::trunc (zreal);
00518   zimag = qucs::trunc (zimag);
00519   return nr_complex_t (zreal, zimag);
00520 }
00521 
00522 
00528 nr_double_t dB (const nr_complex_t z)
00529 {
00530   return 10.0 * std::log10 (std::norm (z));
00531 }
00532 
00533 
00539 nr_complex_t limexp (const nr_complex_t z)
00540 {
00541   nr_double_t mag = qucs::limexp (real (z));
00542   return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
00543 }
00544 
00545 
00551 nr_complex_t polar (const nr_double_t mag, const nr_double_t ang )
00552 {
00553 #ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
00554     return std::polar (mag, ang);
00555 #else
00556     return nr_complex_t (mag * cos (ang), mag * sin (ang));
00557 #endif
00558 }
00559 
00566 nr_complex_t polar (const nr_complex_t a, const nr_complex_t p)
00567 {
00568 #ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
00569     return std::polar (a, p);
00570 #else
00571     return a * exp (nr_complex_t (imag (p),-real (p)));
00572 #endif
00573 }
00574 
00575 
00581 nr_complex_t ztor (const nr_complex_t z, nr_complex_t zref) {
00582     return (z - zref) / (z + zref);
00583 }
00584 
00590 nr_complex_t rtoz (const nr_complex_t r, nr_complex_t zref) {
00591     return zref * (1.0 + r) / (1.0 - r);
00592 }
00593 
00599 nr_complex_t ytor (const nr_complex_t y, nr_complex_t zref) {
00600     return (1.0 - y * zref) / (1.0 + y * zref);
00601 }
00602 
00608 nr_complex_t rtoy (const nr_complex_t r, nr_complex_t zref) {
00609     return (1.0 - r) / (1.0 + r) / zref;
00610 }
00611 
00612 
00613 
00614 
00615 
00623 nr_complex_t   floor (const nr_complex_t z) {
00624   return nr_complex_t (std::floor (real (z)), std::floor (imag (z)));
00625 }
00626 
00627 
00634 nr_complex_t ceil (const nr_complex_t z) {
00635   return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z)));
00636 }
00637 
00645 nr_complex_t fix (const nr_complex_t z) {
00646   nr_double_t x = real (z);
00647   nr_double_t y = imag (z);
00648   x = (x > 0) ? std::floor (x) : std::ceil (x);
00649   y = (y > 0) ? std::floor (y) : std::ceil (y);
00650   return nr_complex_t (x, y);
00651 }
00652 
00653 
00654 
00662 nr_complex_t    fmod (const nr_complex_t x, const nr_complex_t y) {
00663   nr_complex_t n = qucs::floor (x / y);
00664   return x - n * y;
00665 }
00666 
00667 
00673 nr_complex_t sqr (const nr_complex_t z) {
00674   nr_double_t r = real (z);
00675   nr_double_t i = imag (z);
00676   return nr_complex_t (r * r - i * i, 2 * r * i);
00677 }
00678 
00679 
00680 
00681 
00682 
00691 nr_complex_t step (const nr_complex_t z)
00692 {
00693     nr_double_t x = real (z);
00694     nr_double_t y = imag (z);
00695     if (x < 0.0)
00696         x = 0.0;
00697     else if (x > 0.0)
00698         x = 1.0;
00699     else
00700         x = 0.5;
00701     if (y < 0.0)
00702         y = 0.0;
00703     else if (y > 0.0)
00704         y = 1.0;
00705     else
00706         y = 0.5;
00707     return nr_complex_t (x, y);
00708 }
00709 
00710 
00711 //using namespace fspecial;
00712 
00713 
00714 // === bessel ===
00715 
00716 /* FIXME : what about libm jn, yn, isn't that enough? */
00717 
00718 nr_complex_t cbesselj (unsigned int, nr_complex_t);
00719 
00720 #include "cbesselj.cpp"
00721 
00729 nr_complex_t jn (const int n, const nr_complex_t z)
00730 {
00731     return cbesselj (n, z);
00732 }
00733 
00734 
00742 nr_complex_t yn (const int n, const nr_complex_t z)
00743 {
00744     return nr_complex_t (::yn (n, std::real (z)), 0);
00745 }
00746 
00747 
00754 nr_complex_t i0 (const nr_complex_t z)
00755 {
00756     return nr_complex_t (fspecial::i0 (std::real (z)), 0);
00757 }
00758 
00759 
00766 nr_complex_t erf (const nr_complex_t z)
00767 {
00768 #ifdef HAVE_STD_ERF
00769   nr_double_t zerf = std::erf (std::real (z)); // c++11
00770 #elif HAVE_ERF
00771   nr_double_t zerf = ::erf (std::real (z));
00772 #else
00773   nr_double_t zerf = fspecial::erf (std::real (z));
00774 #endif
00775   return nr_complex_t (zerf, 0);
00776 }
00777 
00784 nr_complex_t erfc (const nr_complex_t z)
00785 {
00786 #ifdef HAVE_STD_ERF
00787   nr_double_t zerfc = std::erfc (std::real (z)); // c++11
00788 #elif HAVE_ERFC
00789   nr_double_t zerfc = ::erfc (std::real (z));
00790 #else
00791   nr_double_t zerfc = fspecial::erfc (std::real (z));
00792 #endif
00793   return nr_complex_t (zerfc, 0);
00794 }
00795 
00802 nr_complex_t erfinv (const nr_complex_t z)
00803 {
00804     return nr_complex_t (fspecial::erfinv (std::real (z)), 0);
00805 }
00806 
00813 nr_complex_t erfcinv (const nr_complex_t z)
00814 {
00815     return nr_complex_t (fspecial::erfcinv (std::real (z)), 0);
00816 }
00817 
00818 
00819 
00820 // ========================
00821 
00822 
00826 nr_complex_t operator%(const nr_complex_t z1, const nr_complex_t z2)
00827 {
00828     return z1 - z2 * floor (z1 / z2);
00829 }
00830 
00834 nr_complex_t operator%(const nr_complex_t z1, const nr_double_t r2)
00835 {
00836     return z1 - r2 * floor (z1 / r2);
00837 }
00838 
00842 nr_complex_t operator%(const nr_double_t r1, const nr_complex_t z2)
00843 {
00844     return r1 - z2 * floor (r1 / z2);
00845 }
00846 
00847 
00854 bool operator==(const nr_complex_t z1, const nr_complex_t z2)
00855 {
00856     return (std::real (z1) == std::real (z2)) && (std::imag (z1) == std::imag (z2));
00857 }
00858 
00865 bool operator!=(const nr_complex_t z1, const nr_complex_t z2)
00866 {
00867     return (std::real (z1) != std::real (z2)) || (std::imag (z1) != std::imag (z2));
00868 }
00869 
00873 bool operator>=(const nr_complex_t z1, const nr_complex_t z2)
00874 {
00875     return norm (z1) >= norm (z2);
00876 }
00877 
00881 bool operator<=(const nr_complex_t z1, const nr_complex_t z2)
00882 {
00883     return norm (z1) <= norm (z2);
00884 }
00885 
00889 bool operator>(const nr_complex_t z1, const nr_complex_t z2)
00890 {
00891     return norm (z1) > norm (z2);
00892 }
00893 
00897 bool operator<(const nr_complex_t z1, const nr_complex_t z2)
00898 {
00899     return norm (z1) < norm (z2);
00900 }
00901 
00907 nr_double_t rad2deg (const nr_complex_t x) {
00908   return rad2deg (real(x));
00909 }
00910 
00916 nr_double_t deg2rad (const nr_complex_t x) {
00917   return deg2rad (real(x));
00918 }
00919 
00920 } // namespace qucs
00921