Let the ellipse be in the standard orientation $(x/a)^2+(y/b)^2=1$. This defines the coordinate system in use.
Let the random point be $(x_r,y_r)$. Tangents to the ellipse are defined by the slope/derivative $2x/a^2+2yy'/b^2=0$,
so $yy'/b^2=-x/a^2$, so $y'=-xb^2/(ya^2)$. Let the tangent point be $(x_t,y_t)$, so the tangent line that runs through
the tangent point is $(x,y)= (x_t,y_t)+\alpha (1,-x_tb^2/(y_ta^2))$, where $\alpha$ is a parameter along the points of the tangent.
The normal to the tangent has a slope which is -1 divided by the slope of the tangent, so it has the direction $(1,y_ta^2/(x_tb^2))$. The vector from $(x_r,y_r)$ to the tangent point is $(x_t-x_r,y_t-y_r)$. It is decomposed to a component parallel and a component
perpendicular to the tangent line:
$$
\left( \begin{array}{c}x_t-x_r\\ y_t-y_r\end{array}\right)
= \beta \left(\begin{array}{c}1 \\ -x_tb^2/(y_t a^2)\end{array}\right)
+ \gamma \left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right)
$$
Dot products with the direction vectors give due to the orthogonality
$$
\left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right)
\left( \begin{array}{c}x_t-x_r\\ y_t-y_r\end{array}\right)
=
\gamma \left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right)
\left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right);
$$
$$
x_t-x_r+ \frac{y_ta^2}{x_tb^2}(y_t-y_r)
=
\gamma (1+(y_ta^2/(x_tb^2)^2);
$$
So the factor for the component along the normal of the tangent is known:
$$
\gamma=\frac{
x_t-x_r+ \frac{y_ta^2}{x_tb^2}(y_t-y_r)}
{
1+[y_ta^2/(x_tb^2)]^2}.
$$
The line starting at $(x_t,y_t)$ after the reflection is obtained by flipping the component
of the vector before the reflection that is orthogonal.
$$
\beta \left(\begin{array}{c}1 \\ -x_tb^2/(y_t a^2)\end{array}\right)
- \gamma \left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right)
=
\left( \begin{array}{c}x_t-x_r\\ y_t-y_r\end{array}\right)
-2 \gamma \left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right).
$$
Let $(x_h,y_h)$ with $(x_h/a)^2+(y_h/b)^2=1$
be the point at the next collision with the elliptic billiard wall. so the vector emerging from $(x_t,y_t)$
is
$$
\left(
\begin{array}{c}
x_t \\ y_t
\end{array}
\right)
+\delta
[
\left( \begin{array}{c}x_t-x_r\\ y_t-y_r\end{array}\right)
-2 \gamma \left(\begin{array}{c}1 \\ y_t a^2/(x_tb^2)\end{array}\right)
]
=
\left(
\begin{array}{c}
x_h\\ y_h
\end{array}
\right)
$$
The task is to find the parameter $\delta$ which encodes the distance from $(x_t,y_t)$ to $(x_h,y_h)$
along the (non-normalized) direction.
\begin{eqnarray}
x_t+\delta[x_r-x_r)-2\gamma] &=& x_h\\
y_t+\delta[y_t-y_r)-2\gamma y_ta^2/(x_tb^2)] &=& y_h\\
\end{eqnarray}
Inserting these into $(x_h/a)^2+(y_h/b)^2=1$ gives a quadratic equation for $\delta$, and
the solutions which is $\neq 0$ is to be used to find the next collision point.