Complex.h

Go to the documentation of this file.

#ifndef AMPGEN_COMPLEX_H 
#define AMPGEN_COMPLEX_H 1 
 
#include <complex>
 
namespace AmpGen {
  template <typename real_t>
  struct Complex {
    real_t re;
    real_t im;
    using rt = real_t; 
    Complex() = default;
    Complex( const real_t& re, const real_t& im) : re(re), im(im) {}
    // Complex( const float&   re, const float& im) : re(re), im(im) {}
    Complex( const std::complex<double>& f ) : re( f.real() ), im( f.imag() ) {}
    Complex( const std::complex<float>& f  ) : re( f.real() ), im( f.imag() ) {}
    explicit Complex( const real_t& arg ) : re(arg) {};
    Complex(const real_t& re, const double& im  ) : re(re), im(im) {} 
    // explicit Complex( const double& arg ) : re(arg) {};
    inline Complex operator+=(const Complex& rhs ); 
    inline Complex operator-=(const Complex& rhs ); 
    inline Complex operator*=(const Complex& rhs ); 
    inline Complex operator/=(const Complex& rhs ); 
    real_t real() const { return re; }
    real_t imag() const { return im; }
    real_t norm() const { return re*re + im *im ; }
  };
  namespace detail {
    
    template <typename complex_t, typename real_t> complex_t make_complex( const real_t& re, const real_t& im ){
      return complex_t(re, im ); 
    }
    template <typename complex_t> complex_t make_complex( const typename complex_t::rt& re){
      using real_t = typename complex_t::rt;
      return complex_t(re, real_t(0.) ); 
    }
    template <typename complex_t> complex_t make_complex( const complex_t& cmplx ){
      return cmplx; 
    }
  }
 
  template<typename real_t> inline real_t            real(const Complex<real_t>& arg ){ return arg.re ; }
  template<typename real_t> inline real_t            real(const real_t& arg ){ return arg; }
  template<typename real_t> inline real_t            imag(const Complex<real_t>& arg ){ return arg.im ; }
  template<typename real_t> inline real_t             abs(const Complex<real_t>& v ) { return sqrt( v.re * v.re + v.im * v.im ) ; }
  template<typename real_t> inline real_t            norm(const Complex<real_t>& v ) { return  ( v.re * v.re + v.im * v.im ) ; }
  template<typename real_t> inline Complex<real_t>   conj(const Complex<real_t>& arg ){ return Complex<real_t>(arg.re, -arg.im) ; }
  
 
  template<typename real_t, typename R2_t> inline Complex<real_t> operator+( const Complex<real_t>& lhs, const R2_t& rhs ) { return Complex<real_t>(lhs.re + real_t(rhs), lhs.im); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator-( const Complex<real_t>& lhs, const R2_t& rhs ) { return Complex<real_t>(lhs.re - real_t(rhs), lhs.im); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator*( const Complex<real_t>& lhs, const R2_t& rhs ) { return Complex<real_t>(lhs.re*real_t(rhs), lhs.im*real_t(rhs)); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator/( const Complex<real_t>& lhs, const R2_t& rhs ) { return Complex<real_t>(lhs.re/real_t(rhs), lhs.im/real_t(rhs)); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator+( const R2_t& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(real_t(lhs) + rhs.re,  rhs.im); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator-( const R2_t& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(real_t(lhs) - rhs.re, - rhs.im); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator*( const R2_t& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(real_t(lhs)*rhs.re, lhs*rhs.im); }
  template<typename real_t, typename R2_t> inline Complex<real_t> operator/( const R2_t& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>( real_t(lhs) * rhs.re , -real_t(lhs) *rhs.im) / (rhs.re * rhs.re + rhs.im * rhs.im ); }
  template<typename real_t> inline Complex<real_t> operator+( const Complex<real_t>& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(lhs.re + rhs.re, lhs.im + rhs.im); }
  template<typename real_t> inline Complex<real_t> operator-( const Complex<real_t>& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(lhs.re - rhs.re, lhs.im - rhs.im); }
  template<typename real_t> inline Complex<real_t> operator*( const Complex<real_t>& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(lhs.re*rhs.re - lhs.im*rhs.im, lhs.re*rhs.im  + lhs.im*rhs.re); }
  template<typename real_t> inline Complex<real_t> operator/( const Complex<real_t>& lhs, const Complex<real_t>& rhs ) { return Complex<real_t>(lhs.re*rhs.re + lhs.im*rhs.im, -lhs.re*rhs.im  + lhs.im*rhs.re) / (rhs.re * rhs.re + rhs.im * rhs.im ); }
  template<typename real_t> inline Complex<real_t> operator-( const Complex<real_t>& x ) { return -1.f * x; }
  template<typename real_t> inline Complex<real_t> Complex<real_t>::operator+=(const Complex<real_t>& rhs ){ *this = *this + rhs; return *this; }
  template<typename real_t> inline Complex<real_t> Complex<real_t>::operator-=(const Complex<real_t>& rhs ){ *this = *this - rhs; return *this; }
  template<typename real_t> inline Complex<real_t> Complex<real_t>::operator*=(const Complex<real_t>& rhs ){ *this = *this * rhs; return *this; }
  template<typename real_t> inline Complex<real_t> Complex<real_t>::operator/=(const Complex<real_t>& rhs ){ *this = *this / rhs; return *this; }
  template<typename real_t> inline Complex<real_t> exp( const Complex<real_t>& v ){
    auto [s,c] = sincos(v.im);
    return exp(v.re) * Complex<real_t>(c, s);
  }
  template <typename real_t> 
  inline Complex<real_t> sqrt( const Complex<real_t>& v )
  {
    auto r = abs(v);
    return Complex ( sqrt( 0.5 * (r + v.re) ), sign(v.im) * sqrt( 0.5*( r - v.re ) ) );
  }
  template <typename real_t>
  inline Complex<real_t> log( const Complex<real_t>& v )
  {
    return Complex<real_t>( 0.5 * log( v.norm() ) , atan2(v.im, v.re) );
  }
  template <typename real_t> 
  inline std::ostream& operator<<( std::ostream& os, const Complex<real_t>& obj ) { return os << "( "<< obj.re << ") (" << obj.im << ")"; }
}
 
#endif