Computing the centers of Apollonian circle packings The radii of an Apollonian circle packing are computed from the initial curvatures e.g. (-10, 18, 23, 27) solving Descartes equation $2(a^2+b^2+c^2+d^2)=(a+b+c+d)^2$ and using the four matrices to generate more solutions
$$
 \left[\begin{array}{cccc} -1 & 2 & 2 & 2 \\\\ 0 & 1 & 0 & 0 \\\\ 0 & 0 & 1 & 0 \\\\ 0 & 0 & 0 & 1 \end{array}\right] \hspace{0.25 in} 
\left[\begin{array}{cccc}  1 & 0 & 0 & 0 \\\\ 2 & -1 & 2 & 2 \\\\ 0 & 0 & 1 & 0 \\\\ 0 & 0 & 0 & 1 \end{array}\right] \hspace{0.25in}
\left[\begin{array}{cccc}   1 & 0 & 0 & 0 \\\\ 0 &  1 & 0 & 0\\\\ 2 & 2 & -1 & 2 \\\\ 0 & 0 & 0 & 1 \end{array}\right]
\hspace{0.25in}
\left[\begin{array}{cccc}   1 & 0 & 0 & 0 \\\\ 0 &  1 & 0 & 0\\\\  0 & 0 &  1 & 0 \\\\  2 & 2 &  2&-1  \end{array}\right]
$$
How to compute the centers of circles in the Apollonian circle packing?
The formulas probably simplify if you use complex numbers.
Also in what sense it is the circle packing the limit set of a Kleinian group?

 A: If from the initial circles $a,b,c,d$ with curvatures $k_a, k_b, k_c, k_d$ and centers $c_a, c_b, c_c, c_d$ you generate four new circles with curvatures:
reflection of $a$: $k_e =  2(k_b+k_c+k_d) - k_a$  
reflection of $b$: $k_f =  2(k_c+k_d+k_a) - k_b$ 
reflection of $c$: $k_g =  2(k_d+k_a+k_b) - k_c$ 
reflection of $d$: $k_h =  2(k_a+k_b+k_c) - k_d$
Those new circles have centers (complex plane):
$c_e = \dfrac{2(k_bc_b +k_cc_c +k_dc_d) - k_ac_a}{k_e}$
$c_f = \dfrac{2(k_cc_c +k_dc_d + k_ac_a) - k_bc_b}{k_f}$
$c_g = \dfrac{2(k_dc_d +k_ac_a+k_bc_b) - k_cc_c}{k_g}$
$c_h = \dfrac{2(k_ac_a +k_bc_b +k_cc_c) - k_dc_d}{k_h}$
After the first step, each new circle can generate three more unique circles, e.g., if a circle $e$ is created with $a, b, c, d$ as above then three new circles can be generated with:
$k_i =  2(k_c+k_d+k_e) - k_b$
$k_j =  2(k_b+k_d+k_e) - k_c$
$k_k =  2(k_c+k_b+k_e) - k_d$
$c_i = \dfrac{2(k_cc_c +k_dc_d +k_ec_e) - k_bc_b}{k_i}$
$c_j = \dfrac{2(k_bc_b +k_dc_d + k_ec_e) - k_cc_c}{k_j}$
$c_k = \dfrac{2(k_cc_c +k_bc_b+k_ec_e) - k_dc_d}{k_k}$
...and so on ad infinitum, with the number of circles going up a factor of three for each step.
see: Beyond the Descartes circle theorem
A: To see a link between Apolonnian packings and Kleinian groups, I suggest you to have a look on the work of Hee Oh and her collaborators: 
http://www.math.brown.edu/~heeoh/
A: Have a look at chapter 7 of Indra's Pearls. The equivalence between the matrix description and the Appolonian picture follows from a special isomorphism of Lie groups $PO(3,1)\cong PSL(2,\mathbb{C})\cup c PSL(2,\mathbb{C})=Isom(H^3)$, where $c$ corresponds to complex conjugation. The matrices you've written give reflection generators of the symmetries in $PO(3,1)$, which can be translated via the isomorphism into symmetries of the Appolonian gasket from the action of $Isom(H^3)$ on $\mathbb{CP}^1$. 
As for finding the centers of the circles, I don't know a recursive formula, however if you need to draw a picture of the Appolonian gasket, you need to be able to compute the centers of the circles, so clearly many people have solved this problem. I would suggest looking again at Indra's Pearls to see if they explain how to compute the circles, or look at the source code of McMullen's program lim which generates postcript to draw limit sets of Kleinian groups. 
A: 
"How to compute the centers of circles in the Apollonian circle packing?"

The paper
"Apollonian Circle Packings: Geometry and Group Theory:
I. The Apollonian Group" (arXiv:math/0010298v5), by
Ronald Graham,
Jeffrey Lagarias, 
Colin Mallows,
Allan Wilks,
and C(atherine?) Yan
gives a joint characterization of the circle curvatures (inverse radii) and coordinates.
For four circles $D=\lbrace C_1,C_2,C_3,C_4 \rbrace$, they define a $4 \times 3$ matrix
$$ N_D = \left( 
\begin{array}{cc} 
1/r_1 & (1/r_1) ( x_{1,1}, x_{1,2} )\\
1/r_2 & (1/r_2) ( x_{2,1}, x_{2,2} )\\
1/r_3 & (1/r_3) ( x_{3,1}, x_{3,2} )\\
1/r_4 & (1/r_4) ( x_{4,1}, x_{4,2} )
\end{array} 
\right) 
$$
where $r_i$ are the circle radii and
$(x_{i,1},x_{i,2})$ are the circle coordinates. This matrix satisfies
$$
\frac{1}{2} N^T_D Q_2 N_D =
\left(
\begin{array}{ccc}
0 & 0 & 0 \\
0 & 2 & 0 \\
0 & 0 & 2
\end{array}
\right)
$$
where $Q_2$ is the Descartes quadratic form (a $4 \times 4$ matrix of $\pm 1$s),
and $\sum_i (1/r_i) > 0$.
They call this result the Complex Descartes Theorem.
I believe these "curvature$\times$center coordinates" are the method of choice to compute
the circle centers.  So not only are all curvatures integers (the numbers in your image),
but the circle centers are Gaussian integers rationals (if placed intelligently).
A: $$2(a^2+b^2+c^2+d^2)=(a+b+c+d)^2$$
$$a=4k(k+s)$$
$$b=4s(k+s)$$
$$c=p^2+k^2+s^2+2pk+2ps-2ks$$
$$d=p^2+k^2+s^2-2pk-2ps-2ks$$
