#ifndef __COMPLEX_GSL_H
#define __COMPLEX_GSL_H

/** @file
 * @brief 複素数ライブラリ
 *
 * 複素数を定義したライブラリ。
 * 多分誰が作っても似たような感じになると思われる。
 */

#include "param/complex.h"
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

/**
 * @brief 複素数
 *
 * 複素数をあらわすクラスです。
 *
 * @param FloatT 演算精度、doubleなど
 */
template <>
class Complex<double> : public gsl_complex{
  private:
    typedef double FloatT;
    typedef Complex<FloatT> self_t;
    typedef gsl_complex super_t;

  public:
    Complex(const super_t &c) : super_t(c) {}

  	/**
  	 * コンストラクタ。
  	 *
  	 * @param real 実数部
  	 * @param imaginary 虚数部
  	 */
    Complex(const FloatT &real, const FloatT &imaginary){
      GSL_SET_COMPLEX(this, real, imaginary);
    }
    /**
     * コンストラクタ。
     * 虚数部は0で初期化されます。
     *
     * @param real 実数部
     */
    Complex(const FloatT &real){
      GSL_SET_COMPLEX(this, real, 0);
    }
    /**
     * コンストラクタ。
     * 実数部、虚数部ともに0に初期化されます。
     *
     */
    Complex(){
      GSL_SET_COMPLEX(this, 0, 0);
    }
    /**
     * デストラクタ。
     *
     */
    ~Complex(){}

    /**
     * 実数部を返します。
     *
     * @return (FloatT) 実数部
     */
    FloatT &real(){return GSL_REAL(*this);}
    /**
     * 虚数部を返します。
     *
     * @return (FloatT) 虚数部
     */
    FloatT &imaginary(){return GSL_IMAG(*this);}

    /**
     * 絶対値の二乗を返します。
     * pow(real(), 2) + pow(imaginary(), 2)をしています。
     *
     * @return (FloatT) 絶対値の二乗
     * @see pow(FloatT, FloatT)
     */
    FloatT abs2() const {
      return gsl_complex_abs2(*this);
    }
    /**
     * 絶対値を返します。
     * sqrt(abs2())をしています。
     *
     * @return (FloatT) 絶対値
     * @see abs2()
     * @see sqrt()
     */
    FloatT abs() const {
      return gsl_complex_abs(*this);
    }
    /**
     * 偏角を返します。
     *
     * @return (double) 偏角
     */
    double arg() const {
      return gsl_complex_arg(*this);
    }
    /**
     * 指定乗します。
     * ただし定義として、複素数が
     * @f[
     *    a + b i = r e^{i \theta}
     * @f]
     * とあらわされるとき、
     * @f[
     *    (a + b i)^{n} = r^{n} e^{i \theta n}
     * @f]
     * とします。
     *
     * @param factor 乗数
     * @return (self_t) 結果
     */
    self_t power(const FloatT &factor) const {
    	return Complex(gsl_complex_pow_real(*this, factor));
    }

    /**
     * 平方根を求めます。定義はpowerに従います。
     *
     * @return (self_t) 結果
     * @see power
     */
    self_t sqrt() const {
      return Complex(gsl_complex_sqrt(*this));
    }

    /**
     * 共役複素数を返します。
     *
     * @return (self_t) 共役複素数
     */
    self_t conjugate() const {
      return Complex(gsl_complex_conjugate(*this));
    }

    /**
     * 等しいかどうか調べます。
     *
     * @return (bool) 等しい場合true
     */
    bool operator==(const self_t &complex) const{
      return (GSL_REAL(*this) == GSL_REAL(complex))
          ? (GSL_IMAG(*this) == GSL_IMAG(complex))
          : false;
    }
    /**
     * 等しくないか調べます。
     *
     * @return (bool) 等しくない場合true
     */
    bool operator!=(const self_t &complex) const{
      return !(*this == complex);
    }

    /**
     * 加算を行います。破壊的です。
     *
     * @return (self_t) 加算結果
     */
    self_t &operator+=(const FloatT &scalar){
      GSL_REAL(*this) += scalar;
      return *this;
    }
    /**
     * 加算を行います。
     *
     * @return (self_t) 加算結果
     */
    self_t operator+(const FloatT &scalar) const{
      self_t result = *this;
      return (result += scalar);
    }
    /**
     * 加算を行います。
     *
     * @return (self_t) 加算結果
     */
    friend self_t operator+(const FloatT &scalar, const self_t complex){return (complex + scalar);}

    /**
     * 減算を行います。破壊的です。
     *
     * @return (self_t) 減算結果
     */
    self_t &operator-=(const FloatT &scalar){return (*this) += (-scalar);}
    /**
     * 減算を行います。
     *
     * @return (self_t) 減算結果
     */
    self_t operator-(const FloatT &scalar) const{return (*this) + (-scalar);}
    /**
     * 減算を行います。
     *
     * @return (self_t) 減算結果
     */
    friend self_t operator-(const FloatT &scalar, const self_t complex){return (complex - scalar);}

    /**
     * 乗算を行います。破壊的です。
     *
     * @return (self_t) 乗算結果
     */
    self_t &operator *=(const FloatT &scalar){
      GSL_REAL(*this) *= scalar;
      GSL_IMAG(*this) *= scalar;
      return *this;
    }
    /**
     * 乗算を行います。
     *
     * @return (self_t) 乗算結果
     */
    self_t operator*(const FloatT &scalar) const{
      self_t result(*this);
      return (result *= scalar);
    }
    /**
     * 乗算を行います。
     *
     * @return (self_t) 乗算結果
     */
    friend self_t operator*(const FloatT &scalar, const self_t complex){return (complex * scalar);}

    /**
     * 除算を行います。破壊的です。
     *
     * @return (self_t) 除算結果
     */
    self_t &operator/=(const FloatT &scalar){return (*this) *= (1/scalar);}
    /**
     * 除算を行います。
     *
     * @return (self_t) 除算結果
     */
    self_t operator/(const FloatT &scalar) const{return (*this) * (1/scalar);}

    /**
     * 単項マイナスオペレータ。
     *
     * @return (self_t) 結果
     */
    self_t operator -() const{return ((*this) * -1);}

    /**
     * 加算を行います。破壊的です。
     *
     * @return (self_t) 加算結果
     */
    self_t &operator+=(const self_t &complex){
      GSL_REAL(*this) += (const_cast<Complex &>(complex)).real();
      GSL_IMAG(*this) += (const_cast<Complex &>(complex)).imaginary();
      return *this;
    }
    /**
     * 加算を行います。
     *
     * @return (self_t) 加算結果
     */
    self_t operator+(const self_t &complex) const{
      self_t result = *this;
      return (result += complex);
    }

    /**
     * 減算を行います。破壊的です。
     *
     * @return (self_t) 減算結果
     */
    self_t &operator-=(const self_t &complex){
      return ((*this) += (-complex));
    }
    /**
     * 減算を行います。
     *
     * @return (self_t) 減算結果
     */
    self_t operator-(const self_t &complex) const{
      return ((-complex) += (*this));
    }

    /**
     * 乗算を行います。
     *
     * @return (self_t) 乗算結果
     */
    self_t operator*(const self_t &complex) const{
      return self_t(gsl_complex_mul(*this, complex));
    }
    /**
     * 乗算を行います。破壊的です。
     *
     * @return (self_t) 乗算結果
     */
    self_t &operator*=(const self_t &complex){
      *this = gsl_complex_mul(*this, complex);
      return *this;
    }

    /**
     * 除算を行います。破壊的です。
     *
     * @return (self_t) 除算結果
     */
    self_t &operator/=(const self_t &complex){
      return (*this) *= (complex.conjugate() / complex.abs2());
    }
    /**
     * 除算を行います。
     *
     * @return (self_t) 除算結果
     */
    self_t operator/(const self_t &complex) const{
      self_t copy = (*this);
      return (copy /= complex);
    }

    /**
     * 除算を行います。
     *
     * @return (self_t) 除算結果
     */
    friend self_t operator/(const FloatT &scalar, const self_t complex){
      return self_t(scalar) / complex;
    }

#if ENABLE_IOSTREAM
    /**
     * みやすい形で複素数を出力します。
     *
     * @param out 出力ストリーム
     * @param complex 出力する複素数
     */
    friend std::ostream &operator<<(std::ostream &out, const self_t &complex){
      out << (const_cast<Complex &>(complex)).real() << " + "
          << (const_cast<Complex &>(complex)).imaginary() << "i";
      return out;
    }
#else
    /**
     * 見やすい形で出力します。
     *
     */
    void inspect(char *buffer, int buffer_size) const{
      snprintf(
        buffer, buffer_size,
        "%f + %fi",
        GSL_REAL(*this),
        GSL_IMAG(*this));
    }
#endif

    /**
     * eの指数乗を求めます。
     *
     * @param imaginary 虚部
     * @return (self_t) 結果
     */
    static self_t exp(const FloatT &imaginary){
      return self_t(cos(imaginary), sin(imaginary));
    }

    /**
     * eの指数乗を求めます。
     *
     * @param real 実部
     * @param imaginary 虚部
     * @return (self_t) 結果
     */
    static self_t exp(const FloatT &real, const FloatT &imaginary){
      return self_t::exp(imaginary) *= ::exp(real);
    }

    /**
     * eの指数乗を求めます。
     *
     * @param complex 複素数
     * @return (self_t) 結果
     */
    static self_t exp(const self_t &complex){
      return self_t::exp(
          const_cast<self_t &>(complex).real(),
          const_cast<self_t &>(complex).imaginary());
    }
};

#if 0
/**
 * expの複素数拡張
 *
 * @param real 実部
 * @param imaginary 虚部
 * @param FloatT 演算精度
 * @return (Complex<FloatT>) 結果
 */
template <class FloatT>
inline Complex<FloatT> iexp(const FloatT &real, const FloatT &imaginary){
  return Complex<FloatT>::exp(real, imaginary);
}

/**
 * expの複素数拡張
 * @param imaginary 虚部
 * @param FloatT 演算精度
 * @return (Complex<FloatT>) 結果
 */
template <class FloatT>
inline Complex<FloatT> iexp(const FloatT &imaginary){
  return Complex<FloatT>::exp(imaginary);
}

/**
 * expの複素数拡張
 *
 * @param complex 複素数
 * @param FloatT 演算精度
 * @return (Complex<FloatT>) 結果
 */
template <class FloatT>
inline Complex<FloatT> exp(const Complex<FloatT> &complex){
  return Complex<FloatT>::exp(complex);
}

/**
 * powの複素数拡張
 *
 * @param complex 複素数
 * @param FloatT 演算精度
 * @return (Complex<FloatT>) 結果
 */
template <class FloatT>
inline Complex<FloatT> pow(const Complex<FloatT> &complex, const FloatT &factor){
  return complex.power(factor);
}

/**
 * sqrtの複素数拡張
 *
 * @param complex 複素数
 * @param FloatT 演算精度
 * @return (Complex<FloatT>) 結果
 */
template <class FloatT>
inline Complex<FloatT> sqrt(const Complex<FloatT> &complex){
  return complex.sqrt();
}
#endif

#endif /* __COMPLEX_GSL_H */