Complex.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of BORDERBASIX software. *
*      (C) B. Mourrain,  INRIA                                       *
**********************************************************************/
/*
  Complex.hpp
  definitions of two functions :
  Root_Unity : returns a n- root of unity
  pow : power function to a complex

  complex<Scl<MP_FLOAT> > specialization

*/

#ifndef INCLUDED_COMPLEXH
#define INCLUDED_COMPLEXH
#include <complex.h>
#define COMPLEX complex<double>

#ifdef _GMP
#include "Scl_mpf.hpp"
#endif //

struct Warning_Value_Type_Not_Defined {};

template <class T> struct Real {
  typedef Warning_Value_Type_Not_Defined TYPE; 
};

#ifdef _GMP
template <> struct Real<complex<Scl<MP_FLOAT> > > {
  typedef Scl<MP_FLOAT> TYPE;
  static inline void assign(complex<Scl<MP_FLOAT> > & a, const complex<Scl<MP_FLOAT> > & b) {a = b ;};
  static inline Scl<MP_FLOAT> absolute(const complex<Scl<MP_FLOAT> > & a) 
    {return sqrt(pow(a.real(),2)+pow(a.imag(),2)); }
    };
#endif //

template <> struct Real<double> {
  typedef double TYPE;
  static inline void assign(double & a, const complex<double> & b) { a = b.real() ;};
  static inline double absolute(const double & a) { return abs(a); }
};

#ifdef _GMP
struct MP_FLOAT;
template <> struct Real<Scl<MP_FLOAT> > {
  typedef Scl<MP_FLOAT> TYPE;
  static inline void assign(Scl<MP_FLOAT> & a, const complex<Scl<MP_FLOAT> > & b) {a = b.real() ;};
  static inline Scl<MP_FLOAT> absolute(const Scl<MP_FLOAT> & a) { return abs(a); }
  };
#endif //

template <class T> struct Real<complex<T> > {
  typedef T TYPE; 
  static inline void assign(complex<T> & a, const complex<T> & b) { a = b ;};
  static inline double absolute(const complex<T> & a) { return abs(a); }
};

template <> struct Real<int> {
  typedef int TYPE; 
  static inline void assign(complex<int> & a, const complex<int> & b) { a = b.real() ;};
  //static inline double absolute(const complex<int> & a) { return (double) abs(a); }
};

template <class T> inline void __div_1(T& a, T& b, const T& c, const T& d)
{
  T r(d);r/=c; 
  T s(d); s*=r; s+=c;
  T x(b); x*=r; x+=a;x/=s;
  a*=r;
  b -= a; b /= s;
  a = x;
}
//---------------------------------------------------------------------------------------------------
/* 
   complex<Scl<MP_FLOAT> > specialization
*/

#ifdef _GMP
template <>
inline complex<Scl<MP_FLOAT> >&  complex<Scl<MP_FLOAT> >::operator /=(const complex<Scl<MP_FLOAT> > &y) {
  __div_1(re,im , y.real(), y.imag());
  return *this;
}

template<> inline complex<Scl<MP_FLOAT> > pow(const complex<Scl<MP_FLOAT> >& z, int k) {
    double mod = sqrt(norm(complex<double>(ToDouble(z.real()),ToDouble(z.imag()))));
    if (k == 0)
      return complex<Scl<MP_FLOAT> > (Scl<MP_FLOAT>(1), Scl<MP_FLOAT>(0));
    else
      if (k == -1) 
    return complex<Scl<MP_FLOAT> >(z.real(),Scl<MP_FLOAT>(-1.0*ToDouble(z.imag())/mod));
      else {
    double argument = k * arg(complex<double>(ToDouble(z.real()),ToDouble(z.imag())));
    mod = pow(mod,k);
    return complex<Scl<MP_FLOAT> > (mod*cos(argument), mod*sin(argument));
      }
}
 
inline complex<Scl<MP_FLOAT> >
operator / (const complex<Scl<MP_FLOAT> >& x, double y) {
  Scl<MP_FLOAT> d(y);
  return complex<Scl<MP_FLOAT> >(x.real()/d,x.imag()/d); }

inline complex<Scl<MP_FLOAT> >
operator / (const complex<Scl<MP_FLOAT> >& x,const complex<Scl<MP_FLOAT> > y) {
  complex<Scl<MP_FLOAT> > res(x);
  res /= y;
  return complex<Scl<MP_FLOAT> >(res.real(),res.imag()); }

ostream & operator<<(ostream & os, const complex<Scl<MP_FLOAT> > & z){
  os<<"("<< z.real()<< ","<<z.imag()<<")";
  return(os);
  }
#endif //
//-----------------------------------------------------------------------------------------------------

template<class T> inline void Root_Unity(complex<T> &c, int n) {
  c = complex<T>(cos(8.0*atan(1.0)/n), sin(8.0*atan(1.0)/n));
}

template<class T> inline complex<T> pow(const complex<T>& z, int k) {
    T mod = sqrt(norm(z));
    if (k == 0) return complex<T> (1, T(0));
    else
      if (k == -1) 
      return complex<T>(z.real(), -(z.imag())/mod);
      else {
    double argument = k * arg(z);
    mod = pow(mod,k);
    return complex<T> (mod*cos(argument), mod*sin(argument));
      }
}

//---------------------------------------------------------------------------------------------------
/*template<class T> 
inline boolean operator < (const complex<T>& a,const complex<T> & b)
{
  return norm(a) < norm(b);
  }*/

#endif // // !INCLUDED_COMPLEXH