Branch data Line data Source code
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 : :
26 : : /*!\file complex.cpp
27 : : Implements complex number class and functions
28 : : */
29 : :
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 : :
53 : : /*! \brief Compute complex cosine
54 : : \param[in] z complex angle
55 : : \return cosine of z
56 : : */
57 : 0 : nr_complex_t cos (const nr_complex_t z) {
58 : 0 : return std::cos (z);
59 : : }
60 : :
61 : :
62 : : /*! \brief Compute complex sine
63 : : \param[in] z complex angle
64 : : \return sine of z
65 : : */
66 : 6306 : nr_complex_t sin (const nr_complex_t z) {
67 : 6306 : return std::sin (z);
68 : : }
69 : :
70 : :
71 : : /*! \brief Compute complex tangent
72 : : \param[in] z complex angle
73 : : \return tangent of z
74 : : */
75 : 0 : nr_complex_t tan (const nr_complex_t z) {
76 : 0 : return std::tan (z);
77 : : }
78 : :
79 : :
80 : : /*! \brief Compute complex arc cosine
81 : : \param[in] z complex arc
82 : : \return arc cosine of z
83 : : */
84 : 0 : nr_complex_t acos (const nr_complex_t z) {
85 : : #ifdef HAVE_CXX_COMPLEX_ACOS
86 : 0 : 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 : :
98 : : /*! \brief Compute complex arc sine
99 : : \param[in] z complex arc
100 : : \return arc sine of z
101 : : */
102 : 0 : nr_complex_t asin (const nr_complex_t z) {
103 : : #ifdef HAVE_CXX_COMPLEX_ASIN
104 : 0 : 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 : :
113 : : /*! \brief Compute complex arc tangent
114 : : \param[in] z complex arc
115 : : \return arc tangent of z
116 : : */
117 : 0 : nr_complex_t atan (const nr_complex_t z) {
118 : : #ifdef HAVE_CXX_COMPLEX_ATAN
119 : 0 : 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 : :
131 : : /*! \brief Compute complex hyperbolic cosine
132 : : \param[in] z complex arc
133 : : \return hyperbolic cosine of z
134 : : */
135 : 2053 : nr_complex_t cosh (const nr_complex_t z) {
136 : 2053 : return std::cosh (z);
137 : : }
138 : :
139 : :
140 : : /*! \brief Compute complex hyperbolic sine
141 : : \param[in] z complex arc
142 : : \return hyperbolic sine of z
143 : : */
144 : 4106 : nr_complex_t sinh (const nr_complex_t z) {
145 : 4106 : return std::sinh (z);
146 : : }
147 : :
148 : :
149 : : /*! \brief Compute complex hyperbolic tangent
150 : : \param[in] z complex arc
151 : : \return hyperbolic tangent of z
152 : : */
153 : 0 : nr_complex_t tanh (const nr_complex_t z) {
154 : 0 : return std::tanh (z);
155 : : }
156 : :
157 : :
158 : : /*! \brief Compute complex arc hyperbolic cosine
159 : : \param[in] z complex arc
160 : : \return arc hyperbolic cosine of z
161 : : */
162 : 0 : nr_complex_t acosh (const nr_complex_t z) {
163 : : #ifdef HAVE_CXX_COMPLEX_ACOSH
164 : 0 : return std::acosh (z);
165 : : #else
166 : : return log (z + sqrt (z * z - 1.0));
167 : : #endif
168 : : }
169 : :
170 : :
171 : : /*! \brief Compute complex arc hyperbolic sine
172 : : \param[in] z complex arc
173 : : \return arc hyperbolic sine of z
174 : : */
175 : 0 : nr_complex_t asinh (const nr_complex_t z) {
176 : : #ifdef HAVE_CXX_COMPLEX_ACOSH
177 : 0 : return std::asinh (z);
178 : : #else
179 : : return log (z + sqrt (z * z + 1.0));
180 : : #endif
181 : : }
182 : :
183 : :
184 : : /*! \brief Compute complex arc hyperbolic tangent
185 : : \param[in] z complex arc
186 : : \return arc hyperbolic tangent of z
187 : : */
188 : 0 : nr_complex_t atanh (const nr_complex_t z) {
189 : : #ifdef HAVE_CXX_COMPLEX_ATANH
190 : 0 : 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 : :
201 : : /*! \brief Compute complex exponential
202 : : \param[in] z complex number
203 : : \return exponential of z
204 : : */
205 : 0 : nr_complex_t exp (const nr_complex_t z)
206 : : {
207 : 0 : nr_double_t mag = exp (real (z));
208 : 0 : return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
209 : : }
210 : :
211 : : /*! \brief Compute principal value of natural logarithm of z
212 : : \param[in] z complex number
213 : : \return principal value of natural logarithm of z
214 : : */
215 : 0 : nr_complex_t log (const nr_complex_t z)
216 : : {
217 : 0 : nr_double_t phi = arg (z);
218 : 0 : return nr_complex_t (log (abs (z)), phi);
219 : : }
220 : :
221 : : /*! \brief Compute principal value of decimal logarithm of z
222 : : \param[in] z complex number
223 : : \return principal value of decimal logarithm of z
224 : : */
225 : 140 : nr_complex_t log10 (const nr_complex_t z)
226 : : {
227 : 140 : nr_double_t phi = arg (z);
228 : 140 : return nr_complex_t (log10 (abs (z)), phi * M_LOG10E);
229 : : }
230 : :
231 : :
232 : : /*!\brief Compute power function with real exponent
233 : :
234 : : \param[in] z complex mantisse
235 : : \param[in] d real exponent
236 : : \return z power d (\f$z^d\f$)
237 : : */
238 : 0 : nr_complex_t pow (const nr_complex_t z, const nr_double_t d) {
239 : 0 : return std::pow (z, d);
240 : : }
241 : :
242 : : /*!\brief Compute power function with complex exponent but real mantisse
243 : :
244 : : \param[in] d real mantisse
245 : : \param[in] z complex exponent
246 : : \return d power z (\f$d^z\f$)
247 : : */
248 : 0 : nr_complex_t pow (const nr_double_t d, const nr_complex_t z) {
249 : 0 : return std::pow (d, z);
250 : : }
251 : :
252 : : /*!\brief Compute complex power function
253 : :
254 : : \param[in] z1 complex mantisse
255 : : \param[in] z2 complex exponent
256 : : \return d power z (\f$z_1^{z_2}\f$)
257 : : */
258 : 0 : nr_complex_t pow (const nr_complex_t z1, const nr_complex_t z2) {
259 : 0 : return std::pow (z1, z2);
260 : : }
261 : :
262 : :
263 : : /*!\brief Compute principal value of square root
264 : :
265 : : Compute the square root of a given complex number (except negative
266 : : real), and with a branch cut along the negative real axis.
267 : :
268 : : \param[in] z complex number
269 : : \return principal value of square root z
270 : : */
271 : 3759 : nr_complex_t sqrt (const nr_complex_t z)
272 : : {
273 : 3759 : return std::sqrt (z);
274 : : }
275 : :
276 : :
277 : : /*!\brief Compute euclidian norm of complex number
278 : :
279 : : Compute \f$(\Re\mathrm{e}\;z )^2+ (\Im\mathrm{m}\;z)^2=|z|^2\f$
280 : : \param[in] z Complex number
281 : : \return Euclidian norm of z
282 : : */
283 : 140560 : nr_double_t norm (const nr_complex_t z)
284 : : {
285 : 140560 : return std::norm (z);
286 : : }
287 : :
288 : :
289 : : //
290 : : // Qucs extra math functions
291 : : //
292 : :
293 : : /*!\brief Compute complex cotangent
294 : :
295 : : \param[in] z complex angle
296 : : \return cotangent of z
297 : : */
298 : 0 : nr_complex_t cot (const nr_complex_t z)
299 : : {
300 : 0 : nr_double_t r = 2.0 * std::real (z);
301 : 0 : nr_double_t i = 2.0 * std::imag (z);
302 [ # # ][ # # ]: 0 : 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 : :
305 : : /*!\brief Compute complex arc cotangent
306 : :
307 : : \param[in] z complex arc
308 : : \return arc cotangent of z
309 : : */
310 : 0 : nr_complex_t acot (const nr_complex_t z)
311 : : {
312 [ # # ][ # # ]: 0 : 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 : :
315 : : /*!\brief Compute complex hyperbolic cotangent
316 : :
317 : : \param[in] z complex angle
318 : : \return hyperbolic cotangent of z
319 : : */
320 : 0 : nr_complex_t coth (const nr_complex_t z)
321 : : {
322 : 0 : nr_double_t r = 2.0 * std::real (z);
323 : 0 : nr_double_t i = 2.0 * std::imag (z);
324 [ # # ]: 0 : return 1.0 + 2.0 / (std::polar (std::exp (r), i) - 1.0);
325 : : }
326 : :
327 : : /*!\brief Compute complex argument hyperbolic cotangent
328 : :
329 : : \param[in] z complex arc
330 : : \return argument hyperbolic cotangent of z
331 : : */
332 : 0 : nr_complex_t acoth (const nr_complex_t z)
333 : : {
334 [ # # ]: 0 : return 0.5 * std::log (2.0 / (z - 1.0) + 1.0);
335 : : }
336 : :
337 : :
338 : : /*!\brief Compute complex hyperbolic secant
339 : :
340 : : \param[in] z complex angle
341 : : \return hyperbolic secant of z
342 : : */
343 : 0 : nr_complex_t sech (const nr_complex_t z)
344 : : {
345 [ # # ]: 0 : return (1.0 / std::cosh (z));
346 : : }
347 : :
348 : : /*!\brief Compute complex argument hyperbolic secant
349 : :
350 : : \param[in] z complex arc
351 : : \return argument hyperbolic secant of z
352 : : \todo for symetry reason implement sech
353 : : */
354 : 0 : nr_complex_t asech (const nr_complex_t z)
355 : : {
356 [ # # ]: 0 : return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z);
357 : : }
358 : :
359 : : /*!\brief Compute complex argument hyperbolic cosec
360 : :
361 : : \param[in] z complex arc
362 : : \return argument hyperbolic cosec of z
363 : : */
364 : 0 : nr_complex_t cosech (const nr_complex_t z)
365 : : {
366 [ # # ]: 0 : return (1.0 / std::sinh(z));
367 : : }
368 : :
369 : : /*!\brief Compute complex arc tangent fortran like function
370 : :
371 : : atan2 is a two-argument function that computes the arc tangent of y / x
372 : : given y and x, but with a range of \f$(-\pi;\pi]\f$
373 : :
374 : : \param[in] z complex angle
375 : : \return arc tangent of z
376 : : */
377 : 0 : nr_complex_t atan2 (const nr_complex_t y, const nr_complex_t x)
378 : : {
379 [ # # ][ # # ]: 0 : nr_complex_t a = qucs::atan (y / x); // std::atan missing on OSX 10.6
380 [ # # ]: 0 : return real (x) > 0.0 ? a : -a;
381 : : }
382 : :
383 : :
384 : : //
385 : : // extra math
386 : : //
387 : :
388 : : /*!\brief Compute principal value of binary logarithm of z
389 : :
390 : : \param[in] z complex number
391 : : \return principal value of binary logarithm of z
392 : : */
393 : 0 : nr_complex_t log2 (const nr_complex_t z)
394 : : {
395 : : #ifndef HAVE_CXX_COMPLEX_LOG2
396 : 0 : nr_double_t phi = std::arg (z);
397 : 0 : 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 : :
403 : : /*!\brief complex signum function
404 : :
405 : : compute \f[
406 : : \mathrm{signum}\;z= \mathrm{signum} (re^{i\theta})
407 : : = \begin{cases}
408 : : 0 & \text{if } z=0 \\
409 : : e^{i\theta} & \text{else}
410 : : \end{cases}
411 : : \f]
412 : : \param[in] z complex number
413 : : \return signum of z
414 : : \todo Better implementation z/abs(z) is not really stable
415 : : */
416 : 0 : nr_complex_t signum (const nr_complex_t z)
417 : : {
418 [ # # ]: 0 : if (z == 0.0) return 0;
419 : 0 : return z / abs (z);
420 : : }
421 : :
422 : : /*!\brief complex sign function
423 : :
424 : : compute \f[
425 : : \mathrm{sign}\;z= \mathrm{sign} (re^{i\theta})
426 : : = \begin{cases}
427 : : 1 & \text{if } z=0 \\
428 : : e^{i\theta} & \text{else}
429 : : \end{cases}
430 : : \f]
431 : : \param[in] z complex number
432 : : \return sign of z
433 : : \todo Better implementation z/abs(z) is not really stable
434 : : */
435 : 0 : nr_complex_t sign (const nr_complex_t z)
436 : : {
437 [ # # ]: 0 : if (z == 0.0) return nr_complex_t (1);
438 : 0 : return z / abs (z);
439 : : }
440 : :
441 : :
442 : : /*!\brief Cardinal sine
443 : :
444 : : Compute \f$\mathrm{sinc}\;z=\frac{\sin z}{z}\f$
445 : : \param[in] z complex number
446 : : \return cardianal sine of z
447 : : */
448 : 0 : nr_complex_t sinc (const nr_complex_t z)
449 : : {
450 [ # # ][ # # ]: 0 : if (real(z) == 0.0 && imag(z)) return 1;
[ # # ]
451 [ # # ]: 0 : return std::sin (z) / z;
452 : : }
453 : :
454 : : /*!\brief Euclidean distance function for complex argument
455 : :
456 : : The xhypot() function returns \f$\sqrt{a^2+b^2}\f$.
457 : : This is the length of the hypotenuse of a right-angle triangle with sides
458 : : of length a and b, or the distance
459 : : of the point (a,b) from the origin.
460 : :
461 : : \param[in] a first length
462 : : \param[in] b second length
463 : : \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$
464 : : */
465 : 0 : nr_double_t xhypot (const nr_complex_t a, const nr_complex_t b)
466 : : {
467 : 0 : nr_double_t c = norm (a);
468 : 0 : nr_double_t d = norm (b);
469 [ # # ]: 0 : if (c > d)
470 : 0 : return abs (a) * std::sqrt (1.0 + d / c);
471 [ # # ]: 0 : else if (d == 0.0)
472 : 0 : return 0.0;
473 : : else
474 : 0 : return abs (b) * std::sqrt (1.0 + c / d);
475 : : }
476 : :
477 : : /*!\brief Euclidean distance function for a double b complex */
478 : 0 : nr_double_t xhypot (nr_double_t a, nr_complex_t b)
479 : : {
480 [ # # ]: 0 : return xhypot (nr_complex_t (a), b);
481 : : }
482 : :
483 : : /*!\brief Euclidean distance function for b double a complex */
484 : 0 : nr_double_t xhypot (nr_complex_t a, nr_double_t b)
485 : : {
486 [ # # ]: 0 : return xhypot (a, nr_complex_t (b));
487 : : }
488 : :
489 : :
490 : : /*!\brief Complex round
491 : : Round is the nearest integral value
492 : : Apply round to real and imaginary part
493 : : \param[in] z complex number
494 : : \return rounded complex number
495 : : */
496 : 0 : nr_complex_t round (const nr_complex_t z)
497 : : {
498 : 0 : nr_double_t zreal = real (z);
499 : 0 : nr_double_t zimag = imag (z);
500 : : // qucs::round resolved for double
501 : 0 : zreal = qucs::round (zreal);
502 : 0 : zimag = qucs::round (zimag);
503 : 0 : return nr_complex_t (zreal, zimag);
504 : : }
505 : :
506 : :
507 : : /*!\brief Complex trunc
508 : : Apply round to integer, towards zero to real and imaginary part
509 : : \param[in] z complex number
510 : : \return rounded complex number
511 : : */
512 : 0 : nr_complex_t trunc (const nr_complex_t z)
513 : : {
514 : 0 : nr_double_t zreal = real (z);
515 : 0 : nr_double_t zimag = imag (z);
516 : : // qucs::round resolved for double
517 : 0 : zreal = qucs::trunc (zreal);
518 : 0 : zimag = qucs::trunc (zimag);
519 : 0 : return nr_complex_t (zreal, zimag);
520 : : }
521 : :
522 : :
523 : : /*!\brief Magnitude in dB
524 : : Compute \f$10\log_{10} |z|^2=20\log_{10} |z|\f$
525 : : \param[in] z complex number
526 : : \return Magnitude in dB
527 : : */
528 : 0 : nr_double_t dB (const nr_complex_t z)
529 : : {
530 : 0 : return 10.0 * std::log10 (std::norm (z));
531 : : }
532 : :
533 : :
534 : : /*!\brief Compute limited complex exponential
535 : : \param[in] z complex number
536 : : \return limited exponential of z
537 : : \todo Change limexp(real) limexp(complex) file order
538 : : */
539 : 0 : nr_complex_t limexp (const nr_complex_t z)
540 : : {
541 : 0 : nr_double_t mag = qucs::limexp (real (z));
542 : 0 : return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
543 : : }
544 : :
545 : :
546 : : /*!\brief Construct a complex number using polar notation
547 : : \param[in] mag Magnitude
548 : : \param[in] ang Angle
549 : : \return complex number in rectangular form
550 : : */
551 : 8880 : 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 : 8880 : return nr_complex_t (mag * cos (ang), mag * sin (ang));
557 : : #endif
558 : : }
559 : :
560 : : /*!\brief Extension of polar construction to complex
561 : : \param[in] a Magnitude
562 : : \param[in] p Angle
563 : : \return complex number in rectangular form
564 : : \bug Do not seems holomorph form of real polar
565 : : */
566 : 0 : nr_complex_t polar (const nr_complex_t a, const nr_complex_t p)
567 : : {
568 : : #ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
569 : : return std::polar (a, p);
570 : : #else
571 [ # # ][ # # ]: 0 : return a * exp (nr_complex_t (imag (p),-real (p)));
572 : : #endif
573 : : }
574 : :
575 : :
576 : : /*!\brief Converts impedance to reflexion coefficient
577 : : \param[in] z impedance
578 : : \param[in] zref normalisation impedance
579 : : \return reflexion coefficient
580 : : */
581 : 1500 : nr_complex_t ztor (const nr_complex_t z, nr_complex_t zref) {
582 [ + - ][ + - ]: 1500 : return (z - zref) / (z + zref);
583 : : }
584 : :
585 : : /*!\brief Converts reflexion coefficient to impedance
586 : : \param[in] r reflexion coefficient
587 : : \param[in] zref normalisation impedance
588 : : \return impedance
589 : : */
590 : 761 : nr_complex_t rtoz (const nr_complex_t r, nr_complex_t zref) {
591 [ + - ][ + - ]: 761 : return zref * (1.0 + r) / (1.0 - r);
592 : : }
593 : :
594 : : /*!\brief Converts admittance to reflexion coefficient
595 : : \param[in] y admitance
596 : : \param[in] zref normalisation impedance
597 : : \return reflexion coefficient
598 : : */
599 : 0 : nr_complex_t ytor (const nr_complex_t y, nr_complex_t zref) {
600 [ # # ][ # # ]: 0 : return (1.0 - y * zref) / (1.0 + y * zref);
601 : : }
602 : :
603 : : /*!\brief Converts reflexion coefficient to admittance
604 : : \param[in] r reflexion coefficient
605 : : \param[in] zref normalisation impedance
606 : : \return admittance
607 : : */
608 : 0 : nr_complex_t rtoy (const nr_complex_t r, nr_complex_t zref) {
609 [ # # ][ # # ]: 0 : return (1.0 - r) / (1.0 + r) / zref;
610 : : }
611 : :
612 : :
613 : :
614 : :
615 : :
616 : : /*!\brief Complex floor
617 : :
618 : : floor is the largest integral value not greater than argument
619 : : Apply floor to real and imaginary part
620 : : \param[in] z complex number
621 : : \return floored complex number
622 : : */
623 : 0 : nr_complex_t floor (const nr_complex_t z) {
624 : 0 : return nr_complex_t (std::floor (real (z)), std::floor (imag (z)));
625 : : }
626 : :
627 : :
628 : : /*!\brief Complex ceil
629 : : Ceil is the smallest integral value not less than argument
630 : : Apply ceil to real and imaginary part
631 : : \param[in] z complex number
632 : : \return ceilled complex number
633 : : */
634 : 0 : nr_complex_t ceil (const nr_complex_t z) {
635 : 0 : return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z)));
636 : : }
637 : :
638 : : /*!\brief Complex fix
639 : :
640 : : Apply fix to real and imaginary part
641 : : \param[in] z complex number
642 : : \return fixed complex number
643 : : \todo why not using real fix
644 : : */
645 : 0 : nr_complex_t fix (const nr_complex_t z) {
646 : 0 : nr_double_t x = real (z);
647 : 0 : nr_double_t y = imag (z);
648 [ # # ]: 0 : x = (x > 0) ? std::floor (x) : std::ceil (x);
649 [ # # ]: 0 : y = (y > 0) ? std::floor (y) : std::ceil (y);
650 : 0 : return nr_complex_t (x, y);
651 : : }
652 : :
653 : :
654 : :
655 : : /*!\brief Complex fmod
656 : : Apply fmod to the complex z
657 : : \param[in] x complex number (numerator)
658 : : \param[in] y complex number (denominator)
659 : : \return return \f$x - n * y\f$ where n is the quotient of \f$x / y\f$,
660 : : rounded towards zero to an integer.
661 : : */
662 : 0 : nr_complex_t fmod (const nr_complex_t x, const nr_complex_t y) {
663 [ # # ]: 0 : nr_complex_t n = qucs::floor (x / y);
664 [ # # ][ # # ]: 0 : return x - n * y;
665 : : }
666 : :
667 : :
668 : : /*!\brief Square of complex number
669 : :
670 : : \param[in] z complex number
671 : : \return squared complex number
672 : : */
673 : 0 : nr_complex_t sqr (const nr_complex_t z) {
674 : 0 : nr_double_t r = real (z);
675 : 0 : nr_double_t i = imag (z);
676 : 0 : return nr_complex_t (r * r - i * i, 2 * r * i);
677 : : }
678 : :
679 : :
680 : :
681 : :
682 : :
683 : : /*!\brief Heaviside step function for complex number
684 : :
685 : : Apply Heaviside to real and imaginary part
686 : : \param[in] z Heaviside argument
687 : : \return Heaviside step
688 : : \todo Create Heaviside alias
689 : : \todo Why not using real heaviside
690 : : */
691 : 0 : nr_complex_t step (const nr_complex_t z)
692 : : {
693 : 0 : nr_double_t x = real (z);
694 : 0 : nr_double_t y = imag (z);
695 [ # # ]: 0 : if (x < 0.0)
696 : 0 : x = 0.0;
697 [ # # ]: 0 : else if (x > 0.0)
698 : 0 : x = 1.0;
699 : : else
700 : 0 : x = 0.5;
701 [ # # ]: 0 : if (y < 0.0)
702 : 0 : y = 0.0;
703 [ # # ]: 0 : else if (y > 0.0)
704 : 0 : y = 1.0;
705 : : else
706 : 0 : y = 0.5;
707 : 0 : 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 : :
722 : : /*!\brief Bessel function of first kind
723 : :
724 : : \param[in] n order
725 : : \param[in] z argument
726 : : \return Bessel function of first kind of order n
727 : : \bug Not implemented
728 : : */
729 : 0 : nr_complex_t jn (const int n, const nr_complex_t z)
730 : : {
731 : 0 : return cbesselj (n, z);
732 : : }
733 : :
734 : :
735 : : /*!\brief Bessel function of second kind
736 : :
737 : : \param[in] n order
738 : : \param[in] z argument
739 : : \return Bessel function of second kind of order n
740 : : \bug Not implemented
741 : : */
742 : 0 : nr_complex_t yn (const int n, const nr_complex_t z)
743 : : {
744 : 0 : return nr_complex_t (::yn (n, std::real (z)), 0);
745 : : }
746 : :
747 : :
748 : : /*!\brief Modified Bessel function of first kind
749 : :
750 : : \param[in] z argument
751 : : \return Modified Bessel function of first kind of order 0
752 : : \bug Not implemented
753 : : */
754 : 0 : nr_complex_t i0 (const nr_complex_t z)
755 : : {
756 : 0 : return nr_complex_t (fspecial::i0 (std::real (z)), 0);
757 : : }
758 : :
759 : :
760 : : /*!\brief Error function
761 : :
762 : : \param[in] z argument
763 : : \return Error function
764 : : \bug Not implemented
765 : : */
766 : 0 : nr_complex_t erf (const nr_complex_t z)
767 : : {
768 : : #ifdef HAVE_STD_ERF
769 : 0 : 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 : 0 : return nr_complex_t (zerf, 0);
776 : : }
777 : :
778 : : /*!\brief Complementart error function
779 : :
780 : : \param[in] z argument
781 : : \return Complementary error function
782 : : \bug Not implemented
783 : : */
784 : 0 : nr_complex_t erfc (const nr_complex_t z)
785 : : {
786 : : #ifdef HAVE_STD_ERF
787 : 0 : 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 : 0 : return nr_complex_t (zerfc, 0);
794 : : }
795 : :
796 : : /*!\brief Inverse of error function
797 : :
798 : : \param[in] z argument
799 : : \return Inverse of error function
800 : : \bug Not implemented
801 : : */
802 : 0 : nr_complex_t erfinv (const nr_complex_t z)
803 : : {
804 : 0 : return nr_complex_t (fspecial::erfinv (std::real (z)), 0);
805 : : }
806 : :
807 : : /*!\brief Inverse of complementart error function
808 : :
809 : : \param[in] z argument
810 : : \return Inverse of complementary error function
811 : : \bug Not implemented
812 : : */
813 : 0 : nr_complex_t erfcinv (const nr_complex_t z)
814 : : {
815 : 0 : return nr_complex_t (fspecial::erfcinv (std::real (z)), 0);
816 : : }
817 : :
818 : :
819 : :
820 : : // ========================
821 : :
822 : :
823 : : /*!\brief Modulo
824 : : \todo Why not inline
825 : : */
826 : 0 : nr_complex_t operator%(const nr_complex_t z1, const nr_complex_t z2)
827 : : {
828 [ # # ][ # # ]: 0 : return z1 - z2 * floor (z1 / z2);
829 : : }
830 : :
831 : : /*!\brief Modulo
832 : : \todo Why not inline
833 : : */
834 : 0 : nr_complex_t operator%(const nr_complex_t z1, const nr_double_t r2)
835 : : {
836 [ # # ]: 0 : return z1 - r2 * floor (z1 / r2);
837 : : }
838 : :
839 : : /*!\brief Modulo
840 : : \todo Why not inline
841 : : */
842 : 0 : nr_complex_t operator%(const nr_double_t r1, const nr_complex_t z2)
843 : : {
844 [ # # ]: 0 : return r1 - z2 * floor (r1 / z2);
845 : : }
846 : :
847 : :
848 : : /*!\brief Equality of two complex
849 : : \todo Why not inline
850 : : \note Like equality of double this test
851 : : is meaningless in finite precision
852 : : Use instead fabs(x-x0) < tol
853 : : */
854 : 133746 : bool operator==(const nr_complex_t z1, const nr_complex_t z2)
855 : : {
856 [ - + ][ # # ]: 133746 : return (std::real (z1) == std::real (z2)) && (std::imag (z1) == std::imag (z2));
857 : : }
858 : :
859 : : /*!\brief Inequality of two complex
860 : : \todo Why not inline
861 : : \note Like inequality of double this test
862 : : is meaningless in finite precision
863 : : Use instead fabs(x-x0) > tol
864 : : */
865 : 18 : bool operator!=(const nr_complex_t z1, const nr_complex_t z2)
866 : : {
867 [ - + ][ # # ]: 18 : return (std::real (z1) != std::real (z2)) || (std::imag (z1) != std::imag (z2));
868 : : }
869 : :
870 : : /*!\brief Superior of equal
871 : : \todo Why not inline
872 : : */
873 : 0 : bool operator>=(const nr_complex_t z1, const nr_complex_t z2)
874 : : {
875 : 0 : return norm (z1) >= norm (z2);
876 : : }
877 : :
878 : : /*!\brief Inferior of equal
879 : : \todo Why not inline
880 : : */
881 : 0 : bool operator<=(const nr_complex_t z1, const nr_complex_t z2)
882 : : {
883 : 0 : return norm (z1) <= norm (z2);
884 : : }
885 : :
886 : : /*!\brief Superior
887 : : \todo Why not inline
888 : : */
889 : 0 : bool operator>(const nr_complex_t z1, const nr_complex_t z2)
890 : : {
891 : 0 : return norm (z1) > norm (z2);
892 : : }
893 : :
894 : : /*!\brief Inferior
895 : : \todo Why not inline
896 : : */
897 : 0 : bool operator<(const nr_complex_t z1, const nr_complex_t z2)
898 : : {
899 : 0 : return norm (z1) < norm (z2);
900 : : }
901 : :
902 : : } // namespace qucs
903 : :
|
