🚀 go-pugleaf

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!"