module My

struct D{T<:Number} <: Number
    num::T
    dual::T
end

num(x::D) = x.num
dual(x::D) = x.dual
numdual(x::D) = (num(x), dual(x))

Base.eltype(::Type{D{T}}) where T<:Number = T
Base.zero(::Type{D{T}}) where T<:Number = D(zero(T), zero(T))
Base.one(::Type{D{T}}) where T<:Number = D(one(T), zero(T))

D(x::Number, y::Number) = D(promote(x, y)...)
D(x::Number) = D(x, zero(x))
D(x::D) = x

D{T}(x::Number) where {T<:Number} = D{T}(x, 0)
D{T}(x::D) where {T<:Number} = D{T}(num(x), dual(x))

function (::Type{T})(x::D) where {T<:Number}
    iszero(dual(x)) ? T(num(x))::T : throw(InexactError(nameof(T), T, x))
end

Base.promote_rule(::Type{D{T}}, ::Type{S}) where {T<:Number, S<:Number} = D{promote_type(T, S)}
Base.promote_rule(::Type{D{T}}, ::Type{D{S}}) where {T<:Number, S<:Number} = D{promote_type(T, S)}

const ε = D(false, true)

Base.show(io::IO, x::D) = print(io, "D(", num(x), ", ", dual(x), ')')

function Base.show(io::IO, x::D{T}) where T<:Real
    a, b = numdual(x)
    compact = get(io, :compact, false)
    show(io, a)
    if signbit(b) && !isnan(b)
        print(io, compact ? "-" : " - ")
        if isa(b, Signed) && !isa(b, BigInt) && b == typemin(typeof(b))
            show(io, -widen(b))
        else
            show(io, -b)
        end
    else
        print(io, compact ? "+" : " + ")
        show(io, b)
    end
    if !(isa(b, Integer) && !isa(b, Bool) || isa(b, AbstractFloat) && isfinite(b))
        print(io, "*")
    end
    print(io, "ε")
end

for op in (:+, :-)
    @eval Base.$op(x::D) = D($op(num(x)), $op(dual(x)))  
    @eval Base.$op(x::D, y::D) = D($op(num(x), num(y)), $op(dual(x), dual(y)))  
end
Base.:*(x::D, y::D) = D(num(x)*num(y), num(x)*dual(y) + dual(x)*num(y))
Base.:/(x::D, y::D) = D(num(x)/num(y), (-num(x)*dual(y) + dual(x)*num(y))/num(y)^2)
Base.:\(x::D, y::D) = D(num(x)\num(y), num(x)^2\(num(x)*dual(y) - dual(x)*num(y)))

Base.cos(x::D) = (y0 = cos(num(x)); y1 = -sin(num(x)); D(y0, y1*dual(x)))
Base.sin(x::D) = (y0 = sin(num(x)); y1 =  cos(num(x)); D(y0, y1*dual(x)))
Base.exp(x::D) = (y0 = exp(num(x)); y1 = y0          ; D(y0, y1*dual(x)))
Base.log(x::D) = (y0 = log(num(x)); y1 = 1/num(x)    ; D(y0, y1*dual(x)))

function Base.:^(x::D, y::D)
    z0 = num(x)^num(y)
    z1 = z0*(log(num(x))*dual(y) + dual(x)*num(y)/num(x))
    D(z0, z1)
end

Base.abs(x::D{<:Real}) = D(abs(num(x)), sign(num(x))*dual(x))

function Base.isless(x::D{<:Real}, y::D{<:Real})
    isless(num(x), num(y)) || isless(dual(x), dual(y)) 
end
Base.isless(x::D{<:Real}, y::Real) = isless(promote(x, y)...)
Base.isless(x::Real, y::D{<:Real}) = isless(promote(x, y)...)

function Base.abs(x::D{<:Complex})
    a, b = numdual(x)
    z0 = abs(a)
    z1 = real(conj(a)*b)/z0
    D(z0, z1)
end

end

f(x) = 2x^2 + 3x + 10
f′(x) = My.dual(f(My.D(x, 1)))
f′(1)

@code_warntype f′(1)

using BenchmarkTools, Random; Random.seed!(4649373)
ENV["LINES"] = 10
g(x) = evalpoly(x, (1/1, 1/1, 1/2, 1/6, 1/24, 1/120, 1/720, 1/5040, 1/40320, 1/362880, 1/3628800))
n = 10^6
x = My.D.(rand(n), rand(n))
@btime g.($x)

d = 1My.ε
cos(d), sin(d), exp(d), log(d)

(2 + My.ε)^3

((2 + 1e-10)^3 - 8.0) * 1e10

2^(3 + My.ε)

(2^(3 + 1e-10) - 8.0) * 1e10

(2 + My.ε)^(3 + My.ε)

exp((3 + My.ε)*log(2 + My.ε))

((2 + 1e-10)^(3 + 1e-10) - 8.0) * 1e10

A = [
    1 +  My.ε 2 + 2My.ε
    3 + 3My.ε 4 + 4My.ε
]

inv(A)

A*inv(A)

abs(2 + My.ε)

abs(-2 + My.ε)

(2 + My.ε)^2

((2 + My.ε)^2)^(1/2)

1 < 1 + My.ε < 1 + 2My.ε < 2

B = [
    1im +  My.ε 2im + 2My.ε
    3im + 3My.ε 4im + 4My.ε
]

inv(B)

B*inv(B)

abs(1 + im + My.ε)

abs(1 + im + My.ε) |> typeof

(abs(1 + im + 1e-8) - √2)*1e8