How to calculate $\arcsin$ of dual number: $\arcsin\left( x + y \cdot \varepsilon \right) = \dots$ where $\begin{align*} a, ~b &\in \mathbb{R}\\ \varepsilon^{2} &= 0\\ \varepsilon &\ne 0\\ \end{align*}$?
I recently saw the question on SE where someone asked how to get $\sin$ from a dual number, so I had the idea of reversing it, i.e. getting the $\arcsin$ from a dual number , to investigate. For this I wanted to derive it myself, where several ideas would come to me on how to do it (via the logarithm relationship, the series expansion and calculating backwards). With the individual methods, however, I always came up with other solutions. So I'm wondering what's correct, what I'm doing wrong and if there is a better method for this?
I gave up on the logarithm method, since I would have to use the product of two imaginary units there, which, however, are not commutative with regard to multiplication, which would severely limit the validity of the formulas.
My calculating backwards attempt
I'm using the formula from this SE answer:
$$ \begin{align*} \sin\left( u + v \cdot \varepsilon \right) &= \sin\left( u \right) + v \cdot \cos\left( u \right) \cdot \varepsilon\\ \sin\left( u + v \cdot \varepsilon \right) &= x + y \cdot \varepsilon\\ \\ x = \sin\left( u \right) &\wedge y = v \cdot \cos\left( u \right)\\ \arcsin\left( x \right) = u &\wedge y = v \cdot \cos\left( \arcsin\left( x \right) \right)\\ \arcsin\left( x \right) = u &\wedge y = v \cdot \cos\left( \arcsin\left( x \right) \right)\\ \arcsin\left( x \right) = u &\wedge y = v \cdot \sqrt{1 - x^{2}}\\ \arcsin\left( x \right) = u &\wedge \frac{y}{\sqrt{1 - x^{2}}} = v\\ \\ \sin\left( \arcsin\left( x \right) + \frac{y}{\sqrt{1 - x^{2}}} \cdot \varepsilon \right) &= x + y \cdot \varepsilon\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \arcsin\left( x \right) + \frac{y}{\sqrt{1 - x^{2}}} \cdot \varepsilon\\ \end{align*} $$
And this can only be true for $\left| x \right| < 1$ which I find strange since I didn't expact that $\left| x \right| = 1$ gets lost as a possible argument.
My series expansion attempt
$$ \begin{align*} \arcsin\left( z \right) &= z + \left( \frac{1}{2} \right) \cdot \frac{z^{3}}{3} + \left( \frac{1 \cdot 3}{2 \cdot 4} \right) \cdot \frac{z^{5}}{5} + \left( \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \right) \cdot \frac{z^{7}}{7} + \cdots\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \left( x + y \cdot \varepsilon \right) + \left( \frac{1}{2} \right) \cdot \frac{\left( x + y \cdot \varepsilon \right)^{3}}{3} + \left( \frac{1 \cdot 3}{2 \cdot 4} \right) \cdot \frac{\left( x + y \cdot \varepsilon \right)^{5}}{5} + \left( \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \right) \cdot \frac{\left( x + y \cdot \varepsilon \right)^{7}}{7} + \cdots\\ \arcsin\left( x + y \cdot \varepsilon \right) &= x + y \cdot \varepsilon + \left( \frac{1}{2} \right) \cdot \frac{x^{3} + 3 \cdot x^{2} \cdot y \cdot \varepsilon}{3} + \left( \frac{1 \cdot 3}{2 \cdot 4} \right) \cdot \frac{x^{5} + 5 \cdot x^{4} \cdot y \cdot \varepsilon}{5} + \left( \frac{1 \cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \right) \cdot \frac{x^{7} + 7 \cdot x^{6} \cdot y \cdot \varepsilon}{7} + \cdots\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \arcsin\left( x \right) + \left( 1 + x^{2} + x^{4} + x^{6} + \cdots \right) \cdot y \cdot \varepsilon\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \arcsin\left( x \right) + \left( \sum_{k = 0}^{\infty} x^{2 \cdot k} \right) \cdot y \cdot \varepsilon\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \arcsin\left( x \right) + \frac{1}{1 - x^{2}} \cdot y \cdot \varepsilon, \text{ for } \left| x \right| < 1\\ \arcsin\left( x + y \cdot \varepsilon \right) &= \arcsin\left( x \right) + \frac{y}{1 - x^{2}} \cdot \varepsilon, \text{ for } \left| x \right| < 1\\ \end{align*} $$
But obviously something is wrong: $$ \arcsin\left( x \right) + \frac{y}{1 - x^{2}} \cdot \varepsilon \ne \arcsin\left( x \right) + \frac{y}{\sqrt{1 - x^{2}}} \cdot \varepsilon$$
"Automatic differentiation" is the way to go. Namely, for any real-analytic function $f\colon \Bbb R \to \Bbb R$, its natural extension to $\Bbb R\oplus \Bbb R \varepsilon$ in fact satisfies $$f(x+y\varepsilon) = f(x) + yf'(x)\varepsilon,$$for all $x,y\in \Bbb R$. To see this, write $$f(x) = f(0) + f'(0)x + \frac{f''(0)}{2}x^2 + \frac{f'''(0)}{6} x^3+\cdots$$and replace $x\mapsto x+y\varepsilon$, noting that $$(x+y\varepsilon)^n = x^n + nx^{n-1}y\varepsilon$$for all $x,y\in \Bbb R$, to obtain $$\begin{align} f(x+y\varepsilon) &= f(0) + f'(0)(x+y\varepsilon) + \frac{f''(0)}{2}(x^2+2xy\varepsilon) + \frac{f'''(0)}{6}(x^3+3x^2y\varepsilon) + \cdots \\ &= f(0) + f'(0) x + \frac{f''(0)}{2}x^2 + \frac{f'''(0)}{6}x^3+\cdots \\ &\qquad + y\varepsilon\left(f'(0) + f''(0)x + \frac{f'''(0)}{2}x^2+\cdots\right) \\ &= f(x) + yf'(x)\varepsilon.\end{align}$$ Taking $f = \arcsin$, we have that $$\arcsin(x+y\varepsilon) = \arcsin(x) + \frac{y\varepsilon}{\sqrt{1-x^2}}.$$Trying to do this calculation with $\arcsin$ instead of an abstract $f$ turns out to be a distraction.