/******************************************
 * C++ Complex Numbers
 * Version: 1.0.9
 * Author:  Douglas Wilhelm Harder
 * Date:    2008/03/03
 *
 * Copyright (c) 2006-2008 by Douglas Wilhelm Harder.
 * All rights reserved.
 *
 * Acknowledgments:
 *      Andre Kostro
 ******************************************/

#include "Complex.h"
#include "Support.h"
#include <iostream>
#include <cmath>
#include <string>

template<typename T> inline char Complex<T>::use_symbol( char sym ) {
	char tmp = imaginary_symbol;

	imaginary_symbol = sym;

	return tmp;
}

template<typename T> inline char Complex<T>::get_symbol() {
	return imaginary_symbol;
}

/******************************************
 * Constructors
 ******************************************/

// (a, b) -> a + bi

template <typename T> Complex<T>::Complex( T re, T im ):r(re), i(im) {
	// empty
}

template <typename T> Complex<T>::Complex( T re ):r(re), i(0.0) {
	// empty
}

/******************************************
 * Assignment Operator
 ******************************************/

// z = 3 -> z == 3 + 0i

template <typename T> inline const Complex<T> & Complex<T>::operator = ( T x ) {
	r = x;
	i = 0.0;

	return *this;
}

/******************************************
 * Mutating Arithmetic Operators
 ******************************************/

// (a + bi) + (c + di) = (a + c) + (b + d)i

template <typename T> inline Complex<T> & Complex<T>::operator += ( const Complex<T> & z ) {
	r += z.r;
	i += z.i;

	return *this;
}

// (a + bi) + x = (a + c) + bi

template <typename T> inline Complex<T> & Complex<T>::operator += ( T x ) {
	r += x;

	return *this;
}

// (a + bi) - (c + di) = (a - c) + (b - d)i

template <typename T> inline Complex<T> & Complex<T>::operator -= ( const Complex<T> & z ) {
	r -= z.r;
	i -= z.i;

	return *this;
}

// (a + bi) - x = (a - x) + bi

template <typename T> inline Complex<T> & Complex<T>::operator -= ( T x ) {
	r -= x;

	return *this;
}

// (a + bi) * (c + di) = (a*c - b*d) + (a*c + b*d)i

template <typename T> inline Complex<T> & Complex<T>::operator *= ( const Complex<T> & z ) {
	T RE = r;

	r =  r * z.r - i * z.i;
	i = RE * z.i + z.r * i;

	return *this;
}

// if ( x != oo )
//     (a + bi) * x = ax + bxi
// else
//     (a + 0i) * oo =  oo +   0i
//     (0 + bi) * oo =   0 +  ooi
//     (0 + 0i) * oo = NaN + NaNi

template <typename T> inline Complex<T> & Complex<T>::operator *= ( T x ) {
	if ( Support<T>::is_inf( x ) && !is_zero() ) {
		r *= ( r == 0.0 ) ? Support<T>::sign( x ) : x;
		i *= ( i == 0.0 ) ? Support<T>::sign( x ) : x;
	} else {
		r *= x;
		i *= x;
	}

	return * this;
}

// n = c^2 + d^2
// (a + bi) / (c + di) = (a*c + b*d)/n + ((a*c - b*d)/n)i

template <typename T> inline Complex<T> & Complex<T>::operator /= ( const Complex<T> & z ) {
	T denom = z.norm();
	T RE = r;

	r = ( RE*z.r + i * z.i)/denom;
	i = (-RE*z.i + z.r * i)/denom;

	return * this;
}

// if ( x != 0 )
//     (a + bi) / x = (a/x) + (b/x)i
// else
//     (a + 0i) / 0 =  oo +   0i
//     (0 + bi) / 0 =   0 +  ooi
//     (0 + 0i) / 0 = NaN + NaNi

template <typename T> inline Complex<T> & Complex<T>::operator /= ( T x ) {
	if ( ( x == 0.0 ) && !is_zero() ) {
		r /= ( r == 0.0 ) ? Support<T>::sign( x ) : x;
		i /= ( i == 0.0 ) ? Support<T>::sign( x ) : x;
	} else {
		r /= x;
		i /= x;
	}

	return * this;
}

// ++(a + bi) = (a + 1) + bi
// Returns (a + 1) + bi

template <typename T> inline Complex<T> & Complex<T>::operator ++() {
	++r;

	return *this;
}

// (a + bi)++ = (a + 1) + bi
// Returns a + bi

template <typename T> inline Complex<T> Complex<T>::operator ++( int ) {
	Complex<T> copy = *this;
	++r;

	return copy;
}

// --(a + bi) = (a - 1) + bi
// Returns (a - 1) + bi

template <typename T> inline Complex<T> & Complex<T>::operator --() {
	--r;

	return *this;
}

// (a + bi)-- = (a - 1) + bi
// Returns a + bi

template <typename T> inline Complex<T> Complex<T>::operator --( int ) {
	Complex<T> copy = *this;
	--r;

	return copy;
}

/******************************************
 * Real-valued Functions
 ******************************************/

template <typename T> inline T Complex<T>::real() const {
	return r;
}

template <typename T> inline T Complex<T>::operator []( int n ) const {
	return reinterpret_cast<const T *>( this )[n];
}

template <typename T> inline T & Complex<T>::operator []( int n ) {
	return reinterpret_cast<T *>( this )[n];
}

template <typename T> inline T Complex<T>::imag_i() const {
	return i;
}

template <typename T> inline T Complex<T>::csgn() const {
	return is_zero() ? 0 : Support<T>::sign( r );
}

// if a != oo && b != oo
//     (a + bi) -> sqrt( a^2 + b^2 )
// else
//     ( a + ooi) -> oo
//     (oo +  bi) -> oo
//     (oo + ooi) -> oo

template <typename T> inline T Complex<T>::abs() const {
	return is_inf()? Support<T>::POS_INF : std::sqrt( r*r + i*i );
}

// if a != oo && b != oo
//     (a + bi) -> a^2 + b^2
// else
//     ( a + ooi) -> oo
//     (oo +  bi) -> oo
//     (oo + ooi) -> oo

template <typename T> inline T Complex<T>::norm() const {
	return is_inf() ? Support<T>::POS_INF : r*r + i*i;
}

// a + bi -> |b|

template <typename T> inline T Complex<T>::abs_imag() const {
	return ( i >= 0.0 ) ? i : -i;
}

// a + bi -> b^2

template <typename T> inline T Complex<T>::norm_imag() const {
	return i*i;
}

// (a + bi) -> arctan( b, a )

template <typename T> inline T Complex<T>::arg() const {
	return std::atan2( i, r );
}

/******************************************
 * Complex<T>-valued Functions
 ******************************************/

// (a + bi) -> (0 + bi)

template <typename T> inline Complex<T> Complex<T>::imag() const {
	return Complex<T>( 0.0, i );
}

// (a + bi) -> (a - bi)

template <typename T> inline Complex<T> Complex<T>::conj() const {
	return Complex<T>( r, -i );
}

template <typename T> inline Complex<T> Complex<T>::operator * () const {
	return Complex<T>( r, -i );
}

template <typename T> inline Complex<T> Complex<T>::signum() const {
	T absq = abs();

	if ( absq == 0.0 || Support<T>::is_nan( absq ) || Support<T>::is_inf( absq ) ) {
		return *this;
	} else {
		return Complex<T>( r/absq, i/absq );
	}
}

// (a + bi)^2 = (a^2 - b^2) + 2abi

template <typename T> inline Complex<T> Complex<T>::sqr() const {
	return Complex<T>( r*r - i*i, 2*r*i );
}

// Branch cut:    (-oo, 0)
// Thanks:   Chris Saunders

template <typename T> inline Complex<T> Complex<T>::sqrt() const {
	T srss = std::sqrt( r*r + i*i );

	return Complex<T>(
		0.5*std::sqrt( 2*(srss + r) ),
		0.5*Complex<T>( i, -r ).csgn()*std::sqrt( 2*( srss - r ) )
	);
}

/******************************************
 * Boolean-valued Functions
 ******************************************/

// (0 + bi) -> true
// false otherwise

template <typename T> inline bool Complex<T>::is_imaginary() const {
	return ( r == 0.0 );
}

// if a == +/-oo || b == +/- oo
//     (a + bi) -> true
// false otherwise

template <typename T> inline bool Complex<T>::is_inf() const {
  	return
		( r == Support<T>::POS_INF ) ||
		( r == Support<T>::NEG_INF ) ||
		( i == Support<T>::POS_INF ) ||
		( i == Support<T>::NEG_INF );
}

// if a == NaN or b == NaN
//     (a + bi) -> true
// false otherwise

template <typename T> inline bool Complex<T>::is_nan() const {
  	return ( r != r ) || ( i != i );
}

// -oo + 0i -> true

template <typename T> inline bool Complex<T>::is_neg_inf() const {
  	return ( r == Support<T>::NEG_INF ) && i == 0.0;
}

// oo + 0i -> true

template <typename T> inline bool Complex<T>::is_pos_inf() const {
  	return ( r == Support<T>::POS_INF ) && ( i == 0.0 );
}

// +/-oo + 0i -> true

template <typename T> inline bool Complex<T>::is_real_inf() const {
  	return ( r == Support<T>::POS_INF || r == Support<T>::NEG_INF ) && i == 0.0;
}

// a + 0i -> true

template <typename T> inline bool Complex<T>::is_real() const {
	return ( i == 0.0 );
}

// 0 + 0i -> true

template <typename T> inline bool Complex<T>::is_zero() const {
	return ( r == 0.0 ) && ( i == 0.0 );
}

/******************************************
 * Exponential and Logarithmic Functions
 ******************************************/

/******************************************
 * Exponential Function
 *
 * Checked for errors.
 ******************************************/

template <typename T> inline Complex<T> Complex<T>::exp() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::exp( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( std::cos( i ), std::sin( i ) );
	} else {
		T expr = std::exp( r );

		return Complex<T>(
			expr * std::cos( i ),
			expr * std::sin( i )
		);
	}
}

/******************************************
 * Logarithmic Function
 *
 * Checked for errors.
 ******************************************/

template <typename T> inline Complex<T> Complex<T>::log() const {
	if ( i == 0 ) {
		if ( r > 0 ) {
			return Complex<T>( std::log( r ), i );
		} else if ( r == 0 ) {
			return Complex<T>( Support<T>::NEG_INF, Support<T>::NaN );
		} else {
			return Complex<T>( std::log( -r ), Support<T>::sign( i )*Support<T>::PI );
		}
	} else if ( r == 0 ) {
		if ( i > 0 ) {
			return Complex<T>( std::log( i ), Support<T>::PI2 );
		} else {
			return Complex<T>( std::log( -i ), -Support<T>::PI2 );
		}
	} else {
		return Complex<T>(
			0.5 * std::log( r*r + i*i ),
			std::atan2( i, r )
		);
	}
}

// Branch cut:    (-oo, 0]

template <typename T> inline Complex<T> Complex<T>::log10() const {
	{
		T ln10 = std::log( 10.0 );

		return Complex<T>(
			0.5 * std::log( r*r + i*i ) / ln10,
			std::atan2( i, r ) / ln10
		);
	}
}

// (a + bi)^(c + di) -> exp( log(a + bi)*(c + di) )

template <typename T> inline Complex<T> Complex<T>::pow( const Complex<T> & z ) const {
	if ( is_zero() && z.is_zero() ) {
		return Complex<T>( Support<T>::NaN, Support<T>::NaN );
	} else if ( is_zero() ) {
		if ( Support<T>::sign( r ) == 1 ) {
			return Complex<T>( 0.0, 0.0 );
		} else {
			return Complex<T>( Support<T>::POS_INF, Support<T>::POS_INF );
		}
	} else if ( z.is_zero() ) {
			return Complex<T>( 1.0, 0.0 );
	}

	if ( i == 0 ) {
		if ( z.i == 0 ) {
			if ( r > 0 ) {
				return Complex<T>( std::pow( r, z.r ), 0.0 );
			} else {
				std::cout << "B" << std::endl;

				return Complex<T>(
					std::exp(z.r*std::log(-r))*std::cos(z.r*Support<T>::PI),
					-std::exp(z.r*std::log(-r))*std::sin(z.r*Support<T>::PI)
				);
			}
		} else if ( z.r == 0 ) {
			if ( r > 0 ) {
				return Complex<T>(
					std::cos( z.i*std::log( r ) ),
					std::sin( z.i*std::log( r ) )
				);
			} else {
				std::cout << "D" << std::endl;
				if ( Support<T>::sign( i ) == 1 ) {
					return Complex<T>(
						std::exp( -z.i*Support<T>::PI )*std::cos( z.i*std::log( -r ) ),
						std::exp( -z.i*Support<T>::PI )*std::sin( z.i*std::log( -r ) )
					);
				} else {
					return Complex<T>(
						std::exp( z.i*Support<T>::PI )*std::cos( z.i*std::log( -r ) ),
						std::exp( z.i*Support<T>::PI )*std::sin( z.i*std::log( -r ) )
					);
				}
			}
		} else {
			if ( r > 0 ) {
				return Complex<T>(
					std::exp( z.r*std::log(r) )*std::cos( z.i*std::log(r)),
					std::exp( z.r*std::log(r) )*std::sin( z.i*std::log(r) )
				);
			} else {
				if ( Support<T>::sign( i ) == 1 ) {
					return Complex<T>(
						std::exp( z.r*std::log(-r) - z.i*Support<T>::PI )*
						std::cos( z.i*std::log(-r) + z.r*Support<T>::PI ),
						std::exp( z.r*std::log(-r) - z.i*Support<T>::PI )*
						std::sin( z.i*std::log(-r) + z.r*Support<T>::PI )
					);
				} else {
					return Complex<T>(
						std::pow( -r, z.r )*
						std::exp(z.i*Support<T>::PI)*std::cos(-z.i*std::log(-r) + z.r*Support<T>::PI),
						-std::pow( -r, z.r )*
						std::exp(z.i*Support<T>::PI)*std::sin(-z.i*std::log(-r) + z.r*Support<T>::PI)
					);
				}
			}
		}
	} else if ( r == 0 ) {
		if ( z.i == 0 ) {
			return Complex<T>(
				std::exp(0.5*z.r*std::log(i*i))*std::cos(0.5*z.r*Support<T>::sign(i)*Support<T>::PI),
				std::exp(0.5*z.r*std::log(i*i))*std::sin(0.5*z.r*Support<T>::sign(i)*Support<T>::PI)
			);
		} else if ( z.r == 0 ) {
			return Complex<T>(
				std::exp(-0.5*z.i*Support<T>::sign(i)*Support<T>::PI)*std::cos(0.5*z.i*std::log(i*i)),
				std::exp(-0.5*z.i*Support<T>::sign(i)*Support<T>::PI)*std::sin(0.5*z.i*std::log(i*i))
			);
		} else {
			return Complex<T>(
				std::exp(0.5*z.r*std::log(i*i) - 0.5*z.i*Support<T>::sign(i)*Support<T>::PI)*std::cos(0.5*z.i*std::log(i*i) + 0.5*z.r*Support<T>::sign(i)*Support<T>::PI),
				std::exp(0.5*z.r*std::log(i*i) - 0.5*z.i*Support<T>::sign(i)*Support<T>::PI)*std::sin(0.5*z.i*std::log(i*i) + 0.5*z.r*Support<T>::sign(i)*Support<T>::PI)
			);
		}
	} else {
		if ( z.i == 0 ) {
			return Complex<T>(
				std::exp(0.5*z.r*std::log(r*r + i*i))*std::cos(z.r*std::atan2(i,r)),
				std::exp(0.5*z.r*std::log(r*r + i*i))*std::sin(z.r*std::atan2(i,r))
			);
		} else if ( z.r == 0 ) {
			return Complex<T>(
				std::exp(-z.i*std::atan2(i,r))*std::cos(0.5*z.i*std::log(r*r + i*i)),
				std::exp(-z.i*std::atan2(i,r))*std::sin(0.5*z.i*std::log(r*r + i*i))
			);
		} else {
			return Complex<T>(
				std::exp(0.5*z.r*std::log(r*r + i*i)-z.i*std::atan2(i,r))*std::cos(0.5*z.i*std::log(r*r + i*i) + z.r*std::atan2(i,r)),
				std::exp(0.5*z.r*std::log(r*r + i*i)-z.i*std::atan2(i,r))*std::sin(0.5*z.i*std::log(r*r + i*i) + z.r*std::atan2(i,r))
			);
		}
	}
}

template <typename T> inline Complex<T> Complex<T>::pow(T x) const {
	if ( i == 0 ) {
		return Complex<T>(
			std::exp( x*std::log( std::fabs(r) ) )*std::cos( x*(0.5 - 0.5*Support<T>::sign(r) )*Support<T>::PI ),
			std::exp( x*std::log( std::fabs(r) ) )*std::sin( x*(0.5 - 0.5*Support<T>::sign(r) )*Support<T>::PI )
		);
	} else if ( r == 0 ) {
		return Complex<T>(
			std::exp( 0.5*x*std::log( i*i ) )*std::cos( 0.5*x*Support<T>::sign(i)*Support<T>::PI ),
			std::exp( 0.5*x*std::log( i*i ) )*std::sin( 0.5*x*Support<T>::sign(i)*Support<T>::PI )
		);
	} else {
		return Complex<T>(
			std::exp( 0.5*x*std::log( r*r + i*i ) )*std::cos( x*std::atan2( i, r ) ),
			std::exp( 0.5*x*std::log( r*r + i*i ) )*std::sin( x*std::atan2( i, r ) )
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::inverse() const {
	if ( is_zero() ) {
		return Complex( Support<T>::sign( r )*Support<T>::POS_INF, -Support<T>::sign( i )*Support<T>::POS_INF );
	} else if ( is_inf() ) {
		return Complex( Support<T>::sign( r )*0.0, -Support<T>::sign( i )*0.0 );
	} else if ( is_nan() ) {
		return Complex( Support<T>::NaN, Support<T>::NaN );
	} else {
		T denom = norm();

		return Complex( r/denom, -i/denom );
	}
}

/******************************************
 * Trigonometric and Hyperbolic Functions
 ******************************************/

// a + 0i -> sin(a) + 0i
// 0 + bi ->     0  + sinh(b)i
// a + bi -> [sin(a) cosh(b)] + [cos(a) sinh(b)]i


template <typename T> inline Complex<T> Complex<T>::sin() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::sin( r ), std::cos( r )*i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, std::sinh( i ) );
	} else {
		return Complex<T>(
			std::sin( r ) * std::cosh( i ),
			std::cos( r ) * std::sinh( i )
		);
	}
}

/******************************************
 * Cosine Function
 *
 * a + 0i -> cos(a) + 0i
 * 0 + bi ->     0  + cosh(b)i
 * a + bi -> [cos(a) cosh(b)] - [sin(a) sinh(b)]i
 *
 * Checked for errors.
 ******************************************/

template <typename T> inline Complex<T> Complex<T>::cos() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::cos( r ), -i*std::sin( r ) );
	} else if ( r == 0.0 ) {
		return Complex<T>( std::cosh( i ), -i*r );
	} else {
		return Complex<T>(
			 std::cos( r ) * std::cosh( i ),
			-std::sin( r ) * std::sinh( i )
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::tan() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::tan( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, std::tanh( i ) );
	} else {
		T cosr =  std::cos(r);
		T sinhi = std::sinh(i);

		T denom = cosr*cosr + sinhi*sinhi;

		return Complex<T>(
			std::sin( r )  * cosr  / denom,
			sinhi * std::cosh( i ) / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::sec() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::sec( r ), -i );
	} else if ( r == 0.0 ) {
		return Complex<T>( Support<T>::sech( i ), r*Support<T>::sign( i ) );
	} else {
		T cosr =  std::cos(r);
		T sinhi = std::sinh(i);

		T denom = cosr*cosr + sinhi*sinhi;

		return Complex<T>(
			cosr  * std::cosh( i )  / denom,
			std::sin( r ) * sinhi / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::csc() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::csc( r ), -i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, -Support<T>::csch( i ) );
	} else {
		T sinr =  std::sin(r);
		T sinhi = std::sinh(i);

		T denom = sinr*sinr + sinhi*sinhi;

		return Complex<T>(
			 sinr  * std::cosh( i ) / denom,
			-std::cos( r ) * sinhi  / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::cot() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::cot( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, -Support<T>::coth( i ) );
	} else {
		T sinr =  std::sin(r);
		T sinhi = std::sinh(i);

		T denom = sinr*sinr + sinhi*sinhi;

		return Complex<T>(
			 sinr  * std::cos( r )  / denom,
			-sinhi * std::cosh( i ) / denom
		);
	}
}

// a + 0i ->  sinh(a) + 0i
// 0 + bi ->      0  + sin(b)i
// a + bi -> [sinh(a)  cos(b)] + [cosh(a) sin(b)]i

template <typename T> inline Complex<T> Complex<T>::sinh() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::sinh( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, std::sin( i ) );
	} else {
		return Complex<T>(
			std::sinh( r ) * std::cos( i ),
			std::cosh( r ) * std::sin( i )
		);
	}
}

// a + 0i ->  cosh(a) + 0i
// 0 + bi ->   cos(b) + 0i
// a + bi -> [cosh(a) cos(b)] + [sinh(a) sin(b)]i

template <typename T> inline Complex<T> Complex<T>::cosh() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::cosh( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( std::cos( i ), -r );
	} else {
		return Complex<T>(
			std::cosh( r ) * std::cos( i ),
			std::sinh( r ) * std::sin( i )
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::tanh() const {
	if ( i == 0.0 ) {
		return Complex<T>( std::tanh( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, std::tan( i ) );
	} else {
		T sinhr = std::sinh(r);
		T cosi =  std::cos(i);

		T denom = sinhr*sinhr + cosi*cosi;

		return Complex<T>(
			sinhr  * std::cosh( r )  / denom,
			std::sin( i ) * cosi / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::sech() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::sech( r ), -i );
	} else if ( r == 0.0 ) {
		return Complex<T>( Support<T>::sec( i ), r );
	} else {
		T sinhr = std::sinh(r);
		T cosi =  std::cos(i);

		T denom = sinhr*sinhr + cosi*cosi;

		return Complex<T>(
			 std::cosh( r )  * cosi  / denom,
			-sinhr * std::sin( i ) / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::csch() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::csch( r ), -i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, -Support<T>::csc( i ) );
	} else {
		T sinhr =  std::sinh(r);
		T sini = std::sin(i);

		T denom = sinhr*sinhr + sini*sini;

		return Complex<T>(
			 sinhr  * std::cos( i )  / denom,
			-std::cosh( r ) * sini / denom
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::coth() const {
	if ( i == 0.0 ) {
		return Complex<T>( Support<T>::coth( r ), i );
	} else if ( r == 0.0 ) {
		return Complex<T>( r, -Support<T>::cot( i ) );
	} else {
		T sinhr =  std::sinh(r);
		T sini = std::sin(i);

		T denom = sinhr*sinhr + sini*sini;

		return Complex<T>(
			 sinhr  * std::cosh( r )  / denom,
			-std::cos( i ) * sini / denom
		);
	}
}

// Branch cut:    (-oo, -1) U (1, oo)

template <typename T> inline Complex<T> Complex<T>::asin() const {
	if ( i == 0 ) {
		if ( r >= -1 && r <= 1 ) {
			return Complex<T>( std::asin( r ), i );
		} else {
			if ( r > 1 ) {
				return Complex<T>(
					Support<T>::PI2,
					Support<T>::sign( i )*std::log(
						r + std::sqrt(r*r  - 1)
					)
				);
			} else {
				return Complex<T>(
					-Support<T>::PI2,
					Support<T>::sign( i )*std::log(
						-r + std::sqrt(r*r  - 1)
					)
				);
			}
		}
	} else if ( r == 0 ) {
		return Complex<T>( r, std::log( i + std::sqrt( i*i + 1 ) ) );
	} else {
		T ss = r*r + i*i + 1;
		T ssp2r = std::sqrt( ss + 2*r );
		T ssm2r = std::sqrt( ss - 2*r );
		T sum = 0.5*( ssp2r + ssm2r );

		return Complex<T>(
			std::asin( ssp2r/2 - ssm2r/2 ),
			Complex<T>( i, -r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			)
		);
	}
}

// Branch cut:    (-oo, -1) U (1, oo)

template <typename T> inline Complex<T> Complex<T>::acos() const {
	if ( i == 0 ) {
		if ( r >= -1 && r <= 1 ) {
			return Complex<T>( std::acos( r ), -i );
		} else if ( r > 1 ) {
			return Complex<T>(
				0.0,
				-Support<T>::sign( i )*std::log(
					r + std::sqrt(r*r  - 1)
				)
			);
		} else {
			return Complex<T>(
				Support<T>::PI,
				-Support<T>::sign( i )*std::log(
					-r + std::sqrt(r*r  - 1)
				)
			);
		}
	} else {
		T ss = r*r + i*i + 1;
		T ssp2r = std::sqrt( ss + 2*r );
		T ssm2r = std::sqrt( ss - 2*r );
		T sum = 0.5*( ssp2r + ssm2r );

		return Complex<T>(
			std::acos( 0.5*(ssp2r - ssm2r) ),
			-Complex<T>( i, -r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			)
		);
	}
}

/**********************************************************
 * Arc Tangent Function
 *
 * Branch cut:    (-ooi, -i] U [i, ooi)
 *
 * There do not appear to be any errors in this function
 * for finite arguments.
 * Should simplify for r == 0, however.
 **********************************************************/

template <typename T> inline Complex<T> Complex<T>::atan() const {
	if ( r == 0 ) {
		if ( i < -1 ) {
			T ip = i + 1;
			T im = i - 1;

			return Complex<T>(
				Support<T>::sign(r)*Support<T>::PI2,
				0.25*std::log( ip*ip/(im*im) )
			);
		} else if ( i == -1 ) {
			return Complex<T>( Support<T>::NaN, Support<T>::NEG_INF );
		} else if ( i == 0 ) {
			return *this;
		} else if ( i < 1 ) {
			T ip = i + 1;
			T im = i - 1;

			return Complex<T>(
				r,
				0.25*std::log( ip*Support<T>::PI/(im*im) )
			);
		} else if ( i == 1 ) {
			return Complex<T>( Support<T>::NaN, Support<T>::POS_INF );
		} else {
			T ip = i + 1;
			T im = i - 1;

			return Complex<T>(
				Support<T>::sign(r)*Support<T>::PI2,
				0.25*std::log( ip*ip/(im*im) )
			);
		}
	} else if ( i == 0 ) {
		return Complex<T>( std::atan( r ), i );
	} else {
      T opi = 1.0 + i;
      T omi = 1.0 - i;
		T rr = r*r;

		return Complex<T>(
			0.5*( std::atan2( r, omi ) - std::atan2( -r, opi ) ),
			0.25*std::log( (rr + opi*opi)/(rr + omi*omi) )
		);
	}
}

/**********************************************************
 * Arc Secant Function
 *
 * Branch cut:    (-1, 1)
 *
 * There do not appear to be any errors in this function
 * for finite arguments.
 **********************************************************/

template <typename T> inline Complex<T> Complex<T>::asec() const {
	if ( i == 0 ) {
		if ( ( r <= -1 ) || ( r >= 1 ) ) {
			return Complex<T>( std::acos( 1/r ), i );
		} else if ( r == 0 ) {
			return Complex<T>( Support<T>::POS_INF, Support<T>::POS_INF );
		} else {
			if ( r < 0 ) {
				return Complex<T>(
					Support<T>::PI,
					Support<T>::sign( i )*std::log(
						-1/r + std::sqrt(1/(r*r)  - 1)
					)
				);
			} else {
				return Complex<T>(
					0.0,
					Support<T>::sign( i )*std::log(
						1/r + std::sqrt(1/(r*r)  - 1)
					)
				);
			}
		}
	} else {
		T ss = r*r + i*i;
		T p1 = r/ss;

		T p1p1 = p1 + 1;
		p1p1 = p1p1 * p1p1;

		T p1m1 = p1 - 1;
		p1m1 = p1m1 * p1m1;

		T p2 = i/ss;
		p2 = p2*p2;

		T p1p1p2 = std::sqrt( p1p1 + p2 );
		T p1m1p2 = std::sqrt( p1m1 + p2 );

		T sum = 0.5*( p1p1p2 + p1m1p2 );

		return Complex<T>(
			std::acos( 0.5*(p1p1p2 - p1m1p2 ) ),
			Complex<T>( i, r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			)
		);
	}
}

/**********************************************************
 * Arc Cosecant Function
 *
 * Branch cut:    (-1, 1)
 *
 * There do not appear to be any errors in this function
 * for finite arguments.
 **********************************************************/

template <typename T> inline Complex<T> Complex<T>::acsc() const {
	if ( i == 0 ) {
		if ( ( r <= -1 ) || ( r >= 1 ) ) {
			// (-inf, 1] U [1, inf)

			return Complex<T>( std::asin( 1/r ), -i );
		} else if ( r == 0 ) {
			// 0

			return Complex<T>( Support<T>::POS_INF, Support<T>::POS_INF );
		} else if ( r < 0 ) {
			// (-1, 0)

			return Complex<T>(
				-Support<T>::PI2,
				-Support<T>::sign( i )*std::log(
					-1/r + std::sqrt(1/(r*r)  - 1)
				)
			);
		} else {
			// (0, 1)

			return Complex<T>(
				Support<T>::PI2,
				-Support<T>::sign( i )*std::log(
					1/r + std::sqrt(1/(r*r)  - 1)
				)
			);
		}
	} else if ( r == 0 ) {
		return Complex<T>(
			r,
			-std::log( 1/i + std::sqrt( 1 + 1/(i*i) ) )
		);
	} else {
		T ss = r*r + i*i;
		T p1 = r/ss;

		T p1p1 = p1 + 1;
		p1p1 = p1p1 * p1p1;

		T p1m1 = p1 - 1;
		p1m1 = p1m1 * p1m1;

		T p2 = i/ss;
		p2 = p2*p2;

		T p1p1p2 = std::sqrt( p1p1 + p2 );
		T p1m1p2 = std::sqrt( p1m1 + p2 );

		T sum = 0.5*( p1p1p2 + p1m1p2 );

		return Complex<T>(
			std::asin( 0.5*(p1p1p2 - p1m1p2) ),
			-Complex<T>( i, r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			)
		);
	}
}

// Branch cut:    (-i, i)

template <typename T> inline Complex<T> Complex<T>::acot() const {
	if ( i == 0 ) {
		return Complex<T>( Support<T>::PI2 - std::atan( r ), -i );
	} else if ( r == 0 ) {
		if ( i < -1 ) {
			// (-inf i, -i)

			return Complex<T>(
				Support<T>::sign( r ) == 1 ? 0.0 : Support<T>::PI,
				0.5*std::log( (-i + 1)/(-i - 1) )
			);
		} else if ( i == -1 ) {
			return Complex<T>( Support<T>::NaN, Support<T>::POS_INF );
		} else if ( i < 1 ) {
			// (-i, i)

			return Complex<T>(
				Support<T>::PI2,
				-0.25*std::log( (i + 1)*(i + 1)/((i - 1)*(i - 1)) )
			);
		} else if ( i == 1 ) {
			return Complex<T>( Support<T>::NaN, Support<T>::NEG_INF );
		} else {
			// (i, inf i)

			return Complex<T>(
				Support<T>::sign( r ) == 1 ? 0.0 : Support<T>::PI,
				0.5*std::log( (-i + 1)/(-i - 1) )
			);
		}
	} else {
      T opi = 1.0 + i;
      T omi = 1.0 - i;
		T rr = r*r;

		return Complex<T>(
			Support<T>::PI2 + 0.5*( std::atan2( -r, opi ) - std::atan2( r, omi ) ),
			-0.25*std::log( (rr + opi*opi)/(rr + omi*omi) )
		);
	}
}

/**********************************************************
 * Arc Hyperbolic Sine Function
 *
 * Branch cut:    (-ooi, -i) U (i, ooi)
 *
 * There do not appear to be any errors in this function
 * for finite arguments.
 **********************************************************/

template <typename T> inline Complex<T> Complex<T>::asinh() const {
	if ( r == 0 ) {
		if ( ( i >= -1 ) && ( i <= 1 ) ) {
			// [-i, i]
			return Complex<T>( r, std::asin( i ) );
		} else if ( i > 1 ) {
			// (i, ooi)

			return Complex<T>(
				Support<T>::sign( r )*std::log( i + std::sqrt( i*i - 1 ) ),
				Support<T>::PI2
			);
		} else {
			// (-ooi, -i)

			return Complex<T>(
				Support<T>::sign( r )*std::log( -i + std::sqrt( i*i - 1 ) ),
				-Support<T>::PI2
			);
		}
	} else if ( i == 0 ) {
		return Complex<T>( std::log( r + std::sqrt( r*r + 1 ) ), i );
	} else {
		T ss = r*r + i*i + 1;
		T ssp2i = std::sqrt( ss + 2*i );
		T ssm2i = std::sqrt( ss - 2*i );
		T sum = 0.5*( ssp2i + ssm2i );

		return Complex<T>(
			Complex<T>( r, i ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			),
			std::asin( ssp2i/2 - ssm2i/2 )
		);
	}
}

// Branch cut:    (-oo, 1)

template <typename T> inline Complex<T> Complex<T>::acosh() const {
	if ( i == 0 ) {
		if ( r > 1 ) {
			// (1, oo)

			return Complex<T>( std::log( r + std::sqrt( (r - 1)*(r + 1) ) ), i );
		} else if ( r == 1 ) {
			// 1

			return Complex<T>( 0.0, i );
		} else if ( r == 0 ) {
			// 0

			return Complex<T>( 0.0, Support<T>::sign( i )*Support<T>::PI2 );
		} else if ( r > -1 ) {
			// (-1, 0) U (0, 1)

			return Complex<T>( 0.0, Support<T>::sign( i )*std::acos( r ) );
		} else if ( r == -1 ) {
			// -1

			return Complex<T>( 0.0, Support<T>::sign( i )*Support<T>::PI );
		} else {
			// (-oo, -1)
			return Complex<T>( std::log( -r + std::sqrt( (-1 - r)*(1 - r) ) ), Support<T>::sign( i )*Support<T>::PI );
		}
	} else if ( r == 0 ) {
		return Complex<T>( std::log( std::sqrt( i*i + 1 ) + Support<T>::sign( i )*i ), Support<T>::sign( i )*Support<T>::PI2 );
	} else {
		T ss = r*r + i*i + 1;
		T ssp2r = std::sqrt( ss + 2*r );
		T ssm2r = std::sqrt( ss - 2*r );
		T sum = 0.5*( ssp2r + ssm2r );

		return Complex<T>(
			Complex<T>( -i, r - 1.0 ).csgn()*
			Complex<T>( -i, r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			),
			-Complex<T>( -i, r - 1.0 ).csgn()*std::acos( ssp2r/2 - ssm2r/2 )
		);
	}
}

// Branch cut:    (-oo, -1] U [1, oo)

template <typename T> inline Complex<T> Complex<T>::atanh() const {
	if ( i == 0 ) {
		if ( r > 1 ) {
			// [1, oo)

			return Complex<T>(
				0.25*std::log( (r + 1)*(r + 1)/((r - 1)*(r - 1)) ),
				Support<T>::sign( i )*Support<T>::PI2
			);
		} else if ( r == 1 ) {
			return Complex<T>( Support<T>::POS_INF, Support<T>::NaN );
		} else if ( r > 0 ) {
			return Complex<T>(
				0.5*std::log( (1 + r)/(1 - r) ),
				i
			);
		} else if ( r == 0 ) {
			return *this;
		} else if ( r > -1 ) {
			return Complex<T>(
				-0.5*std::log( (1 - r)/(1 + r) ),
				i
			);
		} else if ( r == -1 ) {
			return Complex<T>( Support<T>::NEG_INF, Support<T>::NaN );
		} else {
			T rp = r + 1;
			T rm = r - 1;
			rp *= rp;
			rm *= rm;
			T logrprm = 0.25*std::log( rp/rm );

			return Complex<T>( logrprm, Support<T>::sign( i )*Support<T>::PI2 );
		}
	} else if ( r == 0 ) {
		return Complex<T>( r, std::atan( i ) );
	} else {
      T opr = 1.0 + r;
      T omr = 1.0 - r;
		T ii = i*i;

		return Complex<T>(
			0.25*std::log( (ii + opr*opr)/(ii + omr*omr) ),
			0.5*( std::atan2( i, opr ) - std::atan2( -i, omr ) )
		);
	}
}

// Branch cut:    (-oo, 0] U (1, oo)

template <typename T> inline Complex<T> Complex<T>::asech() const {
	if ( i == 0 ) {
		if ( r < -1 ) {
			return Complex<T>( 0.0, -Support<T>::sign( i )*std::acos( 1/r ) );
		} else if ( r == -1 ) {
			return Complex<T>( 0.0, -Support<T>::sign( i )*Support<T>::PI );
		} else if ( r < 0 ) {
			return Complex<T>(
				std::log( -1/r + std::sqrt( -(1 - 1/r)*(1 + 1/r) ) ),
				-Support<T>::sign( i )*Support<T>::PI
			);
		} else if ( r == 0 ) {
			return Complex<T>( Support<T>::POS_INF, Support<T>::NaN );
		} else if ( r <= 1 ) {
			return Complex<T>(
				std::log( 1/r + std::sqrt( (1/r + 1)*(1/r - 1) ) ),
				-Support<T>::sign( i )*0.0
			);
		} else {
			return Complex<T>( 0.0, -Support<T>::sign( i )*std::acos( 1/r ) );
		}
	} else if ( r == 0 ) {
		if ( i > 0 ) {
			return Complex<T>(
				std::log( std::sqrt( 1 + 1/(i*i) ) + 1/i ), -Support<T>::PI2
			);
		} else {
			return Complex<T>(
				std::log( std::sqrt( 1 + 1/(i*i) ) - 1/i ), Support<T>::PI2
			);
		}
	} else {
		T ss = r*r + i*i;
		T Rr = r/ss;

		T Rrp1 = Rr + 1;
		Rrp1 = Rrp1 * Rrp1;

		T Rrm1 = Rr - 1;
		Rrm1 = Rrm1 * Rrm1;

		T Ri = i/ss;
		Ri = Ri*Ri;

		T Rrp1pRi = std::sqrt( Rrp1 + Ri );
		T Rrm1pRi = std::sqrt( Rrm1 + Ri );

		T sum = 0.5*( Rrp1pRi + Rrm1pRi );

		return Complex<T>(
			-Complex<T>( -i, ss - r ).csgn()*Complex<T>( i, r ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			),
			Complex<T>( -i, ss - r ).csgn()*(
				std::acos( 0.5*(Rrp1pRi - Rrm1pRi) )
			)
		);
	}
}

// Branch cut:    (-i, i)
// NOT YET FINISHED

template <typename T> inline Complex<T> Complex<T>::acsch() const {
	if ( r == 0 ) {
		if ( i <= -1 ) {
			return Complex<T>( r, -std::asin( 1/i ) );
		} else if ( i == -1 ) {
			return Complex<T>( r, -Support<T>::PI2 );
		} else if ( i < 0 ) {
			return Complex<T>(
				Support<T>::sign( r )*std::log( -1/i + std::sqrt( 1/(i*i) - 1 ) ),
				Support<T>::PI2
			);
		} else if ( i == 0 ) {
			return Complex<T>( Support<T>::POS_INF, Support<T>::NEG_INF );
		} else if ( i < 1 ) {
			return Complex<T>(
				Support<T>::sign( r )*std::log( 1/i + std::sqrt( 1/(i*i) - 1 ) ),
				-Support<T>::PI2
			);
		} else if ( i == 1 ) {
			return Complex<T>( r, -Support<T>::PI2 );
		} else {
			return Complex<T>( r, -std::asin( 1/i ) );
		}
	} else if ( i == 0 ) {
		return Complex<T>(
			std::log( 1/r + std::sqrt( 1/(r*r) + 1 ) ),
			-i
		);
	} else {
		T ss = r*r + i*i;
		T p1 = i/ss;

		T p1p1 = p1 + 1;
		p1p1 = p1p1 * p1p1;

		T p1m1 = p1 - 1;
		p1m1 = p1m1 * p1m1;

		T p2 = r/ss;
		p2 = p2*p2;

		T p1p1p2 = std::sqrt( p1p1 + p2 );
		T p1m1p2 = std::sqrt( p1m1 + p2 );

		T sum = 0.5*( p1p1p2 + p1m1p2 );

		return Complex<T>(
			Complex<T>( r, -i ).csgn()*std::log(
				sum + std::sqrt(sum*sum  - 1)
			),
			std::asin( 0.5*(p1m1p2 - p1p1p2) )
		);
	}
}

// Branch cut:    [-1, 1]
// NOT YET FINISHED

template <typename T> inline Complex<T> Complex<T>::acoth() const {
	if ( i == 0 ) {
		if ( r < -1 ) {
			return Complex<T>(
				0.5*std::log( (r + 1)/(r - 1) ),
				-i
			);
		} else if ( r == -1 ) {
			return Complex<T>( Support<T>::NEG_INF, Support<T>::NaN );
		} else if ( r < 0 ) {
			return Complex<T>(
				0.5*std::log( (r + 1)/(1 - r) ),
				-Support<T>::sign( i )*Support<T>::PI2
			);
		} else if ( r == 0 ) {
			return Complex( r, -Support<T>::sign( i )*Support<T>::PI2 );
		} else if ( r < 1 ) {
			return Complex<T>(
				0.5*std::log( (r + 1)/(1 - r) ),
				-Support<T>::sign( i )*Support<T>::PI2
			);
		} else if ( r == 1 ) {
			return Complex<T>( Support<T>::POS_INF, Support<T>::NaN );
		} else {
			return Complex<T>(
				0.5*std::log( (r + 1)/(r - 1) ),
				-i
			);
		}
	} else if ( r == 0 ) {
		if ( i > 0 ) {
			return Complex<T>( r, -Support<T>::PI2 + std::atan( i ) );
		} else {
			return Complex<T>( r, Support<T>::PI2 + std::atan( i ) );
		}
	} else {
      T rp1 = r + 1.0;
      T rm1 = r - 1.0;
		T ii = i*i;

		return Complex<T>(
			0.25*std::log( (ii + rp1*rp1)/(ii + rm1*rm1) ),
			0.5*( std::atan2( i, rp1 ) - std::atan2( i, rm1 ) )
		);
	}
}

template <typename T> inline Complex<T> Complex<T>::bessel_J( int n ) const {
	if ( n < 0 ) {
		return (n & 1) ? -bessel_J( -n ) : bessel_J( -n );
	}

	Complex<T> z = *this / 2.0;

	Complex<T> term = 1.0;

	for ( int i = 1; i <= n; ++i ) {
		term *= z/i;
	}

	Complex<T> z2 = z.sqr();

	Complex<T> result = term;
	Complex<T> result_old = term;

	for ( T i = 1; ; i += 2.0 ) {
		term *= z2 / (i*(n + i));
		result -= term;

		term *= z2 / ((i + 1.0)*(n + i + 1.0));
		result += term;

		if ( result_old == result ) {
			return result;
		}

		result_old = result;
	}
}

/******************************************
 * Integer Functions
 ******************************************/

// (a + bi) -> ceil(a) + ceil(b)i

template <typename T> inline Complex<T> Complex<T>::ceil() const {
	return Complex<T>(
		std::ceil(r), std::ceil(i)
	);
}

// (a + bi) -> floor(a) + floor(b)i

template <typename T> inline Complex<T> Complex<T>::floor() const {
	return Complex<T>(
		std::floor(r), std::floor(i)
	);
}


/******************************************
 * Horner's Rule
 *
 *   The polynomial is defined by giving the highest
 *   coefficient first:
 *
 *            n - 1         n - 2
 *      v[0]*z      + v[1]*z      + ... + v[n-2]*z + v[n-1]
 *
 *   This is the same as with Matlab.
 ******************************************/

template <typename T> inline Complex<T> Complex<T>::horner( T * v, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = *v;

		for ( unsigned int j = 1; j < n; ++j ) {
			ar = ar*r + *(++v);
		}

		return Complex<T>( ar, 0 );
	} else if ( r == 0 ) {              // imaginary
		T ar = *v;
		T ai = 0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T aR = ar;

			ar = -ai*i + *(++v);
			ai = aR*i;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = *v;
		T ai = 0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T aR = ar;

			ar = ar*r - ai*i + *(++v);
			ai = aR*i + ai*r;
		}

		return Complex<T>( ar, ai );
	}
}

template <typename T> inline Complex<T> Complex<T>::horner( T * v, T * c, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = *v;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - *(++c);

			ar = ar*dr + *(++v);
		}

		return Complex<T>( ar, 0 );
	} else if ( r == 0 ) {              // imaginary
		T ar = *v;
		T ai = 0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = -*(++c);
			T aR = ar;

			ar = ar*dr - ai*i + *(++v);
			ai = aR*i + ai*dr;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = *v;
		T ai = 0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - *(++c);
			T aR = ar;

			ar = ar*dr - ai*i + *(++v);
			ai = aR*i + ai*dr;
		}

		return Complex<T>( ar, ai );
	}

}

template <typename T> inline Complex<T> Complex<T>::horner( Complex<T> * v, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			ar = ar*r + (*(++v)).r;
			ai = ai*r + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else if ( r == 0 ) {              // imaginary
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T aR = ar;

			ar = -ai*i + (*(++v)).r;
			ai = aR*i + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T aR = ar;

			ar = ar*r - ai*i + (*(++v)).r;
			ai = aR*i + ai*r + (*v).i;
		}

		return Complex<T>( ar, ai );
	}
}

template <typename T> inline Complex<T> Complex<T>::horner( Complex<T> * v, T * c, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - *(++c);

			ar = ar*dr + (*(++v)).r;
			ai = ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else if ( r == 0 ) {              // imaginary
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = -*(++c);
			T aR = ar;

			ar = ar*dr - ai*i + (*(++v)).r;
			ai = aR*i + ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - *(++c);
			T aR = ar;

			ar = ar*dr - ai*i + (*(++v)).r;
			ai = aR*i + ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	}
}

template <typename T> inline Complex<T> Complex<T>::horner( T * v, Complex<T> * c, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = *v;
		T ai = 0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - (*(++c)).r;
			T di = -(*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + *(++v);
			ai = aR*di + ai*dr;
		}

		return Complex<T>( ar, ai );
	} else if ( r == 0 ) {              // imaginary
		T ar = *v;
		T ai = 0.0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = -(*(++c)).r;
			T di = i - (*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + *(++v);
			ai = aR*di + ai*dr;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = *v;
		T ai = 0.0;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - (*(++c)).r;
			T di = i - (*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + *(++v);
			ai = aR*di + ai*dr;
		}

		return Complex<T>( ar, ai );
	}
}

template <typename T> inline Complex<T> Complex<T>::horner( Complex<T> * v, Complex<T> * c, unsigned int n ) const {
	if ( n == 0 ) {
		return ZERO;
	}

	if ( i == 0 ) {                     // real
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - (*(++c)).r;
			T di = -(*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + (*(++v)).r;
			ai = aR*di + ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else if ( r == 0 ) {              // imaginary
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = -(*(++c)).r;
			T di = i - (*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + (*(++v)).r;
			ai = aR*di + ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	} else {                            // complex
		T ar = (*v).r;
		T ai = (*v).i;

		for ( unsigned int j = 1; j < n; ++j ) {
			T dr = r - (*(++c)).r;
			T di = i - (*c).i;
			T aR = ar;

			ar = ar*dr - ai*di + (*(++v)).r;
			ai = aR*di + ai*dr + (*v).i;
		}

		return Complex<T>( ar, ai );
	}
}

/******************************************
 * Random Factories
 ******************************************/

// -> a + bi

template <typename T> inline Complex<T> Complex<T>::random() {
	return Complex<T>(
		(static_cast<T>( rand() ))/RAND_MAX,
		(static_cast<T>( rand() ))/RAND_MAX
	);
}

// -> 0 + bi

template <typename T> inline Complex<T> Complex<T>::random_imag() {
	return Complex<T>( 0.0, (static_cast<T>( rand() ))/RAND_MAX );
}

// -> a + 0i

template <typename T> inline Complex<T> Complex<T>::random_real() {
	return Complex<T>( (static_cast<T>( rand() ))/RAND_MAX, 0.0 );
}

/******************************************
 * Binary Arithmetic Operators
 ******************************************/

// (a + bi) + (c + di) -> (a + c) + (b + d)i

template <typename T> inline Complex<T> Complex<T>::operator + ( const Complex<T> & z ) const {
	return Complex<T>( r + z.r, i + z.i );
}

// (a + bi) + x -> (a + x) + bi

template <typename T> inline Complex<T> Complex<T>::operator + ( T x ) const {
	return Complex<T>( r + x, i );
}

// x + (a + bi) -> (a + x) + bi

template <typename T> Complex<T> operator + ( T x, const Complex<T> & z ) {
	return Complex<T>( x + z.real(), z.imag_i() );
}

template <typename T> Complex<T> operator + ( long x, const Complex<T> & z ) {
	return Complex<T>( static_cast<T>( x ) + z.real(), z.imag_i() );
}

// (a + bi) - (c + di) -> (a - c) + (b - d)i

template <typename T> inline Complex<T> Complex<T>::operator - ( const Complex<T> & z ) const {
	return Complex<T>( r - z.r, i - z.i );
}

// (a + bi) - x -> (a - x) + bi

template <typename T> inline Complex<T> Complex<T>::operator - ( T x ) const {
	return Complex<T>( r - x, i );
}

// x - (a + bi) -> (x + a) - bi

template <typename T> Complex<T> operator - ( T x, const Complex<T> & z ) {
	return Complex<T>( x - z.real(), -z.imag_i() );
}

template <typename T> Complex<T> operator - ( long x, const Complex<T> & z ) {
	return Complex<T>( static_cast<T>( x ) - z.real(), -z.imag_i() );
}

template <typename T> inline Complex<T> Complex<T>::operator * ( const Complex<T> & z ) const {
	return Complex<T>( r*z.r - i*z.i, r*z.i + i*z.r );
}

template <typename T> inline Complex<T> Complex<T>::operator * ( T x ) const {
	if ( Support<T>::is_inf( x ) && norm() > 0 ) {
		return Complex<T>(
			( (  r == 0 )?Support<T>::sign(x):x )*r,
			( (  i == 0 )?Support<T>::sign(x):x )*i
		);
	} else {
		return Complex<T>( x*r, x*i );
	}
}

template <typename T> Complex<T> operator * ( T x, const Complex<T> & z ) {
	return Complex<T>( x*z.real(), x*z.imag_i() );
}

template <typename T> Complex<T> operator * ( long x, const Complex<T> & z ) {
	return Complex<T>( static_cast<T>( x )*z.real(), static_cast<T>( x )*z.imag_i() );
}

template <typename T> inline Complex<T> Complex<T>::operator / ( const Complex<T> & z ) const {
	T denom = z.norm();
	return Complex<T>( (r*z.r + i*z.i)/denom, (-r*z.i + i*z.r)/denom );
}

template <typename T> inline Complex<T> Complex<T>::operator / ( T x ) const {
	if ( ( x == 0.0 ) && ( norm() > 0 ) ) {
		return Complex<T>(
			r / ( ( r == 0 )?Support<T>::sign( x ):x ),
			i / ( ( i == 0 )?Support<T>::sign( x ):x )
		);
	} else {
		return Complex<T>( r/x, i/x );
	}
}

template <typename T> Complex<T> operator / ( T x, const Complex<T> & z ) {
	T mltplr = x/z.norm();
	return Complex<T>( mltplr*z.real(), -mltplr*z.imag_i() );
}

template <typename T> Complex<T> operator / ( long x, const Complex<T> & z ) {
	T mltplr = static_cast<T>( x )/z.norm();
	return Complex<T>( mltplr*z.real(), -mltplr*z.imag_i() );
}

/******************************************
 * Unary Arithmetic Operators
 ******************************************/

// -(a + bi) -> -a - bi

template <typename T> inline Complex<T> Complex<T>::operator - () const {
	return Complex<T>( -r, -i );
}

/******************************************
 * Binary Boolean Operators
 ******************************************/

// (a + bi) == (c + di) -> (a == c) && (b == d)

template <typename T> inline bool Complex<T>::operator == ( const Complex<T> & z ) const {
	return ( r == z.r ) && ( i == z.i );
}

// (a + bi) == x -> (a == x) && (b == 0)

template <typename T> inline bool Complex<T>::operator == ( T x ) const {
	return ( r == x ) && ( i == 0.0 );
}

// (a + bi) != (c + di) -> (a != c) || (b != d)

template <typename T> inline bool Complex<T>::operator != ( const Complex<T> & z ) const {
	return ( r != z.r ) || ( i != z.i );
}

// (a + bi) != x -> (a != x) || (b != 0)

template <typename T> inline bool Complex<T>::operator != ( T x ) const {
	return ( r != x ) || ( i != 0.0 );
}

// x == (a + bi) -> (x == a) && (b == 0)

template <typename T> bool operator == ( T x, const Complex<T> & z ) {
	return ( z.real() == x ) && ( z.imag_i() == 0.0 );
}

template <typename T> bool operator == ( long x, const Complex<T> & z ) {
	return ( z.real() == static_cast<T>( x ) ) && ( z.imag_i() == 0.0 );
}

// x != (a + bi) -> (x != a) || (b != 0)

template <typename T> bool operator != ( T x, const Complex<T> & z ) {
	return ( z.real() != x ) || ( z.imag_i() != 0.0 );
}

template <typename T> bool operator != ( long x, const Complex<T> & z ) {
	return ( z.real() != static_cast<T>( x ) ) || ( z.imag_i() != 0.0 );
}

template<typename T> std::ostream & operator << ( std::ostream & out, const Complex<T> & z ) {
	Support<T>::print_real( z.real(), out );
	Support<T>::print_imaginary( z.imag_i(), Complex<T>::get_symbol(), out );

	return out;
}

/******************************************
 * ************************************** *
 * *                                    * *
 * *      Procedural Functions          * *
 * *                                    * *
 * ************************************** *
 ******************************************/

/******************************************
 * Real-valued Functions
 ******************************************/

template <typename T> T real( const Complex<T> & z ) {
	return z.real();
}

template <typename T> T imag_i( const Complex<T> & z ) {
	return z.imag_i();
}

template <typename T> T csgn( const Complex<T> & z ) {
	return z.csgn();
}

template <typename T> T abs( const Complex<T> & z ) {
	return z.abs();
}

template <typename T> T norm( const Complex<T> & z ) {
	return z.norm();
}

template <typename T> T abs_imag( const Complex<T> & z ) {
	return z.abs_imag();
}

template <typename T> T norm_imag( const Complex<T> & z ) {
	return z.norm_imag();
}

/******************************************
 * Complex-valued Functions
 ******************************************/

template <typename T> Complex<T> imag( const Complex<T> & z ) {
	return z.imag();
}

template <typename T> Complex<T> conj( const Complex<T> & z ) {
	return z.conj();
}

template <typename T> Complex<T> signum( const Complex<T> & z ) {
	return z.signum();
}

template <typename T> Complex<T> sqr( const Complex<T> & z ) {
	return z.sqr();
}

template <typename T> Complex<T> sqrt( const Complex<T> & z ) {
	return z.sqrt();
}

/******************************************
 * Exponential and Logarithmic Functions
 ******************************************/

template <typename T> Complex<T> exp( const Complex<T> & z ) {
	return z.exp();
}

template <typename T> Complex<T> log( const Complex<T> & z ) {
	return z.log();
}

template <typename T> Complex<T> log10( const Complex<T> & z ) {
	return z.log10();
}

template <typename T> Complex<T> pow( const Complex<T> & z, const Complex<T> & w ) {
	return z.pow( w );
}

template <typename T> Complex<T> pow( const Complex<T> & z, T x ) {
	return z.pow( x );
}

template <typename T> Complex<T> inverse( const Complex<T> & z ) {
	return z.inverse();
}

/******************************************
 * Trigonometric and Hyperbolic Functions
 ******************************************/

template <typename T> Complex<T> sin( const Complex<T> & z ) {
	return z.sin();
}

template <typename T> Complex<T> cos( const Complex<T> & z ) {
	return z.cos();
}

template <typename T> Complex<T> tan( const Complex<T> & z ) {
	return z.tan();
}

template <typename T> Complex<T> sec( const Complex<T> & z ) {
	return z.sec();
}

template <typename T> Complex<T> csc( const Complex<T> & z ) {
	return z.csc();
}

template <typename T> Complex<T> cot( const Complex<T> & z ) {
	return z.cot();
}

template <typename T> Complex<T> sinh( const Complex<T> & z ) {
	return z.sinh();
}

template <typename T> Complex<T> cosh( const Complex<T> & z ) {
	return z.cosh();
}

template <typename T> Complex<T> tanh( const Complex<T> & z ) {
	return z.tanh();
}

template <typename T> Complex<T> sech( const Complex<T> & z ) {
	return z.sech();
}

template <typename T> Complex<T> csch( const Complex<T> & z ) {
	return z.csch();
}

template <typename T> Complex<T> coth( const Complex<T> & z ) {
	return z.coth();
}

template <typename T> Complex<T> asin( const Complex<T> & z ) {
	return z.asin();
}

template <typename T> Complex<T> acos( const Complex<T> & z ) {
	return z.acos();
}

template <typename T> Complex<T> atan( const Complex<T> & z ) {
	return z.atan();
}

template <typename T> Complex<T> asec( const Complex<T> & z ) {
	return z.asec();
}

template <typename T> Complex<T> acsc( const Complex<T> & z ) {
	return z.acsc();
}

template <typename T> Complex<T> acot( const Complex<T> & z ) {
	return z.acot();
}

template <typename T> Complex<T> asinh( const Complex<T> & z ) {
	return z.asinh();
}

template <typename T> Complex<T> acosh( const Complex<T> & z ) {
	return z.acosh();
}

template <typename T> Complex<T> atanh( const Complex<T> & z ) {
	return z.atanh();
}

template <typename T> Complex<T> asech( const Complex<T> & z ) {
	return z.asech();
}

template <typename T> Complex<T> acsch( const Complex<T> & z ) {
	return z.acsch();
}

template <typename T> Complex<T> acoth( const Complex<T> & z ) {
	return z.acoth();
}

template <typename T> Complex<T> bessel_J( int n, const Complex<T> & z ) {
	return z.bessel_J( n );
}

/******************************************
 * Integer Functions
 ******************************************/

template <typename T> Complex<T> floor( const Complex<T> & z ) {
	return z.floor();
}

template <typename T> Complex<T> ceil( const Complex<T> & z ) {
	return z.ceil();
}

/******************************************
 * Horner's Rule
 ******************************************/

template <typename T> Complex<T> horner( const Complex<T> & z, T * v, unsigned int n ) {
	return z.horner( v, n );
}

template <typename T> Complex<T> horner( const Complex<T> & z, T * v, T * c, unsigned int n ) {
	return z.horner( v, c, n );
}

template <typename T> Complex<T> horner( const Complex<T> & z, Complex<T> * v, unsigned int n ) {
	return z.horner( v, n );
}

template <typename T> Complex<T> horner( const Complex<T> & z, Complex<T> * v, T * c, unsigned int n ) {
	return z.horner( v, c, n );
}

template <typename T> Complex<T> horner( const Complex<T> & z, T * v, Complex<T> * c, unsigned int n ) {
	return z.horner( v, c, n );
}

template <typename T> Complex<T> horner( const Complex<T> & z, Complex<T> * v, Complex<T> * c, unsigned int n ) {
	return z.horner( v, c, n );
}

/**************************************************
 * ********************************************** *
 * *                                            * *
 * *    Double-precision Floating-point         * *
 * *    Instance of Template                    * *
 * *                                            * *
 * ********************************************** *
 **************************************************/

template class Complex<double>;

template <> const Complex<double> Complex<double>::ZERO = Complex<double>( 0, 0 );
template <> const Complex<double> Complex<double>::ONE  = Complex<double>( 1, 0 );
template <> const Complex<double> Complex<double>::I    = Complex<double>( 0, 1 );

template <> const Complex<double> Complex<double>::UNITS[2] = {
	Complex<double>::ONE,
	Complex<double>::I
};

template <> char Complex<double>::imaginary_symbol = 'i';

template std::ostream & operator << ( std::ostream & out, const Complex<double> & );

template Complex<double> operator + ( double, const Complex<double> & );
template Complex<double> operator + ( long, const Complex<double> & );
template Complex<double> operator - ( double, const Complex<double> & );
template Complex<double> operator - ( long, const Complex<double> & );
template Complex<double> operator * ( double, const Complex<double> & );
template Complex<double> operator * ( long, const Complex<double> & );
template Complex<double> operator / ( double, const Complex<double> & );
template Complex<double> operator / ( long, const Complex<double> & );

template bool operator == ( double, const Complex<double> & );
template bool operator == ( long, const Complex<double> & );
template bool operator != ( double, const Complex<double> & );
template bool operator != ( long, const Complex<double> & );

template double real( const Complex<double> & );
template double imag_i( const Complex<double> & );
template double csgn( const Complex<double> & );
template double abs( const Complex<double> & );
template double norm( const Complex<double> & );
template double abs_imag( const Complex<double> & );
template double norm_imag( const Complex<double> & );
template Complex<double> imag( const Complex<double> & );
template Complex<double> conj( const Complex<double> & );
template Complex<double> signum( const Complex<double> & );
template Complex<double> sqr( const Complex<double> & );
template Complex<double> sqrt( const Complex<double> & );
template Complex<double> exp( const Complex<double> & );
template Complex<double> log( const Complex<double> & );
template Complex<double> log10( const Complex<double> & );
template Complex<double> pow( const Complex<double> &, const Complex<double> & );
template Complex<double> pow( const Complex<double> &, double );
template Complex<double> inverse( const Complex<double> & );
template Complex<double> sin( const Complex<double> & );
template Complex<double> cos( const Complex<double> & );
template Complex<double> tan( const Complex<double> & );
template Complex<double> sec( const Complex<double> & );
template Complex<double> csc( const Complex<double> & );
template Complex<double> cot( const Complex<double> & );
template Complex<double> sinh( const Complex<double> & );
template Complex<double> cosh( const Complex<double> & );
template Complex<double> tanh( const Complex<double> & );
template Complex<double> sech( const Complex<double> & );
template Complex<double> csch( const Complex<double> & );
template Complex<double> coth( const Complex<double> & );
template Complex<double> asin( const Complex<double> & );
template Complex<double> acos( const Complex<double> & );
template Complex<double> atan( const Complex<double> & );
template Complex<double> asec( const Complex<double> & );
template Complex<double> acsc( const Complex<double> & );
template Complex<double> acot( const Complex<double> & );
template Complex<double> asinh( const Complex<double> & );
template Complex<double> acosh( const Complex<double> & );
template Complex<double> atanh( const Complex<double> & );
template Complex<double> asech( const Complex<double> & );
template Complex<double> acsch( const Complex<double> & );
template Complex<double> acoth( const Complex<double> & );
template Complex<double> bessel_J( int, const Complex<double> & );
template Complex<double> floor( const Complex<double> & );
template Complex<double> ceil( const Complex<double> & );
template Complex<double> horner( const Complex<double> &, double *, unsigned int );
template Complex<double> horner( const Complex<double> &, double *, double *, unsigned int );
template Complex<double> horner( const Complex<double> &, Complex<double> *, unsigned int );
template Complex<double> horner( const Complex<double> &, Complex<double> *, double *, unsigned int );
template Complex<double> horner( const Complex<double> &, double *, Complex<double> *, unsigned int );
template Complex<double> horner( const Complex<double> &, Complex<double> *, Complex<double> *, unsigned int );

/**************************************************
 * ********************************************** *
 * *                                            * *
 * *    Floating-point Instance of Template     * *
 * *                                            * *
 * ********************************************** *
 **************************************************/

template class Complex<float>;

template <> const Complex<float> Complex<float>::ZERO = Complex<float>( 0, 0 );
template <> const Complex<float> Complex<float>::ONE  = Complex<float>( 1, 0 );
template <> const Complex<float> Complex<float>::I    = Complex<float>( 0, 1 );

template <> const Complex<float> Complex<float>::UNITS[2] = {
	Complex<float>::ONE,
	Complex<float>::I
};

template <> char Complex<float>::imaginary_symbol = 'i';

template std::ostream & operator << ( std::ostream & out, const Complex<float> & );

template Complex<float> operator + ( float, const Complex<float> & );
template Complex<float> operator + ( long, const Complex<float> & );
template Complex<float> operator - ( float, const Complex<float> & );
template Complex<float> operator - ( long, const Complex<float> & );
template Complex<float> operator * ( float, const Complex<float> & );
template Complex<float> operator * ( long, const Complex<float> & );
template Complex<float> operator / ( float, const Complex<float> & );
template Complex<float> operator / ( long, const Complex<float> & );

template bool operator == ( float, const Complex<float> & );
template bool operator == ( long, const Complex<float> & );
template bool operator != ( float, const Complex<float> & );
template bool operator != ( long, const Complex<float> & );

template float real( const Complex<float> & );
template float imag_i( const Complex<float> & );
template float csgn( const Complex<float> & );
template float abs( const Complex<float> & );
template float norm( const Complex<float> & );
template float abs_imag( const Complex<float> & );
template float norm_imag( const Complex<float> & );
template Complex<float> imag( const Complex<float> & );
template Complex<float> conj( const Complex<float> & );
template Complex<float> signum( const Complex<float> & );
template Complex<float> sqr( const Complex<float> & );
template Complex<float> sqrt( const Complex<float> & );
template Complex<float> exp( const Complex<float> & );
template Complex<float> log( const Complex<float> & );
template Complex<float> log10( const Complex<float> & );
template Complex<float> pow( const Complex<float> &, const Complex<float> & );
template Complex<float> pow( const Complex<float> &, float );
template Complex<float> inverse( const Complex<float> & );
template Complex<float> sin( const Complex<float> & );
template Complex<float> cos( const Complex<float> & );
template Complex<float> tan( const Complex<float> & );
template Complex<float> sec( const Complex<float> & );
template Complex<float> csc( const Complex<float> & );
template Complex<float> cot( const Complex<float> & );
template Complex<float> sinh( const Complex<float> & );
template Complex<float> cosh( const Complex<float> & );
template Complex<float> tanh( const Complex<float> & );
template Complex<float> sech( const Complex<float> & );
template Complex<float> csch( const Complex<float> & );
template Complex<float> coth( const Complex<float> & );
template Complex<float> asin( const Complex<float> & );
template Complex<float> acos( const Complex<float> & );
template Complex<float> atan( const Complex<float> & );
template Complex<float> asec( const Complex<float> & );
template Complex<float> acsc( const Complex<float> & );
template Complex<float> acot( const Complex<float> & );
template Complex<float> asinh( const Complex<float> & );
template Complex<float> acosh( const Complex<float> & );
template Complex<float> atanh( const Complex<float> & );
template Complex<float> asech( const Complex<float> & );
template Complex<float> acsch( const Complex<float> & );
template Complex<float> acoth( const Complex<float> & );
template Complex<float> bessel_J( int, const Complex<float> & );
template Complex<float> floor( const Complex<float> & );
template Complex<float> ceil( const Complex<float> & );
template Complex<float> horner( const Complex<float> &, float *, unsigned int );
template Complex<float> horner( const Complex<float> &, float *, float *, unsigned int );
template Complex<float> horner( const Complex<float> &, Complex<float> *, unsigned int );
template Complex<float> horner( const Complex<float> &, Complex<float> *, float *, unsigned int );
template Complex<float> horner( const Complex<float> &, float *, Complex<float> *, unsigned int );
template Complex<float> horner( const Complex<float> &, Complex<float> *, Complex<float> *, unsigned int );

/************************************************************************
 * ******************************************************************** *
 * *                                                                  * *
 * *    Long Double-Precision Floating-point Instance of Template     * *
 * *                                                                  * *
 * ******************************************************************** *
 ************************************************************************/

template class Complex<long double>;

template <> const Complex<long double> Complex<long double>::ZERO = Complex<long double>( 0, 0 );
template <> const Complex<long double> Complex<long double>::ONE  = Complex<long double>( 1, 0 );
template <> const Complex<long double> Complex<long double>::I    = Complex<long double>( 0, 1 );

template <> const Complex<long double> Complex<long double>::UNITS[2] = {
	Complex<long double>::ONE,
	Complex<long double>::I
};

template <> char Complex<long double>::imaginary_symbol = 'i';

template std::ostream & operator << ( std::ostream & out, const Complex<long double> & );

template Complex<long double> operator + ( long double, const Complex<long double> & );
template Complex<long double> operator + ( long, const Complex<long double> & );
template Complex<long double> operator - ( long double, const Complex<long double> & );
template Complex<long double> operator - ( long, const Complex<long double> & );
template Complex<long double> operator * ( long double, const Complex<long double> & );
template Complex<long double> operator * ( long, const Complex<long double> & );
template Complex<long double> operator / ( long double, const Complex<long double> & );
template Complex<long double> operator / ( long, const Complex<long double> & );

template bool operator == ( long double, const Complex<long double> & );
template bool operator == ( long, const Complex<long double> & );
template bool operator != ( long double, const Complex<long double> & );
template bool operator != ( long, const Complex<long double> & );

template long double real( const Complex<long double> & );
template long double imag_i( const Complex<long double> & );
template long double csgn( const Complex<long double> & );
template long double abs( const Complex<long double> & );
template long double norm( const Complex<long double> & );
template long double abs_imag( const Complex<long double> & );
template long double norm_imag( const Complex<long double> & );
template Complex<long double> imag( const Complex<long double> & );
template Complex<long double> conj( const Complex<long double> & );
template Complex<long double> signum( const Complex<long double> & );
template Complex<long double> sqr( const Complex<long double> & );
template Complex<long double> sqrt( const Complex<long double> & );
template Complex<long double> exp( const Complex<long double> & );
template Complex<long double> log( const Complex<long double> & );
template Complex<long double> log10( const Complex<long double> & );
template Complex<long double> pow( const Complex<long double> &, const Complex<long double> & );
template Complex<long double> pow( const Complex<long double> &, long double );
template Complex<long double> inverse( const Complex<long double> & );
template Complex<long double> sin( const Complex<long double> & );
template Complex<long double> cos( const Complex<long double> & );
template Complex<long double> tan( const Complex<long double> & );
template Complex<long double> sec( const Complex<long double> & );
template Complex<long double> csc( const Complex<long double> & );
template Complex<long double> cot( const Complex<long double> & );
template Complex<long double> sinh( const Complex<long double> & );
template Complex<long double> cosh( const Complex<long double> & );
template Complex<long double> tanh( const Complex<long double> & );
template Complex<long double> sech( const Complex<long double> & );
template Complex<long double> csch( const Complex<long double> & );
template Complex<long double> coth( const Complex<long double> & );
template Complex<long double> asin( const Complex<long double> & );
template Complex<long double> acos( const Complex<long double> & );
template Complex<long double> atan( const Complex<long double> & );
template Complex<long double> asec( const Complex<long double> & );
template Complex<long double> acsc( const Complex<long double> & );
template Complex<long double> acot( const Complex<long double> & );
template Complex<long double> asinh( const Complex<long double> & );
template Complex<long double> acosh( const Complex<long double> & );
template Complex<long double> atanh( const Complex<long double> & );
template Complex<long double> asech( const Complex<long double> & );
template Complex<long double> acsch( const Complex<long double> & );
template Complex<long double> acoth( const Complex<long double> & );
template Complex<long double> bessel_J( int, const Complex<long double> & );
template Complex<long double> floor( const Complex<long double> & );
template Complex<long double> ceil( const Complex<long double> & );
template Complex<long double> horner( const Complex<long double> &, long double *, unsigned int );
template Complex<long double> horner( const Complex<long double> &, long double *, long double *, unsigned int );
template Complex<long double> horner( const Complex<long double> &, Complex<long double> *, unsigned int );
template Complex<long double> horner( const Complex<long double> &, Complex<long double> *, long double *, unsigned int );
template Complex<long double> horner( const Complex<long double> &, long double *, Complex<long double> *, unsigned int );
template Complex<long double> horner( const Complex<long double> &, Complex<long double> *, Complex<long double> *, unsigned int );