/*  scomplex.h by K.Tsuru  */
/*********************************************************
SN library
SComplex class
It provides a multi-precision complex number arithmetic.
I referred to two styles.
(a) gcc's "complex"
friend SComplex& ccmult(SComplex* ths, const SComplex& z);

(b) B.Stroustrup "The C++ Programming Language Third Edition"
first operator *=() is defined
second operator *(a, b)

I paied attention to
1.copy constructor is heavy and spends memory.
2.I want to use statements &x == &y in the operator x == y to speed-up.
Almost part has been rewritten in ver.2.18.
*********************************************************/
#ifndef  SCOMPLEX_H
#define  SCOMPLEX_H
#ifndef SN_H
#include "sn.h"
#endif

class SComplex {
private:
    SDouble re, im;
/********* List of unusable functions which have no body.******/
    SComplex operator>(const SComplex&) const;
    SComplex operator>=(const SComplex&) const;
    SComplex operator<(const SComplex&) const;
    SComplex operator<=(const SComplex&) const;
public:
    // default constructor do nothing. For large size array declaration.
    SComplex() {} // "re" and "im" are not initialized and occupy no memory.
    // constructor by single SDouble values
    // No "explicit" causes (SD Op SD) error
    // Cannot use "SComplex c = x;"(NG) instead of "SComplex c(x);"(OK).
    explicit SComplex(const SDouble& rp) : re(rp), im(0.0){}
    // constructor by two SDouble values
    SComplex(const SDouble& rp, const SDouble& ip) : re(rp), im(ip){}

    // constructor by two double values since ver 2.18
    // This enables us to use "return 0.0".
    SComplex(double rp, double ip = 0.0) : re(rp), im(ip){}
    // default copy constructor is used according to refference (b)
//  SComplex(const SComplex& z) :re(z.re), im(z.im){}
    SComplex(const complex<double>& z) :re(z.real()), im(z.imag()){}
    ~SComplex(){}

    SDouble Real() const { return re;}  // real part
    SDouble Imag() const { return im;}  // imaginary part
    SDouble Norm() const;   // an inline function defined below
    SComplex Conj() const  { return SComplex(re, -im); }

    // multiple of i or -i
    SComplex& MultI() { im.ChangeSign(9401); std::swap(re, im); return *this; } // i*z
    SComplex& MultMI() { re.ChangeSign(9402); std::swap(re, im); return *this; } // -i*z = z/i

    uint Head() const {  return min( re.Head(), im.Head() ); }
    uint Tail() const {  return min( re.Tail(), im.Tail() ); }
    uint reHead() const {  return re.Head(); }
    uint reTail() const {  return re.Tail(); }
    uint imHead() const {  return im.Head(); }
    uint imTail() const {  return im.Tail(); }
    void FigureClear(uint from, uint to) {
        re.FigureClear(from, to);   im.FigureClear(from, to);
    }
    void SetInt(int i) { re.SetInt(i);  im.SetZero(); }
    void SetZero(){ re.SetZero();   im.SetZero(); }
    void ChangeSign(int id = 91){   // z -> -z
        re.ChangeSign(id);  im.ChangeSign(id);
    }
// If you want to change the real part only, use such as z.Set(newValue, z.Imag());
// The statement im = z.Imag(); does not copy, because &im == &ip.
    void Set(const SDouble& rp, const SDouble& ip) { re = rp; im = ip; }

    bool IsZero(int id=90) const { return (re.Sign(id) == 0) && (im.Sign(id) == 0); }
//sign operators
    SComplex operator+() const { return *this; }
    SComplex operator-() const;

  //Operators whose rhs is scalar.
// SComplex& operator=(const SDouble& x);       // SD Op SD error
    SComplex& operator+=(const SDouble& x);
    SComplex& operator-=(const SDouble& x);
    SComplex& operator*=(const SDouble& x);
    SComplex& operator/=(const SDouble& x);

/****** mixed-mode operations with double since ver 2.18 ******/
//  SComplex& operator=(double x);
    SComplex& operator+=(double x);
    SComplex& operator-=(double x);
    SComplex& operator*=(double x);
    SComplex& operator/=(double x);

  //Operators whose rhs is complex object.
    SComplex& operator=(const SComplex& z);
    SComplex& operator=(const complex<double>& z);
    SComplex& operator+=(const SComplex& z);
    SComplex& operator-=(const SComplex& z);
    SComplex& operator*=(const SComplex& z);    // 903
    SComplex& operator/=(const SComplex& z);    // 904

  // Multiplication by a small integer.
    SComplex& CsMult(ulong n);
  // Division by small integer.
    SComplex& CsDiv(ulong n);
/*
 At the present time it does not support reading from file.
 Please use SDouble class' Put(s) and write/read the real and imaginary parts, separately.
 Due to the value of "fmt" it outputs in the following form.
 iFMT     : (real part)[+|-]i*(imaginary part)
 BracketFMT : (real part, imaginary part)
 midCR : put crlf between real and imaginary part or not
 endCR : put crlf at line end or not
*/
    enum { iFMT = 64, BracketFMT = 128, MID_CR = 256, END_CR = 512};
    long Put(int fmt= MID_CR|iFMT) const;   // 900 change since ver.2.21
    long Puts(int fmt = MID_CR|iFMT) const{ return Put(fmt|END_CR); }
    friend ostream& operator<<(ostream& os, const SComplex& sl);    // added since version 2.21 902
    static int scLineFormat;  // set by BRACKET  | MID_CR| END_CR 
    static void SetFormat(int fmt);
    static long ioCount;
};
SComplex SetComplex(const char* s); // set SComplex value by (rp,ip).() is not necessary Use delimiter ',' only.
istream& operator>>(istream& s, SComplex& a);   // added since version 2.21  903

///// Cautions : All functions having its body must be "inline".///////
/**********************************
* Implementation of member functions
***********************************/
inline SDouble SComplex::Norm() const {
    return ( re * re + im * im );
}

inline SComplex SComplex::operator-() const {
    SComplex z(*this);
    z.re = -z.re;   z.im = -z.im;
    return z;
}

////// Operators whose rhs is real scalar. //////
inline SComplex& SComplex::operator+=(const SDouble& x) {
    re += x;    return *this;
}

inline SComplex& SComplex::operator-=(const SDouble& x) {
    re -= x;    return *this;
}

inline SComplex& SComplex::operator*=(const SDouble& x) {
    re *= x;    im *= x;    return *this;
}
inline SComplex& SComplex::operator/=(const SDouble& x) {
    re /= x;    im /= x;    return *this;
}

/****** Conversion SDouble to SComplex *******************/
/****************************
Cannot use SC+SD because SC + SD yeilds an error as [Notes]candidates are: operator+() const ... MinGW gcc 4.8.1
Please use SC + SComplex(SD, 0.0); or below
******************************/
inline SComplex SDtoSC(const SDouble& x) {
    return SComplex(x, 0.0);
}
inline SComplex SDtoISC(const SDouble& x) {
    return SComplex(0.0, x); // i*x
}
/****** mixed-mode operations with double since ver 2.18 ******/
/*
inline SComplex& SComplex::operator=(double x) {
    re = x;     im = 0.0;   return *this;
}
*/
inline SComplex& SComplex::operator+=(double x) {
    re += x;    return *this;
}

inline SComplex& SComplex::operator-=(double x) {
    re -= x;    return *this;
}

inline SComplex& SComplex::operator*=(double x) {
    re *= x;    im *= x;    return *this;
}

inline SComplex& SComplex::operator/=(double x) {
    re /= x;    im /= x;    return *this;
}

///// Operators whose rhs is complex object. /////
inline SComplex& SComplex::operator=(const complex<double>& z) {
    re = z.real();  im = z.imag();
    return *this;
}

// substitution operator SC=SC
inline SComplex& SComplex::operator=(const SComplex& z) { // default is not enough
    if(this == &z) return *this;    // z = z;
    re = z.re;  im = z.im;
    return *this;
}

inline SComplex& SComplex::operator+=(const SComplex& z) {
    re += z.re; im += z.im;
    return *this;
}

inline SComplex& SComplex::operator-=(const SComplex& z) {
    re -= z.re; im -= z.im;
    return *this;
}

/******* Non member functions in the following *******/
//// Two term operators /////
//// plus operators x + y
//// for MIXED_OPERATIONS_WITH_SD see below
//// plus operator x + y
inline SComplex operator+(const SComplex& x, const SComplex& y) {
    return SComplex(x.Real() + y.Real(), x.Imag() + y.Imag());
}

//// minus operator x - y
inline SComplex operator-(const SComplex& x, const SComplex& y) {
    return SComplex(x.Real() - y.Real(), x.Imag() - y.Imag());
}

//// multiplication operator x * y
inline SComplex operator*(const SComplex& x, const SComplex& y) {
    SComplex r = x;
    return r *= y;
}

//// division operator x / y
inline SComplex operator/(const SComplex& x, const SComplex& y) {
    SComplex r = x;
    return r /= y;
}

// Multiplication by a small integer.
inline SComplex& SComplex::CsMult(ulong n) {
    if(n) {
        re = DsMult(re, n); im = DsMult(im, n);
    } else SetZero();
    return *this;
}
inline SComplex CsMult(const SComplex& z, ulong n) {
    SComplex r = z;
    return r.CsMult(n);
}
// Division by a small integer.
inline SComplex& SComplex::CsDiv(ulong n) {
    re = DsDiv(re, n);  im = DsDiv(im, n);
    return *this;
}
inline SComplex CsDiv(const SComplex& z, ulong n) {
    SComplex r = z;
    return r.CsDiv(n);
}

/****** mixed-mode operations with double since ver 2.18 ******/
inline SComplex operator+(const SComplex& x, double d) {            // SC + d
    return SComplex(x.Real() + d, x.Imag());
}
inline SComplex operator+(double d, const SComplex& y) {            // d + SC
    return SComplex(d + y.Real(), y.Imag());
}

inline SComplex operator-(const SComplex& x, double d) {            // SC - d
    return SComplex(x.Real() - d, x.Imag());
}
inline SComplex operator-(double d, const SComplex& y) {            // d - SC
    return SComplex(d - y.Real(), -y.Imag());
}

inline SComplex operator*(const SComplex& x, double d) {            // SC * d
    return SComplex(x.Real() * d, x.Imag() * d);
}

inline SComplex operator*(double d, const SComplex& y) {            // d * SC
    return SComplex(d * y.Real(), d * y.Imag());
}

inline SComplex operator/(const SComplex& x, double d) {            // SC / d
    return SComplex(x.Real() / d, x.Imag() / d);
}
inline SComplex operator/(double x, const SComplex& y) {            // d / SC
    SComplex r = x;
    return r /= y;
}

/*****************
 basic functions
******************/
//// relational operators
inline bool operator==(const SComplex& x, const SComplex& y){
    return (x.Real()==y.Real()) && (x.Imag()==y.Imag());
}
inline bool operator!=(const SComplex& x, const SComplex& y){
    return !(x == y);
}

///// using binary splitting method ////////
SComplex CpolarBS(const SDouble& r, const SDouble& x); // r*{cos(x)+i*sin(x)}; 9120
SComplex CexpBS(const SComplex& z); // exp(z) 9121

SComplex CcosBS(const SComplex& z);             // cos z 9121
SComplex CsinBS(const SComplex& z);             // sin z 9122
SComplex CtanBS(const SComplex& z);             // tan z 9123

static const SComplex IU(0.0, 1.0);  // i (imaginary unit)
static const SComplex MI(0.0, -1.0);  // 1/i = -i


////////////////////////////////////////////////////////////////////////////////////
// For reference in the following
#define MIXED_OPERATIONS_WITH_SD 0
#if MIXED_OPERATIONS_WITH_SD
/*
 These operators causes error for 1 + SD etc.
 because candidates are : SD op SD and SC op SC. 
 Then set 0 above.
 Please use "+=", etc. operator in {} below.
*/
inline SComplex& SComplex::operator=(const SDouble& x) {
    re = x;     im = 0.0;   return *this;
}

inline SComplex operator+(const SComplex& x, const SDouble& y) {    // SC + SD
    SComplex r = x;
    return r += y;
}

inline SComplex operator+(const SDouble& x, const SComplex& y) {    // SD + SC
    SComplex r = y;
    return r += x;
}
inline SComplex operator-(const SComplex& x, const SDouble& y) {
    SComplex r = x;
    return r -= y;
}

inline SComplex operator-(const SDouble& x, const SComplex& y) {
    SComplex r(x, 0.0);
    return r -= y;
}

inline SComplex operator*(const SComplex& x, const SDouble& y) {
    SComplex r = x;
    return r *= y;
}
inline SComplex operator*(const SDouble& x, const SComplex& y) {
    SComplex r = y;
    return r *= x;
}

inline SComplex operator/(const SComplex& x, const SDouble& y) {
    SComplex r = x;
    return r /= y;
}
inline SComplex operator/(const SDouble& x, const SComplex& y) {
    SComplex r(x, 0.0);
    return r /= y;
}
#endif  // MIXED_OPERATIONS_WITH_SD
////////////////////////////////////////////////////////////////////////////////////

#endif  //SCOMPLEX_H