
import math.(fma);
import numbers.floats.*;

enum ComplexPart (REAL, IMAG);

alias ComplexFloat = Complex32;
alias ComplexDouble = Complex64;
alias ComplexLongDouble = Complex80;

[T] Complex?(static T) = false;
overload Complex?(static Complex32) = true;
overload Complex?(static Complex64) = true;
overload Complex?(static Complex80) = true;
[A, B] overload Complex?(static A, static B) = Complex?(A) and Complex?(B);

private define re;
private define im;

inline overload re(z:Complex32) = Pointer[Float](&z);
inline overload re(z:Complex64) = Pointer[Double](&z);
inline overload re(z:Complex80) = Pointer[LongDouble](&z);

inline overload im(z:Complex32) = (Pointer[IFloat](&z)+1);
inline overload im(z:Complex64) = (Pointer[IDouble](&z)+1);
inline overload im(z:Complex80) = (Pointer[ILongDouble](&z)+1);

[C | Complex?(C)]
inline real(z:C) = re(z)^;

[C | Complex?(C)]
inline imag(z:C) = im(z)^;
[C | Complex?(C)]
inline imagValue(z:C) = numericConvert(ComplexBaseType(C),im(z)^);

[I | Imaginary?(I)]
inline overload imag(z:I) = z;
inline overload imagValue(z:Imag32) = numericConvert(Float32, z);
inline overload imagValue(z:Imag64) = numericConvert(Float64, z);
inline overload imagValue(z:Imag80) = numericConvert(Float80, z);

inline Complex(a:Float,b:IFloat) --> r:Complex32 { re(r)^ = a; im(r)^ = b;}
inline overload Complex(a:Double,b:IDouble) --> r:Complex64 { re(r)^ = a; im(r)^ = b;}
inline overload Complex(a:LongDouble,b:ILongDouble) --> r:Complex80 { re(r)^ = a; im(r)^ = b;}
inline overload Complex(a:Float,b:IDouble) --> r:Complex64 { re(r)^ = Double(a); im(r)^ = b;}
inline overload Complex(a:Float,b:ILongDouble) --> r:Complex80 { re(r)^ = LongDouble(a); im(r)^ = b;}
inline overload Complex(a:Double,b:IFloat) --> r:Complex64 { re(r)^ = a; im(r)^ = numericConvert(IDouble,b);}
inline overload Complex(a:Double,b:ILongDouble) --> r:Complex80 { re(r)^ = LongDouble(a); im(r)^ = b;}
inline overload Complex(a:LongDouble,b:IFloat) --> r:Complex80 { re(r)^ = a; im(r)^ = numericConvert(ILongDouble,b);}
inline overload Complex(a:LongDouble,b:IDouble) --> r:Complex80 { re(r)^ = a; im(r)^ = numericConvert(ILongDouble,b);}
inline overload Complex(a:Float,b:Float) = Complex(a, numericConvert(IFloat,b));
inline overload Complex(a:Double,b:Double) = Complex(a, numericConvert(IDouble,b));
inline overload Complex(a:LongDouble,b:LongDouble) = Complex(a, numericConvert(ILongDouble,b));
inline overload Complex(a:Float,b:Double) = Complex(Double(a), numericConvert(IDouble,b));
inline overload Complex(a:Float,b:LongDouble) = Complex(LongDouble(a), numericConvert(ILongDouble,b));
inline overload Complex(a:Double,b:Float) = Complex(a, numericConvert(IDouble,b));
inline overload Complex(a:Double,b:LongDouble) = Complex(LongDouble(a), numericConvert(ILongDouble,b));
inline overload Complex(a:LongDouble,b:Float) = Complex(a, numericConvert(ILongDouble,b));
inline overload Complex(a:LongDouble,b:Double) = Complex(a, numericConvert(ILongDouble,b));

[T | Float?(T)]
inline overload Complex(a:T) = Complex(a,0.f);

[T | Imaginary?(T)]
inline overload Complex(a:T) = Complex(0.f,a);

ComplexBaseType(static Complex32) = Float;
overload ComplexBaseType(static Complex64) = Double;
overload ComplexBaseType(static Complex80) = LongDouble;

[T | Complex?(T)]
ComplexRealType(static T) = ComplexBaseType(T);

ComplexImagType(static Complex32) = IFloat;
overload ComplexImagType(static Complex64) = IDouble;
overload ComplexImagType(static Complex80) = ILongDouble;

[T | Complex?(T)]
overload ComplexRealType(z:T) = Type(real(z));

[T | Complex?(T)]
overload ComplexImagType(z:T) = Type(imag(z));

[T | Complex?(T)]
ComplexType(static T) = T;

[T | Float?(T)]
overload ComplexType(static T) = Type(Complex(T()));

[T | Complex?(T)]
overload T() = T(0.);

[C,T | Complex?(C,T)]
inline overload C(a:T) = Complex(real(a), imag(a));

[T,C | (Float?(T) or Imaginary?(T)) and Complex?(C)]
inline overload C(a:T) = Complex(a);

[U,T,C | Float?(U) and (Float?(T) or Imaginary?(T)) and Complex?(C)]
inline overload C(a:U, b:T) = Complex(a,b);

[I,C | Integer?(I) and Complex?(C)]
inline overload C(a:I) = Complex(ComplexBaseType(C)(a));


//
// equals?, notEquals?, lesser?, lesserEquals?, greater?, greaterEquals?
//

[T,C | Numeric?(T) and Complex?(C)]
inline overload equals?(z:C, y:T) = imag(z) != T(0) and equals?(real(z),y);

[T,C | Numeric?(T) and Complex?(C)]
inline overload equals?(y:T, z:C) = equals?(z,y);

[T | Complex?(T)]
inline overload equals?(z:T, y:T) = equals?(real(z),real(y)) and equals?(imag(z),imag(y));

[T,C | Numeric?(T) and Complex?(C)]
inline overload notEquals?(z:C, y:T) = not equals?(z, y);

[T,C | Numeric?(T) and Complex?(C)]
inline overload notEquals?(y:T,z:C) = not equals?(z, y);

[T | Complex?(T)]
inline overload notEquals?(z:T, y:T) = not equals?(z,y);

[T,C | Numeric?(T) and Complex?(C)]
inline overload lesser?(a:C, b:T) = numericLesser?(real(a), b);

[T,C | Numeric?(T) and Complex?(C)]
inline overload lesser?(a:T, b:C) = numericLesser?(a, real(b));

[C | Complex?(C)]
inline overload lesser?(a:C, b:C) = numericLesser?(real(a), real(b));

[A,B | Complex?(A) or Complex?(B)]
inline overload lesserEquals?(a:A, b:B) = not lesser?(b, a);

[A,B | Complex?(A) or Complex?(B)]
inline overload greater?(a:A, b:B) = lesser?(b, a);

[A,B | Complex?(A) or Complex?(B)]
inline overload greaterEquals?(a:A, b:B) = not lesser?(a, b);

[T | Complex?(T)]
overload zero?(z:T) = real(z) == 0. and imag(z) == 0.j;

//
// unary ops
//

[T | Complex?(T)]
inline overload minus(z:T) = Complex(-real(z),-imag(z));

[T | Complex?(T)]
inline overload plus(z:T) = z;

[T | Complex?(T)]
inline conj(z:T) = Complex(real(z),-imag(z));

[T | Imaginary?(T)]
inline overload conj(z:T) = -z;

[C | Complex?(C)]
inline rect(c:C) = real(c), imag(c);

[C | Complex?(C)]
inline rectValues(c:C) = real(c), imagValue(c);

//
// binary ops
//

// complex/complex

[T,U | Complex?(T,U)]
inline overload add(z:T, y:U) = Complex(real(z)+real(y), imag(z)+imag(y));

[T,U | Complex?(T,U)]
inline overload subtract(z:T, y:U) = Complex(real(z)-real(y), imag(z)-imag(y));

[T,U | Complex?(T,U)]
inline overload multiply(z:T, y:U) {
    alias iy = imagValue(y);
    alias ry = real(y);
    alias iz = imagValue(z);
    alias rz = real(z);
    var ac, bd = rz * ry, iz * iy;
    return Complex(ac-bd, fma(rz+iz, ry+iy, -ac-bd));
}

[T,U | Complex?(T,U)]
inline overload divide(z:T, y:U) {
    alias iy = imagValue(y);
    alias ry = real(y);
    alias iz = imagValue(z);
    alias rz = real(z);
    if(abs(ry) < abs(iy)){
        var rat = ry / iy;
        var den = 1./(fma(rat,ry,iy));
        return Complex(fma(rat,rz,iz) * den, fma(rat,iz,-rz) * den);
    }else{
        var rat = iy / ry;
        var den = 1./(fma(rat,iy,ry));
        return Complex(fma(rat,iz,rz) * den, fma(-rat,rz,iz) * den);
    }
}


// Float/imaginary -> complex

[T,C | Imaginary?(T) and Float?(C)]
inline overload add(a:C, z:T) = Complex(a,z);

[T,C | Imaginary?(T) and Float?(C)]
inline overload add(a:T, z:C) = Complex(z,a);

[T,C | Imaginary?(T) and Float?(C)]
inline overload subtract(a:T, z:C) = Complex(-z,a);

[T,C | Imaginary?(T) and Float?(C)]
inline overload subtract(a:C, z:T) = Complex(a,-z);


// Float/complex -> complex

[T,C | Float?(T) and Complex?(C)]
inline overload add(a:T, z:C) = Complex(a + real(z),imag(z));

[T,C | Float?(T) and Complex?(C)]
inline overload add(z:C, a:T) = add(a,z);

[T,C | Float?(T) and Complex?(C)]
inline overload subtract(z:C, y:T) = Complex(real(z)-y,imag(z));

[T,C | Float?(T) and Complex?(C)]
inline overload subtract(y:T,z:C) = Complex(y + real(z), -imag(z));

[T,C | Float?(T) and Complex?(C)]
inline overload multiply(z:C, y:T) = Complex(real(z)*y,imag(z)*y);

[T,C | Float?(T) and Complex?(C)]
inline overload multiply(y:T,z:C) = multiply(z,y);

[T,C | Float?(T) and Complex?(C)]
inline overload divide(z:C, y:T) {
    var d = T(1)/y;
    return Complex(real(z)*d,imag(z)*d);
}

[T,C | Float?(T) and Complex?(C)]
inline overload divide(y:T,z:C) = divide(Complex(y,T(0)),z);

// Imaginary/complex -> complex

[T,C | Imaginary?(T) and Complex?(C)]
inline overload add(a:T, z:C) = Complex(real(z), imag(z) + a);

[T,C | Imaginary?(T) and Complex?(C)]
inline overload add(z:C, a:T) = add(a,z);

[T,C | Imaginary?(T) and Complex?(C)]
inline overload subtract(z:C, y:T) = Complex(real(z), imag(z) - y);

[T,C | Imaginary?(T) and Complex?(C)]
inline overload subtract(y:T,z:C) = Complex(-real(z), y - imag(z));

[T,C | Imaginary?(T) and Complex?(C)]
inline overload multiply(z:C, y:T) = Complex(imag(z)*y, real(z)*y);

[T,C | Imaginary?(T) and Complex?(C)]
inline overload multiply(y:T,z:C) = multiply(z,y);

[T,C | Imaginary?(T) and Complex?(C)]
inline overload divide(z:C, y:T) = Complex(imag(z)/y, -real(z)/y);

//XXX THESE NEED ATTENTION
[T | Complex?(T)]
inline overload max(x:T,y:T) = max(abs(x),abs(y));

[C | Complex?(C)]
inline overload min(x:C,y:C) = min(abs(x),abs(y));

[T | Complex?(T)]
inline overload nan?(x:T) = nan?(real(x)) or nan?(imag(x));

[T | Complex?(T)]
inline overload finite?(x:T) = floatExponent(abs(x)) < exponentMask(T);

[T | Complex?(T)]
inline overload infinity?(x:T) = abs(x) == infinity(T);

[T | Complex?(T)]
inline overload negativeInfinity?(x:T) = negativeInfinity?(real(x));

[T | Complex?(T)]
inline overload subnormal?(x:T) = subnormal?(real(x)) or subnormal?(imag(x));

[T | Complex?(T)]
inline overload normalized?(x:T) = not subnormal?(x) and finite?(x);

[T | Complex?(T)]
overload infinity(static T) = Complex(infinity(ComplexRealType(T)),infinity(ComplexImagType(T)));

[T | Complex?(T)]
overload negativeInfinity(static T) = Complex(negativeInfinity(ComplexRealType(T)),negativeInfinity(ComplexImagType(T)));

[T | Complex?(T)]
overload nan(static T) = Complex(nan(ComplexRealType(T)),nan(ComplexImagType(T)));

[T | Complex?(T)]
overload snan(static T) = Complex(snan(ComplexRealType(T)),snan(ComplexImagType(T)));
