Qucs-core  0.0.18
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
complex.cpp
Go to the documentation of this file.
1 /*
2  * complex.cpp - complex number class implementation
3  *
4  * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org>
5  * Copyright (C) 2006 Gunther Kraut <gn.kraut@t-online.de>
6  *
7  * This is free software; you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation; either version 2, or (at your option)
10  * any later version.
11  *
12  * This software is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this package; see the file COPYING. If not, write to
19  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
20  * Boston, MA 02110-1301, USA.
21  *
22  * $Id$
23  *
24  */
25 
30 #if HAVE_CONFIG_H
31 # include <config.h>
32 #endif
33 
34 #include <cmath>
35 #include <assert.h>
36 #include <errno.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 
40 #include "constants.h"
41 #include "precision.h"
42 #include "complex.h"
43 #include "consts.h"
44 #include "fspecial.h"
45 
46 
47 namespace qucs {
48 
49 //
50 // trigonometric complex
51 //
52 
58  return std::cos (z);
59 }
60 
61 
67  return std::sin (z);
68 }
69 
70 
76  return std::tan (z);
77 }
78 
79 
85 #ifdef HAVE_CXX_COMPLEX_ACOS
86  return std::acos (z);
87 #else
88  // missing on OSX 10.6
89  nr_double_t r = real (z);
90  nr_double_t i = imag (z);
91  nr_complex_t y = sqrt (z * z - 1.0);
92  if (r * i < 0.0) y = -y;
93  return nr_complex_t (0, -1.0) * log (z + y);
94 #endif
95 }
96 
97 
103 #ifdef HAVE_CXX_COMPLEX_ASIN
104  return std::asin (z);
105 #else
106  // missing on OSX 10.6
107  nr_double_t r = real (z);
108  nr_double_t i = imag (z);
109  return nr_complex_t (0.0, -1.0) * log (nr_complex_t (-i, r) + sqrt (1.0 - z * z));
110 #endif
111 }
112 
118 #ifdef HAVE_CXX_COMPLEX_ATAN
119  return std::atan (z);
120 #else
121  // missing on OSX 10.6
122  return nr_complex_t (0.0, -0.5) * log (nr_complex_t (0, 2) / (z + nr_complex_t (0, 1)) - 1.0);
123 #endif
124 }
125 
126 
127 //
128 // hyperbolic complex
129 //
130 
136  return std::cosh (z);
137 }
138 
139 
145  return std::sinh (z);
146 }
147 
148 
154  return std::tanh (z);
155 }
156 
157 
163 #ifdef HAVE_CXX_COMPLEX_ACOSH
164  return std::acosh (z);
165 #else
166  return log (z + sqrt (z * z - 1.0));
167 #endif
168 }
169 
170 
176 #ifdef HAVE_CXX_COMPLEX_ACOSH
177  return std::asinh (z);
178 #else
179  return log (z + sqrt (z * z + 1.0));
180 #endif
181 }
182 
183 
189 #ifdef HAVE_CXX_COMPLEX_ATANH
190  return std::atanh (z);
191 #else
192  return 0.5 * log ( 2.0 / (1.0 - z) - 1.0);
193 #endif
194 }
195 
196 
197 //
198 // transcendentals overloads
199 //
200 
206 {
207  nr_double_t mag = exp (real (z));
208  return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
209 }
210 
216 {
217  nr_double_t phi = arg (z);
218  return nr_complex_t (log (abs (z)), phi);
219 }
220 
226 {
227  nr_double_t phi = arg (z);
228  return nr_complex_t (log10 (abs (z)), phi * M_LOG10E);
229 }
230 
231 
238 nr_complex_t pow (const nr_complex_t z, const nr_double_t d) {
239  return std::pow (z, d);
240 }
241 
248 nr_complex_t pow (const nr_double_t d, const nr_complex_t z) {
249  return std::pow (d, z);
250 }
251 
259  return std::pow (z1, z2);
260 }
261 
262 
272 {
273  return std::sqrt (z);
274 }
275 
276 
283 nr_double_t norm (const nr_complex_t z)
284 {
285  return std::norm (z);
286 }
288 
289 //
290 // Qucs extra math functions
291 //
292 
299 {
300  nr_double_t r = 2.0 * std::real (z);
301  nr_double_t i = 2.0 * std::imag (z);
302  return nr_complex_t (0.0, 1.0) + nr_complex_t (0.0, 2.0) / (std::polar (std::exp (-i), r) - 1.0);
303 }
304 
311 {
312  return nr_complex_t (0.0, -0.5) * std::log (nr_complex_t (0, 2) / (z - nr_complex_t (0, 1)) + 1.0);
313 }
314 
321 {
322  nr_double_t r = 2.0 * std::real (z);
323  nr_double_t i = 2.0 * std::imag (z);
324  return 1.0 + 2.0 / (std::polar (std::exp (r), i) - 1.0);
325 }
326 
333 {
334  return 0.5 * std::log (2.0 / (z - 1.0) + 1.0);
335 }
336 
337 
344 {
345  return (1.0 / std::cosh (z));
346 }
347 
355 {
356  return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z);
357 }
358 
365 {
366  return (1.0 / std::sinh(z));
367 }
368 
378 {
379  nr_complex_t a = qucs::atan (y / x); // std::atan missing on OSX 10.6
380  return real (x) > 0.0 ? a : -a;
381 }
382 
383 
384 //
385 // extra math
386 //
387 
394 {
395 #ifndef HAVE_CXX_COMPLEX_LOG2
396  nr_double_t phi = std::arg (z);
397  return nr_complex_t (std::log (std::abs (z)) * M_LOG2E, phi * M_LOG2E);
398 #else
399  return std::log2 (z);
400 #endif
401 }
402 
417 {
418  if (z == 0.0) return 0;
419  return z / abs (z);
420 }
421 
436 {
437  if (z == 0.0) return nr_complex_t (1);
438  return z / abs (z);
439 }
440 
441 
449 {
450  if (real(z) == 0.0 && imag(z)) return 1;
451  return std::sin (z) / z;
452 }
453 
465 nr_double_t xhypot (const nr_complex_t a, const nr_complex_t b)
466 {
467  nr_double_t c = norm (a);
468  nr_double_t d = norm (b);
469  if (c > d)
470  return abs (a) * std::sqrt (1.0 + d / c);
471  else if (d == 0.0)
472  return 0.0;
473  else
474  return abs (b) * std::sqrt (1.0 + c / d);
475 }
476 
478 nr_double_t xhypot (nr_double_t a, nr_complex_t b)
479 {
480  return xhypot (nr_complex_t (a), b);
481 }
482 
484 nr_double_t xhypot (nr_complex_t a, nr_double_t b)
485 {
486  return xhypot (a, nr_complex_t (b));
487 }
488 
489 
497 {
498  nr_double_t zreal = real (z);
499  nr_double_t zimag = imag (z);
500  // qucs::round resolved for double
501  zreal = qucs::round (zreal);
502  zimag = qucs::round (zimag);
503  return nr_complex_t (zreal, zimag);
504 }
505 
506 
513 {
514  nr_double_t zreal = real (z);
515  nr_double_t zimag = imag (z);
516  // qucs::round resolved for double
517  zreal = qucs::trunc (zreal);
518  zimag = qucs::trunc (zimag);
519  return nr_complex_t (zreal, zimag);
520 }
521 
522 
528 nr_double_t dB (const nr_complex_t z)
529 {
530  return 10.0 * std::log10 (std::norm (z));
531 }
532 
533 
540 {
541  nr_double_t mag = qucs::limexp (real (z));
542  return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
543 }
544 
545 
551 nr_complex_t polar (const nr_double_t mag, const nr_double_t ang )
552 {
553 #ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
554  return std::polar (mag, ang);
555 #else
556  return nr_complex_t (mag * cos (ang), mag * sin (ang));
557 #endif
558 }
559 
567 {
568 #ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
569  return std::polar (a, p);
570 #else
571  return a * exp (nr_complex_t (imag (p),-real (p)));
572 #endif
573 }
574 
575 
582  return (z - zref) / (z + zref);
583 }
584 
591  return zref * (1.0 + r) / (1.0 - r);
592 }
593 
600  return (1.0 - y * zref) / (1.0 + y * zref);
601 }
602 
609  return (1.0 - r) / (1.0 + r) / zref;
610 }
611 
612 
613 
614 
615 
624  return nr_complex_t (std::floor (real (z)), std::floor (imag (z)));
625 }
626 
627 
635  return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z)));
636 }
637 
646  nr_double_t x = real (z);
647  nr_double_t y = imag (z);
648  x = (x > 0) ? std::floor (x) : std::ceil (x);
649  y = (y > 0) ? std::floor (y) : std::ceil (y);
650  return nr_complex_t (x, y);
651 }
652 
653 
654 
663  nr_complex_t n = qucs::floor (x / y);
664  return x - n * y;
665 }
666 
667 
674  nr_double_t r = real (z);
675  nr_double_t i = imag (z);
676  return nr_complex_t (r * r - i * i, 2 * r * i);
677 }
678 
679 
680 
681 
682 
692 {
693  nr_double_t x = real (z);
694  nr_double_t y = imag (z);
695  if (x < 0.0)
696  x = 0.0;
697  else if (x > 0.0)
698  x = 1.0;
699  else
700  x = 0.5;
701  if (y < 0.0)
702  y = 0.0;
703  else if (y > 0.0)
704  y = 1.0;
705  else
706  y = 0.5;
707  return nr_complex_t (x, y);
708 }
709 
710 
711 //using namespace fspecial;
712 
713 
714 // === bessel ===
715 
716 /* FIXME : what about libm jn, yn, isn't that enough? */
717 
718 nr_complex_t cbesselj (unsigned int, nr_complex_t);
719 
720 #include "cbesselj.cpp"
721 
729 nr_complex_t jn (const int n, const nr_complex_t z)
730 {
731  return cbesselj (n, z);
732 }
733 
734 
742 nr_complex_t yn (const int n, const nr_complex_t z)
743 {
744  return nr_complex_t (::yn (n, std::real (z)), 0);
745 }
746 
747 
755 {
756  return nr_complex_t (fspecial::i0 (std::real (z)), 0);
757 }
758 
759 
767 {
768 #ifdef HAVE_STD_ERF
769  nr_double_t zerf = std::erf (std::real (z)); // c++11
770 #elif HAVE_ERF
771  nr_double_t zerf = ::erf (std::real (z));
772 #else
773  nr_double_t zerf = fspecial::erf (std::real (z));
774 #endif
775  return nr_complex_t (zerf, 0);
776 }
777 
785 {
786 #ifdef HAVE_STD_ERF
787  nr_double_t zerfc = std::erfc (std::real (z)); // c++11
788 #elif HAVE_ERFC
789  nr_double_t zerfc = ::erfc (std::real (z));
790 #else
791  nr_double_t zerfc = fspecial::erfc (std::real (z));
792 #endif
793  return nr_complex_t (zerfc, 0);
794 }
795 
803 {
804  return nr_complex_t (fspecial::erfinv (std::real (z)), 0);
805 }
806 
814 {
815  return nr_complex_t (fspecial::erfcinv (std::real (z)), 0);
816 }
817 
818 
819 
820 // ========================
821 
822 
827 {
828  return z1 - z2 * floor (z1 / z2);
829 }
830 
834 nr_complex_t operator%(const nr_complex_t z1, const nr_double_t r2)
835 {
836  return z1 - r2 * floor (z1 / r2);
837 }
838 
842 nr_complex_t operator%(const nr_double_t r1, const nr_complex_t z2)
843 {
844  return r1 - z2 * floor (r1 / z2);
845 }
846 
847 
854 bool operator==(const nr_complex_t z1, const nr_complex_t z2)
855 {
856  return (std::real (z1) == std::real (z2)) && (std::imag (z1) == std::imag (z2));
857 }
858 
865 bool operator!=(const nr_complex_t z1, const nr_complex_t z2)
866 {
867  return (std::real (z1) != std::real (z2)) || (std::imag (z1) != std::imag (z2));
868 }
869 
873 bool operator>=(const nr_complex_t z1, const nr_complex_t z2)
874 {
875  return norm (z1) >= norm (z2);
876 }
877 
881 bool operator<=(const nr_complex_t z1, const nr_complex_t z2)
882 {
883  return norm (z1) <= norm (z2);
884 }
885 
889 bool operator>(const nr_complex_t z1, const nr_complex_t z2)
890 {
891  return norm (z1) > norm (z2);
892 }
893 
897 bool operator<(const nr_complex_t z1, const nr_complex_t z2)
898 {
899  return norm (z1) < norm (z2);
900 }
901 
902 } // namespace qucs
903 
nr_double_t erf(nr_double_t)
Definition: fspecial.cpp:374
std::complex< nr_double_t > nr_complex_t
Definition: complex.h:31
nr_complex_t pow(const nr_complex_t z1, const nr_complex_t z2)
Compute complex power function.
Definition: complex.cpp:258
matrix real(matrix a)
Real part matrix.
Definition: matrix.cpp:568
matrix abs(matrix a)
Computes magnitude of each matrix element.
Definition: matrix.cpp:531
nr_complex_t erf(const nr_complex_t z)
Error function.
Definition: complex.cpp:766
#define M_LOG2E
Binary logartihm of Euler's constant ( )
Definition: consts.h:79
bool operator<=(const nr_complex_t z1, const nr_complex_t z2)
Inferior of equal.
Definition: complex.cpp:881
nr_complex_t ceil(const nr_complex_t z)
Complex ceil Ceil is the smallest integral value not less than argument Apply ceil to real and imagin...
Definition: complex.cpp:634
nr_double_t erfinv(nr_double_t)
Definition: fspecial.cpp:388
nr_complex_t coth(const nr_complex_t z)
Compute complex hyperbolic cotangent.
Definition: complex.cpp:320
nr_complex_t pow(const nr_complex_t z, const nr_double_t d)
Compute power function with real exponent.
Definition: complex.cpp:238
nr_complex_t cos(const nr_complex_t z)
Compute complex cosine.
Definition: complex.cpp:57
nr_complex_t step(const nr_complex_t z)
Heaviside step function for complex number.
Definition: complex.cpp:691
nr_complex_t signum(const nr_complex_t z)
complex signum function
Definition: complex.cpp:416
nr_complex_t atan(const nr_complex_t z)
Compute complex arc tangent.
Definition: complex.cpp:117
nr_complex_t cot(const nr_complex_t z)
Compute complex cotangent.
Definition: complex.cpp:298
Global physical constants header file.
nr_complex_t asech(const nr_complex_t z)
Compute complex argument hyperbolic secant.
Definition: complex.cpp:354
nr_double_t erfcinv(nr_double_t)
Definition: fspecial.cpp:589
nr_complex_t acot(const nr_complex_t z)
Compute complex arc cotangent.
Definition: complex.cpp:310
n
Definition: parse_citi.y:147
r
Definition: parse_mdl.y:515
bool operator>(const nr_complex_t z1, const nr_complex_t z2)
Superior.
Definition: complex.cpp:889
nr_complex_t acosh(const nr_complex_t z)
Compute complex arc hyperbolic cosine.
Definition: complex.cpp:162
nr_complex_t rtoz(const nr_complex_t r, nr_complex_t zref)
Converts reflexion coefficient to impedance.
Definition: complex.cpp:590
nr_complex_t sign(const nr_complex_t z)
complex sign function
Definition: complex.cpp:435
nr_complex_t asinh(const nr_complex_t z)
Compute complex arc hyperbolic sine.
Definition: complex.cpp:175
nr_complex_t cosh(const nr_complex_t z)
Compute complex hyperbolic cosine.
Definition: complex.cpp:135
i
Definition: parse_mdl.y:516
nr_complex_t sqr(const nr_complex_t z)
Square of complex number.
Definition: complex.cpp:673
matrix imag(matrix a)
Imaginary part matrix.
Definition: matrix.cpp:581
nr_complex_t fix(const nr_complex_t z)
Complex fix.
Definition: complex.cpp:645
nr_complex_t sqrt(const nr_complex_t z)
Compute principal value of square root.
Definition: complex.cpp:271
nr_double_t xhypot(const nr_complex_t a, const nr_complex_t b)
Euclidean distance function for complex argument.
Definition: complex.cpp:465
nr_complex_t trunc(const nr_complex_t z)
Complex trunc Apply round to integer, towards zero to real and imaginary part.
Definition: complex.cpp:512
nr_complex_t acos(const nr_complex_t z)
Compute complex arc cosine.
Definition: complex.cpp:84
nr_complex_t log2(const nr_complex_t z)
Compute principal value of binary logarithm of z.
Definition: complex.cpp:393
nr_complex_t fmod(const nr_complex_t x, const nr_complex_t y)
Complex fmod Apply fmod to the complex z.
Definition: complex.cpp:662
nr_complex_t polar(const nr_complex_t a, const nr_complex_t p)
Extension of polar construction to complex.
Definition: complex.cpp:566
nr_complex_t cosech(const nr_complex_t z)
Compute complex argument hyperbolic cosec.
Definition: complex.cpp:364
nr_complex_t yn(const int n, const nr_complex_t z)
Bessel function of second kind.
Definition: complex.cpp:742
nr_complex_t tanh(const nr_complex_t z)
Compute complex hyperbolic tangent.
Definition: complex.cpp:153
nr_complex_t sin(const nr_complex_t z)
Compute complex sine.
Definition: complex.cpp:66
nr_complex_t ztor(const nr_complex_t z, nr_complex_t zref)
Converts impedance to reflexion coefficient.
Definition: complex.cpp:581
nr_double_t dB(const nr_complex_t z)
Magnitude in dB Compute .
Definition: complex.cpp:528
nr_complex_t asin(const nr_complex_t z)
Compute complex arc sine.
Definition: complex.cpp:102
nr_complex_t acoth(const nr_complex_t z)
Compute complex argument hyperbolic cotangent.
Definition: complex.cpp:332
Global math constants header file.
bool operator!=(const nr_complex_t z1, const nr_complex_t z2)
Inequality of two complex.
Definition: complex.cpp:865
x
Definition: parse_mdl.y:498
nr_complex_t log10(const nr_complex_t z)
Compute principal value of decimal logarithm of z.
Definition: complex.cpp:225
nr_complex_t floor(const nr_complex_t z)
Complex floor.
Definition: complex.cpp:623
nr_complex_t sech(const nr_complex_t z)
Compute complex hyperbolic secant.
Definition: complex.cpp:343
nr_complex_t erfc(const nr_complex_t z)
Complementart error function.
Definition: complex.cpp:784
nr_complex_t cbesselj(unsigned int, nr_complex_t)
Main entry point for besselj function.
Definition: complex.cpp:287
nr_complex_t ytor(const nr_complex_t y, nr_complex_t zref)
Converts admittance to reflexion coefficient.
Definition: complex.cpp:599
nr_complex_t erfinv(const nr_complex_t z)
Inverse of error function.
Definition: complex.cpp:802
nr_complex_t rtoy(const nr_complex_t r, nr_complex_t zref)
Converts reflexion coefficient to admittance.
Definition: complex.cpp:608
nr_complex_t i0(const nr_complex_t z)
Modified Bessel function of first kind.
Definition: complex.cpp:754
compute complex bessel J function
nr_complex_t atan2(const nr_complex_t y, const nr_complex_t x)
Compute complex arc tangent fortran like function.
Definition: complex.cpp:377
nr_complex_t sinc(const nr_complex_t z)
Cardinal sine.
Definition: complex.cpp:448
nr_complex_t erfcinv(const nr_complex_t z)
Inverse of complementart error function.
Definition: complex.cpp:813
bool operator==(const nr_complex_t z1, const nr_complex_t z2)
Equality of two complex.
Definition: complex.cpp:854
nr_complex_t limexp(const nr_complex_t z)
Compute limited complex exponential.
Definition: complex.cpp:539
nr_double_t norm(const nr_complex_t z)
Compute euclidian norm of complex number.
Definition: complex.cpp:283
nr_complex_t sinh(const nr_complex_t z)
Compute complex hyperbolic sine.
Definition: complex.cpp:144
y
Definition: parse_mdl.y:499
bool operator>=(const nr_complex_t z1, const nr_complex_t z2)
Superior of equal.
Definition: complex.cpp:873
nr_complex_t jn(const int n, const nr_complex_t z)
Bessel function of first kind.
Definition: complex.cpp:729
nr_double_t erfc(nr_double_t)
Definition: fspecial.cpp:330
nr_complex_t tan(const nr_complex_t z)
Compute complex tangent.
Definition: complex.cpp:75
nr_complex_t exp(const nr_complex_t z)
Compute complex exponential.
Definition: complex.cpp:205
nr_complex_t operator%(const nr_complex_t z1, const nr_complex_t z2)
Modulo.
Definition: complex.cpp:826
#define M_LOG10E
Decimal logartihm of Euler's constant ( )
Definition: consts.h:75
nr_complex_t polar(const nr_double_t mag, const nr_double_t ang)
Construct a complex number using polar notation.
Definition: complex.cpp:551
zref
Definition: parse_zvr.y:130
nr_complex_t log(const nr_complex_t z)
Compute principal value of natural logarithm of z.
Definition: complex.cpp:215
bool operator<(const nr_complex_t z1, const nr_complex_t z2)
Inferior.
Definition: complex.cpp:897
nr_complex_t round(const nr_complex_t z)
Complex round Round is the nearest integral value Apply round to real and imaginary part...
Definition: complex.cpp:496
nr_double_t i0(nr_double_t)
Definition: fspecial.cpp:501
matrix arg(matrix a)
Computes the argument of each matrix element.
Definition: matrix.cpp:555
nr_complex_t atanh(const nr_complex_t z)
Compute complex arc hyperbolic tangent.
Definition: complex.cpp:188