Previously, I mentioned some methods to build a table of logarithms based off of adding or subtracting fractional exponents to attain an arbitrary fixed-point decimal power of a base. In general, this is a fast but possibly inaccurate method for computing arbitrary logarithms and exponents.
Here, I present methods that may have better numerical stability on accuracy, are easy to understand how they are fully general, but may come with the extreme expense of being computationally expensive.
So, first of all, computing nth roots using Newton’s method. It’s rather easy once you setup the right equations.
To determine x = sqrt(n):
n = x^2
0 = x^2 - n
y = x^2 - n
f(x) = x^2 - n
To determine x = n^(1/r):
f(x) = x^r - n
Once you have the principle function/equation setup, you can use
Newton’s root-finding method on it to find the zeros, or roots, of the
function. That is, the values of x that will result in y = 0. In
this case, it can clearly be seen that the zeros of these functions
correspond to the square root or rth root of the number n. One
iteration of Newton’s method is given as follows.
y_(a+1) = y_a - f(x) / f'(x)
The power rule can easily be used to compute the derivative of f(x).
f'(x) = r*x^(r - 1)
To use Newton’s method, first you start out with a good starting guess, it must be reasonably good or else the solution may not converge on the correct value, then you run the iteration multiple times to get successively better approximations of the root. Once you’ve computed to reasonable accuracy, you can stop.
Now, for computing e^x. This is a rather simple application of the
Maclaurin series that results in a power series quite similar to the
factorial series used to compute number e itself.
e^x = sum_{n = 0}^{oo} x^n/n!
Happily, this computational form allows you to easily raise number e
to a non-integer power. The individual terms of the power series
itself all rase x to an integer power, and raising a non-integer to
an integer power is trivial by simply multiplying it by itself an
integer number of times.
We can use Newton’s method to compute natural logarithms with ease.
x = ln(n)
f(x) = e^x - n
f'(x) = e^x
The reason, of course, why this is still easy to compute is because we
already know how to compute e^x. So we can use that to compute the
derivative that is required to be known while using Newton’s method.
In regard to a good starting estimate, we can determine the
approximate base-2 logarithm by counting the number of binary bits of
precision. Base 2 is, after all, approximately equal to a number e
base.
Alternatively, we could also use a Maclaurin series expansion since
the derivative of ln |x| is 1/x. Let’s show some derivatives for
example.
f(a) = ln |a|
f'(a) = 1/a = a^(-1)
f''(a) = -a^(-2)
f'''(a) = 2*a^(-3)
f^(4)(a) = -6*a^(-4)
f^(5)(a) = -24*a^(-5)
f^(6)(a) = 120*a^(-6)
n >= 1, f^(n)(a) = -1^(n - 1) * (n - 1)! * a^(-n)
Now, in a power series, we’re also dividing by a factorial, so we can further simplify this as the factors of the factorials cancel.
-1^(n - 1) * (n - 1)! * a^(-n) * (x - a)^n / n!
= -1^(n - 1) * a^(-n) * (x - a)^n / n
Now unfortunately, we can’t turn this into a Maclaurin series due to
the trouble of evaluating ln(0), but we can get reasonably simple
behavior if we use a = 1.
Finally, to compute arbitrary logarithms and exponents, we can use the logarithm rules to rephrase the equations in terms of number e as the base.
a^b = e^(ln(a^b)) = e^(b * ln(a))
log_a(b) = ln(b) / ln(a)