Is there a good algorithm to divide two integers without using division directly?

Problem. Given positive integers $a$ and $b$, obtain $\frac{a}{b}$ without using division ($/$) directly, though addition ($+$), subtraction ($-$), multiplication ($\times$) and bit-shifts ($\gg$ and $\ll$) are allowed.

Is there a good algorithm to solve the aforementioned problem? This algorithm does not have to be very accurate and $4$ decimal places would be good enough, e.g.,

$$\frac{222}{103} \approx 2.1553$$

Could anyone please give some advice or references?

Newton's method is pretty good for this. Specifically, a good strategy is to first calculate an approximation $x \approx \dfrac{1}{b}$ so that $ax$ is then a good approximation for $\dfrac{a}{b}$. To do this, we will use Newton's method to approximate a root of the function $f(x) = \dfrac{1}{x} - b$.

Newton's method asks us to recursively calculate approximations $x_n$ given by the formula
$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$
A priori, this looks like using division but since $f'(x) = \dfrac{-1}{x^2}$, the formula boils down to
$$x_{n+1} = x_n + x_n^2\left(\frac{1}{x_n} - b\right) = x_n + (x_n - bx_n^2) = (2 - bx_n)x_n,\tag{$*$}$$ which only involves addition and multiplication.

In order to ensure good convergence of this method, it helps to have $b$ very close to $1$. Of course, there are very few integers that are close to $1$. However, if $b$ is an $n$-bit integer then $b_0 = \dfrac{b}{2^{n-1}}$ is between $1$ and $2$. Since the denominator is a power of $2$ this is really just a bit shift and not a division per se. Incorporating this in the formula ($*$) we obtain the iteration
$$
x_{n+1} = (2 - b_0x_n)x_n = \left(2 - \frac{bx_n}{2^{n-1}}\right)x_n = \frac{(2^n - bx_n)x_n}{2^{n-1}}
\tag{$\dagger$}$$
which still involves no division. Note that this gives a method for approximating $1/b_0$ but we can recover a good approximation to $\dfrac{1}{b} = \dfrac{1/b_0}{2^{n-1}}$ by shifting.

Further analysis shows that, starting with $x_0 = 0.75$, the absolute error $\varepsilon_n = \left|x_n - \dfrac{2^{n-1}}{b}\right|$ for $(\dagger)$ satisfies $\varepsilon_{n+1} \leq \varepsilon_n^2$. So the number of correct bits doubles at each step!

Let's illustrate by approximating $22/7$. Then $n=3$ and $b_0 = 7/4 = 1.75$. Then $(\dagger)$ gives the approximations
$$\begin{aligned}
x_0 &= 0.750\,000\,000 \\
x_1 &= 0.515\,625\,000 \\
x_2 &= 0.565\,979\,004 \\
x_3 &= 0.571\,376\,600 \\
x_4 &= 0.571\,428\,567 \\
x_5 &= 0.571\,428\,571 \\
\end{aligned}$$
All the digits shown for $x_5$ are exact for
$$
\frac{1}{b_0} = \frac{4}{7} = 0.571\,428\,571\,428\,571\,428\,571\,428\ldots
$$
Thus
$$
\frac{22x_5}{4} = 3.142\,857\,141
$$
is a good approximation to
$$
\frac{22}{7} = 3.142\,857\,142\,857\,142\,857\,142\,857\ldots
$$

StackExchange Computer Science Q#142623, answer score: 5