/Users/mourrain/Devel/mmx/numerix/include/numerix/complex.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : complex.hpp
* DESCRIPTION: Complexified numbers
* COPYRIGHT  : (C) 2005  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_COMPLEX_HPP
#define __MMX_COMPLEX_HPP
#include <basix/double.hpp>
#include <basix/function.hpp>
#include <numerix/floating.hpp>
#include <numerix/rounded.hpp>
namespace mmx {
#define TMPL_DEF template<typename C>
#define TMPL template<typename C>
#define Complex complex<C>
TMPL class complex;
TMPL inline C Re (const Complex& z);
TMPL inline C Im (const Complex& z);
TMPL inline C& Re (Complex& z);
TMPL inline C& Im (Complex& z);

/******************************************************************************
* The complex type
******************************************************************************/

TMPL_DEF
class complex {
MMX_ALLOCATORS;
private:
  C re;
  C im;
public:
  inline complex () {}
  inline complex (const format<C>& fm):
    re (get_sample (fm)), im (get_sample (fm)) {}
  template<typename T> inline complex (const T& x):
    re (as<C> (x)), im (0) {}
  template<typename T> inline complex (const complex<T>& z):
    re (Re (z)), im (Im (z)) {}
  template<typename XT,typename YT> inline
  complex (const XT& x, const YT& y):
    re (x), im (y) {}
  friend C Re LESSGTR (const Complex& z);
  friend C Im LESSGTR (const Complex& z);
  friend C& Re LESSGTR (Complex& z);
  friend C& Im LESSGTR (Complex& z);

  inline Complex& operator += (const Complex& z);
  inline Complex& operator -= (const Complex& z);
  inline Complex& operator *= (const C& x);
  inline Complex& operator /= (const C& x);
  inline Complex& operator *= (const Complex& z);
  inline Complex& operator /= (const Complex& z);
  inline Complex& operator <<= (const xint& shift);
  inline Complex& operator >>= (const xint& shift);
};
DEFINE_UNARY_FORMAT_1 (complex)

STYPE_TO_TYPE(TMPL,scalar_type,Complex,C);
UNARY_RETURN_TYPE(TMPL,Re_op,Complex,C);
UNARY_RETURN_TYPE(TMPL,abs_op,Complex,C);
UNARY_RETURN_TYPE(TMPL,center_op,Complex,complex<Center_type(C) >);
UNARY_RETURN_TYPE(TMPL,radius_op,Complex,Radius_type(C));
STYPE_TO_TYPE(TMPL,default_radius_type,Complex,Default_radius_type(C));
TMPL struct rounding_helper<Complex >: public rounding_helper<C> {};
TMPL struct projective_helper<Complex > { static const bool val= true; };

TMPL
struct fast_helper<Complex > {
  typedef complex<Fast_type(C) > fast_type;
  static inline fast_type dd (const Complex& z) {
    return complex<Fast_type (C) > (fast (Re (z)), fast (Im (z))); }
  static inline Complex uu (const fast_type& z) {
    return complex<C> (slow<C> (Re (z)), slow<C> (Im (z))); }
};

template<typename T,typename F>
struct as_helper<complex<T>,complex<F> > {
  static inline complex<T> cv (const complex<F>& z) {
    return complex<T> (as<T> (Re (z)), as<T> (Im (z))); }
};

template<typename T, typename F> inline void
set_as (complex<T>& r, const complex<F>& z) {
  set_as (Re(r), Re(z));
  set_as (Im(r), Im(z));
}

template<typename T, typename F> inline void
set_as (complex<T>& r, const F& x) {
  set_as (Re(r), x);
  set_as (Im(r), 0);
}

/******************************************************************************
* Basic routines
******************************************************************************/

TMPL inline C Re (const Complex& z) { return z.re; }
TMPL inline C Im (const Complex& z) { return z.im; }
TMPL inline C& Re (Complex& z) { return z.re; }
TMPL inline C& Im (Complex& z) { return z.im; }
TMPL inline format<C> CF (const Complex& z) { return get_format (Re(z)); }

TMPL inline Complex copy (const Complex& z) {
  return Complex (copy (Re (z)), copy (Im (z))); }
TMPL inline Complex duplicate (const Complex& z) {
  return Complex (duplicate (Re (z)), duplicate (Im (z))); }
TMPL inline bool is_exact_zero (const Complex& z) {
  return is_exact_zero (Re (z)) && is_exact_zero (Im (z)); }

template<typename Op, typename C> inline nat
unary_hash (const Complex& z) {
  nat h= Op::op (Re (z));
  return (h<<1) ^ (h<<5) ^ (h>>29) ^ Op::op (Im (z));
}

template<typename Op, typename C> bool
binary_test (const Complex& z1, const Complex& z2) {
  return Op::op (Re (z1), Re (z2)) && Op::op (Im (z1), Im (z2));
}

template<typename S1, typename D> complex<D>
map (const function_1<D,Argument(S1) >& fun, const complex<S1>& z) {
  return complex<D> (fun (Re (z)), fun (Im (z)));
}

template<typename S1, typename D> complex<D>
map (D (*fun) (const S1&), const complex<S1>& z) {
  return complex<D> (fun (Re (z)), fun (Im (z)));
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Complex)
EXACT_IDENTITY_OP_SUGAR(TMPL,Complex)
HARD_IDENTITY_OP_SUGAR(TMPL,Complex)
EQUAL_INT_SUGAR(TMPL,Complex);

TMPL syntactic
flatten (const Complex& z) {
  if (is_nan (z)) return Nan (syntactic);
  if (is_fuzz (z)) return Fuzz (syntactic);
  if (is_infinite (z)) return Infinity (syntactic);
  if (is_exact_zero (z)) return 0;
  if (Re (z) != 0 && Im (z) == 0)
    return flatten (Re (z));
  if (Re (z) == 0 && Im (z) != 0)
    return flatten (Im (z)) * Imaginary (syntactic);
  return flatten (Re (z)) + flatten (Im (z)) * Imaginary (syntactic);
}

TMPL syntactic
flatten (const Complex& z, xnat dd) {
  if (is_nan (z)) return Nan (syntactic);
  if (is_fuzz (z)) return Fuzz (syntactic);
  if (is_infinite (z)) return Infinity (syntactic);
  if (is_exact_zero (z)) return 0;
  if (Re (z) != 0 && Im (z) == 0)
    return flatten (Re (z), dd);
  if (Re (z) == 0 && Im (z) != 0)
    return flatten (Im (z), dd) * Imaginary (syntactic);
  return flatten (Re (z), dd) + flatten (Im (z), dd) * Imaginary (syntactic);
}

TMPL
struct binary_helper<Complex >: public void_binary_helper<Complex > {
  static inline string short_type_name () {
    return "Cc" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Complex", Full_type_name (C)); }
  static inline generic disassemble (const Complex& z) {
    return gen_vec (as<generic> (Re (z)), as<generic> (Im (z))); }
  static inline Complex assemble (const generic& v) {
    return Complex (as<C> (vector_access (v, 0)),
                    as<C> (vector_access (v, 1))); }
  static inline void write (const port& out, const Complex& z) {
    binary_write<C> (out, Re (z));
    binary_write<C> (out, Im (z)); }
  static inline Complex read (const port& in) {
    C r= binary_read<C> (in);
    C i= binary_read<C> (in);
    return Complex (r, i); }
};

/******************************************************************************
* Polar coordinates and complex conjugation
******************************************************************************/

TMPL inline Complex
gaussian (const C& re, const C& im= 0) {
  return Complex (re, im);
}

template<typename C, typename R> inline void
gaussian_as (C& z, const R& x, const R& y) {
  z= gaussian (x, y);
}

TMPL inline Complex
polar (const C& r, const C& t) {
  return Complex (r * cos (t), r * sin (t));
}

TMPL inline C
abs (const Complex& z) {
  return hypot (Re (z), Im (z));
}

TMPL inline C
arg (const Complex& z) {
  return atan2 (Im (z), Re (z));
}

TMPL inline C
norm (const Complex& z) {
  return square (Re (z)) + square (Im (z));
}

TMPL inline Complex
conj (const Complex& z) {
  return Complex (Re (z), -Im (z));
}

TMPL inline Complex
times_i (const Complex& z) {
  return Complex (-Im (z), Re (z));
}

TMPL inline Complex
over_i (const Complex& z) {
  return Complex (Im (z), -Re (z));
}

/******************************************************************************
* Constants
******************************************************************************/

TMPL inline void set_nan (Complex& z) {
  set_nan (Re (z)); set_nan (Im (z)); }
TMPL inline void set_infinity (Complex& z) {
  set_infinity (Re (z)); set_infinity (Im (z)); }
TMPL inline void set_fuzz (Complex& z) {
  set_fuzz (Re (z)); set_fuzz (Im (z)); }
TMPL inline void set_smallest (Complex& z) {
  set_smallest (Re (z)); set_zero (Im (z)); }
TMPL inline void set_largest (Complex& z) {
  set_largest (Re (z)); set_zero (Im (z)); }
TMPL inline void set_accuracy (Complex& z) {
  set_accuracy (Re (z)); set_zero (Im (z)); }
TMPL inline void set_log2 (Complex& z) {
  set_log2 (Re (z)); set_zero (Im (z)); }
TMPL inline void set_pi (Complex& z) {
  set_pi (Re (z)); set_zero (Im (z)); }
TMPL inline void set_euler (Complex& z) {
  set_euler (Re (z)); set_zero (Im (z)); }
TMPL inline void set_imaginary (Complex& z) {
  set_zero (Re (z)); set_one (Im (z)); }
TMPL inline Complex times_infinity (const Complex& z) {
  return z * infinity_cst<Complex > (); }

/******************************************************************************
* Basic arithmetic
******************************************************************************/

TMPL inline Complex
operator + (const Complex& z1, const Complex& z2) {
  return Complex (Re (z1) + Re (z2), Im (z1) + Im (z2));
}

TMPL inline Complex
operator + (const Complex& z1, const C& z2) {
  return Complex (Re (z1) + z2, Im (z1));
}

TMPL inline Complex
operator + (const C& z1, const Complex& z2) {
  return Complex (z1 + Re (z2), Im (z2));
}

TMPL inline Complex
operator - (const Complex& z) {
  return Complex (-Re (z), -Im (z));
}

TMPL inline Complex
operator - (const Complex& z1, const Complex& z2) {
  return Complex (Re (z1) - Re (z2), Im (z1) - Im (z2));
}

TMPL inline Complex
operator - (const Complex& z1, const C& z2) {
  return Complex (Re (z1) - z2, Im (z1));
}

TMPL inline Complex
operator - (const C& z1, const Complex& z2) {
  return Complex (z1 - Re (z2), -Im (z2));
}

TMPL inline Complex
operator * (const C& x, const Complex& z) {
  return Complex (x * Re (z), x * Im (z));
}

TMPL inline Complex
operator * (const Complex& z, const C& x) {
  return Complex (Re (z) * x, Im (z) * x);
}

TMPL inline Complex
operator * (const Complex& z1, const Complex& z2) {
  return Complex (Re (z1) * Re (z2) - Im (z1) * Im (z2),
                  Re (z1) * Im (z2) + Im (z1) * Re (z2));
}

TMPL inline Complex
square (const Complex& z) {
  C h= Re (z) * Im (z);
  return Complex (square (Re (z)) - square (Im (z)), h + h);
}

TMPL inline Complex
invert (const Complex& z) {
  C a= norm (z);
  return Complex (Re (z) / a, -Im (z) / a);
}

TMPL inline Complex
operator / (const Complex& z, const C& x) {
  return Complex (Re (z) / x, Im (z) / x);
}

TMPL inline Complex
operator / (const C& x, const Complex& z) {
  C a= x / norm (z);
  return Complex (a * Re (z), -a * Im (z));
}

TMPL inline Complex
operator / (const Complex& z1, const Complex& z2) {
  C a= norm (z2);
  return Complex ((Re (z1) * Re (z2) + Im (z1) * Im (z2)) / a,
                  (Im (z1) * Re (z2) - Re (z1) * Im (z2)) / a);
}

ARITH_SCALAR_INT_SUGAR(TMPL,Complex);
GCD_SUGAR(TMPL,Complex);

/******************************************************************************
* In place arithmetic operations
******************************************************************************/

TMPL inline Complex&
Complex::operator += (const Complex& z) {
  re += z.re;
  im += z.im;
  return *this;
}

TMPL inline Complex&
Complex::operator -= (const Complex& z) {
  re -= z.re;
  im -= z.im;
  return *this;
}

TMPL inline Complex&
Complex::operator *= (const C& x) {
  re *= x;
  im *= x;
  return *this;
}

TMPL inline Complex&
Complex::operator /= (const C& x) {
  re /= x;
  im /= x;
  return *this;
}

TMPL inline Complex&
Complex::operator *= (const Complex& z) {
  C RE= re, zRE= z.re;
  re= re * z.re - im * z.im;
  im= RE * z.im + im * zRE;
  return *this;
}

TMPL inline Complex&
Complex::operator /= (const Complex& z) {
  C RE= re;
  C a = norm (z);
  re= (re * z.re + im * z.im) / a;
  im= (im * z.re - RE * z.im) / a;
  return *this;
}

TMPL inline void
neg (Complex& r) {
  neg (Re (r));
  neg (Im (r));
}

TMPL inline void
neg (Complex& r, const Complex& z1) {
  neg (Re (r), Re (z1));
  neg (Im (r), Im (z1));
}

TMPL inline void
add (Complex& r, const Complex& z1, const Complex& z2) {
  add (Re (r), Re (z1), Re (z2));
  add (Im (r), Im (z1), Im (z2));
}

TMPL inline void
sub (Complex& r, const Complex& z1, const Complex& z2) {
  sub (Re (r), Re (z1), Re (z2));
  sub (Im (r), Im (z1), Im (z2));
}

TMPL inline void
mul (Complex& r, const C& c, const Complex& z) {
  mul (Re (r), c, Re (z));
  mul (Im (r), c, Im (z));
}

TMPL inline void
mul (Complex& r, const Complex& z, const C& c) {
  mul (Re (r), Re (z), c);
  mul (Im (r), Im (z), c);
}

TMPL inline void
div (Complex& r, const Complex& z, const C& c) {
  div (Re (r), Re (z), c);
  div (Im (r), Im (z), c);
}

/******************************************************************************
* Exponentiation and logarithms
******************************************************************************/

TMPL Complex
sqrt (const Complex& z) {
  C a= hypot (Re (z), Im (z));
  if (Re (z) > 0) {
    C r= sqrt ((a + Re (z)) / promote (2, a));
    return Complex (r, (Im (z) / r) / promote (2, a));
  }
  else if (a == 0) return promote (0, z);
  else {
    C i = sqrt ((a - Re (z)) / promote (2, a));
    if (Im (z) < 0) i= -i;
    return Complex ((Im (z) / i) / promote (2, a), i);
  }
}

TMPL inline Complex
exp (const Complex& z) {
  return polar (exp (Re (z)), Im (z));
}

TMPL inline Complex
log (const Complex& z) {
  return Complex (log (abs (z)), arg (z));
}

TMPL inline Complex
pow (const Complex& z, const Complex& y) {
  C r= log (abs (z));
  C t= arg (z);
  return polar (exp (r * Re (y) - Im (y) * t),
                Im (y) * r + Re (y) * t);
}

TMPL inline Complex
pow (const Complex& z, const C& y) {
  return exp (y * log (z));
}

TMPL inline Complex
pow (const C& z, const Complex& y) {
  return exp (y * log (z));
}

TMPL inline Complex
pow (const Complex& z, const int& y) {
  return exp (y * log (z));
}

TMPL inline Complex
pow (const int& z, const Complex& y) {
  return exp (y * log (z));
}

TMPL inline Complex
cos (const Complex& z) {
  return Complex (cos (Re (z)) * cosh (Im (z)),
                  -sin (Re (z)) * sinh (Im (z)));
}

TMPL inline Complex
sin (const Complex& z) {
  return Complex (sin (Re (z)) * cosh (Im (z)),
                  cos (Re (z)) * sinh (Im (z)));
}

TMPL inline Complex
tan (const Complex& z) {
  return sin (z) / cos (z);
}

TMPL inline Complex
acos (const C& z) {
  return over_i (acosh (z));
}

TMPL inline Complex
asin (const C& z) {
  return over_i (asinh (times_i (z)));
}

TMPL inline Complex
atan (const C& z) {
  return over_i (atanh (times_i (z)));
}

TMPL inline Complex
acos (const Complex& z) {
  return over_i (acosh (z));
}

TMPL inline Complex
asin (const Complex& z) {
  return over_i (asinh (times_i (z)));
}

TMPL inline Complex
atan (const Complex& z) {
  return over_i (atanh (times_i (z)));
}

TMPL inline Complex
cosh (const Complex& z) {
  return Complex (cosh (Re (z)) * cos (Im (z)),
                  sinh (Re (z)) * sin (Im (z)));
}

TMPL inline Complex
sinh (const Complex& z) {
  return Complex (sinh (Re (z)) * cos (Im (z)),
                  cosh (Re (z)) * sin (Im (z)));
}

TMPL inline Complex
tanh (const Complex& z) {
  return sinh (z) / cosh (z);
}

INV_TRIGO_SUGAR(TMPL,Complex)
INV_HYPER_SUGAR(TMPL,Complex)
ARG_HYPER_SUGAR(TMPL,Complex)

/******************************************************************************
* Floating point related functions
******************************************************************************/

TMPL inline bool is_finite (const Complex& z) {
  return is_finite (Re (z)) && is_finite (Im (z)); }
TMPL inline bool is_nan (const Complex& z) {
  return is_nan (Re (z)) || is_nan (Im (z)); }
TMPL inline bool is_infinite (const Complex& z) {
  return !is_nan (z) && (is_infinite (Re (z)) || is_infinite (Im (z))); }
TMPL inline bool is_fuzz (const Complex& z) {
  return !is_nan (z) && (is_fuzz (Re (z)) || is_fuzz (Im (z))); }
TMPL inline bool is_reliable (const Complex& z) {
  return is_reliable (Re (z)); }

TMPL inline Complex change_precision (const Complex& z, xnat prec) {
  return Complex (change_precision (Re (z), prec),
                  change_precision (Im (z), prec)); }
TMPL inline xnat precision (const Complex& z) {
  return min (precision (Re (z)), precision (Im (z))); }

TMPL inline Abs_type(C) rounding_error (const complex<C>& z) {
  typedef Round_up(Abs_type(C) ) Up;
  return Up::add (rounding_error (Re (z)), rounding_error (Im (z))); }
TMPL inline Abs_type(C) additive_error (const complex<C>& z) {
  typedef Round_up(Abs_type(C) ) Up;
  return Up::add (additive_error (Re (z)), additive_error (Im (z))); }
TMPL inline Abs_type(C) multiplicative_error (const complex<C>& z) {
  typedef Round_up(Abs_type(C) ) Up;
  return incexp2 (Up::add (multiplicative_error (Re (z)),
                           multiplicative_error (Im (z))), 3); }
TMPL inline Abs_type(C) elementary_error (const complex<C>& z) {
  typedef Round_up(Abs_type(C) ) Up;
  return incexp2 (Up::add (elementary_error (Re (z)),
                           elementary_error (Im (z))), 3); }

/*
template<typename V> inline floating<V>
additive_error (const complex<floating<V> >& z) {
  floating<V> r (*Re(z)), s (*Im(z));
  mpfr_abs (*r, *Re(z), GMP_RNDN);
  mpfr_abs (*s, *Im(z), GMP_RNDN);
  mpfr_add (*r, *s, GMP_RNDU);
  mpfr_mul_2si (*r, *r, 1-mpfr_get_prec (*Re(z)), GMP_RNDN);
  mpfr_nextabove (*r);
  return r; }
*/

TMPL inline xint exponent (const Complex& z) {
  return max (exponent (Re (z)), exponent (Im (z))); }
TMPL inline double magnitude (const Complex& z) {
  return magnitude (square (Re(z)) + square (Im(z))) / 2.0; }
TMPL inline Complex operator << (const Complex& z, const xint& shift) {
  return Complex (incexp2 (Re (z), shift), incexp2 (Im (z), shift)); }
TMPL inline Complex operator >> (const Complex& z, const xint& shift) {
  return Complex (decexp2 (Re (z), shift), decexp2 (Im (z), shift)); }

TMPL inline Complex&
Complex::operator <<= (const xint& shift) {
  incexp2_assign (re, shift);
  incexp2_assign (im, shift);
  return *this;
}

TMPL inline Complex&
Complex::operator >>= (const xint& shift) {
  decexp2_assign (re, shift);
  decexp2_assign (im, shift);
  return *this;
}

TMPL inline Complex sharpen (const Complex& z) {
  return Complex (sharpen (Re (z)), sharpen (Im (z))); }
template<typename C, typename R> inline Complex
blur (const Complex& z, const R& r) {
  return Complex (blur (Re (z), r), blur (Im (z), r)); }

#undef TMPL_DEF
#undef TMPL
#undef Complex
} // namespace mmx
#endif // __MMX_COMPLEX_HPP