| There is a book by Cody on calculating the elementary functions;
unfortunately I don't remember the exact title. IBM have published a
manual describing the algorithms that they use in their new math
library (comes with VS Fortran, I believe release 2 or later; the
manual is called, if memory serves, "VS Fortran Language and Library
Reference"). The algorithms used by IBM (but not the constants, of
course) should serve very well on a VAX, because IBM use hex
normalization and so can lose up to three bits of accuracy in a single
floating-point operation; thus, they have to jump through hoops to get
good math routines.
The normal technique for calculating any logarithm (they only differ by
a constant factor, anyway) is as follows:
1. Reduce the argument by repeated divisions (don't do divisions,
subtract from the exponent! Do "conceptual divisions"!) to some
standard interval (normally, 1 to 2). Say that you need to make one
division by A to get the argument between 1 and 2.
2. Use table lookup to find a sub-interval of [1,2) that contains the
(reduced) argument. Say that one of the endpoints of that subinterval
is B.
3. Reduce the argument further by dividing it by B, giving C. Now use the
appropriate (fixed!) number of terms from the series expansion
ln(1+x) = x - x**2/2 + x**3/3 - x**4/4 + ...
to calculate ln(C).
4. Your result is ln(A)+ln(B)+ln(C). If you wanted logarithms other
than base e (=2.718281828...), either multiply your result by an
appropriate scaling factor (in the case of common logarithms, it is
ln(10) or approximately 2.3026) or modify the constants in the series
expansion used to calculate ln(C); in either case, you must perform one
additional multiplication.
Note: sometimes, it is advisable to use a polynomial that is not a
straight Taylor-series expansion. The reason for that is that the
error term in the Taylor series is a *lot* bigger at the end points of
the interval used than it is in the middle (after all, x**n grows
fairly rapidly with x). Popular choices are Chebychev (sp? name is
transliterated variously from the Cyrillic alphabet) polynomials and
minimax polynomials. Chebychev polynomials have the *big* disadvantage
that their coefficients are large, so the calculation of higher-degree
Chebychev polynomials is fraught with numerical problems. Minimax
polynomials are polynomials calculated to minimize ABS(function-polynomial)
over the given interval, so they are -- in a sense -- the best one
can use. There are published algorithms for calculating minimax
polynomials; unfortunately I can't give a reference.
|
| Here is what the documentation of the MATH RTL says about the VAX
implementation for LOG (slightly edited):
Computation of the natural logarithm routine is based on the following:
1 ln(X*Y) = ln(X) + ln(Y)
2 ln(1 + X) = X - x^2/2 + X^3/3 - X^2/4 ... for |X| < 1
3 ln(X) = ln(A) + 2(V + V^3/3 + V^5/5 + V^7/7 ...) =
ln(A) + V*p(V^2), where V = (X-A)/(X+A), A>0, and
p(y) = 2*(1+y/3+y^2/5...)
For x = 2^n*f, where n is an integer and f is in the interval of 0.5
to 1 define the following quantities:
if n>=1, then N=n-1 and F=2f
if n<=0, then N=n and F=f
From (1) above it follows that:
4 ln(X) = N*ln(2) + ln(F)
Based on the above relationships LOG is computed as follows:
1 If |F-1|<2^-5, LOG(X) = N*LOG(2) + W + W*p(W), where W=F-1.
2 Otherwise, LOG(X)=N*LOG(2) + LOG(A) + V*p(V^2), where
V=(F-A)/(F+A) and A and LOG(A) are obtained by table look up.
------------------------------------------------------------------
If this is strictly for internal use, I could probably get the MACRO32
code for you.
Topher
|