Implementation of numbers in this paper: http://people.rit.edu/harkin/research/articles/generalized_complex_numbers.pdf

GeneralComplex.jl

# This file is a naive implementation of numbers of the form x + y*i where i^2 = p

#######################
# GeneralComplex type #
#######################

immutable GeneralComplex{p,T} <: Number
    x::T
    y::T
end

function GeneralComplex{p,T<:Real}(::Type{Val{p}}, x::T, y::T)
    return GeneralComplex{p,T}(x, y)
end

function GeneralComplex{p}(P::Type{Val{p}}, x::Number, y::Number)
    return GeneralComplex(P, promote(x, y)...)
end

typealias GDouble{T} GeneralComplex{1,T}
typealias GDual{T} GeneralComplex{0,T}
typealias GComplex{T} GeneralComplex{-1,T}

########################
# Promotion/Conversion #
########################

function Base.promote_rule{p,A,B}(::Type{GeneralComplex{p,A}}, ::Type{GeneralComplex{p,B}})
    return GeneralComplex{p,promote_type(A, B)}
end

function Base.promote_rule{p,A,B}(::Type{GeneralComplex{p,A}}, ::Type{B})
    return GeneralComplex{p,promote_type(A, B)}
end

function Base.convert{p,T}(::Type{GeneralComplex{p,T}}, g::GeneralComplex{p})
    return GeneralComplex{p,T}(g.x, g.y)
end

function Base.convert{p,T}(::Type{GeneralComplex{p,T}}, n::Number)
    return GeneralComplex{p,T}(n, zero(T))
end

##############
# Arithmetic #
##############

function Base.(:*){p}(a::GeneralComplex{p}, b::GeneralComplex{p})
    return GeneralComplex(Val{p}, a.x*b.x + p*a.y*b.y, a.y*b.x + a.x*b.y)
end

Base.(:*){p}(g::GeneralComplex{p}, n::Number) = GeneralComplex(Val{p}, g.x*n, g.y*n)
Base.(:*){p}(n::Number, g::GeneralComplex{p}) = GeneralComplex(Val{p}, n*g.x, n*g.y)

function Base.(:/){p}(a::GeneralComplex{p}, b::GeneralComplex{p})
    numx = a.x*b.x - p*a.y*b.y
    numy = a.y*b.x - a.x*b.y
    return GeneralComplex(Val{p}, numx, numy) / (b.x^2 - p*b.y^2)
end

function Base.(:/){p}(n::Number, g::GeneralComplex{p})
    numx = n*g.x
    numy = -n*g.y
    return GeneralComplex(Val{p}, numx, numy) / (g.x^2 - p*g.y^2)
end

Base.(:/){p}(g::GeneralComplex{p}, n::Number) = GeneralComplex(Val{p}, g.x/n, g.y/n)

function Base.(:+){p}(a::GeneralComplex{p}, b::GeneralComplex{p})
    return GeneralComplex(Val{p}, a.x+b.x, a.y+b.y)
end

function Base.(:-){p}(a::GeneralComplex{p}, b::GeneralComplex{p})
    return GeneralComplex(Val{p}, a.x-b.x, a.y-b.y)
end

Base.(:-){p}(g::GeneralComplex{p}) = GeneralComplex(Val{p}, -g.x, -g.y)

################
# Trigonometry #
################

tanp{p}(theta, P::Type{Val{p}}) = sinp(theta, P) / cosp(theta, P)

@generated function cosp{p}(theta, ::Type{Val{p}})
    if p == 0
        return 1
    else
        sqrt_p = sqrt(abs(p))
        if p < 0
            return :(cos(theta * $sqrt_p))
        else
            return :(cosh(theta * $sqrt_p))
        end
    end
end

@generated function sinp{p}(theta, ::Type{Val{p}})
    if p == 0
        return :(theta)
    else
        sqrt_p = sqrt(abs(p))
        inv_sqrt_p = inv(sqrt_p)
        if p < 0
            return :($inv_sqrt_p * sin(theta * $sqrt_p))
        else
            return :($inv_sqrt_p * sinh(theta * $sqrt_p))
        end
    end
end

@generated function Base.angle{p}(g::GeneralComplex{p})
    sigma = :(g.y/g.x)

    if p == 0
        return sigma
    else
        sqrt_p = sqrt(abs(p))
        inv_sqrt_p = inv(sqrt_p)
        if p < 0
            return :($inv_sqrt_p * atan($sigma * $sqrt_p))
        else
            return :($inv_sqrt_p * atanh($sigma * $sqrt_p))
        end
    end
end

Base.conj{p}(g::GeneralComplex{p}) = GeneralComplex(Val{p}, g.x, -g.y)

Base.abs{p}(g::GeneralComplex{p}) = sqrt(abs(g.x^2 - p*g.y^2))

function Base.exp{p}(g::GeneralComplex{p})
    P = Val{p}
    return exp(g.x) * GeneralComplex(P, cosp(g.y, P), sinp(g.y, P))
end

###############
# Test values #
###############
using DualNumbers

x, y = rand(), rand()

dg = GDual{Float64}(x, y)
d = dual(x, y)

cg = GComplex{Float64}(x, y)
c = x + y*im