mirror of
git://sourceware.org/git/newlib-cygwin.git
synced 2025-01-16 03:19:54 +08:00
70317d8506
* libm/mathfp/s_logarithm.c: Fix error introduced by previous
fix when handling negative input values. Make function
consistent with math directory and glibc version such that
inf and nan values return inf and nan respectively with no
errno setting.
* libm/mathfp/sf_logarithm.c: Ditto.
* libm/math/w_log.c: Set errno to ERANGE when input is 0.0.
* libm/math/wf_log.c: Ditto.
* libm/math/w_log10.c: Ditto.
* libm/math/wf_log10.c: Ditto.
149 lines
3.4 KiB
C
149 lines
3.4 KiB
C
|
|
/* @(#)z_logarithm.c 1.0 98/08/13 */
|
|
/******************************************************************
|
|
* The following routines are coded directly from the algorithms
|
|
* and coefficients given in "Software Manual for the Elementary
|
|
* Functions" by William J. Cody, Jr. and William Waite, Prentice
|
|
* Hall, 1980.
|
|
******************************************************************/
|
|
|
|
/*
|
|
FUNCTION
|
|
<<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
|
|
|
|
INDEX
|
|
log
|
|
INDEX
|
|
logf
|
|
INDEX
|
|
log10
|
|
INDEX
|
|
log10f
|
|
|
|
ANSI_SYNOPSIS
|
|
#include <math.h>
|
|
double log(double <[x]>);
|
|
float logf(float <[x]>);
|
|
double log10(double <[x]>);
|
|
float log10f(float <[x]>);
|
|
|
|
TRAD_SYNOPSIS
|
|
#include <math.h>
|
|
double log(<[x]>);
|
|
double <[x]>;
|
|
|
|
float logf(<[x]>);
|
|
float <[x]>;
|
|
|
|
double log10(<[x]>);
|
|
double <[x]>;
|
|
|
|
float log10f(<[x]>);
|
|
float <[x]>;
|
|
|
|
DESCRIPTION
|
|
Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
|
|
(where e is the base of the natural system of logarithms, 2.71828@dots{}) or
|
|
base 10.
|
|
<<log>> and <<logf>> are identical save for the return and argument types.
|
|
<<log10>> and <<log10f>> are identical save for the return and argument types.
|
|
|
|
RETURNS
|
|
Normally, returns the calculated value. When <[x]> is zero, the
|
|
returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
|
|
When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
|
|
<<errno>> is set to <<EDOM>>. You can control the error behavior via
|
|
<<matherr>>.
|
|
|
|
PORTABILITY
|
|
<<log>> is ANSI. <<logf>> is an extension.
|
|
|
|
<<log10>> is ANSI. <<log10f>> is an extension.
|
|
*/
|
|
|
|
|
|
/******************************************************************
|
|
* Logarithm
|
|
*
|
|
* Input:
|
|
* x - floating point value
|
|
* ten - indicates base ten numbers
|
|
*
|
|
* Output:
|
|
* logarithm of x
|
|
*
|
|
* Description:
|
|
* This routine calculates logarithms.
|
|
*
|
|
*****************************************************************/
|
|
|
|
#include "fdlibm.h"
|
|
#include "zmath.h"
|
|
|
|
#ifndef _DOUBLE_IS_32BITS
|
|
|
|
static const double a[] = { -0.64124943423745581147e+02,
|
|
0.16383943563021534222e+02,
|
|
-0.78956112887481257267 };
|
|
static const double b[] = { -0.76949932108494879777e+03,
|
|
0.31203222091924532844e+03,
|
|
-0.35667977739034646171e+02 };
|
|
static const double C1 = 22713.0 / 32768.0;
|
|
static const double C2 = 1.428606820309417232e-06;
|
|
static const double C3 = 0.43429448190325182765;
|
|
|
|
double
|
|
_DEFUN (logarithm, (double, int),
|
|
double x _AND
|
|
int ten)
|
|
{
|
|
int N;
|
|
double f, w, z;
|
|
|
|
/* Check for range and domain errors here. */
|
|
if (x == 0.0)
|
|
{
|
|
errno = ERANGE;
|
|
return (-z_infinity.d);
|
|
}
|
|
else if (x < 0.0)
|
|
{
|
|
errno = EDOM;
|
|
return (z_notanum.d);
|
|
}
|
|
else if (!isfinite(x))
|
|
{
|
|
if (isnan(x))
|
|
return (z_notanum.d);
|
|
else
|
|
return (z_infinity.d);
|
|
}
|
|
|
|
/* Get the exponent and mantissa where x = f * 2^N. */
|
|
f = frexp (x, &N);
|
|
|
|
z = f - 0.5;
|
|
|
|
if (f > __SQRT_HALF)
|
|
z = (z - 0.5) / (f * 0.5 + 0.5);
|
|
else
|
|
{
|
|
N--;
|
|
z /= (z * 0.5 + 0.5);
|
|
}
|
|
w = z * z;
|
|
|
|
/* Use Newton's method with 4 terms. */
|
|
z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
|
|
|
|
if (N != 0)
|
|
z = (N * C2 + z) + N * C1;
|
|
|
|
if (ten)
|
|
z *= C3;
|
|
|
|
return (z);
|
|
}
|
|
|
|
#endif /* _DOUBLE_IS_32BITS */
|