137 lines
2.9 KiB
C
137 lines
2.9 KiB
C
|
|
/*++
|
|
|
|
Copyright (c) 1999 Microsoft Corporation
|
|
|
|
Module Name:
|
|
|
|
log.c
|
|
|
|
Abstract:
|
|
|
|
logarithmic functions
|
|
|
|
Author:
|
|
|
|
|
|
|
|
Revision History:
|
|
|
|
29-sept-1999 ATM Shafiqul Khalid [askhalid] copied from rtl library.
|
|
--*/
|
|
|
|
|
|
#include <math.h>
|
|
#include <trans.h>
|
|
|
|
static double _log_hlp( double x, int flag);
|
|
|
|
/* constants */
|
|
static double const c0 = 0.70710678118654752440; /* sqrt(0.5) */
|
|
static double const c1 = 0.69335937500000000000;
|
|
static double const c2 = -2.121944400546905827679e-4;
|
|
static double const c3 = 0.43429448190325182765;
|
|
|
|
/* coefficients for rational approximation */
|
|
static double const a0 = -0.64124943423745581147e2 ;
|
|
static double const a1 = 0.16383943563021534222e2 ;
|
|
static double const a2 = -0.78956112887491257267e0 ;
|
|
static double const b0 = -0.76949932108494879777e3 ;
|
|
static double const b1 = 0.31203222091924532844e3 ;
|
|
static double const b2 = -0.35667977739034646171e2 ;
|
|
/* b3=1.0 is not used -avoid multiplication by 1.0 */
|
|
|
|
#define A(w) (((w) * a2 + a1) * (w) + a0)
|
|
#define B(w) ((((w) + b2) * (w) + b1) * (w) + b0)
|
|
|
|
|
|
/***
|
|
*double log(double x) - natural logarithm
|
|
*double log10(double x) - base-10 logarithm
|
|
*
|
|
*Purpose:
|
|
* Compute the natural and base-10 logarithm of a number.
|
|
* The algorithm (reduction / rational approximation) is
|
|
* taken from Cody & Waite.
|
|
*
|
|
*Entry:
|
|
*
|
|
*Exit:
|
|
*
|
|
*Exceptions:
|
|
* I P Z
|
|
*******************************************************************************/
|
|
|
|
double Proxylog10(double x)
|
|
{
|
|
return(_log_hlp(x,OP_LOG10));
|
|
}
|
|
|
|
double Proxylog(double x)
|
|
{
|
|
return(_log_hlp(x,OP_LOG));
|
|
}
|
|
|
|
static double _log_hlp(double x, int opcode)
|
|
{
|
|
unsigned int savedcw;
|
|
int n;
|
|
double f,result;
|
|
double z,w,znum,zden;
|
|
double rz,rzsq;
|
|
|
|
/* save user fp control word */
|
|
savedcw = _maskfp();
|
|
|
|
/* check for infinity or NAN */
|
|
if (IS_D_SPECIAL(x)){
|
|
switch (_sptype(x)) {
|
|
case T_PINF:
|
|
RETURN(savedcw, x);
|
|
case T_QNAN:
|
|
return _handle_qnan1(opcode, x, savedcw);
|
|
case T_SNAN:
|
|
return _except1(FP_I, opcode, x, _s2qnan(x), savedcw);
|
|
}
|
|
/* NINF will be handled in the x<0 case */
|
|
}
|
|
|
|
if (x <= 0.0) {
|
|
double qnan;
|
|
if (x == 0.0) {
|
|
return _except1(FP_Z,opcode,x,-D_INF,savedcw);
|
|
}
|
|
qnan = (opcode == OP_LOG ? QNAN_LOG : QNAN_LOG10);
|
|
return _except1(FP_I,opcode,x,qnan,savedcw);
|
|
}
|
|
|
|
if (x == 1.0) {
|
|
// no precision ecxeption
|
|
RETURN(savedcw, 0.0);
|
|
}
|
|
|
|
f = _decomp(x, &n);
|
|
|
|
if (f > c0) {
|
|
znum = (f - 0.5) - 0.5;
|
|
zden = f * 0.5 + 0.5;
|
|
}
|
|
else {
|
|
n--;
|
|
znum = f - 0.5;
|
|
zden = znum * 0.5 + 0.5;
|
|
}
|
|
z = znum / zden;
|
|
w = z * z;
|
|
|
|
rzsq = w * A(w)/B(w) ;
|
|
rz = z + z*rzsq;
|
|
|
|
result = (n * c2 + rz) + n * c1;
|
|
if (opcode == OP_LOG10) {
|
|
result *= c3;
|
|
}
|
|
|
|
RETURN_INEXACT1(opcode,x,result,savedcw);
|
|
}
|