newlib-cygwin/newlib/libm/mathfp/s_logarithm.c

137 lines
3.2 KiB
C
Raw Normal View History

2000-02-18 03:39:52 +08:00
/* @(#)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
2005-10-28 Bob Wilson <bob.wilson@acm.org> * libc/misc/unctrl.c: Replace FUNCTION description. * libc/signal/signal.c: Remove documentation for raise and _raise_r. * libc/stdio/getdelim.c: Fix spelling errors. * libc/stdio/getw.c: Put RETURNS on a separate line. Fix punctuation. * libc/stdio/putw.c: Likewise. * libc/stdlib/a64l.c: Fix formatting, spelling and punctuation in documentation. * libc/stdlib/assert.c: Do not capitalize FUNCTION description. * libc/stdlib/efgcvt.c: Add spaces to FUNCTION description. * libc/stdlib/envlock.c: Use em-dash in FUNCTION description. * libc/stdlib/mlock.c: Likewise. * libc/stdlib/mstats.c: Likewise. * libc/time/tzlock.c: Likewise. * libc/stdlib/rand.c: Use "multi-threaded" and "thread-safe" in NOTES. * libc/stdlib/rand48.c: Remove extra space in FUNCTION description and hyphenate "pseudo-random". * libc/string/bcmp.c: Remove extra blank lines in documentation. * libc/string/strncat.c: Likewise. * libc/string/memchr.c: Remove extra ">" character in documentation. * libc/string/strcspn.c: Use "characters" instead of "chars". * libc/string/strpbrk.c: Likewise. * libc/string/strerror_r.c: Capitalize "GNU". * libc/string/strnlen.c: Likewise. * libc/string/strtok.c: Fix formatting, spelling and punctuation in documentation. Use "multi-threaded" and "thread-safe" in NOTES. * libc/string/wcscat.c: Split PORTABILITY into two paragraphs. * libc/string/wcschr.c: Likewise. * libc/string/wcscmp.c: Likewise. * libc/string/wcscpy.c: Likewise. * libc/string/wcscspn.c: Likewise. * libc/string/wcslen.c: Likewise. * libc/string/wcsncat.c: Likewise. * libc/string/wcsncmp.c: Likewise. * libc/string/wcsncpy.c: Likewise. * libc/string/wcsnlen.c: Likewise. * libc/string/wcspbrk.c: Likewise. * libc/string/wcsrchr.c: Likewise. * libc/string/wcsspn.c: Likewise. * libc/string/wmemchr.c: Likewise. * libc/string/wmemcmp.c: Likewise. * libc/string/wmemcpy.c: Likewise. * libc/string/wmemset.c: Likewise. * libc/string/wmemmove.c: Likewise. Also fix FUNCTION description. * libc/string/wcswidth.c: Formatting and punctuation in documentation. * libc/string/wcwidth.c: Likewise. * libm/common/s_modf.c: Remove extra period from documentation. * libm/math/s_isnan.c: Fix formatting, grammar and punctuation in documentation. * libm/mathfp/s_isnan.c: Likewise. * libm/math/s_ldexp.c: Fix punctuation. * libm/mathfp/s_ldexp.c: Likewise. * libm/math/w_log.c: Likewise. * libm/mathfp/s_logarithm.c: Likewise. * libm/math/w_j0.c: Add spaces to FUNCTION description. * libm/mathfp/w_jn.c: Likewise.
2005-10-29 05:21:08 +08:00
<<log>> is ANSI. <<logf>> is an extension.
<<log10>> is ANSI. <<log10f>> is an extension.
2000-02-18 03:39:52 +08:00
*/
/******************************************************************
* 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 domain error here. */
if (x <= 0.0)
{
errno = ERANGE;
return (z_notanum.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 */