complex.cpp

Go to the documentation of this file.

/*
 * complex.cpp - complex number class implementation
 *
 * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org>
 * Copyright (C) 2006 Gunther Kraut <gn.kraut@t-online.de>
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this package; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
 * Boston, MA 02110-1301, USA.
 *
 * $Id$
 *
 */

#if HAVE_CONFIG_H
# include <config.h>
#endif

#include <cmath>
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>

#include "constants.h"
#include "precision.h"
#include "complex.h"
#include "consts.h"
#include "fspecial.h"


namespace qucs {

//
// trigonometric complex
//

nr_complex_t cos (const nr_complex_t z) {
    return std::cos (z);
}


nr_complex_t sin (const nr_complex_t z) {
    return std::sin (z);
}


nr_complex_t tan (const nr_complex_t z) {
    return std::tan (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
}


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
}

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
}


//
// hyperbolic complex
//

nr_complex_t cosh (const nr_complex_t z) {
    return std::cosh (z);
}


nr_complex_t sinh (const nr_complex_t z) {
    return std::sinh (z);
}


nr_complex_t tanh (const nr_complex_t z) {
    return std::tanh (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
}


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
}


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
}


//
// transcendentals overloads
//

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)));
}

nr_complex_t log (const nr_complex_t z)
{
    nr_double_t phi = arg (z);
    return nr_complex_t (log (abs (z)), phi);
}

nr_complex_t log10 (const nr_complex_t z)
{
    nr_double_t phi = arg (z);
    return nr_complex_t (log10 (abs (z)), phi * M_LOG10E);
}


nr_complex_t pow (const nr_complex_t z, const nr_double_t d) {
    return std::pow (z, d);
}

nr_complex_t pow (const nr_double_t d, const nr_complex_t z) {
    return std::pow (d, z);
}

nr_complex_t pow (const nr_complex_t z1, const nr_complex_t z2) {
    return std::pow (z1, z2);
}


nr_complex_t sqrt (const nr_complex_t z)
{
    return std::sqrt (z);
}


nr_double_t norm (const nr_complex_t z)
{
    return std::norm (z);
}


//
//  Qucs extra math functions
//

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);
}

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);
}

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);
}

nr_complex_t acoth (const nr_complex_t z)
{
    return 0.5 * std::log (2.0 / (z - 1.0) + 1.0);
}


nr_complex_t sech (const nr_complex_t z)
{
    return (1.0 / std::cosh (z));
}

nr_complex_t asech (const nr_complex_t z)
{
    return std::log ((1.0 + std::sqrt (1.0 - z * z)) / z);
}

nr_complex_t     cosech (const nr_complex_t z)
{
    return (1.0 / std::sinh(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
//

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)) * M_LOG2E, phi * M_LOG2E);
#else
    return std::log2 (z);
#endif
}

nr_complex_t signum (const nr_complex_t z)
{
    if (z == 0.0) return 0;
    return z / abs (z);
}

nr_complex_t sign (const nr_complex_t z)
{
    if (z == 0.0) return nr_complex_t (1);
    return z / abs (z);
}


nr_complex_t    sinc (const nr_complex_t z)
{
    if (real(z) == 0.0 && imag(z)) return 1;
    return std::sin (z) / z;
}

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);
}

nr_double_t xhypot (nr_double_t a, nr_complex_t b)
{
    return xhypot (nr_complex_t (a), b);
}

nr_double_t xhypot (nr_complex_t a, nr_double_t b)
{
    return xhypot (a, nr_complex_t (b));
}


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);
}


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);
}


nr_double_t dB (const nr_complex_t z)
{
  return 10.0 * std::log10 (std::norm (z));
}


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)));
}


nr_complex_t polar (const nr_double_t mag, const nr_double_t ang )
{
#ifdef HAVE_CXX_COMPLEX_POLAR_COMPLEX
    return std::polar (mag, ang);
#else
    return nr_complex_t (mag * cos (ang), mag * sin (ang));
#endif
}

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 (imag (p),-real (p)));
#endif
}


nr_complex_t ztor (const nr_complex_t z, nr_complex_t zref) {
    return (z - zref) / (z + zref);
}

nr_complex_t rtoz (const nr_complex_t r, nr_complex_t zref) {
    return zref * (1.0 + r) / (1.0 - r);
}

nr_complex_t ytor (const nr_complex_t y, nr_complex_t zref) {
    return (1.0 - y * zref) / (1.0 + y * zref);
}

nr_complex_t rtoy (const nr_complex_t r, nr_complex_t zref) {
    return (1.0 - r) / (1.0 + r) / zref;
}





nr_complex_t   floor (const nr_complex_t z) {
  return nr_complex_t (std::floor (real (z)), std::floor (imag (z)));
}


nr_complex_t ceil (const nr_complex_t z) {
  return nr_complex_t (std::ceil (real (z)), std::ceil (imag (z)));
}

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);
}



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;
}


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);
}





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;


// === bessel ===

/* FIXME : what about libm jn, yn, isn't that enough? */

nr_complex_t cbesselj (unsigned int, nr_complex_t);

#include "cbesselj.cpp"

nr_complex_t jn (const int n, const nr_complex_t z)
{
    return cbesselj (n, z);
}


nr_complex_t yn (const int n, const nr_complex_t z)
{
    return nr_complex_t (::yn (n, std::real (z)), 0);
}


nr_complex_t i0 (const nr_complex_t z)
{
    return nr_complex_t (fspecial::i0 (std::real (z)), 0);
}


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);
}

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);
}

nr_complex_t erfinv (const nr_complex_t z)
{
    return nr_complex_t (fspecial::erfinv (std::real (z)), 0);
}

nr_complex_t erfcinv (const nr_complex_t z)
{
    return nr_complex_t (fspecial::erfcinv (std::real (z)), 0);
}



// ========================


nr_complex_t operator%(const nr_complex_t z1, const nr_complex_t z2)
{
    return z1 - z2 * floor (z1 / z2);
}

nr_complex_t operator%(const nr_complex_t z1, const nr_double_t r2)
{
    return z1 - r2 * floor (z1 / r2);
}

nr_complex_t operator%(const nr_double_t r1, const nr_complex_t z2)
{
    return r1 - z2 * floor (r1 / z2);
}


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));
}

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));
}

bool operator>=(const nr_complex_t z1, const nr_complex_t z2)
{
    return norm (z1) >= norm (z2);
}

bool operator<=(const nr_complex_t z1, const nr_complex_t z2)
{
    return norm (z1) <= norm (z2);
}

bool operator>(const nr_complex_t z1, const nr_complex_t z2)
{
    return norm (z1) > norm (z2);
}

bool operator<(const nr_complex_t z1, const nr_complex_t z2)
{
    return norm (z1) < norm (z2);
}

} // namespace qucs