General solution for intersection of line and circle If the equation for a circle is $|c-x|^2 = r^2$ and the equation for the line is $n \cdot x=d $,  and assuming that the circle and line intersect in two points, how can I find these points? 
Also as kind of a side note, is there also a general formula for obtaining the equation of a circle of the intersection of a plane and sphere in $\mathbb R^3$, using similar equations?

Here I present a method, I don't think it was not used yet, but don't usually find it in mathematics textbooks :
$$\begin{cases}
ax + by + c = 0 \\
(x-\alpha)² + (y-\beta)² = R_{0}^2
\end{cases}
$$
As it shows, the first equation describes the line $\Delta$ with a normal vector $\left( \begin{matrix}
a \\
b
\end{matrix} \right) $ and the second equation describes every point lying on the circumference of the circle with radius $R_0^2$ and a center point lying at $(\alpha,\beta)$.
An easy but very helpful transformation is used here : let $X = x-\alpha $ and $Y = y-\beta$, thus reducing the system of equation to :
$$\begin{cases}
aX + bY = C \\
X² + Y² = R_0^2
\end{cases}
$$
with $ C = -c - a\alpha - b \beta$. (this is equivalent to dealing with the same equations with the center of the circle being the origin of our plane).
let $\left( \begin{matrix}
X \\ Y \end{matrix} \right) \left( \begin{matrix} a & b \\ -b & a \end{matrix} \right) = \left( \begin{matrix} C \\ \gamma \end{matrix} \right)$ with $\gamma \in \mathbb R, \gamma = -bX + aY$.
This expression is equivalent to $z' = z(a-ib)$ with $z' = C + i\gamma$ and $z = X + iY$. Using a basic yet a powerful rule of complex numbers we can safely say that :
$$|z'|² = |z|²|a-ib|² \\
\implies C² + \gamma² = (X² + Y²)(a²+b²)$$
and according to our initial system of equations : $X²+Y²=R_0^2$ and $\gamma = (-bX + aY) $.
$$ \implies (-bX+aY)² = R_0^2(a²+b²)-C²$$
hence :
$$\implies (-bX+aY) = \pm \sqrt{R_0^2(a²+b²)-C²}$$
which explains why we will find
$$\begin{cases}
0 \text{ solutions} & \text{if }   R_0^2(a²+b²)-C²<0 & \implies R_0<\frac{|a\alpha+b\beta+c|}{\sqrt{a²+b²}} \\
1 \text{ solution} & \text{if }   R_0^2(a²+b²)-C²=0 & \implies R_0=\frac{|a\alpha+b\beta+c|}{\sqrt{a²+b²}} \\
2 \text{ solutions} & \text{if }   R_0^2(a²+b²)-C²>0 & \implies R_0>\frac{|a\alpha+b\beta+c|}{\sqrt{a²+b²}}
\end{cases}$$
Back to our linear transformation we have now :
$$\left( \begin{matrix}
X \\ Y \end{matrix} \right) \left( \begin{matrix} a & b \\ -b & a \end{matrix} \right) = \left( \begin{matrix} C \\ \pm \sqrt{R_0^2(a²+b²)-C²} \end{matrix} \right)$$
or equivalently expressed as :
$$(X+iY)(a-ib)=(C \pm i \ \sqrt{R_0^2(a²+b²)-C²})$$
Dealing with complex division is way too easy comparing to matrix inversion, and that's why we keep using them, and now we can directly deduce that
$$X+iY = \frac{(C \pm i \ \sqrt{R_0^2(a²+b²)-C²})(a+ib)}{(a²+b²)} \\
\implies \begin{cases}
X = \frac{1}{a²+b²}(aC \mp b\sqrt{R_0^2(a²+b²)-C²}) \\
Y = \frac{1}{a²+b²}(\pm a\sqrt{R_0^2(a²+b²)-C²} +bC)
\end{cases}\\
\implies \begin{cases}
x = \frac{1}{a²+b²}(aC \mp b\sqrt{R_0^2(a²+b²)-C²}) + \alpha \\
y = \frac{1}{a²+b²}(\pm a\sqrt{R_0^2(a²+b²)-C²} +bC) + \beta
\end{cases} $$
