gslpp_complex.cpp

Go to the documentation of this file.

/* 
 * Copyright (C) 2012 HEPfit Collaboration
 * All rights reserved.
 *
 * For the licensing terms see doc/COPYING.
 */

#ifndef GSLPP_COMPLEX_H
#include "gslpp_complex.h"
#include <gsl/gsl_sf.h>
#endif
#include <math.h>

#define GSLEPS 1.e-15

namespace gslpp
{
  complex::complex()
  {
    assign(0., 0., false);
  }

  complex::complex(const double& real, const double& imag, bool polar)
  {
    assign(real, imag, polar);
  }

  complex::complex(const complex& z)
  {
    assign(z.real(), z.imag(), false);
  }
  complex::complex(const double& a)
  {
    assign(a, 0., false);
  }

  complex::~complex()
  {
  }

  bool complex::is_real() const
  {
      return ::fabs(imag()/real()) < GSLEPS;
  }


  bool complex::is_imag() const
  {
      return ::fabs(real()/imag()) < GSLEPS;
  }


  const double& complex::real() const
  {
      return GSL_REAL(_complex);
  }


  const double& complex::imag() const
  {
      return GSL_IMAG(_complex);
  }


  double& complex::real()
  {
    return GSL_REAL(_complex);
  }


  double& complex::imag()
  {
    return GSL_IMAG(_complex);
  }

  double complex::arg() const
  {
      return gsl_complex_arg(_complex);
  }

  double complex::abs() const
  {
      return gsl_complex_abs(_complex);
  }

  double complex::abs2() const
  {
      return gsl_complex_abs2(_complex);
  }

  double complex::log_of_abs() const
  {
      return gsl_complex_logabs(_complex);
  }

  complex& complex::operator=(const complex& z)
  {
    GSL_SET_COMPLEX(&_complex, z.real(), z.imag());
    return *this;
  }

  complex& complex::operator=(const double& x)
  {
    GSL_SET_COMPLEX(&_complex, x, 0);
    return *this;
  }

  bool complex::operator==(const complex& z1) const
  {
      return GSL_COMPLEX_EQ(_complex, z1._complex);
  }

  bool complex::operator!=(const complex& z1) const
  {
      return !(GSL_COMPLEX_EQ(_complex, z1._complex));
  }

  gsl_complex* complex::as_gsl_type_ptr() const
  {
      return const_cast<gsl_complex*>(&_complex);
  }

  gsl_complex& complex::as_gsl_type()
  {
    return _complex;
  }

  const gsl_complex& complex::as_gsl_type() const
  {
      return _complex;
  }

  complex::operator gsl_complex& ()
  {
    return _complex;
  }

  complex::operator const gsl_complex& () const
  {
      return _complex;
  }

  std::ostream& operator<<(std::ostream& output, const complex& z)
  {
    double absim = ::fabs(z.imag());
    output << "(" <<  z.real() << (z.imag() < 0.? "-" : "+");
    if (absim != 1.)
      output << fabs(z.imag()) << "*";
    output << "i)";
    return output;
  }

  // const complex* complex::_i = 0;
  const complex& complex::i()
  {
    static complex _i(0.,1.);
    return _i;
  }

  complex::complex(const gsl_complex* z)
  {
//    _complex = (gsl_complex*)malloc(sizeof(gsl_complex));
    GSL_SET_COMPLEX(&_complex, GSL_REAL(*z), GSL_IMAG(*z));
  }

  complex::complex(const gsl_complex& z)
  {
//    _complex = (gsl_complex*)malloc(sizeof(gsl_complex));
    GSL_SET_COMPLEX(&_complex, GSL_REAL(z), GSL_IMAG(z));
  }

  void
  complex::assign(const double& real=0., const double& imag=0.,
                  bool polar=false)
  {
    if (polar)
      _complex = gsl_complex_polar(real,imag);
    else
      {
        GSL_SET_COMPLEX(&_complex, real, imag);
        // *_complex = gsl_complex_rect(real,imag);
      }
  }

  complex complex::operator-() const
    {
      gsl_complex t = gsl_complex_negative(_complex);
      return complex(&t);
    }

  complex complex::operator+(const complex& z1) const
  {
    gsl_complex rl = gsl_complex_add(_complex, z1._complex);
    return complex(rl);
  }

  complex complex::operator-(const complex& z1) const
  {
    gsl_complex rl = gsl_complex_sub(_complex, z1._complex);
    return complex(rl);
  }

  complex complex::operator*(const complex& z1) const
  {
    gsl_complex rl = gsl_complex_mul(_complex, z1._complex);
    return complex(rl);
  }

  complex complex::operator/(const complex& z1) const
  {
    gsl_complex rl = gsl_complex_div(_complex, z1._complex);
    return complex(rl);
  }

  complex& complex::operator+=(const complex& z1)
  {
    _complex = gsl_complex_add(_complex, z1._complex);
    return *this;
  }

  complex& complex::operator-=(const complex& z1)
  {
    _complex = gsl_complex_sub(_complex, z1._complex);
    return *this;
  }

  complex& complex::operator*=(const complex& z1)
  {
    _complex = gsl_complex_mul(_complex, z1._complex);
    return *this;
  }

  complex& complex::operator/=(const complex& z1)
  {
    _complex = gsl_complex_div(_complex, z1._complex);
    return *this;
  }

  complex complex::operator+(const double& a) const
  {
    gsl_complex rl = gsl_complex_add_real(_complex,a);
    return complex(rl);
  }

  complex complex::operator-(const double& a) const
  {
    gsl_complex rl = gsl_complex_sub_real(_complex,a);
    return complex(rl);
  }

  complex complex::operator*(const double& a) const
  {
    gsl_complex rl = gsl_complex_mul_real(_complex,a);
    return complex(rl);
  }

  complex complex::operator/(const double& a) const
  {
    gsl_complex rl = gsl_complex_div_real(_complex,a);
    return complex(rl);
  }

  complex& complex::operator+=(const double& a)
  {
    _complex = gsl_complex_add_real(_complex,a);
    return *this;
  }

  complex& complex::operator-=(const double& a)
  {
    _complex = gsl_complex_sub_real(_complex,a);
    return *this;
  }

  complex&
  complex::operator*=(const double& a)
  {
    _complex = gsl_complex_mul_real(_complex,a);
    return *this;
  }

  complex& complex::operator/=(const double& a)
  {
    _complex = gsl_complex_div_real(_complex,a);
    return *this;
  }

  complex complex::conjugate() const
  {
      gsl_complex t = gsl_complex_conjugate(_complex);
      return complex(t);
  }

  complex complex::inverse() const
  {
      gsl_complex t = gsl_complex_inverse(_complex);
      return complex(t);
  }

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

  complex operator+(const double& x1, const complex& z2)
  {
    complex z1(x1, 0.);
    return z1 + z2;
  }

  complex operator-(const double& x1, const complex& z2)
  {
    complex z1(x1, 0.);
    return z1 - z2;
  }

  complex operator*(const double& x1, const complex& z2)
  {
    complex z1(x1, 0.);
    return z1 * z2;
  }

  complex operator/(const double& x1, const complex& z2)
  {
    complex z1(x1, 0);
    return z1 / z2;
  }

  //==================================================================
  complex exp(const complex& z)
  {
    return complex(gsl_complex_exp(z.as_gsl_type()));
  }

  complex log(const complex& z)
  {
    return complex(gsl_complex_log(z.as_gsl_type()));
  }

  complex log10(const complex& z)
  {
    return complex(gsl_complex_log10(z.as_gsl_type()));
  }

  complex log(const complex& z,
              const complex& b)
  {
    return complex(gsl_complex_log_b(z.as_gsl_type(),b.as_gsl_type()));
  }
  complex dilog(const complex& z)
  {
    gsl_sf_result re, im;
    gsl_sf_complex_dilog_xy_e(z.real(), z.imag(), &re, &im);
    return complex(re.val, im.val, false);
  }
  //==================================================================
  complex sqrt(const complex& z)
  {
    return complex(gsl_complex_sqrt(z.as_gsl_type()));
  }

  complex pow(const complex& z1,
              const complex& z2)
  {
    return complex(gsl_complex_pow(z1.as_gsl_type(),
                                   z2.as_gsl_type()));
  }

  complex pow(const complex& z, const double  x)
  {
    return complex(gsl_complex_pow_real(z.as_gsl_type(), x));
  }
  //==================================================================
  complex sin(const complex& z)
  {
    return complex(gsl_complex_sin(z.as_gsl_type()));
  }

  complex cos(const complex& z)
  {
    return complex(gsl_complex_cos(z.as_gsl_type()));
  }

  complex tan(const complex& z)
  {
    return complex(gsl_complex_tan(z.as_gsl_type()));
  }

  complex sec(const complex& z)
  {
    return complex(gsl_complex_sec(z.as_gsl_type()));
  }

  complex csc(const complex& z)
  {
    return complex(gsl_complex_csc(z.as_gsl_type()));
  }

  complex cot(const complex& z)
  {
    return complex(gsl_complex_cot(z.as_gsl_type()));
  }

  complex arcsin(const complex& z)
  {
    return complex(gsl_complex_arcsin(z.as_gsl_type()));
  }

  complex arccos(const complex& z)
  {
    return complex(gsl_complex_arccos(z.as_gsl_type()));
  }

  complex arctan(const complex& z)
  {
    return complex(gsl_complex_arctan(z.as_gsl_type()));
  }

  complex arcsec(const complex& z)
  {
    return complex(gsl_complex_arcsec(z.as_gsl_type()));
  }

  complex arccsc(const complex& z)
  {
    return complex(gsl_complex_arccsc(z.as_gsl_type()));
  }

  complex arccot(const complex& z)
  {
    return complex(gsl_complex_arccot(z.as_gsl_type()));
  }

  //==================================================================
  complex sinh(const complex& z)
  {
    return complex(gsl_complex_sinh(z.as_gsl_type()));
  }

  complex cosh(const complex& z)
  {
    return complex(gsl_complex_cosh(z.as_gsl_type()));
  }

  complex tanh(const complex& z)
  {
    return complex(gsl_complex_tanh(z.as_gsl_type()));
  }

  complex sech(const complex& z)
  {
    return complex(gsl_complex_sech(z.as_gsl_type()));
  }

  complex csch(const complex& z)
  {
    return complex(gsl_complex_csch(z.as_gsl_type()));
  }

  complex coth(const complex& z)
  {
    return complex(gsl_complex_coth(z.as_gsl_type()));
  }

  complex arcsinh(const complex& z)
  {
    return complex(gsl_complex_arcsinh(z.as_gsl_type()));
  }

  complex arccosh(const complex& z)
  {
    return complex(gsl_complex_arccosh(z.as_gsl_type()));
  }

  complex arctanh(const complex& z)
  {
    return complex(gsl_complex_arctanh(z.as_gsl_type()));
  }

  complex arcsech(const complex& z)
  {
    return complex(gsl_complex_arcsech(z.as_gsl_type()));
  }

  complex arccsch(const complex& z)
  {
    return complex(gsl_complex_arccsch(z.as_gsl_type()));
  }

  complex arccoth(const complex& z)
  {
    return complex(gsl_complex_arccoth(z.as_gsl_type()));
  }

}
//
// EOF
//