Wikipedia has a characterization of the Cantor normal form of the sum and product of ordinals given by their normal forms, which yields an algorithm for calculation of sums and products on such expressions. But what about exponentiation?
It is easy to see that if $\alpha,\beta<\varepsilon_0$ then $\alpha^\beta<\varepsilon_0$, so there is a function on ordinal notations that matches the exponential, but of course this simple argument doesn't show it's computable. Nevertheless I highly doubt it is much more complicated than the algorithms for sum and product.
Obviously to answer this we need to have a good handle on multiplying a Cantor normal form by itself. Let $\alpha=\alpha'+m=\omega^{\alpha_0}a_0+\dots+\omega^{\alpha_n}a_n+m>1$ where $m<\omega$ is possibly zero and $\alpha_n>0$.
If $m=0$, then $\alpha^2=\omega^{\alpha_0+\alpha_0}a_0+\dots+\omega^{\alpha_0+\alpha_n}a_n=\omega^{\alpha_0}\alpha$. It follows that $\alpha^{k+1}=\omega^{\alpha_0k}\alpha$, and $\alpha^\omega=\omega^{\alpha_0\omega}$ (after reducing $\alpha_0\omega$) because the leading exponent $\alpha_0k$ of $\alpha^k$ is growing, so each term is bounded between $\omega^{\alpha_0k}$ and $\omega^{\alpha_0k+1}$ and hence the limit is $\large\omega^{\sup_{k<\omega}\alpha_0k}=\omega^{\alpha_0\omega}$.
Then $\alpha^{\omega+k+1}=\omega^{\alpha_0(\omega+k)}\alpha$, so $\alpha^{\omega2}=\omega^{\alpha_0\omega2}$, and at this point the pattern is clear: at limit ordinals $\beta$, $\alpha^\beta=\omega^{\alpha_0\beta}$, and at successors, $\alpha^{\beta+k+1}=\omega^{\alpha_0(\beta+k)}\alpha$.
If $m\ne0$, then let $\alpha=\alpha'+m$ where $\alpha'$ is the remainder of the Cantor normal form. Then $\alpha^2=\omega^{\alpha_0}\alpha'+\alpha' m+m$, where $\alpha' m=\omega^{\alpha_0}a_0m+\dots+\omega^{\alpha_n}a_n$ only adds a coefficient to the leading term. Working out another power, $\alpha^3=\omega^{\alpha_02}\alpha'+\omega^{\alpha_0}\alpha'm+\alpha' m+m$, and now the general pattern is apparent: $\alpha^{k+1}=\omega^{\alpha_0k}\alpha'+\sum_{i=k-1}^0\omega^{\alpha_0i}\alpha'm+m$, which is a CNF expansion with $(k+1)(n-1)+1$ terms.
At $\omega$, since the leading exponent is still $\alpha_0k$, for the same reasons as above $\alpha^\omega=\omega^{\alpha_0\omega}$. Then $\alpha^{\omega+k+1}=\omega^{\alpha_0(\omega+k)}\alpha'+\sum_{i=k-1}^0\omega^{\alpha_0(\omega+i)}\alpha'm+\omega^{\alpha_0\omega}m$, so $\alpha^{\omega2}=\omega^{\alpha_0\omega2}$ and so on. In general we have $\alpha^\beta=\omega^{\alpha_0\beta}$ and $$(\alpha'+m)^{\beta+k+1}=\omega^{\alpha_0(\beta+k)}\alpha'+\sum_{i=k-1}^0\omega^{\alpha_0(\beta+i)}\alpha'm+\omega^{\alpha_0\beta}m.$$
Another possibility is that $\alpha_0$ doesn't exist, because $\alpha=m>1$ is finite. In this case, of course $m^k=m^k$, and $m^\omega=\omega$. $m^{\omega+k}=\omega m^k$, and $m^{\omega2}=\omega^2$. This generalizes to $m^{\omega i+k}=\omega^im^k$, so $m^{\omega^2}=\omega^\omega$, and so the pattern is $m^{\omega\beta+k}=\omega^\beta m^k$ (where $\beta$ is obtained by factoring out $\omega$ from the infinite part of the exponent).