EDIT: See the necro-update below-- it provides a much more efficient algorithm, even 1 multiply better than most suggested implementations of Newton-Raphson division
In the study of numerical algorithms, quadratic convergence refers to an iterative algorithm that approximates a function and each iteration doubles the digits of accuracy. In this post, I will expose an algorithm for division that can double, triple, quadruple, etc. the number of digits of accuracy at each iteration. I will then show that the optimal algorithm is that which offers cubic convergence (digits of accuracy tripled).
What does it mean to multiply the number of correct digits by a?Suppose that we have a sequence
##x_{0},x_{1},x_{2},x_{3},...## and suppose that
##x_{n}\rightarrow c## Then by the definition of convergence,
##\lim\limits_{n\rightarrow\infty}{|x_{n}-c|}=0##. You should always think of |x-y| as a "distance" between two points, and in this case, that distance is the error of our approximations,
##x_{n}##.
If the number of correct digits is multiplied by
a, then that means
##\log(|x_{n+1}-c|)=a\cdot\log|x_{n}-c|##, or rewritten using logarithm tricks,
##|x_{n+1}-c|=|x_{n}-c|^{a}##. This assumes that the error is initially less than 1.
The AlgorithmThe algorithm is very simple, but to break it down we need to look at "obvious" bits of math. Firstly,
##\frac{N}{D}=\frac{N\cdot c}{D\cdot c}##, and secondly, all real numbers can be written in the form of
##x=y\cdot 2^{m}, y\in(.5,1]##. You may have said "duh" at the first statement, but if you need convincing of the second, it basically says that if you divide or multiply a number by 2 enough times, then it will eventually get between 1/2 and 1. Or in another way, shown by example: 2
11<3767<2
12, so dividing it all by 2
12, .5<3767/4096<1. All real numbers can be bounded by consecutive power of 2, therefore they can all be written in the above form.
The final piece of the puzzle is to recognize that
##\frac{N}{1}={N}##. What we will do is recursively multiply D by some constant
c at each iteration in order to drive it closer to 1. If we use range reduction techniques to get D on (.5,1] (multiply both N and D by the power of 2 that gets D in that range). Then what we want is to choose
c so that
##|1-D_{n}\cdot c|= |1-D_{n}|^{a}##. If
##0\leq D_{n}\leq 1##, then we have
##1-D_{n}\cdot c = (1-D_{n})^{a}##. When
a is a natural number greater than 1 (if it is 1, then there is no convergence
), we have the following:
##1-D_{n}\cdot c = \sum\limits_{k=0}^{a}{{a \choose k}(-D_{n})^{k}}####D_{n}\cdot c = 1-\sum\limits_{k=0}^{a}{{a \choose k}(-D_{n})^{k}}####D_{n}\cdot c = 1-(1+\sum\limits_{k=1}^{a}{{a \choose k}(-D_{n})^{k}})####D_{n}\cdot c = -\sum\limits_{k=1}^{a}{{a \choose k}(-D_{n})^{k}}####D_{n}\cdot c = {a \choose 1}D_{n}-{a \choose 2}D_{n}^{2}+{a \choose 3}D_{n}^{3}+\cdots####c = {a \choose 1}-{a \choose 2}D_{n}+{a \choose 3}D_{n}^{2}+\cdots##Using a technique similar to Horner's Method, we can obtain:
##c = {a \choose 1}-D_{n}({a \choose 2}-D_{n}({a \choose 3}-D_{n}({a \choose 4}-D_{n}\cdots##Then the selection for
c requires
a-2 multiplications and
a-1 subtractions. If we take
##N_{n+1}=N_{n}\cdot c## and
##D_{n+1}=D_{n}\cdot c##, then the number of digits of accuracy is multiplied by
a using a total of
a multiplications and
a-1.
ExampleLets say we want quadratic convergence. That is, a=2. Then
##c_{n+1}=2-D_{n}## That's right-- all you have to do is multiply by 2 minus the denominator at each iteration and if you started with 1 digit of accuracy, you next get 2 digits, then 4, then 8, 16, 32,.... Let's perform 49/39. Divide numerator and denominator by 2
6=64 and we get .765625/.609375.
First iteration:
c=2-D=1.390625
Nc=1.064697266
Dc=.8474121094
Second iteration:
c=1.152587891
Nc=1.227157176
Dc=.9767169356
next iterations in the form {c,N,D}:
{1.023283064,1.255729155,.9994578989}
{1.000542101,1.256409887,.9999997061}
{1.000000294,1.256410256,.9999999999}
{1.000000294,1.256410256,1.000000000}
And in fact, 49/39 is 1.256410256 according to my calc.
Can we find an optimal choice for a?The only variable over which we have control is
a-- how quickly the algorithm converges. Therefore we should find an optimal choice of
a. In this case, "optimal" means getting the most accuracy for the least work (computation time).
First,
##1-D_{0}<.5=2^{-1}## So then after n iterations, we have at least a
n bits of accuracy and since
##D_{0}## may get arbitrarily close to .5, this is the best lower bound on the number of bits of accuracy achieved for n iterations of a randomly chosen
D. As well, at
n iterations, we need to use
##a\cdot n## multiplications and the subtractions are treated as trivial. Then in order to get
m bits of accuracy, we need roughly
##a\log_{a}(m)## multiplications. We want to minimize this function of
a for
a is an integer greater than 1. To do that, we find when
##0=\frac{d}{da}(a\log_{a}(m))##:
##0=\frac{d}{da}(a\log_{a}(m))####0=\frac{d}{da}(a\frac{\log(m)}{\log(a)})####0=\frac{d}{da}(\frac{a}{\log(a)})####0=\frac{1}{\log(a)}-\frac{1}{\log^{2}(a)}##Since a>1:
##0=1-\frac{1}{\log(a)}####1=\frac{1}{\log(a)}####\log(a)=1####a=e##However,
e is not an integer, but it is bounded by 2 and 3, so check which is smaller: 2/log(2) or 3/log(3) and we find that a=3 is the optimal value that achieves the most accuracy for the least work.
Example:Our algorithm body is:
c=3-D*(3-D)
N=N*c
D=D*c
Then 49/39=.765625/.609375 and {c,N,D}=:
{1.543212891,1.181522369,.9403953552}
{1.063157358,1.256144201,.9997882418}
{1.000211803,1.256410256,1.000000000}
ConclusionIn reality, their is a specific M so that for all m>M, a=3 always requires less computational power than a=2. Before that cutoff, there are values for which they require the same amount of work or even the cost is in favor of a=2. For example, 32768 bits of accuracy can be achieved with the same amount of work. The cutoff appears to be m=2^24 which is pretty gosh darn huge (16777216 bits of accuracy is over 5 million digits).
