panon@cheddar.cc.ubc.ca (Paul-Andre Panon for Jonathan Thornburg) (05/03/90)
I'm posting this for Jonathan Thornburg in our Astronomy dept. who seems
to have a read-only news connection.
---------------------------------
There's a bug in log(3) in Microport's System V/AT 2.3.0-U. This is on an
AT clone with no hardware floating point. The problem is that, for arguments
between 0.5 and 1.0, the answers are wrong. In particular, log(0.8) gives a
*positive* result!.
The (my) fix is to rewrite log(3). I'm not brave enough to change libm, so I
put my new log(3) in a new library, "libfixm"; I do "-lfixm -lm" whenever I
used to do "-lm". (Yes, this is a kludge...)
Here's the code for a new log(3) I cooked up:
---------- file "log.c" ----------
/* log - my log(), log10() functions for IEEE double precision floating pt.
*/
/*
** <<<useful constants>>>
* log10 - base 10 logarithm function (calls log())
* log - natural logarithm function
*/
#include <math.h>
#define then
#define EDOM DOMAIN /* from <math.h> */
/***************************************************************************
**/
/*
* This software written by
* Jonathan Thornburg
* Dept of Geophysics & Astronomy
* The University of British Columbia thornburg@mtsg.ubc.ca
* Vancouver BC V6T 1W5 userbkis@ubcmtsg.bitnet
* Canada ...!ubc-vision!ubcmtsg.bitnet!thornburg
* on 5 April 1990. It's in the public domain.
*
* The constants and coefficients are from
* John F. Hart, E. W. Cheney, Charles L. Lawson, Hans J. Maehly,
* Charles K. Mesztenyi, John R. Rice, Henry G. Thacher, Jr.,
* and Christoph Witzgall,
* "Computer Approximations",
* Krieger (1978)
*/
/* useful constants, from appendix C */
#define LOG10_OF_E 0.43429448190325182765
#define SQRT_OF_HALF 0.70710678118654752440
#define LOGE_OF_2 0.69314718055994530942
/***************************************************************************
**/
/*
* This is the base 10 logarithm function.
*/
public double log10(x)
double x;
{
public double log();
return( LOG10_OF_E * log(x) );
}
/***************************************************************************
**/
/*
* This is the natural logarithm function.
*
* The algorithm used is:
* 1. Check for domain error.
* 2. Range reduce to [1/2, 1).
* 3. Range adjust to [1/sqrt(2), sqrt(2))
* 4. Compute rational approximation on this range (#2704 in Hart et al).
* 5. Compute result.
*
* Bugs:
* - The System V matherr() error handling system is *not* supported. Instead,
* we return -HUGE (defined in <math.h> and set errno to EDOM (defined above)
* on a domain error (argument <= 0).
* - This function was put together in a hurry (when I found that my system's
* log() returns garbage for some arguments), with only minimal attention to
* minimizing cancellation and roundoff error. My rather limited tests of it
* have found no errors larger than 2 units in the least significant bit, but
* I don't guarantee this.
* - For the same reason, this function is somewhat slower (perhaps up to a
* factor of 2 on machines with software floating point) than necessary. It
* uses floating point throughout, where a "tuned" implementation could use
* fixed point. In stage 3, a more extensive range reduction might well be
* faster. In stage 4, an "economized" (eg Pan form, c.f. Hart or Knuth v.2)
* polynomial evaluation scheme might be faster. In stage 5, it does an
* integer to floating point conversion and a floating point multiplication
* where a "tuned" implementation could do a table lookup for most arguments.
*/
public double log(x)
double x;
{
double f, z, z2;
int n;
double P_of_z2, Q_of_z2, log_of_f;
double log_of_x;
static double P[] = {
-0.90174691662040536329e+02,
0.93463900642858538247e+02,
-0.18327870372215593212e+02
};
static double Q[] = {
-0.45087345831020305748e+02,
0.61761065598471302843e+02,
-0.20733487895513939345e+02,
1.0 /* not actually used */
};
/*
* Stage 1: Check for domain error
*/
if (x <= 0.0)
then {
errno = EDOM;
return(- HUGE);
}
/*
* Stage 2: Range reduce to [1/2, 1).
*
* We compute (f,n) such that x = f * 2**n.
*/
f = frexp(x, &n);
/*
* Stage 3: Range adjust to [1/sqrt(2), sqrt(2)).
*
* We maintain the invariant that x = f * 2**n.
*
* As mentioned under "Bugs:" in the function comments, it might be worthwhile
* doing a extensive range reduction, perhaps to [1/(2**.25), 2**.25). This
* would allow a lower-degree rational approximation to be used in stage 4.
*/
if (f < SQRT_OF_HALF)
then {
f = ldexp(f, 1);
--n;
}
/*
* Stage 4: Compute rational approximation to log(f) on [1/sqrt(2), sqrt(2)).
*
* The approximation is of the form
* log(f) = z * P(z^2) / Q(z^2)
* It's #2704 in the Hart et al reference given above (note the typo in the
* heading of the index listing).
*
* As mentioned under "Bugs:" in the function comments, it might be worthwhile
* using an "economized" polynomial evaluation technique (eg Pan form, c.f.
* Hart or Knuth v.2) here. Another possibility is to use the J-fraction form
* for the rational approximation.
*/
z = (f - 1.0) / (f + 1.0);
z2 = z * z;
P_of_z2 = P[0] + z2 * (P[1] + z2 * P[2]);
Q_of_z2 = Q[0] + z2 * (Q[1] + z2 * (Q[2] + z2 /* * Q[3] */));
log_of_f = z * P_of_z2 / Q_of_z2;
/*
* Stage 5: Compute result log(x).
*
* As mentioned under "Bugs:" in the function comments, we could replace the
* integer to floating point conversion and floating point multiply
* LOGE_OF_2 * ((double) n)
* by a table lookup of values of this expression. Since IEEE double precision
* floating point format has an exponent range of about +/- 1024, a table with
* entries for all possible exponents would be quite large (2K entries @ 8
* bytes = 16K bytes). Thus the best thing would probably be to have the table
* only cover (say) the IEEE single precison floating point exponent range of
* about +/- 128 (==> 1K bytes), and do the integer to floating point
* conversion and floating point multiply for the remaining (very large or very
* small) exponents.
*/
log_of_x = LOGE_OF_2 * ((double) n) + log_of_f;
return(log_of_x);
}
---------- end of file "log.c" ----------
- Jonathan Thornburg
Dept of Geophysics & Astronomy
The University of British Columbia thornburg@mtsg.ubc.ca
Vancouver BC V6T 1W5 userbkis@ubcmtsg.bitnet
Canada ...!ubc-vision!ubcmtsg.bitnet!thornburg
--
Paul-Andre_Panon@staff.ucs.ubc.ca or USERPAP1@UBCMTSG
or Paul-Andre_Panon@undergrad.cs.ubc.ca or USERPAP1@mtsg.ubc.ca
Looking for a .signature? "We've already got one. It is ver-ry ni-sce!"