DJ3EI's blog

What is this about? 🔗

You want to divide two numbers with your computer.

No problem: Just let your computer do it!

But in rare cases, your computer may not know how to divide.

This could be the case if said computer is a mere µcontroller. Even more so if you do your own arithmetic, e.g., BCD or fixed point.

So this post assumes your computer knows how to do addition $+$, subtraction $-$, multiplication $*$, and comparison with $<$, but it does not yet know division $/$. It offers a software solution that teaches your computer division.

I understand one normal expert advice for teaching computers division is to use the CORDIC algorithm family (which can even be used to teach your computer multiplication, in case it has not learned that yet, either). But that CORDIC stuff is quite a beast to wrap your head around (which I, honestly, have not really tried yet). Apparently, it also needs a bit of ROM space for lookup tables. Finally, from my present partial understanding, CORDIC seems to be intended for binary numeric representation; the approach offered here can be used for BCD numbers unchanged.

So this toot dares to suggest a quite simple home-grown approach to the division problem.

That approach has a reasonably decent runtime performance, though it may not be optimal. It results in short program code which is easy to get right. And it doesn’t need lookup tables.

And, of course: If you just want to read this for the fun of it, without any particular application in mind, you are warmly invited to do so.

Full disclosure: I wrote this for the fun of it, without any particular application in mind.

In a nutshell 🔗

You are given a number $a>0$ and desire to calculate $\frac{1}{a}$.

(That is actually all that is missing to do general division: Instead of dividing something by $a$, multiply it with $\frac{1}{a}$.)

To do so, start with any value $x_0$ that is a reasonable approximation of $\frac{1}{a}$.

Here, the requirements for “reasonable” are lax. Any $x_0$ with $0<x_0<\frac{2}{a}$ will work. (But the better your initial approximation is, the faster the calculation completes.)

From your start value $x_0$, define a sequence $x_0,x_1,x_2,\dots$ of better and better approximations for $\frac{1}{a}$ via the simple formula

Here is a numerical example with $a=3$, using $x_0=0.5$ as an initial approximation for $\frac{1}{3}$.

The sequence stabilizes with all subsequent values being identically to $x_6$. That value incidentally is the best that can be achieved, as the next higher value that can be represented by the floating point number scheme used is 0.333333333333333370, which is less precise.

In all numerical calculations presented here, I used whatever floating point number scheme my PC’s libraries offer. I believed this to be the ubiquitous IEEE 754 binary64, but I did nothing to verify this.

Once $x_2$ has found the first correct digit, each of the few next subsequent values has twice as many correct digits as its predecessor. This is typical of the quadratic convergence afforded by the method presented here.

When to stop? 🔗

One might be tempted to stop when $x_n=x_{n+1}$. This will work most of the time, but not always. In rare cases, the sequence will oscillate.

One example found by brute-force search is $a=338$ with start value $x_0=0.002$:

The value returned by the build-in division in this case prints as 0.0029585798816568047.

It is conceivable that in some cases, the sequence might end up in a merry-go-round of three or even more numbers. With a brief brute-force search, I did not find any occurrences of that, but it should not be ruled out.

In the realm of pure mathematics and precise real numbers, the algorithm employed guarantees an increasing sequence $x_1<x_2<x_3<\dots$. (This statement purposefully omits the possibly too-high start-value $x_0$.)

In the practical realm of computers with finite-precision arithmetic, I suggest as a straightforward termination criterion $x_n\ge x_{n+1}$. When that triggers, a good finite-precision approximation of the desired result $\frac{1}{a}$ is given by $x_{n+1}$.

For a posh twist, you can use that only if $x_n=x_{n+1}$, and $0.5(x_{n+1}+x_{n+2})$ if not. This results in a marginal increase in typical precision, at the expense of a marginal increase in runtime.

Choice of start value 🔗

One can try to find a good start value $x_0$ by analyzing the given number $a$ to be inverted.

Where floating point numbers are used, analysis of the representation of $a$ yields information about the magnitude of $a$, from which a rough approximation of its inverse can be obtained that is good enough to serve as $x_0$.

A more pedestrian approach starts with the value $\mathrm{x}=1$ and calculates the product $\mathrm{x}\mathrm{a}$. If that product is close enough to $1$, then $\mathrm{x}$ can be used as is. If not, $\mathrm{x}$ can be either increased or decreased as pertinent, and $\mathrm{x}\mathrm{a}$ calculated anew.

This pedestrian approach is used in the following sample implementation:

Sample implementation 🔗

Here is a sample implementation for this algorithm in Python:

def sample_implementation(a):
    assert 0 < a

    # Find decent start value x
    # between 0.17/a and 1.8/a:
    x = 1
    xa = x * a
    while 1.8 < xa:
        x *= 0.1
        xa = x * a
    while xa < 0.17:
        x *= 10
        xa = x * a
    # The multiplication factors
    # of 0.1 and 10 in the above
    # have been pulled
    # out of thin air.

    # They'll work, but better
    # choices may be possible.

    # If representation is binary,
    # using the values 0.125 and 8
    # will most likely
    # speed up calculations.

    # Omit the possibly too high
    # start value x, proceed to the
    # next value in the sequence
    # right away:
    x = x * (2 - xa)

    while True:
        x_new = x * (2 - x * a)
        if x_new <= x:
            return x_new
        x = x_new

To verify, the results of this implementation where compared with that of the system’s floating point division for 1,000,000 values in the range 1e-300 through 1e+300.

The results from this implementation where either the same as that computed by the system’s 1/a, or it was the next lower floating point value that can be represented, or the next higher. In no case was the result from this implementation further off the system’s, it was always one of those three.

In percentages:

For a further experiment, the posh implementation of calculating the return value was tried:

if x_new <= x:
    if x_new == x:
        return x_new
    else:
        return 0.5 * (
          x_new +
          x_new * (2 - x_new * a)
        )

This was only slightly better:

Suggested exercises 🔗

• Implement this algorithm in a programming language of your choice and play around with it.
• Extend the algorithm so it can deal with negative $a$ and also does whatever the appropriate error handling for your platform is when $a=0$.

The math behind this: a theorem. 🔗

Theorem’s assumptions 🔗

Assume a real-valued function $f$ is defined over the open non-empty interval of real numbers between two numbers $u$ and $v$ and is twice continuously differentiable over that interval.

Assume a point $z$ exists in that interval, that is, $u<z<v$, which is a fixpoint of $f$, that is, $f\left(z\right)=z$.

Assume further $f^{\prime }\left(z\right)=0$, and also $0\le f^{\prime }\left(x\right)$ for all $x$ such that , and finally, assume that no fixpoint of $f$ exists larger than $u$ and smaller than $z$.

Theorem’s claims 🔗

• Then, any $x$ such that $u<x<z$ fulfills $x<f\left(x\right)\le z$.

Assuming that first claim to be true, given any $x_0$ such that $u<x_0\le z$, we can define a sequence recursively via $x_{n+1}=f\left(x_n\right)$ for all indices (natural numbers) $n$.

This sequence has the following properties:

• $u<x_n\le z$ will hold for any index $n$. (This helps the sequence definition to make sense.)
• This is an increasing sequence, that is, $x_n\le x_{n+1}$ for any index $n$.
• This sequence converges to $z$.
• Finally, it does so in quadratic fashion, that is, there exist some number $m$ and some constant $k>0$ such that $\left|x_{n+1}-z\right|\le k{|x_n-z|}^2$ for any index $n$ such that $n>m$.

Application 🔗

Assuming for a moment the theorem given were true. It can be applied to our inversion problem as follows:

• $u=0$, $v=\frac{2}{a}$.
• $f\left(x\right)=x(2-ax)$ is a quadratic parable with zeros at $0$ and $\frac{2}{a}$, opening downwards.
• $z=\frac{1}{a}$ is clearly a fixpoint of $f$.
• A straightforward calculation shows that the only other fixpoint of $f$ is $0$.
• The maximal value of $f$ is in the middle between its zeros, so at $z$.
• The first derivative of $f$ at $z$ is $0$, as $f$ has its maximum there.
• The first derivative of $f$ has its only zero at $z$.
• As $f$ opens downwards, the derivative has to be positive before $z$ and negative thereafter.
• The theorem gives us that our sequence is increasing and we have quadratic convergence, for any start value $x_0$ such that $0<x_0\le \frac{1}{a}$.
• Now $f$ is symmetrical around $\frac{1}{a}$, so we also have quadratic convergence for any start value $x_0$ such that $\frac{1}{a}\le x_0<\frac{2}{a}$.

Suggested exercise 🔗

Convince yourself (in a “pedestrian” way) that for any $x$ with $0<x<\frac{1}{a}$ the assertion $x<f\left(x\right)\le z$ holds in our particular $\frac{1}{a}$ situation.

Suggested exercise 🔗

Convince yourself that the theorem allows the use of the function $f\left(x\right)=x+\sin \left(x\right)$ to calculate $\pi$.

(There are much more convenient ways to calculate $\pi$. But this way is at least legitimate, as the $\sin$ function can be calculated via its power series without prior knowledge of $\pi$.)

Suggested exercise 🔗

In music, the “equal temperament” demands that the quotient of the frequencies of two adjacent (half-)tones is $z=\sqrt[12]{2}$.

• Convince yourself that, for any real $\lambda$, this number $z$ is a fixpoint of $f_{\lambda }\left(x\right)=(1-\lambda )x+\lambda (x^{13}-x)$.
• Find a value $\lambda$ that might be suitable to make the theorem applicable.
• Verify that the theorem indeed is applicable.
• Calculate $z$ via $f_{\lambda }$.

Proof of the theorem 🔗

We first need to prove: Given $x$ such that $u<x<z$, it is true that $x<f\left(x\right)\le z$.

For $f\left(x\right)\le z$: According to the mean value theorem, there must exist some $\xi$ with $x<\xi <z$ such that $f^{\prime }\left(\xi \right)=\frac{f\left(z\right)-f\left(x\right)}{z-x}$. According to our assumptions, $f^{\prime }\left(\xi \right)\ge 0$. As $x<z$, we see $f\left(x\right)\le f\left(z\right)=z$.

For $x<f\left(x\right)$, consider the function $p\left(x\right)=f\left(x\right)-x$. This function is zero precisely where $f$ has a fixpoint. In particular, $p\left(z\right)=0$. Now $f$ has no fixpoint between $u$ and $z$ (exclusively), so $p$ becomes zero nowhere in that interval. Being a continuous function, $p$ cannot change sign over an interval without assuming a zero value in that interval. So $p$ cannot change sign between $u$ and $z$ (exclusively), in other words, must either be always negative or always positive there. Now $p$ is continuously differentiable. Its derivative evaluates to $-1$ at $z$. As the derivative is continuous, it must be negative for some interval of points slightly smaller than $z$. By another application of the mean value theorem, $p$ must be positive on that same interval. We established $p$ does not change sign in the interval strictly between $u$ and $z$, and we have found points in that interval where $p$ is positive, so we now know it must be positive everywhere between $u$ and $z$ (exclusively), which in particular proves $p\left(x\right)>0$. Now $p\left(x\right)>0$ is equivalent to $x<f\left(x\right)$.

These two together prove that the sequence increases and is bounded from above by $z$.

Any increasing bounded sequence of real numbers converges. (How to prove this fact, which is both profound and basic, depends on your definition of real numbers.)

Let us call the limit $\zeta =\underset{n\to \infty }{\lim }{x}_n$. As $f$ is continuous:

So $\zeta$ is a fixpoint of $f$, which forces $\zeta =z$. So the sequence $x_0,x_1,x_2,\dots$ converges with limit $z$.

The only claim not yet proven is that this convergence is quadratic.

This is a direct consequence of Taylor’s theorem. It allows us to rewrite $f\left(x_n\right)=f\left(z\right)+f^{\prime }\left(z\right)(x_n-z)+R_1(x_n)$

In our case, $f^{\prime }\left(z\right)=0$, so this reduces to $f\left(x_n\right)=f\left(z\right)+R_1(x_n)$ or $f\left(x_n\right)-f\left(z\right)=R_1(x_n)$ or $x_{n+1}-z=R_1(x_n)$. The Lagrange form of the remainder boils down to $R_1(x_n)=\frac{f^{″}\left(\xi \right)}{2}{(x_n-z)}^2$ for some suitable value $\xi$ with $x_n\le \xi \le z$. Now $f^{″}$ is continuous, hence bounded over any closed interval containing $z$. If we have an upper bound $k$ of $\left|f\right|$ over some such interval and let $n$ become large enough so all subsequent $x_n$ are contained in that interval, then indeed $|x_{n+1}-z|\le \frac{k}{2}{|x_n-z|}^2$, which completes the proof.

Suggested exercise 🔗

Formulate the mirror-image theorem, where convergence ends up not being from below, but from above, and prove it.

Use your theorem to prove that Heron’s method works, which uses the recursion formula $x_{n+1}=\frac{1}{2}(x_n+\frac{a}{x_n})$ to calculate $\sqrt{a}$.

Use your theorem to prove that recursion via the function $f\left(x\right)=x+\sin \left(x\right)$ converges to $\pi$ for any start value $x_0$ with $\pi \le x_0<2\pi$.

Suggested exercise 🔗

Above, there is a definition of “quadratic convergence” that loosely means that after a certain number of digits of the final result have been calculated correctly, the next iteration step tends to produce twice as many correct digits. There are approximation sequences that exhibit cubic convergence: The next iteration step tends to produce three times as many correct digits as the previous one has established.

• Give a formal definition of cubic convergence.
• Show that, given any start value $x_0$ with $0<x_0<2\pi$, the sequence defined by $x_{n+1}=x_{\mathrm{n}}+\sin x_{\mathrm{n}}$ exhibits cubic convergence behaviour, converging towards $\pi$.
• Show that, given any $\mathrm{a}>0$ and any start value $x_0$ with $0<x_0<\sqrt{\frac{7}{3}a}$, the sequence defined by $x_{n+1}=\left(\right(\left(\frac{3}{8a^2}{x}_n^2\right)-\frac{5}{4a})x_n^2+\frac{15}{8})x_n$ exhibits cubic convergence behaviour, converging to $\sqrt{a}$. – When all subterms that do not involve $x_n$ are computed only once beforehand, four multiplications and two additions/subtractions need to be done in each iteration of the loop. It depends on circumstances whether this is faster or slower, compared with the ancient Heron’s method as mentioned above. 🔗

Author’s remark 🔗

This was also an experiment regarding typing math in via raw mathml. Result: Can be done, but is not entirely convenient.

I did not bother installing web fonts, but rely on your browser’s support of mathml. Depending what browser you use, your mileage may vary.

Discussion opportunity 🔗

If you want to comment or discuss this piece and have a Fediverse account, feel invited to answer my pertinent toot.