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: complex.cpp 1825 2011-03-11 20:42:14Z ela $
 *
 */

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

#include <math.h>

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

using namespace fspecial;

nr_complex_t rect (const nr_double_t x, const nr_double_t y) {
  return nr_complex_t (x, y);
}

#ifndef HAVE_CXX_COMPLEX_NORM

nr_double_t norm (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return r * r + i * i;
}
#endif

#ifndef HAVE_CXX_COMPLEX_POLAR

nr_complex_t polar (const nr_double_t mag, const nr_double_t ang) {
  return rect (mag * cos (ang), mag * sin (ang));
}
#endif

#ifndef HAVE_CXX_COMPLEX_POLAR_COMPLEX

nr_complex_t polar (const nr_complex_t a, const nr_complex_t p) {
  return a * exp (rect (imag (p), -real (p)));
}
#endif

#ifndef HAVE_CXX_COMPLEX_COS 

nr_complex_t cos (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return (polar (exp (-i), r) + polar (exp (i), -r)) / 2.0;
}
#endif

#ifndef HAVE_CXX_COMPLEX_ACOS

nr_complex_t acos (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
#if 0
  return rect (0.0, -2.0) * log (M_SQRT1_2 * (sqrt (z + 1) + sqrt (z - 1)));
#else
  nr_complex_t y = sqrt (z * z - 1.0);
  if (r * i < 0.0) y = -y;
  return rect (0, -1.0) * log (z + y);
#endif
}
#endif

#ifndef HAVE_CXX_COMPLEX_COSH

nr_complex_t cosh (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return (polar (exp (r), i) + polar (exp (-r), -i)) / 2.0;
}
#endif

#ifndef HAVE_CXX_COMPLEX_ACOSH

nr_complex_t acosh (const nr_complex_t z) {
  return log (z + sqrt (z * z - 1.0));
}
#endif

#ifndef HAVE_CXX_COMPLEX_EXP

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

nr_complex_t limexp (const nr_complex_t z) {
  nr_double_t mag = limexp (real (z));
  return nr_complex_t (mag * cos (imag (z)), mag * sin (imag (z)));
}

#ifndef HAVE_CXX_COMPLEX_LOG 

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

#ifndef HAVE_CXX_COMPLEX_LOG10 

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

#ifndef HAVE_CXX_COMPLEX_LOG2

nr_complex_t log2 (const nr_complex_t z) {
  nr_double_t phi = arg (z);
  return nr_complex_t (log (abs (z)) * M_LOG2E, phi * M_LOG2E);
}
#endif

#ifndef HAVE_CXX_COMPLEX_POW

nr_complex_t pow (const nr_complex_t z, const nr_double_t d) {
  return polar (pow (abs (z), d), arg (z) * d);
}

nr_complex_t pow (const nr_double_t d, const nr_complex_t z) {
  return exp (z * log (d));
}

nr_complex_t pow (const nr_complex_t z1, const nr_complex_t z2) {
  return exp (z2 * log (z1));
}
#endif

#ifndef HAVE_CXX_COMPLEX_SIN 

nr_complex_t sin (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return (polar (exp (-i), r - M_PI_2) - polar (exp (i), -r - M_PI_2)) / 2.0;
}
#endif

#ifndef HAVE_CXX_COMPLEX_ASIN

nr_complex_t asin (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return nr_complex_t (0.0, -1.0) * log (rect (-i, r) + sqrt (1.0 - z * z));
}
#endif

#ifndef HAVE_CXX_COMPLEX_SINH

nr_complex_t sinh (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return (polar (exp (r), i) - polar (exp (-r), -i)) / 2.0;
}
#endif

#ifndef HAVE_CXX_COMPLEX_ASINH

nr_complex_t asinh (const nr_complex_t z) {
  return log (z + sqrt (z * z + 1.0));
}
#endif 

#ifndef HAVE_CXX_COMPLEX_SQRT

nr_complex_t sqrt (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
#if 0
  if (i == 0.0) {
    return r < 0.0 ? nr_complex_t (0.0, sqrt (-r)) : nr_complex_t (sqrt (r));
  } else {
    nr_double_t phi = arg (z);
    return polar (sqrt (fabs (z)), phi / 2.0);
  }
#else /* better numerical behaviour, avoiding arg() and polar() */
  if (r == 0.0 && i == 0.0) {
    return nr_complex_t (0.0, 0.0);
  } else {
    nr_double_t x = fabs (r);
    nr_double_t y = fabs (i);
    nr_double_t k;
    if (x >= y)  {
      nr_double_t t = y / x;
      k = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));
    } else {
      nr_double_t t = x / y;
      k = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));
    }
    if (r >= 0.0) {
      return nr_complex_t (k, i / (2.0 * k));
    } else {
      nr_double_t vi = (i >= 0) ? k : -k;
      return nr_complex_t (i / (2.0 * vi), vi);
    }
  }
#endif
}
#endif

#ifndef HAVE_CXX_COMPLEX_TAN

nr_complex_t tan (const nr_complex_t z) {
  nr_double_t r = 2.0 * real (z);
  nr_double_t i = 2.0 * imag (z);
  return rect (0.0, -1.0) + rect (0.0, 2.0) / (polar (exp (-i), r) + 1.0);
}
#endif

#ifndef HAVE_CXX_COMPLEX_ATAN

nr_complex_t atan (const nr_complex_t z) {
  return rect (0.0, -0.5) * log (rect (0, 2) / (z + rect (0, 1)) - 1.0);
}
#endif

#ifndef HAVE_CXX_COMPLEX_ATAN2

nr_complex_t atan2 (const nr_complex_t y, const nr_complex_t x) {
  nr_complex_t a = atan (y / x);
  return real (x) > 0.0 ? a : -a;
}
#endif

#ifndef HAVE_CXX_COMPLEX_TANH

nr_complex_t tanh (const nr_complex_t z) {
  nr_double_t r = 2.0 * real (z);
  nr_double_t i = 2.0 * imag (z);
  return 1.0 - 2.0 / (polar (exp (r), i) + 1.0);
}
#endif

#ifndef HAVE_CXX_COMPLEX_ATANH

nr_complex_t atanh (const nr_complex_t z) {
  return 0.5 * log ( 2.0 / (1.0 - z) - 1.0);
}
#endif

nr_complex_t cot (const nr_complex_t z) {
  nr_double_t r = 2.0 * real (z);
  nr_double_t i = 2.0 * imag (z);
  return rect (0.0, 1.0) + rect (0.0, 2.0) / (polar (exp (-i), r) - 1.0);
}

nr_complex_t acot (const nr_complex_t z) {
  return rect (0.0, -0.5) * log (rect (0, 2) / (z - rect (0, 1)) + 1.0);
}

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

nr_complex_t coth (const nr_complex_t z) {
  nr_double_t r = 2.0 * real (z);
  nr_double_t i = 2.0 * imag (z);
  return 1.0 + 2.0 / (polar (exp (r), i) - 1.0);
}

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

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

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

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 rtoz (const nr_complex_t r, nr_complex_t zref) {
  return zref * (1.0 + r) / (1.0 - r);
}

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

#ifndef HAVE_CXX_COMPLEX_FLOOR 

nr_complex_t floor (const nr_complex_t z) {
  return rect (floor (real (z)), floor (imag (z)));
}
#endif

#ifndef HAVE_CXX_COMPLEX_FMOD

nr_complex_t fmod (const nr_complex_t x, const nr_complex_t y) {
  nr_complex_t n = floor (x / y);
  return x - n * y;
}

nr_complex_t fmod (const nr_complex_t x, const nr_double_t y) {
  nr_complex_t n = floor (x / y);
  return x - n * y;
}

nr_complex_t fmod (const nr_double_t x, const nr_complex_t y) {
  nr_complex_t n = floor (x / y);
  return x - n * y;
}
#endif

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

nr_complex_t sign (const nr_complex_t z) {
  if (z == 0) return nr_complex_t (1);
  return z / abs (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) * sqrt (1 + d / c);
  else if (d == 0)
    return 0;
  else
    return abs (b) * sqrt (1 + c / d);
}

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

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

nr_complex_t sinc (const nr_complex_t z) {
  if (z == 0) return 1;
  return sin (z) / z;
}

nr_complex_t ceil (const nr_complex_t z) {
  return rect (ceil (real (z)), 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) ? floor (x) : ceil (x);
  y = (y > 0) ? floor (y) : ceil (y);
  return rect (x, y);
}

nr_complex_t round (const nr_complex_t z) {
  return rect (round (real (z)), round (imag (z)));
}

nr_complex_t trunc (const nr_complex_t z) {
  return rect (trunc (real (z)), trunc (imag (z)));
}

nr_complex_t sqr (const nr_complex_t z) {
  nr_double_t r = real (z);
  nr_double_t i = imag (z);
  return rect (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 rect (x, y);
}

nr_complex_t cbesselj (unsigned int, nr_complex_t);

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 rect (yn (n, real (z)), 0);
}

nr_complex_t i0 (const nr_complex_t z) {
  return rect (i0 (real (z)), 0);
}

nr_complex_t erf (const nr_complex_t z) {
  return rect (erf (real (z)), 0);
}

nr_complex_t erfc (const nr_complex_t z) {
  return rect (erfc (real (z)), 0);
}

nr_complex_t erfinv (const nr_complex_t z) {
  return rect (erfinv (real (z)), 0);
}

nr_complex_t erfcinv (const nr_complex_t z) {
  return rect (erfcinv (real (z)), 0);
}

bool operator==(const nr_complex_t z1, const nr_complex_t z2) {
  return (real (z1) == real (z2)) && (imag (z1) == imag (z2));
}

bool operator!=(const nr_complex_t z1, const nr_complex_t z2) {
  return (real (z1) != real (z2)) || (imag (z1) != 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);
}

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