title: Division without division and other computational tricks slug: inverse date: 2026-05-22 12:22:29 modified: 2026-05-25 00:17:36 UTC+02:00 type: text tags: nerd_fun ## What is this about? [🔗](#intro){.small} {: #intro} 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](https://en.wikipedia.org/wiki/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 [🔗](#summary){.small} {: #summary} 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}, ...

of better and better approximations for

\frac{1}{a}

via the simple formula

x_{n + 1} = x_{n} (2 - a x_{n})

Here is a numerical example with

a = 3

, using

x_{0} = 0.5

as an initial approximation for

\frac{1}{3}

. | value --|------

x_{0}

| 0.500000000000000000

x_{1}

| 0.250000000000000000

x_{2}

| 0.312500000000000000

x_{3}

| 0.332031250000000000

x_{4}

| 0.333328247070312500

x_{5}

| 0.333333333255723119

x_{6}

| 0.333333333333333315

x_{7}

| 0.333333333333333315

x_{8}

| 0.333333333333333315 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](https://en.wikipedia.org/wiki/Double-precision_floating-point_format), 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? [🔗](#stop_condition){.small} {: #stop_condition} 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

: | value -|-

x_{0}

| 0.0020000000000000000

x_{1}

| 0.0026479999999999997

x_{2}

| 0.0029259764480000002

x_{3}

| 0.0029582205931032645

x_{4}

| 0.0029585798380249708

x_{5}

| 0.0029585798816568042

x_{6}

| 0.0029585798816568051

x_{7}

| 0.0029585798816568042

x_{8}

| 0.0029585798816568051 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} < ...

. (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} ​​​\geq 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 [🔗](#x0){.small} {: #x0} 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

x = 1

and calculates the product

x a

. If that product is close enough to

1

, then

x

can be used as is. If not,

x

can be either increased or decreased as pertinent, and

x a

calculated anew. This pedestrian approach is used in the following sample implementation: ## Sample implementation [🔗](#sample_impl){.small} {: #sample_impl} Here is a sample implementation for this algorithm in Python: ``` 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: what | frequency -----|------ same value | 67.8 % next lower value | 23.2 % next larger value | 9.0 % For a further experiment, the posh implementation of calculating the return value was tried: ``` python 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: what | frequency -----|------ same value | 68.1 % next lower value | 22.9 % next larger value | 9.0 % ### Suggested exercises [🔗](#codeit){.small} {: #codeit} - 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){.small} {: #theorem} ### Theorem's assumptions [🔗](#assumptions){.small} {: #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 (z) = z

. Assume further

f^{'} (z) = 0

, and also

0 \leq f^{'} (x)

for all

x

such that

u <x\leqz,
and finally, assume that no fixpoint offexists
larger thanuand smaller thanz.

### Theorem's claims  [🔗](#claims){.small} {: #claims}

- Then, anyxsuch thatu < x < zfulfillsx < f (x) \leq z.

Assuming that first claim to be true, given
anyx_{0}such
thatu < x_{0} \leq z,
we can define a sequence recursively viax_{n + 1} = f (x_{n})for all indices (natural numbers)n.

This sequence has the following properties:

-u < x_{n} \leq zwill hold for any indexn.
  (This helps the sequence definition to make sense.)
- This is an increasing sequence, that
  is,x_{n} \leq x_{n + 1}for any indexn.
- This sequence converges toz.
- Finally, it does so in quadratic fashion,
that is, there exist some numbermand some
constantk > 0such
that| x_{n + 1} - z | \leq k {| x_{n} - z |}^{2}for any indexnsuch thatn > m.
{: #quadratic}

## Application [🔗](#inv_from_theorem){.small} {: #inv_from_theorem}

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 (x) = x (2 - a x)is a quadratic parable with zeros at0and\frac{2}{a}, opening downwards.
*z = \frac{1}{a}is clearly a fixpoint
  off.
* A straightforward calculation shows that the only other fixpoint
  offis0.
* The maximal value offis in the middle between its zeros,
  so atz.
* The first derivative offatzis0, asfhas its maximum there.
* The first derivative offhas its _only_ zero
  atz.
* Asfopens downwards,
  the derivative has to be positive beforezand negative thereafter.
* The theorem gives us that our sequence is increasing and we have quadratic convergence, for any start
  valuex_{0}such that0 < x_{0} \leq \frac{1}{a}.
* Nowfis symmetrical around\frac{1}{a},
so we also have quadratic convergence for any
start valuex_{0}such that\frac{1}{a} \leq x_{0} < \frac{2}{a}.

## Suggested exercise [🔗](#check_inequality){.small} {: #check_inequality}

Convince yourself (in a "pedestrian" way) that for anyxwith0 < x < \frac{1}{a}the assertionx < f (x) \leq zholds in our particular\frac{1}{a}situation.

## Suggested exercise [🔗](#pi){.small} {: #pi}

Convince yourself that the theorem allows the use of the
functionf (x) = x + \sin (x)to calculateπ.

(There are much more convenient ways to calculateπ.
But this way is at least legitimate, as the\sinfunction
can be calculated via its [power series](https://en.wikipedia.org/wiki/Sine_and_cosine#Series_and_polynomials)
without prior knowledge ofπ.)

## Suggested exercise [🔗](#temperament){.small} {: #temperament}

In music, the "[equal
temperament](https://en.wikipedia.org/wiki/Equal_temperament)" demands
that the quotient of the frequencies of two adjacent (half-)tones
isz = \sqrt[12]{2}.

* Convince yourself that, for any realλ, this
numberzis a fixpoint off_{λ} (x) = (1 - λ) x + λ (x^{13} - x).
* Find a valueλthat might be suitable to make the theorem
applicable.
* Verify that the theorem indeed is applicable.
* Calculatezviaf_{λ}.

## Proof of the theorem [🔗](#proof){.small} {: #proof}

We first need to prove: Givenxsuch thatu < x < z,
it is true thatx < f (x) \leq z.

Forf (x) \leq z:
According to the [mean value theorem](https://en.wikipedia.org/wiki/Mean_value_theorem), there
must exist someξwithx < ξ < zsuch thatf^{'} (ξ) = \frac{f (z) - f (x)}{z - x}.  According to our assumptions,f^{'} (ξ) \geq 0.
Asx < z,
we seef (x) \leq f (z) = z.


Forx < f (x), consider the
functionp (x) = f (x) - x.
This function is zero precisely wherefhas a fixpoint.
In particular,p (z) = 0.
Nowfhas no fixpoint betweenuandz(exclusively),
sopbecomes zero nowhere in that interval.  Being a continuous function,pcannot
change sign over an interval without assuming a zero value in that interval.
Sopcannot change sign betweenuandz(exclusively), in other words,
must either be always negative or always positive there.
Nowpis continuously differentiable.
Its derivative evaluates to-1atz.
As the derivative is continuous, it must be negative for some interval of points
slightly smaller thanz.
By another application of the mean value theorem,pmust be positive
on that same interval.  We establishedpdoes not change sign
in the interval strictly betweenuandz,
and we have found points in that interval
wherepis positive, so we now know
it must be positive everywhere betweenuandz(exclusively), which in particular
provesp (x) > 0.
Nowp (x) > 0is equivalent
tox < f (x).

These two together prove that the sequencex_{0}, x_{1}, x_{2},...increases and is bounded from above byz.

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ζ = \lim_{n \to \infty} x_{n}.
Asfis continuous:f (ζ) = f (\lim_{n \to \infty} x_{n}) = \lim_{n \to \infty} f (x_{n}) = \lim_{n \to \infty} x_{n + 1} = \lim_{n \to \infty} x_{n} = ζSoζis a fixpoint off, which
forcesζ = z.  So the sequencex_{0}, x_{1}, x_{2}, ...converges with limitz.

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

This is a direct consequence of
[Taylor's theorem](https://en.wikipedia.org/wiki/Taylor%27s_theorem).
It allows us to rewritef (x_{n}) = f (z) + f^{'} (z) (x_{n} - z) + R_{1} (x_{n})In our case,f^{'} (z) = 0,
so this reduces tof (x_{n}) = f (z) + R_{1} (x_{n})orf (x_{n}) - f (z) = R_{1} (x_{n})

x_{n + 1} - z = R_{1} (x_{n})

. The Lagrange form of the remainder boils down to

R_{1} (x_{n}) = \frac{f^{″} (ξ)}{2} {(x_{n} - z)}^{2}

for some suitable value

ξ

with

x_{n} \leq ξ \leq z

. Now

f^{″}

is continuous, hence bounded over any closed interval containing

z

. If we have an upper bound

k

| f |

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 | ≤ \frac{k}{2} {| x_{n} - z |}^{2}

, which completes the proof. ## Suggested exercise [🔗](#mirrortheorem){.small} {: #mirrortheorem} 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](https://en.wikipedia.org/wiki/Square_root_algorithms#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}

. {: #herons} Use your theorem to prove that recursion via the function

f (x) = x + \sin (x)

converges to

π

for any start value

x_{0}

with

π \leq x_{0} < 2 π

. ## Suggested exercise [🔗](#cubic){.small} {: #cubic} Above, there is a [definition of "quadratic convergence"](#quadratic) 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 π

, the sequence defined by

x_{n + 1} = x_{n} + \sin x_{n}

exhibits cubic convergence behaviour, converging towards

π

. * Show that, given any

a > 0

and any start value

x_{0}

with

0 < x_{0} < \sqrt{\frac{7}{3} a}

, the sequence defined by

x_{n + 1} = (((\frac{3}{8 a^{2}} x_{n}^{2}) - \frac{5}{4 a}) 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](#herons). [🔗](#fast_sqrt){.small} {: #fast_sqrt} ## Author's remark [🔗](#mathmlpain){.small} {: #mathmlpain} 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 [🔗](#comments){.small} {: #comments} If you want to comment or discuss this piece and have a Fediverse account, feel invited to answer [my pertinent toot](https://mastodon.radio/@dj3ei/116618396589431125).