TComplex.h

Go to the documentation of this file.

//                                                                      //
// TComplex                                                             //
//                                                                      //
// Implement complex numbers in ROOT                                    //
//                                                                      //

#ifndef LITR_TComplex
#define LITR_TComplex


#undef   TCOMPLEX_USE_RTYPES_MACROS
#define  TCOMPLEX_USE_DOUBLE

#ifdef TCOMPLEX_USE_DOUBLE
# define TComplexBase  Double_t
#else
# define TComplexBase  Real_t
#endif

#include <Riostream.h>
#include <TBuffer.h>
#include <TMath.h>

class TCplxRhoPhi {
protected:

  Double_t fRho;    //module of complex number
  Double_t fPhi;    //phase of complex number
public:
  TCplxRhoPhi(Double_t x,Double_t y) { fRho = TMath::Sqrt(x*x+y*y);
                                       fPhi = atan2(y,x); }
  Double_t GetRho() { return fRho; }
  Double_t GetPhi() { return fPhi; }
  void Sqrt() { fRho = TMath::Sqrt(fRho); fPhi = fPhi/2.0; }
};

class TComplex {

protected:
  friend TComplex operator + (const Double_t,const TComplex&);
  friend TComplex operator - (const Double_t,const TComplex&);
  friend TComplex operator * (const Int_t,const TComplex&);
  friend TComplex operator * (const Double_t,const TComplex&);
  friend TComplex operator * (const TComplex&,const TComplex&);
  friend TComplex operator * (TComplex&, TComplex&);
  friend TComplex operator % (const TComplex&,const TComplex&);
  friend TComplex operator / (const Double_t,const TComplex&);
  friend Bool_t   operator == (const TComplex&, const TComplex&);
  friend Bool_t   operator != (const TComplex&, const TComplex&);
  friend Bool_t   operator == (const TComplex&, const Double_t&);
  friend Bool_t   operator != (const TComplex&, const Double_t&);
  friend Double_t Abs(const TComplex&);
  friend Double_t Abs2(const TComplex&);
  friend TComplex Sqrt(const TComplex&);
  friend TComplex Sqrt3(const TComplex&, Int_t);
  friend TComplex Exp(const TComplex&);
  friend TComplex Log(const TComplex&);
  friend TComplex Log10(const TComplex&);
  friend TComplex Log2(const TComplex&);
  friend TComplex Sin(const TComplex&);
  friend TComplex Cos(const TComplex&);
  friend TComplex Power(const TComplex &z,Int_t n) { return Exp(n*Log(z)); }
  friend TComplex Power(const TComplex &z,Double_t x) { return Exp(x*Log(z)); }
  friend Double_t Phase(TComplex  z) { return TMath::ATan2(z.fIm,z.fRe);};
  friend istream &operator>>(istream&,TComplex&);
  friend ostream &operator<<(ostream&,TComplex&);

public:

  TComplexBase fRe;     //real part
  TComplexBase fIm;     //imaginary part
  TComplex() { fRe=0.0; fIm =0.0; }
  TComplex(Double_t x,Double_t y=0) { fRe = x; fIm=y; }
  TComplex(const TComplex &z) { fRe = z.fRe; fIm = z.fIm; }
  Double_t Re()  const  { return fRe; }
  Double_t Im()  const  { return fIm; }
  Double_t Phase() { return TMath::ATan2(fIm,fRe);};
  void Set(Double_t x,Double_t y) { fRe = x; fIm = y; }
  void SetReal(Double_t x) { fRe = x; }
  void SetIm(Double_t y)   { fIm = y; }
  void C()                 { fIm = - fIm; }
  void PureReal()          { fIm = 0.0; }
  void PureIm()            { fRe = 0.0; }
  void SetNeg()            { fRe = -fRe; fIm = -fIm; }
  void Dump();

  operator Bool_t() { return fRe||fIm; };
  operator Double_t() { return fRe;};
  operator Int_t() { return Int_t(fRe);};
  
  

  Bool_t    operator!() { return (!fRe)&&(!fIm);};
  Bool_t    operator&&(TComplex &z) {
    return Bool_t()&&Bool_t(z);
  }

  TComplex  operator+(const TComplex& z) {
    //+ with other complex number
    TComplex z4;
    z4.fRe = fRe + z.fRe;
    z4.fIm = fIm + z.fIm;
    return z4;
  }

  TComplex  operator+(const Double_t x) {
    TComplex z4;
    z4.fRe = fRe + x;
    z4.fIm = fIm;
    return z4;
  }
  TComplex  operator+(const Int_t x) {
    TComplex z4;
    z4.fRe = fRe + x;
    z4.fIm = fIm;
    return z4;
  }

  TComplex operator -() {
    TComplex z;
    z.fRe = -fRe;
    z.fIm  = -fIm;
    return z;
  }

  TComplex  operator-(const TComplex& z) {
    //+ with other complex number
    TComplex z4;
    z4.fRe = fRe - z.fRe;
    z4.fIm = fIm - z.fIm;
    return z4;
  }
  
  TComplex  operator-(const Double_t x) {
    //- with other complex number
    TComplex z4;
    z4.fRe = fRe - x;
    z4.fIm = fIm;
    return z4;
  }
  TComplex  operator-(const Int_t x) {
    //- with other complex number
    TComplex z4;
    z4.fRe = fRe - x;
    z4.fIm = fIm;
    return z4;
  }





  TComplex& operator=(const TComplex &z) {
    //defines = for complex
    if (this != &z) {
      fRe = z.fRe;
      fIm = z.fIm;
    }
    return *this;
  }
  TComplex& operator=(const Double_t x) {
    //defines = for complex
    fRe = x;
    fIm = 0.0;
    return *this;
  }

  
  TComplex& operator+=(const Double_t& x) {
    fRe += x;
    return *this;
  }
  TComplex& operator+=(const TComplex& z) {
    //Add itself with a complex numbers
    fRe += z.fRe;
    fIm += z.fIm;
    return *this;
  }

  TComplex& operator-=(const Double_t& x) {
    fRe -= x;
    return *this;
  }
  TComplex& operator-=(const TComplex& z) {
    //Add itself with a complex numbers
    fRe -= z.fRe;
    fIm -= z.fIm;
    return *this;
  }


  TComplex& operator/=(const Double_t& x) {
    fRe /= x;
    fIm /= x;
    return *this;
  }
  
  TComplex& operator/=(const TComplex& z) {
    Double_t tmp,d;
    d = z.fRe*z.fRe + z.fIm*z.fIm;
    tmp = (fRe*z.fRe + fIm*z.fIm)/d;
    fIm = (fIm*z.fRe - fRe*z.fIm)/d;
    fRe = tmp;
    
    return *this;
  }
  
  
  TComplex& operator*=(const Double_t & x) {
    //Multiply itself with a float numbers
    fIm *=x;
    fRe *=x;
    return *this;
  }

  TComplex& operator*=(const TComplex& z) {
    //Multiply itself with a complex numbers
    register Double_t tmp =  fRe*z.fRe - fIm*z.fIm;
    fIm = fRe*z.fIm + fIm*z.fRe;
    fRe = tmp;
    return *this;
  }


  
  
  void FromRhoPhi(TCplxRhoPhi &r) {
    //returns from phase representation
    register Double_t rho;
    register Double_t phi;
    rho = r.GetRho();
    phi = r.GetPhi();
    fRe = rho*TMath::Cos(phi);
    fIm = rho*TMath::Sin(phi);
  }
  Bool_t IsComplex(Double_t limreal) {
    //if the abs. value of the imaginary part is smaller than limreal, make the
    //number real and returns kFALSE else returns kTRUE.
    const Double_t zero = 0.0;
    Bool_t cplx = kTRUE;
    if (TMath::Abs(fIm)<TMath::Abs(limreal)) {
      fIm = zero;
      cplx = kFALSE;
    }
    return cplx;
  }
  void BetterConj(TComplex &z2) {
    //Make this and z2 better complex conjugate
    const Double_t deux = 2.0;
    Double_t a,b;
    a = (fRe + z2.fRe)/deux;
    b = (fIm - z2.fIm)/deux;
    z2.fRe =  a;
    z2.fIm = -b;
    fRe    = a;
    fIm    = b;
  }
  void CosFromSin() {
    //  Given a complex number representing the sinus of a complex angle, returns
    //the complex value of the cosinus of the same complex angle. The value chosen
    //is the one with a positive real part, if the imaginary part is 0.
    //If the imaginary part is non-zero, the value chosen is the one with a
    //negative imaginary part. THE IMAGINARY PART OF THE RESULT IS ALWAYS CHOSEN
    //NEGATIVE IF IT IS NON-ZERO. This is an arbitrary choice motivated by the
    //fact that in Litrani, the phases of waves have to be negative when a wave
    //is absorbed, never positive which would correspond to an unphysical explosion
    //of the wave.
    const Double_t zero   = 0.0;
    const Double_t un     = 1.0;
    const Double_t vsmall = 1.0e-12;
    TComplex z;
    z = *this;z*=z;
    z = Sqrt(un - z);
    if (z.fIm>vsmall) z = zero - z;
    if (TMath::Abs(z.fIm)<=vsmall) z.fIm = zero;
    *this = z;
  }
  void RPhi(Double_t &r, Double_t &phi) {
    //  Calculates module (with sign !) and phase of a complex number.
    //If the real part of the complex number is negativ, returns a negative module.
    //Advantage : in case of a negative real number, the module stays negative and
    //the phase 0. A negative real number is not considered as a non-real, complex
    //number with a phase of pi.
    const Double_t zero   = 0.0;
    const Double_t wsmall = 1.0e-300;
    Double_t axr,axi;
    axr = TMath::Abs(fRe);
    axi = TMath::Abs(fIm);
    if ((axr<wsmall) && (axi<wsmall)) {
      r   = zero;
      phi = zero;
    }
    else {
      r   = TMath::Sqrt(fRe*fRe+fIm*fIm);
      if (fRe<zero) r = -r;
      if (axi<wsmall) {
    phi = zero;
      }
      else {
    phi = atan2(fIm,fRe);
    if (fRe<zero) {
      if (fIm>=zero) phi = phi - TMath::Pi();
      else phi = phi + TMath::Pi();
    }
      }
    }
  }
  
  TComplex operator*(const Double_t x) {
    TComplex z4;
    z4.fRe = fRe*x;
    z4.fIm = fIm*x;
    return z4;
  }
  TComplex operator*(const Int_t x) {
    TComplex z4;
    Double_t xx = x;
    z4.fRe = fRe*xx;
    z4.fIm = fIm*xx;
    return z4;
  }

  TComplex operator/(const TComplex &z) {
    //division with other complex number
    TComplex z4;
    Double_t d;
    d = z.fRe*z.fRe + z.fIm*z.fIm;
    z4.fRe = (fRe*z.fRe + fIm*z.fIm)/d;
    z4.fIm = (fIm*z.fRe - fRe*z.fIm)/d;
    return z4;
  }
  TComplex operator/(const Double_t x) {
    //division with other complex number
    TComplex z4;
    z4.fRe = fRe/x;
    z4.fIm = fIm/x;
    return z4;
  }
  TComplex operator/(const Int_t x) {
    //division with other complex number
    TComplex z4;
    Double_t xx = x;
    z4.fRe = fRe/xx;
    z4.fIm = fIm/xx;
    return z4;
  }



  void Print();
  void Test(Double_t,const TComplex&);

#ifdef TCOMPLEX_USE_RTYPES_MACROS
  
  // The is the standard version using Rtypes.h-Macro
  //
  // ==> this results in a sizeof(TComplex)==20

  ClassDef(TComplex,0)
    
#else

  // This Part I got from $(ROOTSYS)/include/Rtypes.h
  // and it correspond to the ClassDef-Macro. But because
  // we do not want to inherits from this base type
  //
  // ==> this results in a sizeof(TComplex)==16

 private: 
  static TClass *fgIsA; 
 public: 
  static TClass *Class(); 
  static const char *Class_Name(); 
  static Version_t Class_Version() { return 0; } 
  static void Dictionary(); 
  TClass *IsA() const { return TComplex::Class(); } 
  void ShowMembers(TMemberInspector &insp, char *parent); 
  void Streamer(TBuffer &b); 
  void StreamerNVirtual(TBuffer &b) { TComplex::Streamer(b); } 
  static const char *DeclFileName() { return __FILE__; } 
  static int DeclFileLine() { return __LINE__; } 
  static const char *ImplFileName(); 

# if !defined(R__ACCESS_IN_SYMBOL) || defined(__CINT__)
  
  friend void ROOT__ShowMembersFunc(TComplex *obj, TMemberInspector &R__insp, 
                    char *R__parent); 
  
# endif
  static int ImplFileLine();

#endif

};

// Global Functions

inline Double_t Dist(TComplex& z1, TComplex& z2) {
  return Abs(z1-z2);
}

inline TComplex operator + (const Double_t x, const TComplex &z) {
//addition of a real and a complex
  return TComplex(x + z.fRe,z.fIm);
}

inline TComplex operator - (const Double_t x, const TComplex &z) {
//substraction of a complex from a real
  return TComplex(x - z.fRe,-z.fIm);
}
inline TComplex operator * (const Int_t x, const TComplex &z) {
//multiplication of a real by a complex
  return TComplex(x*z.fRe,x*z.fIm);
}
inline TComplex operator * (const Double_t x, const TComplex &z) {
//multiplication of a real by a complex
  return TComplex(x*z.fRe,x*z.fIm);
}
inline TComplex operator * (const TComplex &x, const TComplex &y) {
//multiplication of a complex by a complex
  return TComplex(x.fRe*y.fRe - x.fIm*y.fIm,x.fRe*y.fIm + x.fIm*y.fRe);
}
inline TComplex operator * (TComplex &x, TComplex &y) {
//multiplication of a complex by a complex
  return TComplex(x.fRe*y.fRe - x.fIm*y.fIm,x.fRe*y.fIm + x.fIm*y.fRe);
}
inline TComplex operator % (const TComplex &x, const TComplex &y) {
//multiplication of a complex-transposed by a complex
  return TComplex(x.fRe*y.fRe + x.fIm*y.fIm,x.fRe*y.fIm - x.fIm*y.fRe);
}
inline TComplex operator / (const Double_t x, const TComplex &z) {
//division of a real by a complex
  TComplex z4;
  Double_t d;
  d = z.fRe*z.fRe + z.fIm*z.fIm;
  z4.fRe =  (x*z.fRe)/d;
  z4.fIm = -(x*z.fIm)/d;
  return z4;
}

inline Bool_t operator ==(const TComplex &z1, const TComplex &z2) {
  return ((z1.fRe==z2.fRe)&&(z1.fIm==z2.fIm));
}

inline Bool_t operator !=(const TComplex &z1, const TComplex &z2) {
  return ((z1.fRe!=z2.fRe)||(z1.fIm!=z2.fIm));
}
inline Bool_t operator ==(const TComplex &z1, const Double_t &d) {
  return ((z1.fRe==d)&&(z1.fIm==0.0));
}

inline Bool_t operator !=(const TComplex &z1, const Double_t &d) {
  return ((z1.fRe!=d)||(z1.fIm!=0.0));
}
inline Double_t Abs(const TComplex &z) {
//Calculates the absolute value of a complex number
  return TMath::Sqrt(Abs2(z));
}
inline Double_t Abs2(const TComplex &z) {
//Calculates the absolute value of a complex number
  return z.fRe*z.fRe + z.fIm*z.fIm;
}

inline TComplex Sqrt(const TComplex &z) {
//Calculates the square root of a complex number
  const Double_t zero = 0.0;
  TComplex zz;
  if (z.fRe >= zero) {
    TCplxRhoPhi r(z.fRe,z.fIm);
    r.Sqrt();
    zz.FromRhoPhi(r);
  }
  else {
    TCplxRhoPhi r(-z.fRe,-z.fIm);
    r.Sqrt();
    zz.FromRhoPhi(r);
    TComplex j(0,1);
    //zz = zz*j;
    zz*=j;
  }
  return zz;
}

inline TComplex Sqrt3(const TComplex &z,Int_t k) {
//
//  Calculates the cubic root of a complex number.
//  If the real part is negative, calculates the cubic root of minus the complex
//number and then changes the sign of the solution.
//  There are 3 possible solutions. k decides which one is taken :
//    - k = 0  ==>  phi --> phi/3
//    - k = 1  ==>  phi --> (phi+2pi)/3
//    - k = 2  ==>  phi --> (phi-2pi)/3
//
  const Double_t zero   = 0.0;
  const Double_t deux   = 2.0;
  const Double_t trois  = 3.0;
  const Double_t small  = 1.0e-7;
  const Double_t wsmall = 1.0e-300;
  Double_t rho,phi;
  Bool_t isreal,isneg;
  TComplex z1,zz;
  isreal = (TMath::Abs(z.fIm)<wsmall);
  isneg  = (z.fRe<zero);
  if (isneg) z1 = zero - z;
  else       z1 = z;
  k      = k % 3;
  TCplxRhoPhi r(z1.fRe,z1.fIm);
  rho = r.GetRho();
  phi = r.GetPhi();
  rho = TMath::Exp(TMath::Log(rho)/trois);
  switch (k) {
  case 0:
    phi = phi/trois;
    break;
  case 1:
    phi = (phi + deux*TMath::Pi())/trois;
    break;
  case 2:
    phi = (phi - deux*TMath::Pi())/trois;
    break;
  }
  zz.fRe = rho*TMath::Cos(phi);
  zz.fIm = rho*TMath::Sin(phi);
  if (isreal && (TMath::Abs(zz.fIm)<small)) zz.fIm = zero;
  if (isneg) zz = zero - zz;
  return zz;
}

inline TComplex Exp(const TComplex &z) {
//Calculates the exponential of a complex number
  register Double_t ex;
  TComplex zz;
  ex = TMath::Exp(z.fRe);
  zz.fRe = ex*TMath::Cos(z.fIm);
  zz.fIm = ex*TMath::Sin(z.fIm);
  return zz;
}

inline TComplex Log(const TComplex &z) {
//Calculates the logarithm of a complex number
  TComplex zz;
  TCplxRhoPhi r(z.fRe,z.fIm);
  zz.fRe = TMath::Log(r.GetRho());
  zz.fIm = r.GetPhi();
  return zz;
}
inline TComplex Log2(const TComplex &z) {
//Calculates the logarithm of a complex number
  TComplex zz;
  TCplxRhoPhi r(z.fRe,z.fIm);
  zz.fRe = TMath::Log2(r.GetRho());
  zz.fIm = r.GetPhi();
  return zz;
}
inline TComplex Log10(const TComplex &z) {
//Calculates the logarithm of a complex number
  TComplex zz;
  TCplxRhoPhi r(z.fRe,z.fIm);
  zz.fRe = TMath::Log10(r.GetRho());
  zz.fIm = r.GetPhi();
  return zz;
}

inline TComplex Sin(const TComplex &z) {
//Calculates the sinus of a complex number
  const Double_t un   = 1.0;
  const Double_t deux = 2.0;
  TComplex zz;
  Double_t a,b;
  a = TMath::Exp(z.fIm);
  b = un/a;
  zz.fRe = ((a+b)*TMath::Sin(z.fRe))/deux;
  zz.fIm = ((a-b)*TMath::Cos(z.fRe))/deux;
  return zz;
}

inline TComplex Cos(const TComplex &z) {
//Calculates the cosinus of a complex number
  const Double_t un   = 1.0;
  const Double_t deux = 2.0;
  TComplex zz;
  Double_t a,b;
  a = TMath::Exp(z.fIm);
  b = un/a;
  zz.fRe = ((a+b)*TMath::Cos(z.fRe))/deux;
  zz.fIm = ((b-a)*TMath::Sin(z.fRe))/deux;
  return zz;
}


extern       TComplex ComplexZero;
extern       TComplex ComplexOne;
extern       TComplex I;
#endif

---