I’m actually not using a symbolic package, I wanted total control over evaluation and simplification so I rolled my own. Its still in the early stages but I’ll share some snippets here to show you whats going on but I think I’ll put the whole thing up on a github repo soon.

The idea is that I’ve added methods to all the basic mathematical functions such that

```
julia> :x + :y
:(x + y)
julia> :x - 5
:(x + -5) #I'm not sure I'm super happy with this choice but it makes simplification much easier.
julia> sin(:x)
:(sin(x))
```

etc.

And then I made a simplification algorithm which looks for common mathematical patterns in the expression tree and applies replacement rules to them and then keeps on looping until the expression stops changing ie.

```
julia> :(x^1) |> expansion_loop
x
julia> :(0 * x) |> expansion_loop
0
julia> :(1 * x) |> expansion_loop
x
julia> :(0 + x) |> expansion_loop
x
```

etc.

Finally, to take derivatives of symbolic functions I made a struct

```
Mathy = Union{Number, Symbol, Expr}
struct Dual_Number
real::Mathy
infinitesimal::Mathy
end
function dual_number(a::Mathy, b::Mathy)
if b == 0
a
else
Dual_Number(a, b)
end
end
dual_number(a) = dual_number(0, a)
```

and then screwed around with how dual numbers are displayed, ie.

```
function Base.show(io::IO, dual::Dual_Number)
if (dual.infinitesimal isa Number) && (dual.infinitesimal < 0)
op = "-"
else
op = "+"
end
real_string = dual.real == 0 ? "" : "$(dual.real) $op "
infinitesimal_string = dual.infinitesimal == 0 ? "" :
abs(dual.infinitesimal) == 1 ? "ϵ" :
(dual.infinitesimal isa Expr) && (dual.infinitesimal.args[1] == :+ || dual.infinitesimal.args[1] == :-) ? "($(dual.infinitesimal))ϵ" :
"$(dual.infinitesimal)ϵ"
print(io, "$real_string$infinitesimal_string")
end
```

Next, I’ve made methods for some standard mathematical functions for working with dual numbers that are agnostic about numbers versus symbols and expressions (because I can compose all of them). The idea here is that if I define

```
conj(a::Dual_Number) = dual_number(a.real, -(a.infinitesimal))
+(a::Dual_Number, b::Dual_Number) = dual_number(a.real + b.real, a.infinitesimal + b.infinitesimal)
+(a::Dual_Number, b::Mathy) = dual_number(a.real + b, a.infinitesimal)
+(a::Mathy, b::Dual_Number) = dual_number(a + b.real, b.infinitesimal)
-(a::Dual_Number, b::Dual_Number) = dual_number(a.real - b.real, a.infinitesimal - b.infinitesimal)
-(a::Dual_Number, b::Mathy) = dual_number(a.real - b, a.infinitesimal)
-(a::Mathy, b::Dual_Number) = dual_number(a - b.real, b.infinitesimal)
*(a::Dual_Number, b::Dual_Number) = dual_number(a.real * b.real, a.infinitesimal * b.real + a.real * b.infinitesimal)
*(a::Dual_Number, b::Mathy) = dual_number(a.real * b, a.infinitesimal * b)
*(a::Mathy, b::Dual_Number) = dual_number(a * b.real, a * b.infinitesimal)
/(a::Dual_Number, b::Dual_Number) = b.real != 0 ? (a * conj(b))/(b.real)^2 : Inf
/(a::Dual_Number, b::Mathy) = dual_number(a.real / b, a.infinitesimal / b)
/(a::Mathy, b::Dual_Number) = b.real != 0 ? (a * conj(b))/(b.real)^2 : Inf
^(a::Dual_Number, b::Mathy) = dual_number(a.real^b, b * a.real^(b-1) * a.infinitesimal)
^(a::Dual_Number, b::Integer) = dual_number(a.real^b, b * a.real^(b-1) * a.infinitesimal)
^(a::Mathy, b::Dual_Number) = dual_number(b^a.real, log(b) * a^b.real * a.infinitesimal)
log(a::Dual_Number) = dual_number(log(a.real), 1/a.real * a.infinitesimal)
sin(a::Dual_Number) = dual_number(sin(a.real), cos(a.real)*a.infinitesimal)
cos(a::Dual_Number) = dual_number(cos(a.real), -sin(a.real)*a.infinitesimal)
```

The idea here is that if I define `ϵ = dual_number(0, 1)`

then `ϵ^2 == 0`

so `ϵ`

is like a very small (infinitesimal) number and for any function `f(x)`

, `f(x + ϵ) = f(x) + ϵ f'(x)`

I then made a derivative function which takes a function `f`

and returns a function which when evaluated will evaluate `f(t + ϵ)`

and extracts out the infinitesimal part which will be the automatic derivative of `f`

and then simplifies that result.

```
ϵ = dual_number(0, 1)
function D(f::Function)
t::Mathy -> f(t + ϵ).infinitesimal |> expansion_loop
end
```

With this I can do stuff like

```
julia> f(x) = log(sin((x^2)));
julia> D(f)(:x)
:((1 / sin(x ^ 2)) * (cos(x ^ 2) * (2x)))
```

which is correct(!!!) but could be simplified a bit more.

Sorry that turned into a lot of code but I was excited enough that I wanted to share what I had done. Once this is a bit more polished I’ll maybe make some blog posts about it if people are interested.