All points \(\bm{p} = (x, y, z)\) are on a sphere with radius \(r\) centered around \(\bm{o} = (0, 0, 0)\) if the following equation holds:
$$
x^2 + y^2 + z^2 = r^2
$$

If the center of the sphere \(\bm{c} = (c_x, c_y, c_z)\) is not the origin of our coordinate system, we need to subtract the coordinates of the center from the coordinates of the points in the previous equation:
$$
(x - c_x)^2 + (y - c_y)^2 + (z - c_z)^2 = r^2
$$

More compactly, we can express this equation using the dot product of the difference between vectors \(\bm{p}\) and \(\bm{c}.\)
$$
(\bm{p} - \bm{c}) \cdot (\bm{p} - \bm{c}) = r^2
$$

Rays

We represent our rays with parametric line equations of the form
$$
\bm{r}(t) = \bm{o} + t \bm{d}
$$
where \(\bm{o}\) is the origin of our coordinate system as above and \(\bm{d}\) is a directional vector, which may or may not be a unit vector.

Intersection

Since we want to find intersection points between sphere and rays, we insert \(\bm{r}(t)\) for \(\bm{p}\) in the sphere equation.
$$
(\bm{o} + t \bm{d} - \bm{c}) \cdot (\bm{o} + t \bm{d} - \bm{c}) = r^2
$$
The only unknown here is our parameter \(t.\)

Now we regroup the terms a bit and bring everything to the left side.

We expand the product of the difference, which is fine because the dot product is distributive over vector addition.
$$
\begin{align*}
& (\bm{o} + t \bm{d} - \bm{c}) \cdot (\bm{o} + t \bm{d} - \bm{c})\\
=\hphantom{.} & \bm{o} \cdot \bm{o} + \bm{o} \cdot t\bm{d} - \bm{o} \cdot \bm{c} \\
&+ t\bm{d} \cdot \bm{o} + t\bm{d} \cdot t\bm{d} - t\bm{d} \cdot \bm{c} \\
&- \bm{c} \cdot \bm{o} - \bm{c} \cdot t\bm{d} + \bm{c} \cdot \bm{c}
\end{align*}
$$

Group these nine terms according to the parameter: one term contains \(t^2\), four terms contain \(t\) and the remaining four do not contain \(t.\)
We simplify further and bring all terms to the left side.

This is just a quadratic equation in standard form \(at^2 + bt + c = 0\) with
$$
\begin{align*}
a &= \bm{d} \cdot \bm{d}\\
b &= 2 \bm{d} \cdot (\bm{o} - \bm{c})\\
c &= (\bm{o} - \bm{c}) \cdot (\bm{o} - \bm{c}) - r^2,
\end{align*}
$$
which we solve with the quadratic formula
$$
t = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}.
$$

Since \(b\) is of the form \(b = 2 b'\), we can simplify the quadratic formula.
$$
\begin{align*}
&\frac{-b \pm \sqrt{b^2 - 4ac}}{2a}\\
=& \frac{-2b' \pm \sqrt{(2b')^2 - 4ac}}{2a}\\
=& \frac{-2b' \pm \sqrt{4b'^2 - 4ac}}{2a}\\
\end{align*}
$$
$$
t = \frac{-b' \pm \sqrt{b'^2 - ac}}{a}
$$

The term under the square root is known as the discriminant.
If it is negative, we know that there are no (real) solutions.
If it's zero, there is one solution.
If it's positive, there are two.