|| ||           || Doxygen's  

---

HxComplex.h

/*
 *  Copyright (c) 2001, University of Amsterdam, The Netherlands.
 *  All rights reserved.
 *
 *
 *  Author(s):
 *  Jan-Mark Geusebroek (mark@science.uva.nl)
 *  Dennis Koelma (koelma@wins.uva.nl)
 *  Edo Poll (poll@wins.uva.nl)
 */

#ifndef HxComplex_h
#define HxComplex_h

//#include <stdlib>

#include "HxIoFwd.h"
#include "HxMath.h"
#include "HxString.h"

#undef min
#undef max

class HxScalarInt;
#include  "HxScalarDouble.h"
class HxVec2Int;
class HxVec2Double;
class HxVec3Int;
class HxVec3Double;


class L_HXBASIS HxComplex
{
public:


                            HxComplex();

                            HxComplex(double re, double im);
                            
                            HxComplex(const HxComplex& rhs);

    void*                   operator new(size_t, void * = 0);

    HxComplex&              operator=(const HxComplex& rhs);


    int                     dim() const;
    
    double                  x() const;
    
    double                  y() const;

    double                  getValue(int dimension) const; // dimension: 1 or 2
    void                    setValue(int dimension, double value); // dimension: 1 or 2

/** \name Conversion*/
                            operator HxScalarInt() const;

                            operator HxScalarDouble() const;

                            operator HxVec2Int() const;

                            operator HxVec2Double() const;

                            operator HxVec3Int() const;

                            operator HxVec3Double() const;


    int                     operator==(const HxComplex& v) const;

    int                     operator!=(const HxComplex& v) const;

    int                     operator< (const HxComplex& v) const;

    int                     operator<=(const HxComplex& v) const;

    int                     operator> (const HxComplex& v) const;

    int                     operator>=(const HxComplex& v) const;

    static const HxComplex SMALL_VAL;

    static const HxComplex LARGE_VAL;


    HxComplex               operator-() const;

    HxComplex               complement() const;

    HxComplex               conjugate() const;

    HxScalarDouble          abs() const;

    HxScalarDouble          arg() const;

    HxComplex               ceil() const;

    HxComplex               floor() const;

    HxComplex               round() const;

    HxComplex               sum() const;
    
    HxComplex               product() const;

    HxScalarDouble          min() const;

    HxScalarDouble          max() const;

    HxScalarDouble          norm1() const;

    HxScalarDouble          norm2() const;

    HxScalarDouble          normInf() const;

    HxComplex               sqrt() const;

    HxComplex               sin() const;

    HxComplex               cos() const;

    HxComplex               tan() const;

    HxComplex               asin() const;

    HxComplex               acos() const;

    HxComplex               atan() const;

    HxComplex               atan2() const;

    HxComplex               sinh() const;

    HxComplex               cosh() const;

    HxComplex               tanh() const;

    HxComplex               exp() const;

    HxComplex               log() const;

    HxComplex               log10() const;


    HxComplex&              operator+=(const HxComplex& v);

    HxComplex&              operator-=(const HxComplex& v);

    HxComplex&              operator*=(const HxComplex& v);

    HxComplex&              operator/=(const HxComplex& v);


    friend HxComplex        operator+(const HxComplex& v1,
                                      const HxComplex& v2);

    friend HxComplex        operator-(const HxComplex& v1,
                                      const HxComplex& v2);

    friend HxComplex        operator*(const HxComplex& v1,
                                      const HxComplex& v2);

    friend HxComplex        operator/(const HxComplex& v1,
                                      const HxComplex& v2);


    HxComplex               min(const HxComplex& v) const;

    HxComplex&              minAssign(const HxComplex& v);

    HxComplex               max(const HxComplex& v) const;

    HxComplex&              maxAssign(const HxComplex& v);

    HxComplex               inf(const HxComplex& v) const;

    HxComplex&              infAssign(const HxComplex& v);

    HxComplex               sup(const HxComplex& v) const;

    HxComplex&              supAssign(const HxComplex& v);

    HxComplex               pow(const HxComplex& v) const;

    HxComplex               mod(const HxComplex& v) const;

    HxComplex               and(const HxComplex& v) const;

    HxComplex               or(const HxComplex& v) const;

    HxComplex               xor(const HxComplex& v) const;

    HxComplex               leftShift(const HxComplex& v) const;

    HxComplex               rightShift(const HxComplex& v) const;

    HxScalarDouble          dot(const HxComplex& v) const;

    HxComplex               cross(const HxComplex& v) const;


    STD_OSTREAM&            put(STD_OSTREAM& os) const;

    HxString                toString() const;

private:
    double                  _values[2];
};

typedef HxComplex (*HxUpoComplex)(const HxComplex& a);
typedef HxComplex (*HxBpoComplex)(const HxComplex& a1, const HxComplex& a2);

inline STD_OSTREAM&
operator<<(STD_OSTREAM& os, const HxComplex v)
{
    return v.put(os);
}

inline
HxComplex::HxComplex() 
{ 
}

inline
HxComplex::HxComplex(double re, double im)
{
    _values[0] = re;
    _values[1] = im;
}

inline
HxComplex::HxComplex(const HxComplex& v)
{
    _values[0] = v._values[0];
    _values[1] = v._values[1];
}

inline HxComplex&
HxComplex::operator=(const HxComplex& rhs)
{
    _values[0] = rhs._values[0];
    _values[1] = rhs._values[1];
    return *this;
}

inline void*
HxComplex::operator new(size_t size, void *m) 
{
    return m ? m : new char[size];
}

inline int
HxComplex::dim() const 
{ 
    return 2;
}

inline double
HxComplex::x() const 
{ 
    return _values[0];
}

inline double
HxComplex::y() const 
{ 
    return _values[1]; 
}

inline double
HxComplex::getValue(int dim) const
{ 
    return _values[dim - 1]; 
}

inline void
HxComplex::setValue(int dim, double val)
{ 
    _values[dim - 1] = val; 
}

inline int
HxComplex::operator==(const HxComplex& v) const
{
    return (_values[0] == v._values[0]) && (_values[1] == v._values[1]);
}

inline int
HxComplex::operator!=(const HxComplex& v) const
{
    return (_values[0] != v._values[0]) || (_values[1] != v._values[1]);
}

inline int
HxComplex::operator<(const HxComplex& v) const
{
    return (fabs(_values[0]) + fabs(_values[1])) <
           (fabs(v._values[0]) + fabs(v._values[1]));
}

inline int
HxComplex::operator<=(const HxComplex& v) const
{
    return (fabs(_values[0]) + fabs(_values[1])) <=
           (fabs(v._values[0]) + fabs(v._values[1]));
}

inline int
HxComplex::operator>(const HxComplex& v) const
{
    return (fabs(_values[0]) + fabs(_values[1])) >
           (fabs(v._values[0]) + fabs(v._values[1]));
}

inline int
HxComplex::operator>=(const HxComplex& v) const
{
    return (fabs(_values[0]) + fabs(_values[1])) >=
           (fabs(v._values[0]) + fabs(v._values[1]));
}

inline HxComplex
HxComplex::operator-() const
{
    return HxComplex(-_values[0], -_values[1]);
}

inline HxComplex
HxComplex::complement() const
{
    return HxComplex(-_values[0], -_values[1]);
}

inline HxComplex
HxComplex::conjugate() const
{
    return HxComplex(_values[0], -_values[1]);
}

inline HxScalarDouble
HxComplex::abs() const
{
    return HxScalarDouble(::sqrt(_values[0]*_values[0]+_values[1]*_values[1]));
}

inline HxScalarDouble
HxComplex::arg() const
{
    return HxScalarDouble(::atan2(_values[1], _values[0]));
}

inline HxComplex
HxComplex::ceil() const
{
    return HxComplex(::ceil(_values[0]), ::ceil(_values[1]));
}

inline HxComplex
HxComplex::floor() const
{
    return HxComplex(::floor(_values[0]), ::floor(_values[1]));
}

inline HxComplex
HxComplex::round() const
{
    return HxComplex((int) (_values[0] + ((_values[0] >= 0) ? 0.5 : -0.5)),
                        (int) (_values[1] + ((_values[1] >= 0) ? 0.5 : -0.5)));
}

inline HxComplex
HxComplex::sqrt() const
{
    double a = _values[0];
    double b = _values[1];
    double sq = a*a+b*b;
    double arg = ::atan(b/a)*0.5;
    double mul = ::pow(sq, 0.25)*::exp(a*0.5);

    return HxComplex(mul*::cos(arg), mul*::sin(arg));
}

inline HxComplex
HxComplex::sin() const
{
    return HxComplex(::sin(_values[0])*::cosh(_values[1]),
        ::cos(_values[0])*::sinh(_values[1]));
}

inline HxComplex
HxComplex::cos() const
{
    return HxComplex(::cos(_values[0])*::cosh(_values[1]),
        -::sin(_values[0])*::sinh(_values[1]));
}

inline HxComplex
HxComplex::tan() const
{
    double den = ::cos(2*_values[0])+::cosh(2*_values[1]);

    return HxComplex(::sin(2*_values[0])/den, ::sinh(2*_values[1])/den);
}

inline HxComplex
HxComplex::asin() const
{
    double a = _values[0];
    double b = _values[1];

    HxComplex c = HxComplex(1-a*a+b*b, 2*a*b).sqrt() + HxComplex(-b,a);

    return HxComplex(c.arg().x(), -::log(c.abs().x()));
}

inline HxComplex
HxComplex::acos() const
{
    double a = _values[0];
    double b = _values[1];

    HxComplex c = HxComplex(1-a*a+b*b, 2*a*b).sqrt() + HxComplex(-b,a);

    return HxComplex(M_PI/2.0 - c.arg().x(), ::log(c.abs().x()));
}

inline HxComplex
HxComplex::atan() const
{
    double a = _values[0];
    double b = _values[1];

    HxComplex cp = HxComplex(1-b,a);
    HxComplex cm = HxComplex(1+b,-a);

    return HxComplex(0.5*(cp.arg().x()-cm.arg().x()),
        0.5*(::log(cm.abs().x())-::log(cp.abs().x())));
}

inline HxComplex
HxComplex::sinh() const
{
    return HxComplex(::sinh(_values[0])*::cos(_values[1]),
        ::cosh(_values[0])*::sin(_values[1]));
}

inline HxComplex
HxComplex::cosh() const
{
    return HxComplex(::cosh(_values[0])*::cos(_values[1]),
        ::sinh(_values[0])*::sin(_values[1]));
}

inline HxComplex
HxComplex::tanh() const
{
    double den = ::cosh(2*_values[0])+::cos(2*_values[1]);

    return HxComplex(::sinh(2*_values[0])/den, ::sin(2*_values[1])/den);
}

inline HxComplex
HxComplex::exp() const
{
    double ea = ::exp(_values[0]);
    return HxComplex(::cos(_values[1])*ea, ::sin(_values[1])*ea);
}

inline HxComplex
HxComplex::log() const
{
    double re = _values[0];
    double im = _values[1];
    double mag = ::sqrt(re*re+im*im);

    return HxComplex(::log(mag), ::atan(re/im));
}

const double theLog10 = 1./log(10.0);

inline HxComplex
HxComplex::log10() const
{
    double re = _values[0];
    double im = _values[1];
    double mag = ::sqrt(re*re+im*im);

    return HxComplex(::log(mag)*theLog10, ::atan(re/im)*theLog10);
}

inline HxComplex&
HxComplex::operator+=(const HxComplex& v)
{
    _values[0] += v._values[0];
    _values[1] += v._values[1];
    return *this;
}

inline HxComplex&
HxComplex::operator-=(const HxComplex& v)
{
    _values[0] -= v._values[0];
    _values[1] -= v._values[1];
    return *this;
}

inline HxComplex&
HxComplex::operator*=(const HxComplex& v)
{
    double re = _values[0]*v._values[0] - _values[1]*v._values[1];
    double im = _values[0]*v._values[1] + _values[1]*v._values[0];
    _values[0] = re;
    _values[1] = im;
    return *this;
}

inline HxComplex&
HxComplex::operator/=(const HxComplex& v)
{
    double re = v._values[0];
    double im = v._values[1];
    double sq = re*re+im*im;

    double mulre = _values[0]*re + _values[1]*im;
    double mulim = _values[1]*re - _values[0]*im;
    _values[0] = mulre / sq;
    _values[1] = mulim / sq;

    return *this;
}

inline HxComplex
operator+(const HxComplex& v1, const HxComplex& v2)
{
    return HxComplex(v1._values[0] + v2._values[0],
                        v1._values[1] + v2._values[1]);
}

inline HxComplex
operator-(const HxComplex& v1, const HxComplex& v2)
{
    return HxComplex(v1._values[0] - v2._values[0],
                        v1._values[1] - v2._values[1]);
}

inline HxComplex
operator*(const HxComplex& v1, const HxComplex& v2)
{
    double re = v1._values[0]*v2._values[0] - v1._values[1]*v2._values[1];
    double im = v1._values[0]*v2._values[1] + v1._values[1]*v2._values[0];

    return HxComplex(re, im);
}

inline HxComplex
operator/(const HxComplex& v1, const HxComplex& v2)
{
    double re = v2._values[0];
    double im = v2._values[1];
    double sq = re*re+im*im;

    double mulre = v1._values[0]*re + v1._values[1]*im;
    double mulim = v1._values[1]*re - v1._values[0]*im;

    return HxComplex(mulre / sq, mulim / sq);
}

inline HxComplex
HxComplex::min(const HxComplex& v) const
{
    return (operator<(v)) ? (*this) : v;
}

inline HxComplex&
HxComplex::minAssign(const HxComplex& v)
{
    if (operator<(v))
        return *this;
    operator=(v);
    return *this;
}

inline HxComplex
HxComplex::max(const HxComplex& v) const
{
    return (operator>(v)) ? (*this) : v;
}

inline HxComplex&
HxComplex::maxAssign(const HxComplex& v)
{
    if (operator>(v))
        return *this;
    operator=(v);
    return *this;
}

inline HxComplex
HxComplex::inf(const HxComplex& v) const
{
    return HxComplex((_values[0] < v._values[0]) ? _values[0] : v._values[0],
                       (_values[1] < v._values[1]) ? _values[1] : v._values[1]);
}

inline HxComplex&
HxComplex::infAssign(const HxComplex& v)
{
    _values[0] = (_values[0] < v._values[0]) ? _values[0] : v._values[0];
    _values[1] = (_values[1] < v._values[1]) ? _values[1] : v._values[1];
    return *this;
}

inline HxComplex
HxComplex::sup(const HxComplex& v) const
{
    return HxComplex((_values[0] > v._values[0]) ? _values[0] : v._values[0],
                       (_values[1] > v._values[1]) ? _values[1] : v._values[1]);
}

inline HxComplex&
HxComplex::supAssign(const HxComplex& v)
{
    _values[0] = (_values[0] > v._values[0]) ? _values[0] : v._values[0];
    _values[1] = (_values[1] > v._values[1]) ? _values[1] : v._values[1];
    return *this;
}

inline HxComplex
HxComplex::pow(const HxComplex& v) const
{
    if (v._values[1] == 0) {
        double a = _values[0];
        double b = _values[1];
        double c = v._values[0];
        double sq = a*a+b*b;
        double arg = ::atan(b/a)*c;
        double mul = ::pow(sq, c/2)*::exp(a*c);

        return HxComplex(mul*::cos(arg), mul*::sin(arg));
    }

    return (v * (*this).log()).exp();
}

inline HxComplex
HxComplex::mod(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::and(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::or(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::xor(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::leftShift(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::rightShift(const HxComplex&) const
{
    return (*this);
}

inline HxComplex
HxComplex::cross(const HxComplex&) const
{
    return HxComplex(0, 0);
}

inline HxComplex
HxComplex::sum() const
{
    return *this;
}

inline HxComplex
HxComplex::product() const
{
    return *this;
}

inline HxScalarDouble
HxComplex::min() const
{
    return (_values[0] < _values[1]) ? _values[0] : _values[1];
}

inline HxScalarDouble
HxComplex::max() const
{
    return (_values[0] > _values[1]) ? _values[0] : _values[1];
}

#endif

---