---

CTBcomplex.cxx

Go to the documentation of this file.

/*----------------------------------------------------------------------------*/
/* C Tool Box: Designed and implemented by:                                   */
/*    Walter F.J. Mueller   Gesellschaft fuer Schwerionenforschung (GSI)      */
/*                          Postfach 110552, D-64220 Darmstadt, Germany       */
/*                  Email:  W.F.J.Mueller@gsi.de                              */
/*                  WWW:    http://www-kp3.gsi.de/www/kp3/people/mueller.html */
/*------------------------------------------+---------------------------------*/

#include <math.h>

#include "CTB.hxx"
#include "CTBnum.hxx"
#include "CTBtraits.hxx"
#include "CTBcomplex.hxx"

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Sqrt() const
{
  T f_re = 0.;
  T f_im = 0.;

  if (IsReal()) {
    if (mf_re >= 0.) {
      f_re = sqrt(mf_re);
    } else {
      f_im = sqrt(-mf_re);
    }

  } else {
    T f_abs = Abs();
    if (f_abs > 0.) {
      if (mf_re > 0.) {
        f_re = sqrt(0.5 * (f_abs + mf_re));
        f_im = 0.5 * mf_im / f_re;
      } else {
        f_im = sqrt(0.5 * (f_abs - mf_re));
        if (mf_im < 0.) f_im = -f_im;
        f_re = 0.5 * mf_im / f_im;
      }
    }
  }

  return CTBcomplex<T>(f_re, f_im);
}

//------------------------------------------+-----------------------------------

template <class T>
static inline void sincos(T& f_sin, T& f_cos, T f_x)
{
  f_sin = CTBsin(f_x);
  f_cos = CTBcos(f_x);
  return;
}

//------------------------------------------+-----------------------------------

template <class T>
static inline void sinhcosh(T& f_sinh, T& f_cosh, T f_x)
{
  T f_expp = CTBexp(f_x);
  T f_expm = 1./f_expp;
  f_sinh = 0.5 * (f_expp - f_expm);
  f_cosh = 0.5 * (f_expp + f_expm);
  return;
}

//------------------------------------------+-----------------------------------

template <class T>
static inline void sinhcosh_red(T& f_sinh, T& f_cosh, T f_x)
{
  T f_xabs = CTBabs(f_x);
  T f_expm = (f_xabs < 25.) ? CTBexp(-2.*f_xabs) : 0.;
  f_sinh = CTBsign(1. - f_expm, f_x);
  f_cosh = 1. + f_expm;
  return;
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Exp() const
{
  T f_mag = 1.;
  T f_re  = 1.;
  T f_im  = 0.;

  if (mf_re != 0.) f_mag = CTBexp(mf_re);
  if (mf_im != 0.) sincos(f_im, f_re, mf_im);

  return CTBcomplex<T>(f_mag*f_re, f_mag*f_im);
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Log() const
{
  return CTBcomplex<T>(CTBlog(Abs()),Arg());
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Sin() const
{
  T f_sin_re  = 0.;
  T f_cos_re  = 1.;
  T f_sinh_im = 0.;
  T f_cosh_im = 1.;
  
  if (mf_re != 0.) sincos(f_sin_re, f_cos_re, mf_re);
  if (mf_im != 0.) sinhcosh(f_sinh_im,f_cosh_im, mf_im);

  return CTBcomplex<T>(f_sin_re*f_cosh_im, f_cos_re*f_sinh_im);
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Cos() const
{
  T f_sin_re  = 0.;
  T f_cos_re  = 1.;
  T f_sinh_im = 0.;
  T f_cosh_im = 1.;
  
  if (mf_re != 0.) sincos(f_sin_re, f_cos_re, mf_re);
  if (mf_im != 0.) sinhcosh(f_sinh_im, f_cosh_im, mf_im);

  return CTBcomplex<T>(f_cos_re*f_cosh_im, -f_sin_re*f_sinh_im);
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Tan() const
{
  T f_sin_re  = 0.;
  T f_cos_re  = 1.;
  T f_sinh_im = 0.;
  T f_cosh_im = 1.;
  
  if (mf_re != 0.) sincos(f_sin_re, f_cos_re, mf_re);
  if (mf_im != 0.) sinhcosh_red(f_sinh_im, f_cosh_im, mf_im);

  CTBcomplex<T> z_sin(f_sin_re*f_cosh_im,  f_cos_re*f_sinh_im);
  CTBcomplex<T> z_cos(f_cos_re*f_cosh_im, -f_sin_re*f_sinh_im);

  return z_sin/z_cos;
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Asin() const
{
  return CTBtimesI(CTBlog(-CTBtimesI(*this) + CTBsqrt((T)1.-CTBpow2(*this))));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Acos() const
{
  return - CTBtimesI(CTBlog(*this + CTBtimesI(CTBsqrt((T)1.-CTBpow2(*this)))));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Atan() const
{
  CTBcomplex<T> i(0.,1.);
  
  return (T)0.5 * CTBtimesI( CTBlog((i + *this) / (i - *this) ));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Sinh() const
{
  T f_sinh_re = 0.;
  T f_cosh_re = 1.;
  T f_sin_im  = 0.;
  T f_cos_im  = 1.;
  
  if (mf_re != 0.) sinhcosh(f_sinh_re,f_cosh_re, mf_re);
  if (mf_im != 0.) sincos(f_sin_im, f_cos_im, mf_im);

  return CTBcomplex<T>(f_sinh_re*f_cos_im, f_cosh_re*f_sin_im);
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Cosh() const
{
  T f_sinh_re = 0.;
  T f_cosh_re = 1.;
  T f_sin_im  = 0.;
  T f_cos_im  = 1.;
  
  if (mf_re != 0.) sinhcosh(f_sinh_re, f_cosh_re, mf_re);
  if (mf_im != 0.) sincos(f_sin_im, f_cos_im, mf_im);

  return CTBcomplex<T>(f_cosh_re*f_cos_im, f_sinh_re*f_sin_im);
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Tanh() const
{
  T f_sinh_re = 0.;
  T f_cosh_re = 1.;
  T f_sin_im  = 0.;
  T f_cos_im  = 1.;
  
  if (mf_re != 0.) sinhcosh_red(f_sinh_re, f_cosh_re, mf_re);
  if (mf_im != 0.) sincos(f_sin_im, f_cos_im, mf_im);

  CTBcomplex<T> z_sinh(f_sinh_re*f_cos_im, f_cosh_re*f_sin_im);
  CTBcomplex<T> z_cosh(f_cosh_re*f_cos_im, f_sinh_re*f_sin_im);

  return z_sinh/z_cosh;
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Asinh() const
{
  return CTBlog(*this + CTBsqrt(CTBpow2(*this) + (T)1.));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Acosh() const
{
  return CTBlog(*this + CTBsqrt(CTBpow2(*this) - (T)1.));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Atanh() const
{
  return (T)0.5 * CTBlog( ((T)1. + *this) / ((T)1. - *this));
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Pow(int i_y) const
{
  CTBcomplex<T> z_f(1.);
  int           i_nmul = CTBabs(i_y);

  for (int i = 0; i < i_nmul; i++) z_f *= *this;
  if (i_y < 0) z_f = (T)1./z_f;

  return z_f;
}

//------------------------------------------+-----------------------------------

template <class T>
CTBcomplex<T> CTBcomplex<T>::Pow(const CTBcomplex<T>& z_y) const
{
  return CTBexp(z_y * CTBlog(*this));
}

//------------------------------------------+-----------------------------------
template <class T>
  CTBcomplex<T> CTBcomplex<T>::Divide(const CTBcomplex<T>& z_den) const
{
  T f_re;
  T f_im;
  
  if (fabs(z_den.mf_re) >= fabs(z_den.mf_im)) {
    T f_t1 = z_den.mf_im / z_den.mf_re;
    T f_t2 = 1./ (z_den.mf_re + f_t1 * z_den.mf_im);
    f_re   = (mf_re + f_t1 * mf_im) * f_t2;
    f_im   = (mf_im - f_t1 * mf_re) * f_t2;

  } else {
    T f_t1 = z_den.mf_re / z_den.mf_im;
    T f_t2 = 1. / (f_t1 * z_den.mf_re + z_den.mf_im);
    f_re   = (f_t1 * mf_re + mf_im) * f_t2;
    f_im   = (f_t1 * mf_im - mf_re) * f_t2;
  }

  return CTBcomplex<T>(f_re, f_im);
}

//------------------------------------------+-----------------------------------
template <class T>
  CTBcomplex<T> CTBcomplex<T>::Invert(T f_nom) const
{  
  T f_re;
  T f_im;
  
  if (fabs(mf_re) >= fabs(mf_im)) {
    T f_t1 = mf_im / mf_re;
    T f_t2 = 1./ (mf_re + f_t1 * mf_im);
    f_re   =   f_nom * f_t2;
    f_im   = - f_nom * f_t1 * f_t2;

  } else {
    T f_t1 = mf_re / mf_im;
    T f_t2 = 1./ (f_t1 * mf_re + mf_im);
    f_re   =   f_nom * f_t1 * f_t2;
    f_im   = - f_nom * f_t2;
  }

  return CTBcomplex<T>(f_re, f_im);
}

//------------------------------------------+-----------------------------------

template <class T>
void CTBcomplex<T>::ToStream(ostream& os) const
{
  int i_width = os.width();

  os << setw(0) << "(" << setw(i_width) << mf_re 
     << "," << setw(i_width) << mf_im << ")";
  return;
}

//------------------------------------------+-----------------------------------

template <class T>
void CTBcomplex<T>::FromStream(istream& is)
{
  T    f_re;
  T    f_im;
  char c_char;
  
  if (!is) return;                          // quit if stream not good

  is >> c_char;                             // get '('
  if (c_char != '(') { 
    is.setstate(ios::failbit);
    return;
  }

  is >> f_re;                               // read real part
  if (!is) return;

  is >> c_char;                             // get ',' or ')'

  if (c_char == ',') {
    is >> f_im;                             // read imaginary part
    if (!is) return;
    is >> c_char;                           // get trailing ')'
  }

  if (c_char != ')') { 
    is.setstate(ios::failbit);
    return;
  }

  mf_re = f_re;
  mf_im = f_im;

  return;
}

---