I know there is a well known algorithm which uses the circle equation to approximate it with pixels.
However, I wanted to approach this problem from the most basic principles.
So we start with a square with the side $a_0$ (all lengths are measured in pixels).
To make it more 'circular' we add several pixel rows on each side. The number of rows $k$ is defined to make the distance between the opposite sides as close as possible to the length of the diagonal.
$$(a_0+2k)^2=2a_0^2$$
$$\left[\frac{\sqrt{2}-1}{2} a_0 \right] \leq k \leq \left[\frac{\sqrt{2}-1}{2} a_0 \right]+1$$
Here $[]$ is the floor function.
Now we need to find the length of each row. Let's consider a rectangle formed by two opposite rows. Its diagonal should be as close as possible to the diagonal of the initial square.
$$a_n^2+(a_0+2n)^2=2a_0^2$$
$$\left[ \sqrt{a_0^2-4na_0-4n^2} \right] \leq a_n \leq \left[ \sqrt{a_0^2-4na_0-4n^2} \right]+1$$
Now if we use the 'exact' value of $k$, we get:
$$k=\frac{\sqrt{2}-1}{2} a_0~~~~\rightarrow~~a_k=0$$
Which means that the correct value for $k$ is:
$$k=\left[\frac{\sqrt{2}-1}{2} a_0 \right]$$
Value of $a_n$ is defined from symmetry considerations:
If $a_0$ is even, all $a_n$ should be even.
If $a_0$ is odd, all $a_n$ should be odd.
If we count the number of pixels in our circle (its discrete area), we get:
$$S_k=a_0^2+4 \sum_{n=1}^k a_n$$
$$\pi=\lim_{k \to \infty} \frac{2 S_k}{a_0^2}$$
And here are two examples:
$$a_0=100~~~~~~~~~k=20~~~~~~~~~S=15496~~~~~~~~~\pi = 3.0992$$
$$a_0=10000~~~~~~~~~k=2071~~~~~~~~~S=157059544~~~~~~~~~\pi = 3.14119$$
This sequence converges to $\pi$ from below. The convergence is stochastic and very slow. See the plot below for first $1000$ terms:
Is this algorithm correct (does it create a 'best' circle approximation for a given $a_0$)? Is there some other way to approximate $\pi$ using this algorithm?
If we just want to compute $\pi$ (without creating the circle) we can simplify the algorithm and get the following:
$$\pi=2+4 \lim_{p \to \infty} \frac{1}{p} \sum_{n=1}^k \sqrt{1-\frac{2n}{p}-\frac{n^2}{p^2}}$$
$$k=\left[(\sqrt{2}-1) p \right]$$
Have you seen this limit before? Maybe it's related to some known formula for $\pi$?
You state that
$\pi=2+4 \lim_{p \to \infty} \frac{1}{p} \sum_{n=1}^k \sqrt{1-\frac{2n}{p}-\frac{n^2}{p^2}} $
This is a standard Riemann sum.
$\lim_{p \to \infty} \frac{1}{p} \sum_{n=1}^k \sqrt{1-\frac{2n}{p}-\frac{n^2}{p^2}} \to \int_0^{\sqrt{2}-1} \sqrt{1-2x-x^2} dx $
Since $ \int \sqrt{1-2 x-x^2} dx = \frac12 \sqrt{-x^2-2 x+1} (x+1)+\sin^{-1}\frac{x+1}{\sqrt{2}} $, and $1-2x-x^2 = 0 $ for $x = \sqrt{2}-1$ ($x^2+2x =3-2\sqrt{2}+2\sqrt{2}-2 = 1 $),
$\begin{array}\\ \int_0^{\sqrt{2}-1} \sqrt{1-2 x-x^2} dx &= \left(\frac12 \sqrt{-x^2-2 x+1} (x+1)+\sin^{-1}\frac{x+1}{\sqrt{2}}\right)\big|_0^{\sqrt{2}-1}\\ &=(0+ \sin^{-1}\frac{(\sqrt{2}-1)+1}{\sqrt{2}})-(\frac12+\sin^{-1}\frac1{\sqrt{2}}) \\ &=\sin^{-1}(1)-(\frac12+\sin^{-1}\frac1{\sqrt{2}})\\ &=\frac{\pi}{2}-(\frac12+\frac{\pi}{4})\\ &=\frac{\pi}{4}-\frac12\\ \end{array} $
which verifies your limit.
The standard digital circle drawing algorithm is Bresenhams's.
For example, see here:
https://en.wikipedia.org/wiki/Midpoint_circle_algorithm