Eine Lösungsmöglichkeit stellt die Lösung des Gleichungssystems \[\begin{align} ||P5 - M|| - r = 0 \\ ||g(s)-M|| - r = 0 \\ ||h(t) - M|| - r = 0 \\ \langle M-g(s), v\rangle = 0 \\ \langle M-h(t)), w\rangle = 0 \end{align}\] dar, 5 unabhängige Gleichungen für 5 Variablen liefern eine eindeutige Lösung.
Die Bestimmung ist nur numerisch möglich, wie du das letztendlich implementierst hängt wohl stark von der verwendeten Programmiersprache ab.
Für mathematische Berechnungen ausgelegte Sprachen wie Matlab, R, Julia, usw. sollten aber einfach Funktionen zum Lösen solcher nichtlinearen Gleichungssysteme bieten.

Hier eine mögliche Implementierung in Julia (ist recht häßlich, da ich das ursprünglich bloß zur Überprüfung verwendet habe). Ich hab mal mit ein paar Werten rumgespielt, es scheint soweit sinnvolle Ergebnisse zu liefern.

Wenn du es Bare-Metal implementieren willst, könntest du das Newton-Verfahren verwenden: \[x_{k+1} = x_k - J_f^{-1}(x_k) f(x_k)\] wobei \(f\) hier die Funktionen auf der jeweils linken Seite des obigen Gleichungssystems sind, und \(J_f\) ihre Jakobimatrix ist, \(x_k\) kann man in beliebiger Reihenfolge definieren, z.B. mit \(x_1\) = x-Koord und \(x_2\) = y-Koordinate des Mittelpunkts, \(x_3\) = Radius des Kreises und \(x_{4,5}\) = s und t.

Mit etwas Umstellung muss man nun für jede Iteration \(x_k \rightsquigarrow x_{k+1}\) das lineare Gleichungssystem \(J_f(x_k) \Delta x_k = f(x_k)\) für\(\Delta x_k\) lösen (z.B. Gauß-Algorithmus - man kann natürlich auch einfach die Inverse berechnen), wobei dann \(x_{k+1} = x_k - \Delta x_k\).
Man benötigt noch einen Startwert \(x_0\), bei dessen Wahl man etwas aufpassen muss, normalerweise sollte es jedoch kein Problem darstellen. Üblicherweise beendet man das Verfahren nach Unterschreiten einer gewissen Toleranz (= Abbruchbedingung) wenn \(||x_{k+1} - x_k|| < TOL\)

Für die Jacobimatrix erhält man (Erste 3 Gleichungen quadriert) sowas wie  \[J_f(x) = \begin{pmatrix} -2(E-M)^T & -2r & 0&0 \\ -2(g(s)-M)^T & -2r & 2\langle g(s) - M, v\rangle & 0\\ -2(h(t) - M)^T & -2r & 0 & 2\langle h(t)-M,w\rangle\\ v^T & 0 & -||v||^2 & 0\\w^T & 0 & 0 & -||w||^2\end{pmatrix}\]