[Mplapack-devel] Operator<< for MPFR complex

[Fletcher, J. P <j.p...@as...>]

/*
 * Copyright (c) 2010
 *	Nakata, Maho
 * 	All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 */
/*
Complex class declare for the MPFR
*/

#ifndef __MP_COMPLEX_H__
#define __MP_COMPLEX_H__

#include "mpreal.h"
#include <complex>
#include "mpc.h"

#if defined ___MPACK_BUILD_WITH_GMP___
#include "gmpxx.h"
#include "mpc_class.h"
#endif
#if defined ___MPACK_BUILD_WITH_QD___
#include "qd_complex.h"
#endif
#if defined ___MPACK_BUILD_WITH_DD___
#include "dd_complex.h"
#endif

namespace mpfr {
class mpcomplex {
  private:
   mpc_t mpc;

  public:
  static mpc_rnd_t default_rnd;
  static mp_prec_t default_real_prec;
  static mp_prec_t default_imag_prec;
  static int default_base;
  static int double_bits;

//constructor & deconstructor
    mpcomplex();
    mpcomplex(const mpc_t a);
    mpcomplex(const mpfr_t a, const mpfr_t b);
    mpcomplex(const mpf_t a, const mpf_t b);
    mpcomplex(const char *s, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const mpcomplex& a);
    mpcomplex(const std::complex<double>& a, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const mpreal& a, const mpreal& b);
    mpcomplex(const double& a, const double& b, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const char *s, const char *t, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);

    mpcomplex(const mpreal & a);
    mpcomplex(const mpfr_t a);
    mpcomplex(const mpf_t a);
    mpcomplex(const double a, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    ~mpcomplex(); 

    mpcomplex& operator=(const mpcomplex& a);
    mpcomplex& operator=(const mpc_t a);
    mpcomplex& operator=(const std::complex<double> a);
    mpcomplex& operator=(const char* s);

    mpcomplex& operator=(const mpreal& a);
    mpcomplex& operator=(const mpfr_t a);
    mpcomplex& operator=(const double a);

    //+
    mpcomplex& operator+=(const mpcomplex& a);
    mpcomplex& operator+=(const mpc_t a);
    mpcomplex& operator+=(const std::complex<double> a);

    mpcomplex& operator+=(const mpreal& a);
    mpcomplex& operator+=(const mpfr_t a);
    mpcomplex& operator+=(const double a);
    const mpcomplex operator+() const;
    mpcomplex& operator++();
    const mpcomplex operator ++ (int);

    //-
    mpcomplex& operator-=(const mpcomplex& a);
    mpcomplex& operator-=(const mpc_t a);
    mpcomplex& operator-=(const std::complex<double> a);

    mpcomplex& operator-=(const mpreal& a);
    mpcomplex& operator-=(const mpfr_t a);
    mpcomplex& operator-=(const double a);
    const mpcomplex operator-() const;
    mpcomplex& operator--();
    const mpcomplex operator -- (int);

    //*
    mpcomplex& operator*=(const mpcomplex& a);
    mpcomplex& operator*=(const mpc_t a);
    mpcomplex& operator*=(const std::complex<double> a);
    mpcomplex& operator*=(const mpreal& a);
    mpcomplex& operator*=(const mpfr_t a);
    mpcomplex& operator*=(const double a);

    // /
    mpcomplex& operator/=(const mpcomplex& a);
    mpcomplex& operator/=(const mpc_t a);
    mpcomplex& operator/=(const std::complex<double> a);
    mpcomplex& operator/=(const mpreal& a);
    mpcomplex& operator/=(const mpfr_t a);
    mpcomplex& operator/=(const double a);

    // Output/ Input
    friend std::ostream& operator<<(std::ostream& os, const mpcomplex& v);
    //    friend std::istream& operator>>(std::istream& is, mpcomplex& v);

    //comparison
    friend bool operator == (const mpcomplex& a, const mpcomplex& b);
    friend bool operator == (const mpcomplex& a, const mpreal& b);
    friend bool operator == (const mpreal& a, const mpcomplex& b);
    friend bool operator != (const mpcomplex& a, const mpcomplex& b);
    friend bool operator != (const mpreal& a, const mpcomplex& b);
    friend bool operator != (const mpcomplex& a, const mpreal& b);

    //random
    friend const mpcomplex urandom_c (gmp_randstate_t& state);

//extraction of real and imaginary parts
    const mpreal real() const
    {
        return mpreal(mpc_realref(mpc));
    }
    const mpreal imag() const
    {
        return mpreal(mpc_imagref(mpc));
    }
    mpreal real()
    {
        return mpreal(mpc_realref(mpc));
    }
    mpreal imag()
    {
        return mpreal(mpc_imagref(mpc));
    }
    //other functions
    friend const mpreal abs(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpreal norm(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpcomplex conj(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpreal arg(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpcomplex proj(const mpcomplex& a, mpc_rnd_t mode = mpreal::default_rnd);

    //Powerfunction and Logarithm    
    friend const mpcomplex sqrt(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex pow(const mpcomplex& a, const mpcomplex& b,mpc_rnd_t mode = default_rnd);
    friend const mpcomplex exp(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex log(const mpcomplex& a, mpc_rnd_t mode = default_rnd);

    //trigonometric functions
    friend const mpcomplex sin(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex cos(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex tan(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex sinh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex cosh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex tanh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex asin(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex acos(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex atan(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex asinh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex acosh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex atanh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);

// Type Conversion operators
    operator std::complex<double>() const;
    operator std::string() const;
    std::string to_string(size_t n = 0, int b = default_base, mpc_rnd_t mode = default_rnd) const;

// Set/Get instance properties
   mp_prec_t get_prec() const;
   mp_prec_t get_prec_re() const;
   mp_prec_t get_prec_im() const;
   void set_prec(mp_prec_t prec, mpc_rnd_t rnd_mode = default_rnd);
   void set_prec2(mp_prec_t pr, mp_prec_t pi, mpc_rnd_t rnd_mode = default_rnd);

#if defined ___MPACK_BUILD_WITH_GMP___
mpcomplex(const mpc_class& a);
operator mpc_class() const;
#endif
#if defined ___MPACK_BUILD_WITH_QD___
mpcomplex(const qd_complex& a);
operator qd_complex() const;
mpcomplex& operator=(const qd_complex& a);
#endif
#if defined ___MPACK_BUILD_WITH_DD___
mpcomplex(const dd_complex& a);
operator dd_complex() const;
mpcomplex& operator=(const dd_complex& a);
#endif
};

//+ addition
const mpcomplex operator+(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator+(const mpcomplex& a, const std::complex<double> b);
const mpcomplex operator+(const mpcomplex& a, const char* b);

const mpcomplex operator+(const std::complex<double> a, const mpcomplex& b);
const mpcomplex operator+(const char* a, const mpcomplex& b);

const mpcomplex operator+(const mpcomplex& a, const mpreal b);
const mpcomplex operator+(const mpcomplex& a, const double b);
const mpcomplex operator+(const mpcomplex& a, const int b);

const mpcomplex operator+(const mpreal a, const mpcomplex& b);
const mpcomplex operator+(const double a, const mpcomplex& b);
const mpcomplex operator+(const int a, const mpcomplex& b);
//- subtraction
const mpcomplex operator-(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator-(const mpcomplex& a, const std::complex<double>& b);
const mpcomplex operator-(const mpcomplex& a, const char* b);

const mpcomplex operator-(const std::complex<double>& a, const mpcomplex& b);
const mpcomplex operator-(const char* a, const mpcomplex& b);

const mpcomplex operator-(const mpcomplex& a, const mpreal b);
const mpcomplex operator-(const mpcomplex& a, const double b);
const mpcomplex operator-(const mpcomplex& a, const int b);

const mpcomplex operator-(const mpreal a, const mpcomplex& b);
const mpcomplex operator-(const double a, const mpcomplex& b);
const mpcomplex operator-(const int a, const mpcomplex& b);

//* multiplication
const mpcomplex operator*(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator*(const mpcomplex& a, const mpreal& b);
const mpcomplex operator*(const mpreal& a, const mpcomplex& b);

/// division
const mpcomplex operator/(const mpcomplex& a, const mpcomplex& b);

inline mp_prec_t mpcomplex::get_prec() const
{
	return mpc_get_prec(mpc);
}

inline mp_prec_t mpcomplex::get_prec_re() const
{
      mp_prec_t pr, pi;
      mpc_get_prec2(&pr, &pi, mpc);
      return pr;
}

inline mp_prec_t mpcomplex::get_prec_im() const
{
      mp_prec_t pr, pi;
      mpc_get_prec2(&pr, &pi, mpc);
      return pi;
}

inline void mpcomplex::set_prec(mp_prec_t prec, mpc_rnd_t rnd_mode)
{
    mpcomplex tmp (*this);
    mpc_init2(mpc, prec);
    mpc_set(mpc, tmp.mpc, default_rnd);
}

inline void mpcomplex::set_prec2(mp_prec_t pr, mp_prec_t pi, mpc_rnd_t rnd_mode)
{
    mpcomplex tmp (*this);
    mpc_init3(mpc, pr, pi);
    mpc_set(mpc, tmp.mpc, default_rnd);
}

//constructor and deconstructor
inline mpcomplex::mpcomplex()
{
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_ui(mpc,0UL,default_rnd);
}

inline mpcomplex::mpcomplex(const mpc_t a) 
{
   mp_prec_t pr, pi;
   mpc_get_prec2(&pr, &pi, a);
   mpc_init3(mpc, pr, pi);
   mpc_set(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const mpfr_t a, const mpfr_t b)
{
   mp_prec_t pr, pi;
   pr = mpfr_get_prec(a);  pi = mpfr_get_prec(b);
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, a, b, default_rnd);
}

inline mpcomplex::mpcomplex(const mpf_t a, const mpf_t b)
{
   mp_prec_t pr, pi;
   pr = mpf_get_prec(a);  pi = mpf_get_prec(b);
   mpc_init3(mpc, pr, pi);
   mpc_set_f_f(mpc, a, b, default_rnd);
}

inline mpcomplex::mpcomplex(const char *s, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
    mpc_init3(mpc, pr, pi);
    mpc_set_str(mpc, (char *)s, default_base, mode); 
}

inline mpcomplex::mpcomplex(const mpcomplex& a) 
{
   mp_prec_t pr, pi;
   mpc_get_prec2(&pr, &pi, a.mpc);
   mpc_init3(mpc, pr, pi);
   mpc_set(mpc, a.mpc, default_rnd);
}

inline mpcomplex::mpcomplex(const std::complex<double>& a, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, pr, pi);
   mpc_set_d_d(mpc, a.real(), a.imag(), default_rnd);
}

inline mpcomplex::mpcomplex(const mpreal& a, const mpreal& b)
{
   mp_prec_t pr, pi;
   mpreal tmp1(a), tmp2(b);
   pr = a.get_prec();  pi = b.get_prec();
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, (mpfr_ptr)tmp1, (mpfr_ptr)tmp2, default_rnd);
}

inline mpcomplex::mpcomplex(const double& a, const double& b, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_d_d(mpc, a, b,default_rnd);
}

inline mpcomplex::mpcomplex(const char *s, const char *t, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpfr_t a, b;
   mpfr_init2(a, pr); mpfr_init2(b, pi);
   mpfr_set_str(a, s, default_base, MPC_RND_RE(mode));
   mpfr_set_str(b, t, default_base, MPC_RND_IM(mode));
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, a, b, default_rnd);
   mpfr_clear(a); mpfr_clear(b);
}

inline mpcomplex::mpcomplex(const mpreal& a)
{
   mpreal tmp(a);
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_fr(mpc, (mpfr_ptr)(tmp), default_rnd);
}

inline mpcomplex::mpcomplex(const mpfr_t a)
{
   mp_prec_t pr;
   pr = mpfr_get_prec(a);
   mpc_init2(mpc, pr);
   mpc_set_fr(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const mpf_t a)
{
   mp_prec_t pr;
   pr = mpf_get_prec(a);
   mpc_init2(mpc, pr);
   mpc_set_f(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const double a, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, pr, pi);
   mpc_set_d(mpc, a, default_rnd);
}

inline mpcomplex::~mpcomplex() 
{ 
   mpc_clear(mpc);
}                           

inline mpcomplex& mpcomplex::operator=(const mpcomplex& a)
{
	if (this!= &a) mpc_set(mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpc_t a)
{
	mpc_set(mpc,a,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator=(const std::complex<double> a)
{
    mpc_set_d_d(mpc, a.real(), a.imag(), default_rnd);
    return *this;
}

inline mpcomplex& mpcomplex::operator=(const char* s)
{
     mpc_init3(mpc, default_real_prec, default_imag_prec);
     mpc_set_str(mpc, (char *)s,default_base,default_rnd);				
     return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpreal& a)
{
   mpc_set_fr(mpc, (mpfr_ptr)(&a), default_rnd);
   return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpfr_t a)
{
   mpc_set_fr(mpc, a, default_rnd);
   return *this;
}

inline mpcomplex& mpcomplex::operator=(const double a)
{
   mpc_set_d(mpc, a, default_rnd);
   return *this;
}
// + Addition
inline mpcomplex& mpcomplex::operator+=(const mpcomplex& a)
{
	mpc_add(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator+=(const mpc_t a)
{
	*this += mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator+=(const std::complex<double> a)
{
	return *this += mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator+=(const mpreal& a)
{
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator+=(const mpfr_t a)
{
     mpc_add_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator+=(const double a)
{
     mpreal p = a;
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex mpcomplex::operator+()const
{
	return mpcomplex(*this);
}

inline mpcomplex& mpcomplex::operator++() 
{
	*this += 1.0;
	return *this;
}

inline const mpcomplex mpcomplex::operator++ (int)
{
	mpcomplex x(*this);
	*this += 1.0;
	return x;
}

// + Addition again
inline const mpcomplex operator+(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) += b;
    else return mpcomplex(b) += a;
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp += b;
  }
}

inline const mpcomplex operator+(const mpcomplex& a, const std::complex<double> b)
{
  return mpcomplex(a) += b;
}

inline const mpcomplex operator+(const mpcomplex& a, const char *b)
{
  return mpcomplex(b) += a;
}

inline const mpcomplex operator+(const std::complex<double> a, const mpcomplex& b)
{
  return mpcomplex(b) += a;
}

inline const mpcomplex operator+(const char *a, const mpcomplex& b)
{
  return mpcomplex(a) += b;
}


inline const mpcomplex operator+(const mpcomplex& a, const mpreal b)
{
  mpcomplex tmp(a);
  tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
  return tmp += b;
}

inline const mpcomplex operator+(const mpreal a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  tmp.set_prec2(std::max(b.get_prec_re(), a.get_prec()), b.get_prec_im());
  return tmp += a;
}

inline const mpcomplex operator+(const double a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  return tmp += a;
}

inline const mpcomplex operator+(const int a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  return tmp += a;
}

// - Subtraction
inline mpcomplex& mpcomplex::operator-=(const mpcomplex& a)
{
        mpc_sub(mpc,mpc,a.mpc,default_rnd);
        return *this;
}

inline mpcomplex& mpcomplex::operator-=(const mpc_t a)
{
	*this -= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator-=(const std::complex<double> a)
{
	return *this -= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator-=(const mpreal& a)
{
     mpc_sub_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator-=(const mpfr_t a)
{
     mpc_sub_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator-=(const double a)
{
     mpreal p = a;
     mpc_sub_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex mpcomplex::operator-()const
{
     mpcomplex u(*this);
     mpc_neg (u.mpc, u.mpc, default_rnd);
     return u;
}

inline mpcomplex& mpcomplex::operator--() 
{
	*this -= 1.0;
	return *this;
}

inline const mpcomplex mpcomplex::operator-- (int)
{
	mpcomplex x(*this);
	*this -= 1.0;
	return x;
}

inline const mpcomplex operator-(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) -= b;
    else return -(mpcomplex(b) -= a);
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp -= b;
  }
}

inline const mpcomplex operator-(const mpcomplex& a, const std::complex<double>& b)
{
    return -(mpcomplex(b) -= a);
}

inline const mpcomplex operator-(const mpcomplex& a, const char* b)
{
       return a-mpcomplex(b);
}

inline const mpcomplex operator-(const std::complex<double>& a, const mpcomplex& b)
{
    return mpcomplex(a) -= b;
}

inline const mpcomplex operator-(const char* a, const mpcomplex& b)
{
        return mpcomplex(a)-b;
}

inline const mpcomplex operator-(const mpcomplex& a, const mpreal b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) -= b;
    else return -(mpcomplex(b) -= a);
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
    return tmp -= b;
  }
}

inline const mpcomplex operator-(const mpcomplex& a, const double b)
{
        return mpcomplex(a)-mpreal(b);
}

inline const mpcomplex operator-(const mpcomplex& a, const int b)
{
        return mpcomplex(a)-mpreal(b);
}

inline const mpcomplex operator-(const mpreal a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

inline const mpcomplex operator-(const double a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

inline const mpcomplex operator-(const int a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

// * multiplication
inline mpcomplex& mpcomplex::operator*=(const mpcomplex& a)
{
	mpc_mul(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator*=(const mpc_t a)
{
	*this *= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator*=(const std::complex<double> a)
{
	return *this *= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator*=(const mpreal& a)
{
     mpc_mul_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator*=(const mpfr_t a)
{
     mpc_mul_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator*=(const double a)
{
     mpreal p(a);
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex operator*(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) *= b;
    else return mpcomplex(b) *= a;
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp *= b;
  }
}

inline const mpcomplex operator*(const mpcomplex& a, const mpreal& b)
{
   mpcomplex tmp(a);
   tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
   return tmp *= mpcomplex(b);
}

inline const mpcomplex operator*(const mpreal& a, const mpcomplex& b)
{
   mpcomplex tmp(a);
   tmp.set_prec2(std::max(b.get_prec_re(), a.get_prec()), b.get_prec_im());
   return tmp *= b;
}

// / division
inline mpcomplex& mpcomplex::operator/=(const mpcomplex& a)
{
	mpc_div(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator/=(const mpc_t a)
{
	*this /= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator/=(const std::complex<double> a)
{
	return *this /= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator/=(const mpreal& a)
{
     mpc_div_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator/=(const mpfr_t a)
{
     mpc_div_fr(mpc, mpc, a, default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator/=(const double a)
{
     mpreal p(a);
     mpc_div_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline bool operator == (const mpcomplex& a, const mpcomplex& b)
{
	return (mpc_cmp(a.mpc,b.mpc)==0);
}

inline bool operator == (const mpcomplex& a, const mpreal& b)
{
        mpcomplex c(b);
	return (mpc_cmp(a.mpc, c.mpc)==0);
}

inline bool operator == (const mpreal& a, const mpcomplex& b)
{
        mpcomplex c(a);
	return (mpc_cmp(c.mpc, b.mpc)==0);
}

inline bool operator != (const mpcomplex& a, const mpcomplex& b)
{
	return (!mpc_cmp(a.mpc, b.mpc)==0);
}

inline bool operator != (const mpcomplex& a, const mpreal& b)
{
        mpcomplex c(b);
	return (!mpc_cmp(a.mpc, c.mpc)==0);
}

inline bool operator != (const mpreal& a, const mpcomplex& b)
{
        mpcomplex c(a);
	return (!mpc_cmp(b.mpc, c.mpc)==0);
}

inline mpcomplex::operator std::complex<double>() const
{
  mpreal re, im;
  std::complex<double> tmp;
  re = (*this).real();
  im = (*this).imag();
  tmp.real() = mpfr_get_d(re,MPC_RND_RE(default_rnd));
  tmp.imag() = mpfr_get_d(im,MPC_RND_IM(default_rnd));
  return tmp;
}

inline const mpcomplex operator/(const mpcomplex& a, const mpcomplex& b)
{
  mpcomplex tmp(a);
  tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
  return tmp /= b;
}

inline const mpreal abs(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x;
        mpcomplex y(a);
	mpc_abs((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpreal norm(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x;
        mpcomplex y(a);
	mpc_norm((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpcomplex conj(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
        mpcomplex y(a);
	mpc_conj(y.mpc, y.mpc, rnd_mode);
	return y;
}

inline const mpreal arg(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x(a.real());
        mpcomplex y(a);
	mpc_arg((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpcomplex proj(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
        mpcomplex y(a);
	mpc_proj(y.mpc, y.mpc, rnd_mode);
	return y;
}

inline const mpcomplex sqrt(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sqrt(x.mpc, x.mpc, rnd_mode);
	return x;
}

inline const mpcomplex pow(const mpcomplex& a, const mpcomplex& b, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_pow(x.mpc, x.mpc, b.mpc, rnd_mode);
	return x;
}

inline const mpcomplex exp(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_exp(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex log(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_log(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex sin(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sin(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex cos(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_cos(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex tan(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_tan(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex sinh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sinh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex cosh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_cosh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex tanh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_tanh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex asin(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_asin(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex acos(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_acos(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex atan(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_atan(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex asinh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_asinh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex acosh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_acosh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex atanh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_atanh(x.mpc, a.mpc, rnd_mode);
	return x;
}
inline const mpcomplex urandom_c (gmp_randstate_t& state)
{
	mpcomplex x;
	mpc_urandom(x.mpc,state);
	return x;
}

#if defined ___MPACK_BUILD_WITH_GMP___
inline mpcomplex::mpcomplex(const mpc_class& a)
{
   mp_prec_t pr, pi;
   pr = mpf_get_prec(a.real().get_mpf_t());  pi = mpf_get_prec(a.imag().get_mpf_t());
   mpc_init3(mpc, pr, pi);
   mpc_set_f_f(mpc, a.real().get_mpf_t(), a.imag().get_mpf_t(), default_rnd);
}

inline const mpcomplex operator-(const mpcomplex& a, const mpc_class& b)
{
   mpcomplex tmp(b);
   return mpcomplex(a) -= b;
}

inline const mpcomplex operator-(const mpc_class& a, const mpcomplex& b)
{
   mpcomplex tmp(a);
   mpcomplex tmp2(b);
   return -(mpcomplex(a) -= b);
}

inline mpc_class cast2mpc_class(const mpcomplex &b)
{
//mpcomplex -> mpfr, mpfr -> mpf, mpf -> mpc_class
  mpf_t  tmpre, tmpim;
  mpfr_t tmpre2, tmpim2;
  mpcomplex a(b);
  mp_prec_t pr, pi;

  pr = a.get_prec_re();
  pi = a.get_prec_im();
  
  mpfr_init2(tmpre2, pr); mpfr_init2(tmpim2, pi);
  mpfr_set(tmpre2, (mpfr_ptr)a.real(), MPC_RND_RE(mpcomplex::default_rnd));
  mpfr_set(tmpim2, (mpfr_ptr)a.imag(), MPC_RND_IM(mpcomplex::default_rnd));

  mpf_init2 (tmpre, pr); mpf_init2 (tmpim, pr);
  mpfr_get_f(tmpre, tmpre2, MPC_RND_RE(mpcomplex::default_rnd));  mpfr_get_f(tmpim, tmpim2, MPC_RND_IM(mpcomplex::default_rnd));

  mpc_class tmp(tmpre, tmpim);

  mpf_clear(tmpre);
  mpf_clear(tmpim);

  mpfr_clear(tmpre2);
  mpfr_clear(tmpim2);
  return tmp;
}     

#endif

#if defined ___MPACK_BUILD_WITH_QD___
inline mpcomplex::mpcomplex(const qd_complex& a)
{
   mpfr_t mp_real, mp_imag;
   mpc_init3(mpc, default_real_prec, default_imag_prec);

   mpfr_init2(mp_real, default_real_prec);
   mpfr_set_d (mp_real, a.real().x[0],  MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[1], MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[2], MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[3], MPC_RND_RE(default_rnd));

   mpfr_init2(mp_imag, default_imag_prec);
   mpfr_set_d (mp_imag, a.imag().x[0], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[1], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[2], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[3], MPC_RND_IM(default_rnd));
   
   mpc_set_fr_fr(mpc, mp_real, mp_imag, default_rnd);
   mpfr_clear(mp_imag);
   mpfr_clear(mp_real);
}
inline qd_complex cast2qd_complex(const mpcomplex &b)
{
  std::complex<double> p[4];
  mpcomplex c (b);
  qd_complex q;

  p[0] = c;
  c = c - p[0];
  p[1] = c;
  c = c - p[1];
  p[2] = c;
  c = c - p[2];
  p[3] = c;
//do not optimize
  q = p[0];
  q = q + p[1];
  q = q + p[2];
  q = q + p[3];
  return q;
}     
inline const mpcomplex operator-(const mpcomplex& a, const qd_complex& b)
{
   return mpcomplex(b) -= a;
}

inline const mpcomplex operator-(const qd_complex& a, const mpcomplex& b)
{
   return -(mpcomplex(a) -= b);
}

inline mpcomplex& mpcomplex::operator=(const qd_complex& a)
{
   mpcomplex tmp(a);
   *this = tmp; 
   return *this;
}

#endif

#if defined ___MPACK_BUILD_WITH_DD___
inline mpcomplex::mpcomplex(const dd_complex& a)
{
   mpfr_t mp_real, mp_imag;
   mpc_init3(mpc, default_real_prec, default_imag_prec);

   mpfr_init2(mp_real, default_real_prec);
   mpfr_set_d (mp_real, a.real().x[0],  MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[1], MPC_RND_RE(default_rnd));

   mpfr_init2(mp_imag, default_imag_prec);
   mpfr_set_d (mp_imag, a.imag().x[0], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[1], MPC_RND_IM(default_rnd));
   
   mpc_set_fr_fr(mpc, mp_real, mp_imag, default_rnd);
   mpfr_clear(mp_imag);
   mpfr_clear(mp_real);
}
inline dd_complex cast2dd_complex(const mpcomplex &b)
{
  std::complex<double> p[2];
  mpcomplex c (b);
  dd_complex q;

  p[0] = c;
  c = c - p[0];
  p[1] = c;
  c = c - p[1];
//do not optimize
  q = p[0];
  q = q + p[1];
  return q;
}     
inline const mpcomplex operator-(const mpcomplex& a, const dd_complex& b)
{
   return mpcomplex(b) -= a;
}

inline const mpcomplex operator-(const dd_complex& a, const mpcomplex& b)
{
   return -(mpcomplex(a) -= b);
}

inline mpcomplex& mpcomplex::operator=(const dd_complex& a)
{
   mpcomplex tmp(a);
   *this = tmp; 
   return *this;
}

#endif

}

#endif /* __MP_COMPLEX_H__ */