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? 
 A: $\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*}
A: 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.$$
A: The other answer is not very consistent with the set of Dual Numbers. You don't need any approximation to prove this.
First suppose that $a\neq 0\\$.
The second step is to prove that for a generic Dual Number $z=a+\epsilon b\,\,\text{with } a\neq 0$, holds :
\begin{align}
z=e^{\log(z)}
\end{align}
and holds since for the Dual Numbers set we have the following:
\begin{align*}
\log(z)= \log(a+\epsilon b)=\log(a)+\frac{b}{a}\epsilon
\end{align*}
\begin{align*}
e^{z}= e^{a+\epsilon b}=e^a+e^ab\epsilon
\end{align*}
Combining them we have
\begin{align*}
e^{\log(z)}= e^{\log(a)+\epsilon \frac{b}{a} }=e^{\log(a)}+e^{\log(a)}\frac{b}{a}\epsilon=a+\epsilon b=z
\end{align*}
Coming back to the original problem we have:
\begin{align*}
(a+b\epsilon)^{c+d\epsilon} &= e^{ (c+d\epsilon)\log(a+b\epsilon) } 
  & z=e^{\log z} \\
&=e^{ (c+d\epsilon)(\log a+b\epsilon/a) } 
  & \log(a+\epsilon b)=\log(a)+\frac{b}{a}\epsilon\\
&=e^{ c\log a+\epsilon(d\log a+b c/a) }& \epsilon^2=0 \\
&=e^{c\log a}\,e^{ \epsilon\,d\log a }\,e^{ \epsilon \frac{b\,c}{a} }\\
&=a^{c+\epsilon d}\,e^{\epsilon \frac{b\,c}{a}}
\end{align*}
or, if you prefer
\begin{align*}
(a+b\epsilon)^{c+d\epsilon} &=a^{c+\epsilon d}\,e^{\epsilon \frac{b\,c}{a}}=a^{c+\epsilon d}\,(1+\epsilon \frac{b\,c}{a})
\end{align*}
A: In Mathematica use this code:
$Pre = MatrixFunction[
      Function[ε, #], {{0, 0}, {1, 0}}] /. {{a_, 
        b_}, {0, a_}} -> a + ε b /. {{a_, 0}, {b_, a_}} -> a + ε b &;

(a + b ε)^(c + d ε)

The result:
a^c + a^(-1 + c) (b c + a d Log[a]) ε

