cmplx.cpp

Go to the documentation of this file.

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

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

#include <math.h>

#include "real.h"
#include "cmplx.h"

cmplx::cmplx (nr_double_t real, nr_double_t imag) {
  r = real;
  i = imag;
}

cmplx::cmplx () {
  r = 0.0;
  i = 0.0;
}

cmplx::cmplx (const cmplx& z) {
  r = z.r;
  i = z.i;
}

nr_double_t cmplx::abs (void) {
  return xhypot (r, i);
}

nr_double_t abs (const cmplx z) {
  return xhypot (z.r, z.i);
}

nr_double_t cmplx::norm (void) {
  return r * r + i * i;
}

nr_double_t cmplx::arg (void) {
  return atan2 (i, r);
}

nr_double_t arg (const cmplx z) {
  return atan2 (z.i, z.r);
}

nr_double_t cmplx::real (void) {
  return r;
}

nr_double_t real (const cmplx z) {
  return z.r;
}

nr_double_t cmplx::imag (void) {
  return i;
}

nr_double_t imag (const cmplx z) {
  return z.i;
}

cmplx cmplx::conj (void) {
  return cmplx (r, -i);
}

cmplx conj (const cmplx z) {
  return cmplx (z.r, -z.i);
}

cmplx cmplx::operator+() {
  return cmplx (r, i);
}

cmplx cmplx::operator-() {
  return cmplx (-r, -i);
}

cmplx& cmplx::operator+=(const cmplx z2) {
  r += z2.r;
  i += z2.i;
  return *this;
}

cmplx& cmplx::operator+=(const nr_double_t r2) {
  r += r2;
  return *this;
}

cmplx& cmplx::operator-=(const cmplx z2) {
  r -= z2.r;
  i -= z2.i;
  return *this;
}

cmplx& cmplx::operator-=(const nr_double_t r2) {
  r -= r2;
  return *this;
}

cmplx& cmplx::operator*=(const nr_double_t r2) {
  r *= r2;
  i *= r2;
  return *this;
}

cmplx& cmplx::operator/=(const nr_double_t r2) {
  r /= r2;
  i /= r2;
  return *this;
}

cmplx operator+(const cmplx z1, const cmplx z2) {
  return cmplx (z1.r + z2.r, z1.i + z2.i);
}

cmplx operator+(const nr_double_t r1, const cmplx z2) {
  return cmplx (r1 + z2.r, z2.i);
}

cmplx operator+(const cmplx z1, const nr_double_t r2) {
  return cmplx (z1.r + r2, z1.i);
}

cmplx operator-(const cmplx z1, const cmplx z2) {
  return cmplx (z1.r - z2.r, z1.i - z2.i);
}

cmplx operator-(const nr_double_t r1, const cmplx z2) {
  return cmplx (r1 - z2.r, -z2.i);
}

cmplx operator-(const cmplx z1, const nr_double_t r2) {
  return cmplx (z1.r - r2, z1.i);
}

cmplx operator*(const cmplx z1, const nr_double_t r2) {
  return cmplx (z1.r * r2, z1.i * r2);
}

cmplx operator*(const nr_double_t r1, const cmplx z2) {
  return cmplx (z2.r * r1, z2.i * r1);
}

cmplx operator*(const cmplx z1, const cmplx z2) {
  return cmplx (z1.r * z2.r - z1.i * z2.i, z1.i * z2.r + z1.r * z2.i);
}

cmplx operator/(const cmplx z1, const nr_double_t r2) {
  return cmplx (z1.r / r2, z1.i / r2);
}

cmplx operator/(const cmplx z1, const cmplx z2) {
#if 0
  nr_double_t n = norm (z2);
  return cmplx ((z1.r * z2.r + z1.i * z2.i) / n,
        (z1.i * z2.r - z1.r * z2.i) / n);
#else /* avoid numerical overflow and underrun */
  nr_double_t r, i, n, d;
  if (fabs (z2.r) > fabs (z2.i)) {
    n = z2.i / z2.r;
    d = z2.r + z2.i * n;
    r = (z1.r + z1.i * n) / d;
    i = (z1.i - z1.r * n) / d;
  }
  else {
    n = z2.r / z2.i;
    d = z2.r * n + z2.i;
    r = (z1.r * n + z1.i) / d;
    i = (z1.i * n - z1.r) / d;
  }
  return cmplx (r, i);
#endif
}

cmplx& cmplx::operator/=(const cmplx z) {
#if 0
  nr_double_t n1, n2;
  n1 = norm (z);
  n2 = (r * z.r + i * z.i) / n1;
  i  = (i * z.r - r * z.i) / n1;
  r  = n2;
#else /* avoid numerical overflow and underrun */
  nr_double_t n, d, t;
  if (fabs (z.r) > fabs (z.i)) {
    n = z.i / z.r;
    d = z.r + z.i * n;
    t = (r + i * n) / d;
    i = (i - r * n) / d;
    r = t;
  }
  else {
    n = z.r / z.i;
    d = z.r * n + z.i;
    t = (r * n + i) / d;
    i = (i * n - r) / d;
    r = t;
  }
#endif
  return *this;
}

cmplx operator/(const nr_double_t r1, const cmplx z2) {
#if 0
  nr_double_t n = norm (z2);
  return cmplx (r1 * z2.r / n, -r1 * z2.i / n);
#else /* avoid numerical overflow and underrun */
  nr_double_t r, i, n, d;
  if (fabs (z2.r) > fabs (z2.i)) {
    n = z2.i / z2.r;
    d = z2.r + z2.i * n;
    r = r1 / d;
    i = -n * r1 / d;
  }
  else {
    n = z2.r / z2.i;
    d = z2.r * n + z2.i;
    r = r1 * n / d;
    i = -r1 / d;
  }
  return cmplx (r, i);
#endif
}

cmplx& cmplx::operator=(const cmplx z) {
  r = z.r;
  i = z.i;
  return *this;
}

cmplx& cmplx::operator=(const nr_double_t x) {
  r = x;
  i = 0.0;
  return *this;
}

cmplx& cmplx::operator*=(const cmplx z) {
  nr_double_t n;
  n = r * z.r - i * z.i;
  i = i * z.r + r * z.i;
  r = n;
  return *this;
}

#ifdef DEBUG
#include <stdio.h>
void cmplx::print (void) {
  fprintf (stderr, "%+.2e,%+.2e ", (double) r, (double) i);
}
#endif /* DEBUG */