I have been trying to extend faithfully extend the dual numbers, recreationally, of course.
Background: As a quick summary, our domain of discourse is $\mathbb{R}_\varepsilon=\{a+b\varepsilon:(a,b)\in\mathbb{R}^2\}$, where $\varepsilon\neq 0$, yet $\varepsilon^2=0$. In this way, $\varepsilon$ is effectively an infinitesimal. This makes it a incredibly useful alternative in many analytical situations to limits.
In short, we may formalize its behavior with the isomorphism $\phi(a+b\varepsilon)=\begin{bmatrix} a & b \\ 0 & a \end{bmatrix}$.
Problem / observation: By default, one cannot divide by $\varepsilon$, in the same way that $\begin{bmatrix} 0 & 1 \\ 0 & 0\end{bmatrix}$ has no nice inverse (as far as I'm aware; my linear algebra knowledge is limited however). However, there are certain situations where I have found that definition $\omega=\frac{1}{\varepsilon}$ is consistent. In a similar way as to $\varepsilon$ effectively being an infinitesimal, $\omega$ plays the role of an infinity.
For example, it is not so difficult to show that $$\ln(x)=\omega(x^\varepsilon-1),$$ which can be rearranged to find $$x^\varepsilon=1+\varepsilon\ln(x),$$ which is useful, for example, in efficiently evaluating the following limit. $$\begin{align*} \lim_{x\to 0}\frac{2^x-3^x}{x} &= \omega(2^\varepsilon-3^\varepsilon)\\ &= \omega(1+\varepsilon\ln(2)-1-\varepsilon\ln(3))\\ &= \ln\left(\frac{2}{3}\right). \end{align*}$$
I assume this effectiveness is no coincidence, thus I presume there is some consistent way to define $\omega=\frac{1}{\varepsilon}$. In particular, how do we define $\omega^2$?
My attempt: My initial idea, was to treat $\omega$ as some sort of a hyperinteger. However, this idea seems to conflict with the fact that since $\varepsilon^2=0$, it seems only reasonable that $\omega^2=\infty$. Thus, I thought that a better domain of discourse would be $\{a\omega+b+c\varepsilon:(a,b,c)\in\mathbb{R}^2\}\cup\{-\infty,\infty\}$, where we define the behavior of $\infty$ similar to extended real number line or the Real projective line. This is where it gets way over my head, trying to define how $\omega$ relates to / interacts with / is distinct from $\infty$, and whether or not to introduce a further symbol$\bot$ from Wheel theory, and how might the isomorphism with $2\times 2$ matrices may help or not, etc.
So I'm hoping someone can straighten out this idea into something much more simple than all the conflicting ideas I have. I just want to have some extension to the dual numbers, in which division is defined for infinitesimals. Hence, any advice or ideas would be greatly appreciated.
Your calculations
For $a>0$, the Maclaurin series (Taylor series centered at $0$) for $a^{x}$ is $1+x\ln a+x^{2}\dfrac{\left(\ln a\right)^{2}}{2!}+\cdots$. So if we replace $x$ with $\varepsilon$, it makes sense to define $a^{\varepsilon}=1+\varepsilon\ln a$ (defining functions using their Maclaurin series is common when working with dual numbers). Renaming $a$ to “$x$”, we recover (for positive $x$) “$x^{\varepsilon}=1+\varepsilon\ln x$” as in the question.
Now, since every number $\gamma$ in the dual numbers has the form $a+b\varepsilon$, which one might write in terms of the “real and dual parts” of $\gamma$ as $\mathrm{Re}\left(\gamma\right)+\mathrm{Du}\left(\gamma\right)\varepsilon$. Then certainly $\gamma-\mathrm{Re}\left(\gamma\right)=\mathrm{Du}\left(\gamma\right)\varepsilon$. And so we might want to write a formal inverse $\omega\cdot$ for $\varepsilon$ which either is undefined for dual numbers that have a nonzero real part (this is probably best), or just returns the dual part (ala $\begin{bmatrix}a & b\\ 0 & a \end{bmatrix}\mapsto b$), and then we have $\omega\cdot\left(\gamma-\mathrm{Re}\left(\gamma\right)\right)=\mathrm{Du}\left(\gamma\right)$ and $\ln\left(x\right)=\omega\cdot\left(x^{\varepsilon}-1\right)$, etc.
Now, when faced with a limit like $\lim_{x\to0}\dfrac{2^{x}-3^{x}}{x}$: Since the numerator has limit $0$, we know that replacing $x$ with $\varepsilon$ in the numerator will have the form $\mathrm{Du}\left(\gamma\right)\varepsilon$, which will be amenable to application of $\omega$, which is exactly suggested by the $x$ in the denominator. This will work out fine because L'Hospital's rule says we can take the derivative of numerator and denominator, so that we basically want to take the ratio of the dual parts, which is indeed $\dfrac{\ln2-\ln3}{1}=\ln\left(\dfrac{2}{3}\right)\checkmark$.
What about $\omega^2$
However, if you want to consider something like $\omega^{2}$, it would seem that would be like applying this “removal of $\varepsilon$” twice, which either doesn't make sense or always gives $0$, since $a\varepsilon^{2}=0$ for real $a$. So I think $\omega^2$ at face value is a dead-end.
Extending things
Let's try this $\omega$ method on another limit from Evaluating the limit using Taylor Series: $$\lim_{x\to0}\frac{e^{3x}-\sin(x)-\cos(x)+\ln(1-2x)}{-1+\cos(5x)}$$
$$\substack{?\\ = } \dfrac{\left(1+3\varepsilon\right)-\varepsilon-1-2\varepsilon}{-1+1}=\dfrac{0}{0}\cdots$$
Here, there is no $\varepsilon$ to divide by, and so no help from $\omega$, and we are stuck.
However, if we give up “the simplicity of $\varepsilon^{2}=0$” but keep $\varepsilon^{3}=0$, we can then solve the limit by some intuitive reasoning/approximations:
$$\lim_{x\to0}\frac{e^{3x}-\sin(x)-\cos(x)+\ln(1-2x)}{-1+\cos(5x)}$$$$\substack{?\\ = } \dfrac{\left(1+3\varepsilon+\dfrac{9\varepsilon^{2}}{2}\right)-\varepsilon-\left(1-\dfrac{\varepsilon^{2}}{2}\right)+\left(-2\varepsilon-2\varepsilon^{2}\right)}{-1+\left(1-\dfrac{25}{2}\varepsilon^{2}\right)}$$$$=\dfrac{\dfrac{9\varepsilon^{2}}{2}+\dfrac{\varepsilon^{2}}{2}-2\varepsilon^{2}}{-\dfrac{25}{2}\varepsilon^{2}}\substack{?\\ = } \omega^{2}\cdot\dfrac{\dfrac{9}{2}+\dfrac{1}{2}-2}{-\dfrac{25}{2}}\varepsilon^{2}=-\dfrac{6}{25}\checkmark$$
But in some other more subtle limit problem, maybe we'd have to give up on $\varepsilon^{3}=0$ and only be able to keep $\varepsilon^{17}=0$, if the problem were set up so that it could naturally be done by applying L'Hospital's rule $16$ times.
A fix
Instead of fixing a level at which we have perfect equality with $0$, it would be most generally useful if we said that no power of an infinitesimal like $\varepsilon$ actually equaled zero, but instead that (if we had a limit involving nice functions) we could turn any extra terms with $\varepsilon$s like $-5\varepsilon$ or $+7\varepsilon^{3}$ to zero at the end of the calculation, and get the right answer.
This is pretty close to what what is done in Robinson's non-standard analysis, which you could learn about at an intuitive level from a textbook like “Elementary Calculus: An Infinitesimal Approach” by Keisler.