builtin_cplx.h

Go to the documentation of this file.

// Written: KG  98/08/13
// Changed: KG  98/08/14
// $Id: builtin_cplx.h,v 1.13.2.11 2019/05/28 11:13:02 garloff Exp $

#include "tbci/constants.h"
//#include <iostream>
#include "tbci/basics.h"

#ifndef HAVE_BUILTIN_CPLX
# include "tbci/cplx.h"
#else

#ifndef TBCI_BUILTIN_CPLX_H
#define TBCI_BUILTIN_CPLX_H  

#ifndef CTYP
# define CTYP double
#endif
 
NAMESPACE_TBCI
 
typedef CTYP cplx_base;

#ifndef CPLX
 typedef __complex__ CTYP cplx_t; //ALIGN(MIN_ALIGN2) 
# define CPLX(t) cplx_t
# define CPLXR(t,a,r) cplx_t a (r)
# define CPLXRI(t,a,r,i) cplx_t a; __real__ a = (r); __imag__ a = (i)
# define CPLXC(t,a,c) cplx_t a (c)
# define CPLXVRI(t,r,m) ((cplx_t) ((r) + (m) * 1.0i))
//# define CPLXVRI(t,r,i) (hcplx(r, i))
#endif

#define USE_BUILTIN_CPLX 1

//template <typename CTYP>
class hcplx
{
 protected:
        cplx_t cx ALIGN(MIN_ALIGN2);
 
 public:
        hcplx () {}
        hcplx (const hcplx & c) : cx(c.cx)  {}
        //hcplx (const hcplx & c)  { __real__ cx = __real__ c.cx; __imag__ cx = __imag__ c.cx; }
        hcplx (const cplx_base r, const cplx_base i = 0) { __real__ cx = r; __imag__ cx = i; }
        ~hcplx () {}

        // The most important ones: Provide type conversion
        hcplx (const cplx_t& c) : cx(c)  {}
        //hcplx (const cplx_t& c)  { __real__ cx = __real__ c; __imag__ cx = __imag__ c; }
        //operator cplx_t () const  { return cx; }
        //operator cplx_t& ()       { return cx; }
        operator const cplx_t& () const  { return cx; }
        //hcplx (const cplx_t & c) : cx(c) {}
        //operator cplx_t & () const { return cx; }
        
        /*cplx_base*/
        double fabssqr () const 
        { return (__real__ cx) * (__real__ cx) + (__imag__ cx) * (__imag__ cx); }
        double fabs () const
        { return MATH__ sqrt (fabssqr()); }
        //friend double fabssqr (const hcplx& c)  { return c.fabssqr(); }
        //friend double fabs (const hcplx& c)  { return c.fabs(); }
        
#ifdef PREFER_REF
        const cplx_base& real () const  { return __real__ cx; }
        const cplx_base& imag () const  { return __imag__ cx; }
#else
        cplx_base real () const  { return __real__ cx; }
        cplx_base imag () const  { return __imag__ cx; }
#endif
        cplx_base& real ()  { return __real__ cx; }
        cplx_base& imag ()  { return __imag__ cx; }
        
        cplx_t operator - () const 
        { cplx_t r; __real__ r = - __real__ cx; __imag__ r = - __imag__ cx; return r; }
};

//compat
inline double fabssqr (const cplx_t& a)
{ return (__real__ a) * (__real__ a) + (__imag__ a) * (__imag__ a); }

inline double fabssqr (const TBCI__ hcplx& c)  { return c.fabssqr(); }

#if 0
inline cplx_t operator - (const cplx_t& c) 
{ hcplx cx(c); return -cx; }
#endif

NAMESPACE_END
NAMESPACE_CSTD
inline double fabs (const TBCI__ cplx_t& a)
{ return MATH__ sqrt (TBCI__ fabssqr (a)); }

inline TBCI__ cplx_base abs (const TBCI__ cplx_t& a)
{ return MATH__ sqrt ((__real__ a) * (__real__ a) + (__imag__ a) * (__imag__ a)); }

inline double fabs (const TBCI__ hcplx& c)  
{ return c.fabs(); }

NAMESPACE_CSTD_END

#if 0 //def USE_NS
inline double MATH__ fabs (const TBCI__ cplx_t& c)
{ return TBCI__ fabs (c); }

//inline T MATH__ abs (const TBCI__ cplx_t& c)
//{ return TBCI__ abs (c); }
#endif

NAMESPACE_CPLX
#ifdef PREFER_REF
inline const TBCI__ cplx_base& real (const TBCI__ hcplx& c)  { return c.real(); }
inline const TBCI__ cplx_base& imag (const TBCI__ hcplx& c)  { return c.imag(); }
#else
inline TBCI__ cplx_base real (const TBCI__ hcplx& c)  { return c.real(); }
inline TBCI__ cplx_base imag (const TBCI__ hcplx& c)  { return c.imag(); }
#endif
inline TBCI__ cplx_base& real (const TBCI__ cplx_t& c) { return __real__ c; }
inline TBCI__ cplx_base& imag (const TBCI__ cplx_t& c) { return __imag__ c; }

inline TBCI__ cplx_t conj (const TBCI__ cplx_t& c)
{ return ~c; }

NAMESPACE_CPLX_END
NAMESPACE_TBCI

inline cplx_t dot (const cplx_t a, const cplx_t b)
{ return hcplx ((__real__ a) * (__real__ b) + (__imag__ a) * (__imag__ b),
                (__imag__ a) * (__real__ b) - (__real__ a) * (__imag__ b)); }

inline bool operator <= (const hcplx & a, const hcplx & b)
{ return (fabssqr (a) <= fabssqr(b)); }
inline bool operator <  (const hcplx & a, const hcplx & b)
{ return (fabssqr (a) < fabssqr(b)); }
inline bool operator >= (const hcplx & a, const hcplx & b)
{ return (fabssqr (a) >= fabssqr(b)); }
inline bool operator >  (const hcplx & a, const hcplx & b)
{ return (fabssqr (a) > fabssqr(b)); }

inline bool operator <= (const hcplx & a, const cplx_t & b)
{ return (fabssqr (a) <= fabssqr(b)); }
inline bool operator <  (const hcplx & a, const cplx_t & b)
{ return (fabssqr (a) < fabssqr(b)); }
inline bool operator >= (const hcplx & a, const cplx_t & b)
{ return (fabssqr (a) >= fabssqr(b)); }
inline bool operator >  (const hcplx & a, const cplx_t & b)
{ return (fabssqr (a) > fabssqr(b)); }

inline bool operator <= (const cplx_t & a, const hcplx & b)
{ return (fabssqr (a) <= fabssqr(b)); }
inline bool operator <  (const cplx_t & a, const hcplx & b)
{ return (fabssqr (a) < fabssqr(b)); }
inline bool operator >= (const cplx_t & a, const hcplx & b)
{ return (fabssqr (a) >= fabssqr(b)); }
inline bool operator >  (const cplx_t & a, const hcplx & b)
{ return (fabssqr (a) > fabssqr(b)); }

#ifdef NEED_BCPLX_COMP_OP
// Does not compile ...
inline bool operator <= (const cplx_t & a, const cplx_t & b)
{ return (fabssqr (a) <= fabssqr(b)); }
inline bool operator <  (const cplx_t & a, const cplx_t & b)
{ return (fabssqr (a) < fabssqr(b)); }
inline bool operator >= (const cplx_t & a, const cplx_t & b)
{ return (fabssqr (a) >= fabssqr(b)); }
inline bool operator >  (const cplx_t & a, const cplx_t & b)
{ return (fabssqr (a) > fabssqr(b)); }
#endif

inline STD__ ostream& operator << (STD__ ostream& os, const cplx_t & c)
{ os << "( " << __real__ c;
  if ((__imag__ c) < 0) return os << " - i " << - __imag__ c << " )";
  else return os << " + i " << __imag__ c << " )";
}

inline STD__ istream& operator >> (STD__ istream& in, cplx_t& c)
{
  cplx_base r, i;
  char s = '+', b = 'i', bra = '(';
  
  if (! (in >> bra))
    return in;
  if (bra != '(') in.putback(bra);
  if (! (in >> r))
    return in;
  in >> s;
  if ((s != '+') && (s != '-')) in.putback(s);
  in >> b;
  if ((b != 'i')) in.putback(b);
  if (! (in >> i))
    return in;
  __real__ c = r;
  if (s == '-') __imag__ c =- i;
  else __imag__ c = i;
  if (bra == '(') 
    {
      if (! (in >> bra))
        return in;
      if (bra != ')') in.putback(bra);
    }
  return in;
}

inline cplx_base polar (const cplx_t& c)
{
  if (__real__ c == 0)
    {
      if (__imag__ c > 0) return 0.5*pi;
      else if (__imag__ c < 0) return 1.5*pi;
           else return 0;
    }
  
  if (__imag__ c >= 0) return MATH__ atan2 (__imag__ c, __real__ c);
  else return 2*pi + MATH__ atan2 (__imag__ c, __real__ c);
}
        
inline cplx_t expi (const cplx_base phi)
{
  cplx_base tmp = fmod (phi, pi2);
  return cplx_t (MATH__ cos(tmp) + MATH__ sin(tmp) * 1.0i);
}

inline cplx_t sqr (const cplx_t &c)
{
  return fabssqr (c) * expi (fmod ((cplx_base)2 * polar(c), pi2));
}


#ifdef NEED_IMAG_UNIT
const cplx_t Imag_unit = CPLXVRI (double, 0.0, 1.0);
//const CPLXRI(double,Imag_unit,0.0,1.0);
#endif

NAMESPACE_END

#ifdef ENABLE_BCPLX
# include "obsolete/b_cplx_old.h"
#endif

NAMESPACE_CSTD  

inline TBCI__ cplx_t exp (const TBCI__ cplx_t& z)
{
  return MATH__ exp (__real__ z) * TBCI__ expi (__imag__ z);
}

inline TBCI__ cplx_t sqrt (const TBCI__ cplx_t& c)
{
  return MATH__ sqrt (MATH__ fabs(c)) * TBCI__ expi ((TBCI__ cplx_base)0.5 * TBCI__ polar(c));
}

NAMESPACE_CSTD_END
#endif /* TBCI_BUILTIN_CPLX_H */
#endif /* __GNUC__ */