Finding square root without division and initial guess
OK, let's go with floating point. Any positive floating point number is of the form $m \times 2^e$ where $m$ is in the range $[1,2)$ and $e$ is an integer.

If $e$ is even, then:

$$\sqrt{m \times 2^{2p}} = \sqrt{m} \times 2^p$$

And if $e$ is odd, then:

$$\sqrt{m \times 2^{2p+1}} = \sqrt{m} \times \sqrt{2} \times 2^p$$

We can precompute $\sqrt{2}$. So we have reduced the problem to finding a square root of a number between 1 and 2. So far, so good.

We could use the Newton-Raphson method for square roots, but you specifically didn't want division, and it's a bad idea anyway; division is more expensive than multiplication on a modern FPU.

But we could use the Newton-Raphson method for reciprocal square roots. The iteration looks like this. If $x_0$ is an initial approximation to $\frac{1}{\sqrt{m}}$, then iterate:

$$x_{i+1} = \frac{1}{2} x_i \left( 3 - m x_i^2 \right)$$

When this has converged, $m x_k$ is an approximation to $\sqrt{m}$. No divisions here, and the problem is now just to come up with an initial estimate.

But wait, we can do better.

Given $b_0 = m$, let's suppose we could come up with a sequence of numbers $Y_i$ such that $b_i = b_0 Y_0^2 Y_1^2 \cdots Y_n^2$ approaches $1$. Then $x_i = b_0 Y_0 \cdots Y_{i-1}$ approaches $\sqrt{m}$ and $y_i = Y_0 \cdots Y_{i-1}$ approaches $\frac{1}{\sqrt{m}}$. Please pause for a moment and convince yourself that this is the case.

We already know such a sequence, because we can use the Newton-Raphson iteration for reciprocal square roots:

$$b_i = b_{i-1} Y_{i-1}^2$$
$$Y_i = \frac{1}{2}(3 - b_i)$$

And so:

$$x_0 = m Y_0$$
$$x_{i} = x_{i-1} Y_i$$

and:

$$y_0 = Y_0$$
$$y_{i} = y_{i-1} Y_i$$

But observe that $b_i = x_{i-1} y_{i-1}$. So we don't actually have to retain $b_i$:

$$Y_i = \frac{1}{2}(3 - x_{i-1} y_{i-1}) = 1 + \frac{1}{2}(1 - x_{i-1} y_{i-1})$$

Substituting this into the equations for $x_i$ and $y_i$ gives:

$$x_{i} = x_{i-1} \left( 1 + \frac{1}{2}(1 - x_{i-1} y_{i-1}) \right) = x_{i-1} + \frac{1}{2}(1 - x_{i-1} y_{i-1}) x_{i-1}$$
$$y_{i} = y_{i-1} \left( 1 + \frac{1}{2}(1 - x_{i-1} y_{i-1}) \right) = y_{i-1} + \frac{1}{2}(1 - x_{i-1} y_{i-1})y_{i-1}$$

Removing the redundant computation gives:

$$r_i = \frac{1}{2}(1 - x_{i-1} y_{i-1})$$
$$x_{i} = x_{i-1} + r_i x_{i-1}$$
$$y_{i} = y_{i-1} + r_i y_{i-1}$$

But we're not done! We can save the multiplication by a half by substituting $h_i = \frac{1}{2} y_i$.

Given $Y$ as an approxmiation to $\frac{1}{\sqrt{m}}$, set:

$$x_0 = m Y$$
$$h_0 = \frac{1}{2} Y$$

Iterate:

$$r_i = \frac{1}{2} - x_{i-1} h_{i-1}$$
$$x_{i} = x_{i-1} + r_i x_{i-1}$$
$$h_{i} = h_{i-1} + r_i h_{i-1}$$

Then $x_i$ converges to $\sqrt{m}$ and $y_i$ converges to $\frac{1}{2\sqrt{m}}$. Not only does the inner loop have no divisions, it can be implemented on modern FPUs using three fused multiply-adds and nothing else. Perhaps it's unsurprising that this algorithm, known as Goldschmidt's algorithm, is the one used in your FPU.

The final question: How do we come up with $Y$? Well, you could use $\frac{1}{m}$ or something if you want to allow one division. But you said no divisions. So since $m$ is known to be between 1 and 2, you could use its high-order bits to consult a lookup table, or use an interpolating polynomial approximation, such as a Chebyshev polynomial. I'll leave the details of that to you.