RealOf.hpp

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of BORDERBASIX software. *
*      (C) B. Mourrain,  INRIA                                       *
**********************************************************************
History: 
$Id: RealOf.hpp,v 1.1.1.1 2006/10/06 08:01:40 trebuche Exp $
**********************************************************************/
/*
  RealOf.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>
#define COMPLEX complex<double>
//----------------------------------------------------------------------------
struct Real;
//----------------------------------------------------------------------------
#ifdef _SCL_H_
//----------------------------------------------------------------------------
template<class T> struct Scl;
//----------------------------------------------------------------------------
#endif //



struct Warning_Value_Type_Not_Defined {};

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

template <> struct RealOf<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); }
};

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

template <class T> struct RealOf<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); }
};

#if 0//def _SCL_MPF_H_
//extern "C++" 
typedef __mpf_struct MP_FLOAT; //struct MP_FLOAT;
template<class T> struct Scl;


//extern struct Scl<MP_FLOAT>;
/*
template <> struct RealOf<complex<Scl<MP_FLOAT> > > {
  typedef typename 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)); }
    };

template <> struct RealOf<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); }
  };
*/

//----------------------------------------------------------------------------
//template<class T> struct Scl;
//struct MP_FLOAT;

/* 
   complex<Scl<MP_FLOAT> > specialization
*/

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;
}

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 // //_Scl_H_
//-----------------------------------------------------------------------------------------------------

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