Newton-Raphson Fixed Point Division

Sometimes you need to perform fractional arithmetics on platforms without Floating-Point Unit. And such as certain embedded platforms go, you may even need to implement fractional division without any division operation support. That is where Newton-Raphson division algorithm comes handy.

One such "embedded" platform is the uxn virtual processor. It is a 16-bit stack machine which supports both integer multiplication and division, however only in 16 bit integer domain. For certain applications we needed to implement fixed-point arithmetic with at least 16-bit precision. We opted for 16 bits for integer part and 16 bits for the fractional part.

The engineering challenges of implementing addition, subtraction and multiplication are not that interesting. The real fun start with multiplication. The first thing to realize is that a division can be transformed into a multiplication:

$$\frac{N}{D}=N\cdot D^{-1}$$

So what we actually need is to find the reciprocal of the denominator. We chose the Newton-Raphson algorithm, because of its performance even in such restricted computational setting. Firstly the denominator should be normalized to satisfy the following:

$$0.5\le D\le 1.0$$

In order to do this, a simple algorithm can be used that ensures the denominator is in the required range and we know the bit-shift that got it there. Yes, the normalization is basically about getting the integer part to zero (the value of 1.0 is handled specially), ensuring the fractional part's most significant bit is 1 and knowing how many bits the denominator had to be shifted left or right.

We start with $D$ being the original denominator value and the denominator shift $d=0$.

While $D>1.0$ set $D\leftarrow D÷2$ and $d\leftarrow d-1$.

While $D<0.5$ set $D\leftarrow D\cdot 2$ and $d\leftarrow d+1$.

If $d<0$ we will be shifting to the right at the end and if $d>0$ then the shift will have to be to the left.

Secondly an initial estimate has to be calculated. There are two well-known estimates for the values in the normalized range. First one being a linear estimate:

$$x_0=\frac{48}{17}-\frac{32}{17}D$$

But this did not produce as precise results as we hoped for. Therefore we went with the second option - a quadratic estimate:

$$x_0=\frac{140}{33}-\frac{64}{11}D+\frac{256}{99}{D}^2$$

And thirdly a single iteraction of the Newton-Raphson estimation algorithm has to be computed as following series:

$$x_{n+1}=x_n(2-D\cdot x_n)$$

With a good enough estimate, only two iterations are needed to achieve the required precision in our case. But how do we compute the estimate, without division? It looks almost like catch 22, right? Not really, as we precompute the values. Given the fact the most significant fractional bit is always one, we opted for precomputing the values for all 8 combinations of the next 3 bits. And surprisingly enough, with this three-bit lookup table and two iterations, the computational error is always at most $\frac{1}{65536}$, which is inevitable anyway.

Hope you enjoyed slightly different topic than usual and see you next time for some other weirdness!