Calculating square roots of operators using power series for $\sqrt{1 - z}$ The following is a theorem from Reed & Simon's Methods of Modern Mathematical Physics, Volume I. Here, we are working in a complex Hilbert space $(\mathcal{H}, (\cdot, \cdot))$, and $\mathscr{L}(\mathcal{H})$ is the space of bounded linear operators $\mathcal{H} \to \mathcal{H}$. 
Statement of theorem:
Let $A \in \mathscr{L}(\mathcal{H})$ be positive (i.e., $(x, Ax) \ge 0$, all $x \in \mathcal{H}$). Then there is a unique positive $B \in \mathscr{L}(\mathcal{H})$ such that $B^2 = A.$ 
The proof in the text begins as follows:
"It is sufficient to consider the case where $||A|| \le 1$. First observe that 
$$||I- A|| = \operatorname{sup}_{||\phi|| = 1}|((I-A)\phi, \phi)| \le 1.$$
Next we use that fact that $\sqrt{1-z} = 1 + \sum_{k=1}^\infty c_kz^k$ converges absolutely for complex $z$ satisfying $|z| \le 1$, where the constants $c_k$ are known explicitly. This fact implies that the series $1 + \sum_{k=1}^\infty c_k(I- A)^k$ converges in norm to an operator $B$. Since the convergence is absolute, we can square the series and rearrange terms, which proves that $B^2 = A \dots$ 
I am able to verify most of the statements in this passage. However:   

My question is, how does one successfully work out the calculation for squaring the series?

I have been trying for awhile to do the formal multiplication
 $$(1 + c_1(I - A) + c_2(I-A)^2 + \cdots )(1 + c_1(I - A) + c_2(I-A)^2 + \cdots ),$$
using the fact that $0 = 1 + \sum_{k=1}^\infty c_k$. But I am just getting a bunch of messy terms. Is there a good trick to use, or an easier way to realize that $B^2 = A$?
Hints or solutions are greatly appreciated. 
 A: From $\sqrt{1-z} = 1 + \sum_{k=1}^\infty c_kz^k$ for $|z|\leq1$, you get
$$
1-z=\left(1 + \sum_{k=1}^\infty c_kz^k\right)\,\left(1 + \sum_{k=1}^\infty c_kz^k\right)
=\sum_{k,j=0}^\infty c_kc_jz^{k+j}\\ =\sum_{\ell=0}^\infty\left(\sum_{k=0}^\ell c_kc_{\ell-k}\right)\,z^\ell.
$$
This shows that $c_0=1$, $c_0c_1+c_1c_0=-1$, and $\sum_{k=0}^\ell c_kc_{\ell-k}=0$ for all $\ell\geq2$. 
Now, when you square $B$, the same sums between the coefficients $c_k$ arise, and so 
$$
B^2=I-(I-A)=A.
$$
A: There is more known about this special case of Newtons binomial series,
$$\sqrt{1-z}=\sum_{k=0}^\infty \binom{1/2}{k} (-z)^k$$
for $|z|<1$, and since for $k\ge 1$
$$\binom{1/2}{k}=\frac{\tfrac12(\tfrac12-1)(\tfrac12-2)...(\frac12-k+1)}{k!}
=(-1)^{k-1} \frac{(2k-3)(2k-5)...3\cdot 1}{2^k k!}=\frac{(-1)^{k-1}}{4^k(2k-1)}\binom{2k}{k}$$
this is
$$\sqrt{1-z}=1-\sum_{k=1}^\infty \frac1{2k-1}\binom{2k}{k} \left(\frac z4\right)^k$$
The point is that in the first form, squaring in the sense of the Cauchy product reproduces $1-z$, and this stay true if you replace $z$ by a matrix of norm smaller $1$
