Have you considered finding the intersections using an implicit form for the circles, $$\frac{x^2}{r^2} + \frac{y^2}{r^2} + ax + by + c = 0?$$ This representation doesn't have any coefficients that diverge as the circle approaches a straight line. To find intersections, you'll have to solve a quadratic equation whose leading coefficient could be zero or arbitrarily close to it, but the alternative form of the quadratic formula should be able to deal with that robustly.
You'll then have to do some jiggery-pokery to figure out whether the intersection points lie within the arcs. If the arc's bending angle is smaller than $\pi$, a projection onto the line joining the endpoints will suffice.
(Disclaimer: While all of this feels like it should work, I haven't analyzed it in any detail. Also, there could still be a problem when the circle is close to a line and you want the longer arc. But I can't imagine that's a case that would turn up in any practical application.)
Update: For a concrete example, here is the equation for a circular arc passing through the three points $(0,0)$, $(0.5, h)$, and $(1,0)$: $$\kappa^2 x^2 + \kappa^2 y^2 - \kappa^2 x - 2\eta y = 0,$$ where $$\begin{align}\kappa &= \frac{8h}{4h^2 + 1}, \\ \eta &= \frac{8h(4h^2-1)}{(4h^2+1)^2}.\end{align}$$ As you can see, the coefficients remain bounded as $h \to 0$.
Update 2: Wait, that equation becomes trivial if $h = 0$, which is bad. We really want something like $x^2/r + y^2/r + ax + by + c,$ i.e. multiply the previous expression through by $r$. Then for the same example, our equation becomes $$\kappa x^2 + \kappa y^2 - \kappa x - 2\eta' y = 0,$$ where $\eta' = (4h^2-1)/(4h^2+1)$. Here are some explicit values.
$h = 1/2$: $$2 x^2 + 2 y^2 - 2 x = 0,$$ $h = 0.01$: $$0.07997 x^2 + 0.07997 y^2 - 0.07997 x + 1.998 y = 0,$$ $h = 0$: $$2 y = 0.$$
By the way, in this format, the linear terms will always be simply $-2(x_0/r)x$ and $-2(y_0/r)y$, where the center of the circle is at $(x_0,y_0)$. As the center goes to infinity but the endpoints remain fixed, these coefficients remain bounded and nonzero (i.e. not both zero).
I will give a solution whose starting point is as follows: I know two points $P_1$ and $P_2$ on line 1 and points $Q_1$ and $Q_2$ on line 2.
Form the vectors ${\bf p} = P_2-P_1$ and ${\bf q}=Q_2-Q_1$.
Check for an intersection of the two lines: ${\bf p}t+P_1={\bf q}s+Q_1$ for some pair of real numbers $s$ and $t$. If no solution is found, the lines do not intersect.
Next, check to see if ${\bf p}$ is a scalar multiple of ${\bf q}$ (i.e. ${\bf p}=s{\bf q}$ for some $s$). If this is true, the directions of the lines are parallel.
There are 4 cases:
1) Intersecting + Parallel Directions = They are the same line (FAIL: One line will not determine a plane.)
2) Not Intersecting + Non-parallel Directions = They are skew lines (FAIL: Skew lines do not lie in a common plane.)
3) Not Intersecting + Parallel Directions = Parallel lines. In this case any normal line for one line will be normal for the other. So any point in the plane will lie on a common normal line. If you would like a point "between" the two parallel lines: Form a vector ${\bf r}=Q_1-P_1$ (this points from $P_1$ on line 1 to $Q_1$ on line 2). Project this vector ${\bf r}$ onto ${\bf p}$ (the direction of the first line): ${\bf w}=\mathrm{proj}_{\bf p}({\bf r}) = \frac{{\bf p\cdot r}}{|{\bf p}|^2}{\bf p}$. Then ${\bf r}-{\bf w}$ is a vector which lies in the plane and is orthogonal to the first line and points from the first line to the second. So $P_1+(1/2)({\bf r}-{\bf w})$ is a point exactly half-way between the two lines.
4) Intersecting + Non-parallel Directions = Distinct intersecting lines. In this case compute the cross product ${\bf n}={\bf p} \times {\bf q}$ to get a vector normal to the plane. Then ${\bf v} = {\bf n} \times {\bf p}$ will lie on the plane and be perpendicular to line 1. Likewise ${\bf w}= {\bf n} \times {\bf q}$ for line 2. Then the lines parametrized by ${\bf v}t+P_1$ and ${\bf w}s+Q_1$ will be normal to lines 1 and 2 respectively. Solving ${\bf v}t+P_1={\bf w}s+Q_1$ will give you the point of intersection of these normals.
Edit: Oops! I forgot to mention. If the lines intersect and you create normal lines at that point of intersection, then those normal lines intersect at that point of intersection!
Best Answer
I suggest turning this into a 2D problem and then find the circle from three points on the plane.
Start by finding the normal vector to the plane defined by the three points
$$ \boldsymbol{n} = \mathrm{unitvector}( \boldsymbol{A} \times \boldsymbol{B} + \boldsymbol{B} \times \boldsymbol{C} + \boldsymbol{C} \times \boldsymbol{A}) \tag{1}$$
The find the arbitrary mutually orthogonal directions
$$ \boldsymbol{u} = (\boldsymbol{C}-\boldsymbol{A})\times \boldsymbol{n} \tag{2}$$ $$ \boldsymbol{v} = \boldsymbol{n} \times \boldsymbol{u} \tag{3}$$
and construct a 3×2 rotation matrix from the three direction vectors as columns.
$$ \mathbf{R} = \left[ \boldsymbol{u} \; \boldsymbol{v} \right] \tag{4}$$
You also need the distance of the plane ABC to the origin
$$ d = \boldsymbol{n} \cdot \boldsymbol{A} \tag{5}$$
Now convert the problem into a 2D problem with
$$ \begin{aligned} \boldsymbol{a} = \pmatrix{a_x\\a_y} & = \mathbf{R}^\top \boldsymbol{A}\\ \boldsymbol{b} = \pmatrix{b_x\\b_y} & = \mathbf{R}^\top \boldsymbol{B}\\ \boldsymbol{c} = \pmatrix{b_x\\b_y} & = \mathbf{R}^\top \boldsymbol{B} \end{aligned} \tag{6}$$
Now solve the 2D problem by finding the center point $\boldsymbol{q}=\pmatrix{x\\y}$ using a 2×2 system of equations
$$ \begin{bmatrix} b_x-a_x & b_y -a_y \\ c_x -a_x & c_y - a_y \end{bmatrix} \pmatrix{x\\y} = \pmatrix{ \frac{ (b_x^2+b_y^2)-(a_x^2+a_y^2)}{2} \\ \frac{ (c_x^2+c_y^2)-(a_x^2+a_y^2)}{2} } \tag{7} $$
Lastly construct the 3D point for the circle center $\boldsymbol{Q}$
$$ \boldsymbol{Q} = \boldsymbol{n}\, d + \mathbf{R} \pmatrix{x\\y} \tag{8}$$
If the radius is needed just take the distance from the center to any point
$$ \mathrm{radius} = \| \boldsymbol{A}-\boldsymbol{Q} \| \tag{9} $$
You might also need the sweep angle which you get from the angle between the vectors $\boldsymbol{C}-\boldsymbol{Q}$ and $\boldsymbol{A}-\boldsymbol{Q}$
$$ \theta = \cos^{-1} \left( \frac{ (\boldsymbol{C}-\boldsymbol{Q}) \cdot ( \boldsymbol{A}-\boldsymbol{Q}) }{ \| \boldsymbol{C}-\boldsymbol{Q}\| \,\|\boldsymbol{A}-\boldsymbol{Q} \| } \right) \tag{10} $$
Example
Points:
$$\begin{aligned} \boldsymbol{A} &= \pmatrix{1\\0\\0} & \boldsymbol{B} & = \pmatrix{3\\1\\0} & \boldsymbol{C} &= \pmatrix{2 \\0 \\ -1} \end{aligned} $$
Normal:
$$ \boldsymbol{n} = \pmatrix{ -\tfrac{\sqrt{6}}{6} \\ \tfrac{\sqrt{6}}{3} \\ -\tfrac{\sqrt{6}}{6} } $$
Distance from Origin
$$ d = -\tfrac{\sqrt 6}{6} $$
Rotation:
$$ \mathbf{R} = \begin{bmatrix} \tfrac{\sqrt 3}{3} & \tfrac{\sqrt 2}{2} \\ \tfrac{\sqrt 3}{3} & 0 \\ \tfrac{\sqrt 3}{3} & -\tfrac{\sqrt 2}{2} \end{bmatrix} $$
2D Points:
$$ \begin{aligned} \boldsymbol{a} & = \pmatrix{ \tfrac{\sqrt 3}{3} \\ \tfrac{\sqrt 2}{2} } & \boldsymbol{b} & = \pmatrix{ \tfrac{4 \sqrt 3}{3} \\ \tfrac{3 \sqrt 2}{2} } & \boldsymbol{c} & = \pmatrix{ \tfrac{\sqrt 3}{3} \\ \tfrac{3 \sqrt 2}{2} } \end{aligned} $$
System of equations:
$$\left. \begin{bmatrix} \sqrt{3} & \sqrt{2} \\ 0 & \sqrt{2} \end{bmatrix} \pmatrix{x \\ y} = \pmatrix{\tfrac{9}{2} \\ 2 } \;\right\} \; \pmatrix{x = \tfrac{5 \sqrt 3}{6}\\y = \sqrt 2} $$
Circle Center:
$$ \boldsymbol{Q} = \left( -\tfrac{\sqrt 6}{6}\right) \pmatrix{ -\tfrac{\sqrt{6}}{6} \\ \tfrac{\sqrt{6}}{3} \\ -\tfrac{\sqrt{6}}{6} } + \begin{bmatrix} \tfrac{\sqrt 3}{3} & \tfrac{\sqrt 2}{2} \\ \tfrac{\sqrt 3}{3} & 0 \\ \tfrac{\sqrt 3}{3} & -\tfrac{\sqrt 2}{2} \end{bmatrix} \pmatrix{ \tfrac{5 \sqrt 3}{6}\\ \sqrt 2} = \pmatrix{2 \\ \tfrac{1}{2} \\ 0} $$
Prove that points A, B, and C are equidistant to Q
$$ \mathrm{radius} = \| \boldsymbol{A}-\boldsymbol{Q} \| = \| \boldsymbol{B}-\boldsymbol{Q} \| = \| \boldsymbol{C}-\boldsymbol{Q} \| = \tfrac{\sqrt 5}{2} \;\;\checkmark $$
PS. I am also a fellow roboticist. You do need a basic understanding of linear algebra operations as well as some vector calculus to understand the how and why of things in robotics.