synaps/arithm/CC.h

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*   Author(s): B. Mourrain, GALAAD, INRIA                            *
**********************************************************************
History: 
$Id: CC.h,v 1.4 2005/09/28 06:40:42 mourrain Exp $
**********************************************************************/
#ifndef synaps_arithm_CC_H
#define synaps_arithm_CC_H
/********************************************************************/
#include <assert.h>
#include <synaps/init.h>
#include <synaps/arithm/Scl.h>
#include <mpsolve/mpc.h>
#include <iostream>
/********************************************************************/
__BEGIN_NAMESPACE_SYNAPS
//--------------------------------------------------------------------
typedef  __mpc_struct MPC;
//typedef  mpc_t MPC;
typedef Scl<MPC> CC;
//--------------------------------------------------------------------
template<> inline 
void CC::init() {mpc_init(&rep());}
//-------------------------- Scl<T> (...) ----------------------------
template<> inline
Scl<MPC>::Scl(void)
{
  mpc_init(&rep());
}

template<> inline
Scl<MPC>::Scl(const Scl<MPC>& b)
{
  mpc_init_set(&rep(), &b.rep());
}

template<> inline
Scl<MPC>::Scl(const MPC & b)
{
  mpc_init_set(&rep(), &b);
}

template<> inline
Scl<MPC>::Scl(double d)
{
  mpc_init_set_d(&rep(), d, 0);
}

template<> inline
Scl<MPC>::Scl(signed long int sl)
{
  mpc_init_set_si(&rep(), sl, 0);
}

template<> inline
Scl<MPC>::Scl(unsigned long int ul)
{
  mpc_init_set_ui(&rep(), ul, 0);
}

template<> inline
Scl<MPC>::Scl(int si)
{
  mpc_init_set_si(&rep(), (long) si,0);
}

template<> inline
Scl<MPC>::Scl(const char *string,unsigned int base)
{
  if (mpc_init_set_str(&rep(), string,"0", base))
    std::cerr << "Scl<MPC>::Scl<MPC>: The string " << string 
              << " is not a valid number in base " << base << std::endl;
}

//-------------------------- Scl<MPC> =, +=, -=, *=, /=, %= ---------------------

template<> inline 
Scl<MPC>& Scl<MPC>::operator =  (const Scl<MPC>& rhs)
{
  if (this != &rhs)
    mpc_set(&rep(), &rhs.rep()); return *this;
}

template<> inline 
Scl<MPC>& Scl<MPC>::operator += (const Scl<MPC>& rhs)
{
  mpc_add(&rep(), &rep(), &rhs.rep());  return *this;
}

template<> inline 
Scl<MPC>& Scl<MPC>::operator -= (const Scl<MPC>& rhs)
{
  mpc_sub(&rep(), &rep(), &rhs.rep()); return *this;
}

template<> inline 
Scl<MPC>& Scl<MPC>::operator *= (const Scl<MPC>& rhs)
{
  mpc_mul(&rep(), &rep(), &rhs.rep()); return *this;
}

template<> inline 
Scl<MPC>& Scl<MPC>::operator /= (const Scl<MPC>& rhs)
{
  mpc_div(&rep(), &rep(), &rhs.rep()); return *this;
}

//-------------------------- unsigned long =, +=, -=, *=, /=, %= --------------

template<> inline 
Scl<MPC>&  Scl<MPC>::operator =  (unsigned long ul)
{
  mpc_set_ui(&rep(), ul,0); return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator += (unsigned long ul)
{
  mpc_add_ui(&rep(), &rep(), ul,0);return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator -= (unsigned long ul)
{
  mpc_sub_ui(&rep(), &rep(), ul,0);return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator *= (int ul)
{
   if (ul >= 0)
     mpc_mul_ui(&rep(), &rep(), (unsigned long) ul);
   else {
     mpc_mul_ui(&rep(), &rep(), (unsigned long) (-ul));
     mpc_neg(&rep(), &rep());
   }
   return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator *= (unsigned long ul)
{
  mpc_mul_ui(&rep(), &rep(), ul);return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator /= (unsigned long ul)
{
  mpc_div_ui(&rep(), &rep(), ul);return *this;
}

//-------------------------- long =, +=, -=, *=, /=, %= -----------------------
template<> inline Scl<MPC>& 
Scl<MPC>::operator =  (long sl)
{
  mpc_set_si(&rep(), sl,0);return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator += (long sl)
{
  if (sl >= 0)
    mpc_add_ui(&rep(), &rep(), (unsigned long) sl,0);
  else
    mpc_sub_ui(&rep(), &rep(), ((unsigned long) -sl),0);
  return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator -= (long sl)
{
  if (sl >= 0)
    mpc_sub_ui(&rep(), &rep(), (unsigned long) sl,0);
  else
    mpc_add_ui(&rep(), &rep(), (unsigned long) sl,0);
  return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator *= (long sl)
{
  if (sl >= 0)
    mpc_mul_ui(&rep(), &rep(), (unsigned long)sl);
  else
    {
      mpc_mul_ui(&rep(), &rep(), ((unsigned long) -sl));
      mpc_neg(&rep(), &rep());
    }
  return *this;
}

template<> inline Scl<MPC>& 
Scl<MPC>::operator /= (long sl)
{
  if (sl >= 0)
    mpc_div_ui(&rep(), &rep(), (unsigned long) sl);
  else
    {
      mpc_div_ui(&rep(), &rep(), ((unsigned long) -sl));
      mpc_neg(&rep(), &rep());
    }
  return *this;
}

//----------------------------------------------------------------------
template<> inline Scl<MPC>& Scl<MPC>::operator = (int si)
{
  mpc_init_set_si(&rep(), (long) si,0); return *this;
}
//----------------------------------------------------------------------
template<> inline Scl<MPC>& 
Scl<MPC>::operator += (int i)
{
  *this+= long(i);
  return *this;
}

template<> inline 
Scl<MPC>& Scl<MPC>::pow (unsigned long exp)
{
  //  double r(mpc_get_d(&rep())),x(r); 
  Scl<MPC> r(*this);
  for (unsigned long i = 2; i <= exp; i++) mpc_mul(&rep(),&rep(),&r.rep()); // x *=r;
  // mpc_set_d(&rep(),x);
  
  return *this;
}

inline Scl<MPC> operator * (double d, const Scl<MPC>& b)
{
  /*Scl<MPC> r(b);
    mpc_set_d(&r.rep(),mpc_get_d(&b.rep()) * d);*/
  Scl<MPC> r(d);
  mpc_mul(&r.rep(),&r.rep(),&b.rep());
  return Scl<MPC>(r);
}
//----------------------------------------------------------------------
template<> inline bool
Scl<MPC>::operator == (const Scl<MPC>& x) const
{
  return mpc_eq(&rep(), &x.rep(), 0);
}
//----------------------------------------------------------------------
template<> inline bool
Scl<MPC>::operator != (int si) const
{
  if(si==0)
    return mpc_eq_zero(& this->rep());
  else if(si==1)
    return mpc_eq_one(& this->rep());
  else
    {
      Scl<MPC> tmp(si);
      return mpc_eq(&rep(), &tmp.rep(), 0); //xxx
    }
}
//----------------------------------------------------------------------
inline Scl<MPC>
Norm(const Scl<MPC>& x)
{
 if(x>0)
  return Scl<MPC>(x);
 else
  return Scl<MPC>(-x);
}

//======================================================================
// ASSIGNEMENTS  
//======================================================================

inline void add(MPC& r, const MPC& a, const MPC& b) {mpc_add(&r,&a,&b);}
inline void sub(MPC& r, const MPC& a, const MPC& b) {mpc_sub(&r,&a,&b);}
inline void mul(MPC& r, const MPC& a, const MPC& b) {mpc_mul(&r,&a,&b);}
inline void div(MPC& r, const MPC& a, const MPC& b) {mpc_div(&r,&a,&b);}


// namespace let {

//   inline void assign(Scl<MPC> & p, 
//                   const OP<'+',Scl<MPC>, Scl<MPC> >& M) 
//   {
//     mpc_add(&p.rep(),&M.op1.rep(),&M.op2.rep());
//   }
//   //----------------------------------------------------------------------
//   inline
//   void assign(Scl<MPC> & p, 
//            const OP<'-',Scl<MPC>, Scl<MPC> >& M) 
//   { 
//     mpc_sub(&p.rep(),&M.op1.rep(),&M.op2.rep());
//   }
//   //----------------------------------------------------------------------
//   inline void assign(Scl<MPC> & p, 
//            const OP<'*',Scl<MPC>, Scl<MPC> >& M) 
//   { 
//     mpc_mul(&p.rep(),&M.op1.rep(),&M.op2.rep());
//   }
//   //----------------------------------------------------------------------
//   inline void assign(Scl<MPC> & p, 
//            const OP<'/',Scl<MPC>,Scl<MPC> >& M) 
//   { 
//     assert(M.op2!=0);
//     mpc_div(&p.rep(),&M.op1.rep(),&M.op2.rep());
//   }
// }
//======================================================================
// INPUT OUTPUT
//======================================================================
inline std::ostream&  operator << (std::ostream& os, const Scl<MPC>& b)
{
  mpc_out_str(stdout, 10, mpc_get_prec(&b.rep()), &b.rep()); //mpc_get_prec(&b.rep())
  return os;
}
//----------------------------------------------------------------------
inline std::istream&  operator >> (std::istream& is, Scl<MPC>& x)
{
  std::string a,b; is >> a>> b;
  mpc_init_set_str(&x.rep(),a.c_str(),b.c_str(),10);
  return is;
}
//--------------------------------------------------------------------
template<class T>
CC Complex(const T& a, const T& b)
{
  RR r,i;
  let::assign(r,a);
  let::assign(i,b);
  CC z;
  mpc_set_f(&z.rep(),r.get_mpf_t(),i.get_mpf_t());
  return z; 
}

inline RR real(const CC& z)
{
  return RR(z.rep().r); 
} 

inline RR img(const CC& z)
{
  return RR(z.rep().i); 
} 
/********************************************************************/
__END_NAMESPACE_SYNAPS
#endif //synaps_arithm_CC_H