Qucs-core
0.0.19
|
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