Dual number $(a+b\varepsilon)$ raised to a dual power, e.g. $(a+b\varepsilon)^{(c+d\varepsilon)}$

I'm working on some code which utilizes Newton's method, and I would like to take advantage of dual numbers to simplify taking the derivative. I've worked out a class definition Dual which works great for polynomials, pow (where either base or exponent is real) exp, and log.

However, I am a bit stymied at raising a dual number to the power of another dual number, e.g.

$(a+b\varepsilon)^{(c+d\varepsilon)}$

I used the general form derived from the Taylor series $f(a+b\varepsilon) = f(a) + bf'(a)\varepsilon$ to derive the rules for $x^n$ and $n^x$ where $n$ is real and $x$ is dual. For those, I got:

$x^n = a^n + bna^{n-1}$ , where $ x = a+b\varepsilon$

$n^y = n^c + d\ln(n)n^{c}$ , where $ y = c+d\varepsilon$

Substituting and simplifying, I end up with:

$x^y = a^c + (d\ln(a)a^c + bca^{c-1})\varepsilon$

This seems to work right in my code, but I do not know if this is actually correct. The implementations of exp and pow work properly with this more general code, and $\varepsilon^\varepsilon = 1 + NaN\varepsilon$, which is a good-ish sign (if I got something well-formed, that would be more troubling).

Is this formula valid?

$\def\e{\varepsilon}$\begin{align*} (a+b\e)^{c+d\e} &= \exp((c+d\e)\ln(a+b\e)) & x=e^{\ln x} \\ &=\exp((c+d\e)(\ln a+\ln(1+b\e/a))) \\ &=\exp((c+d\e)(\ln a+b\e/a)) & \ln(1+x)=x+\mathrm{h.o.} \\ &=\exp(c\ln a+\e(d\ln a+bc/a))\\ &=a^c e^{\e(d\ln a+bc/a)} \\ &=a^c (1+\e(d\ln a+bc/a)). & e^x=1+x+\mathrm{h.o.} \end{align*}

The accepted answer is the correct one. The derivation is cool but there is an easier way to get to this result. Notice that the result is just the real function + $\epsilon$ times the sum of the results you got for $x^n$ and $n^y$. You already know $$f(a+b\epsilon) = f(a) + f'(a)b\epsilon.$$ In $x^y$ you have a function with two arguments, and this works in a similar way: $$f(a+b\epsilon,c+d\epsilon) = f(a,c) + \frac{\partial}{\partial a}f(a,c)b\epsilon + \frac{\partial}{\partial c}f(a,c)d\epsilon. $$ Applying this to $x^y$ you get your result: $$\left(a+b\epsilon\right)^{c+d\epsilon} = a^c + ca^{c-1}b\epsilon +\ln(a)a^cd\epsilon.$$