Finding a square root of a power series in terms of another power series I want to find a square root of a power series in the following way:
Let $$A(x)=a_0+a_1x+\cdots=\sum_{j\ge 0}a_jx^j\qquad \mathrm{and}\qquad B(x)=b_0+b_1x+\ldots=\sum_{i\ge 0}b_ix^i$$
such that $a_0>0$ and $b_0>0$ and $A=B^2$.
How can I express the coefficients $a_j$ in terms of $b_i$? I used the Cauchy product to get $a_j=\sum_{k=0}^jb_kb_{j-k}$, but have no idea how to continue.
Can someone help me out with finding a recursive relation between $a_j$ and $b_i$ with initial value for the sequence $(b_i)$?
 A: From the Cauchy product, identify the coefficients by increasing powers.
$$b_0^2=a_0,\\
2b_0b_1=a_1,\\
2b_0b_2+b_1^2=a_2,\\
2b_0b_3+2b_1b_2=a_3,\\
2b_0b_4+2b_1b_3+b_2^2=a_4,\\
\cdots$$
Then
$$b_0=\sqrt{a_0},\\
b_1=\frac{a_1}{2b_0},\\
b_2=\frac{a_2-b_1^2}{2b_0},\\
b_3=\frac{a_3-2b_1b_2}{2b_0},\\
b_4=\frac{a_4-b_2^2-2b_1b_3}{2b_0},\\
\cdots$$
A: There is no known closed form of $\{b_k\}_{k \in \mathbb{N}}$ in terms of $\{a_k\}_{k \in \mathbb{N}}$ in general.
Some particular cases are remarkable. 


*

*For $x \in \mathbb{R}$, we have $A=B^2$, with $$
   A=e^x=\sum_{k=0}^{\infty} \frac{x^k}{k!},\quad
   B=\sqrt{e^x}=e^{\large\frac{x}2}
   = \sum_{k=0}^{\infty} \frac{x^k}{2^k k!}. \tag1$$

*For $|x|<1$, we have $A=B^2$, with $$ A=1+x,\quad
   B=\sqrt{1+x}=(1+x)^{1/2}
   = \sum_{k=0}^{\infty} \binom{1/2}{k}x^k \tag2$$ and generally, for $|x|<1$,  $$ A=(1+x)^{2\alpha},\quad B=(1+x)^{\alpha}
   = \sum_{k=0}^{\infty} \binom{\alpha}{k}x^k \tag3$$ where   $$ \binom{\alpha}{k}:=\frac{\alpha(\alpha-1)\ldots(\alpha-k+1)}{k!},\quad
   \alpha \in \mathbb{R}. \tag4$$

A: We can also use the convolution theorem in Fourier analysis $$\widehat{f*g} = \hat f \cdot \hat g,$$ where $*$ denotes  convolution and $\cdot$ is the ordinary  product.
We can do this because the product of any polynomial or power series corresponds to a discrete convolution of it's coefficients so it becomes a product once we step over into the Fourier domain.
Now what remains to solve $x^2 = c$ in the Fourier domain (for each Fourier coefficient). The trivial solution would just be to choose the same $x = \sqrt{c}$ for some suitable branch of the complex root function. And then lastly to perform the inverse-transform.
A: What you are asking for is a special case of Faa di Bruno's formula, which is a formula expressing the Taylor coefficients of the composition $f(g(x))$ of two functions in terms of the coefficients of the individual functions $f(t)$ and of $g(x)$. Your question is simply about the case where $f(t)=\sqrt{t}$, so just plug in that function into the formula.
Note that Faa di Bruno's formula is fairly complicated -- the $n$th coefficient of $f(g(x))$ is expressed as a sum over the $p(n)$ integer partitions of $n$. Since the number of partitions grows exponentially in $\sqrt{n}$, that can mean a lot of terms to process even for $n$ that's moderately large.
