
import math.*;
import numbers.floats.*;

[U | Complex?(U)]
inline overload abs(x:U) = hypot(real(x),imagValue(x));

[T,C | Numeric?(T) and Complex?(C)]
inline overload max(x:T,y:C) = max(x,abs(y));

[T,C | Numeric?(T) and Complex?(C)]
inline overload max(x:C,y:T) = max(abs(x),y);

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

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

[C | Complex?(C)]
overload sqrt(x:C) --> res:C {
    alias T = ComplexBaseType(C);
    alias ZERO = T(0);
    alias HALF = T(0.5);
    alias HUGE_VAL = GreatestPositiveFinite(T);
    alias xr = real(x);
    alias xi = imagValue(x);
    var rcls,icls = fpclassify(xr), fpclassify (xi);
    if (rcls <= FP_INFINITE or icls <= FP_INFINITE){
        if (icls == FP_INFINITE){
            res <-- Complex(HUGE_VAL,xi);
        }
        else if (rcls == FP_INFINITE){
            if (xr < ZERO){
                res <-- Complex(if(icls == FP_NAN) nan(T) else ZERO, copysign (HUGE_VAL, xi));
            }
            else{
                res <-- Complex(xr,if(icls == FP_NAN) nan(T) else copysign (ZERO, xi));
            }
        }
        else {
                res <-- Complex(nan(T),nan(T));
        }
    }
    else{
        if (icls == FP_ZERO){
            if (xr < ZERO){
                res <-- Complex(ZERO, copysign(sqrt(-xr),xi));
            }
            else{
                res <-- Complex(abs(sqrt(xr)), copysign (ZERO, xi));
            }
        }
        else if (rcls == FP_ZERO){
            var r = sqrt(HALF * abs (xi));
            res <-- Complex(r, copysign(r, xi));
        }
        else {
            var d,r,s = hypot(xr, xi),ZERO,ZERO;
            /* Use the identity   2  Re res  Im res = Im x
            to avoid cancellation error in  d +/- Re x.  */
            if (xr > ZERO){
                r = sqrt(HALF * d + HALF * xr);
                s = (HALF * xi) / r;
            }
            else{
                s = sqrt(HALF * d - HALF * xr);
                r = abs((HALF * xi) / s);
            }
            res <-- Complex(r,copysign(s, xi));
        }
    }
}


[I | Imaginary?(I)]
inline overload sqrt(z:I) = sqrt(Complex(z));

[T | Float?(T)]
inline fromPolar(m:T,ph:T) = Complex(m*sin(ph),m*cos(ph));

[T | Complex?(T)]
inline overload arg(x:T) = atan2(imagValue(x),real(x));

