See https://github.com/JuliaLang/julia/blob/master/base/complex.jl

In [1]:

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

Out[1]:

Main.My

In [2]:

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

Out[2]:

In [3]:

@code_warntype f′(1)

MethodInstance for f′(::Int64)
  from f′(x) in Main at In[2]:2
Arguments
  #self#::Core.Const(f′)
  x::Int64
Body::Int64
1 ─ %1 = Main.My.dual::Core.Const(Main.My.dual)
│   %2 = Main.My.D::Core.Const(Main.My.D)
│   %3 = (%2)(x, 1)::Core.PartialStruct(Main.My.D{Int64}, Any[Int64, Core.Const(1)])
│   %4 = Main.f(%3)::Main.My.D{Int64}
│   %5 = (%1)(%4)::Int64
└──      return %5

In [4]:

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)

  22.527 ms (2 allocations: 15.26 MiB)

Out[4]:

1000000-element Vector{Main.My.D{Float64}}:
 1.4872473757882296 + 1.2279017958415617ε
  2.1741933564723377 + 1.039486994972546ε
  2.208976491862254 + 0.5511949973102477ε
                                        ⋮
 2.3452847947934723 + 1.3505676136686622ε
 2.6183423982243044 + 1.6468869986017245ε

In [5]:

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

Out[5]:

(1.0 - 0.0ε, 0.0 + 1.0ε, 1.0 + 1.0ε, -Inf + Inf*ε)

In [6]:

(2 + My.ε)^3

Out[6]:

8 + 12ε

In [7]:

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

Out[7]:

12.000000992884452

In [8]:

2^(3 + My.ε)

Out[8]:

8.0 + 5.545177444479562ε

In [9]:

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

Out[9]:

5.545182091282186

In [10]:

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

Out[10]:

8.0 + 17.545177444479563ε

In [11]:

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

Out[11]:

7.999999999999998 + 17.54517744447956ε

In [12]:

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

Out[12]:

17.545183084166638

In [13]:

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

Out[13]:

2×2 Matrix{Main.My.D{Int64}}:
 1+1ε  2+2ε
 3+3ε  4+4ε

In [14]:

inv(A)

Out[14]:

2×2 Matrix{Main.My.D{Float64}}:
 -2.0+2.0ε   1.0-1.0ε
  1.5-1.5ε  -0.5+0.5ε

In [15]:

A*inv(A)

Out[15]:

2×2 Matrix{Main.My.D{Float64}}:
 1.0-2.22045e-16ε  1.11022e-16+0.0ε
 0.0-8.88178e-16ε          1.0+0.0ε

In [16]:

abs(2 + My.ε)

Out[16]:

2 + 1ε

In [17]:

abs(-2 + My.ε)

Out[17]:

2 - 1ε

In [18]:

(2 + My.ε)^2

Out[18]:

4 + 4ε

In [19]:

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

Out[19]:

2.0 + 1.0ε

In [20]:

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

Out[20]:

true

In [21]:

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

Out[21]:

2×2 Matrix{Main.My.D{Complex{Int64}}}:
 D(0+1im, 1+0im)  D(0+2im, 2+0im)
 D(0+3im, 3+0im)  D(0+4im, 4+0im)

In [22]:

inv(B)

Out[22]:

2×2 Matrix{Main.My.D{ComplexF64}}:
 D(-0.0+2.0im, -2.0-0.0im)    D(0.0-1.0im, 1.0+0.0im)
   D(0.0-1.5im, 1.5+0.0im)  D(-0.0+0.5im, -0.5-0.0im)

In [23]:

B*inv(B)

Out[23]:

2×2 Matrix{Main.My.D{ComplexF64}}:
         D(1.0+0.0im, 0.0+0.0im)  D(0.0+0.0im, 0.0-1.11022e-16im)
 D(8.88178e-16+0.0im, 0.0+0.0im)  D(1.0+0.0im, 0.0-2.22045e-16im)

In [24]:

abs(1 + im + My.ε)

Out[24]:

1.4142135623730951 + 0.7071067811865475ε

In [25]:

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

Out[25]:

Main.My.D{Float64}

In [26]:

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

Out[26]:

0.7071067731345693

In [ ]: