source: trunk/libs/newlib/src/newlib/libm/mathfp/s_logarithm.c @ 503

Last change on this file since 503 was 444, checked in by satin@…, 6 years ago

add newlib,libalmos-mkh, restructure shared_syscalls.h and mini-libc
File size: 3.2 KB
Line 
1
2/* @(#)z_logarithm.c 1.0 98/08/13 */
3/******************************************************************
4 * The following routines are coded directly from the algorithms
5 * and coefficients given in "Software Manual for the Elementary
6 * Functions" by William J. Cody, Jr. and William Waite, Prentice
7 * Hall, 1980.
8 ******************************************************************/
9
10/*
11FUNCTION
12       <<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
13
14INDEX
15    log
16INDEX
17    logf
18INDEX
19    log10
20INDEX
21    log10f
22
23SYNOPSIS
24       #include <math.h>
25       double log(double <[x]>);
26       float logf(float <[x]>);
27       double log10(double <[x]>);
28       float log10f(float <[x]>);
29
30DESCRIPTION
31Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
32(where e is the base of the natural system of logarithms, 2.71828@dots{}) or
33base 10.
34<<log>> and <<logf>> are identical save for the return and argument types.
35<<log10>> and <<log10f>> are identical save for the return and argument types.
36
37RETURNS
38Normally, returns the calculated value.  When <[x]> is zero, the
39returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
40When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
41<<errno>> is set to <<EDOM>>.  You can control the error behavior via
42<<matherr>>.
43
44PORTABILITY
45<<log>> is ANSI. <<logf>> is an extension.
46
47<<log10>> is ANSI. <<log10f>> is an extension.
48*/
49
50
51/******************************************************************
52 * Logarithm
53 *
54 * Input:
55 *   x - floating point value
56 *   ten - indicates base ten numbers
57 *
58 * Output:
59 *   logarithm of x
60 *
61 * Description:
62 *   This routine calculates logarithms.
63 *
64 *****************************************************************/
65
66#include "fdlibm.h"
67#include "zmath.h"
68
69#ifndef _DOUBLE_IS_32BITS
70
71static const double a[] = { -0.64124943423745581147e+02,
72                             0.16383943563021534222e+02,
73                            -0.78956112887481257267 };
74static const double b[] = { -0.76949932108494879777e+03,
75                             0.31203222091924532844e+03,
76                            -0.35667977739034646171e+02 };
77static const double C1 =  22713.0 / 32768.0;
78static const double C2 =  1.428606820309417232e-06;
79static const double C3 =  0.43429448190325182765;
80
81double
82logarithm (double x,
83        int ten)
84{
85  int N;
86  double f, w, z;
87
88  /* Check for range and domain errors here. */
89  if (x == 0.0)
90    {
91      errno = ERANGE;
92      return (-z_infinity.d);
93    }
94  else if (x < 0.0)
95    {
96      errno = EDOM;
97      return (z_notanum.d);
98    }
99  else if (!isfinite(x))
100    {
101      if (isnan(x))
102        return (z_notanum.d);
103      else
104        return (z_infinity.d);
105    }
106
107  /* Get the exponent and mantissa where x = f * 2^N. */
108  f = frexp (x, &N);
109
110  z = f - 0.5;
111
112  if (f > __SQRT_HALF)
113    z = (z - 0.5) / (f * 0.5 + 0.5);
114  else
115    {
116      N--;
117      z /= (z * 0.5 + 0.5);
118    }
119  w = z * z;
120
121  /* Use Newton's method with 4 terms. */
122  z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
123
124  if (N != 0)
125    z = (N * C2 + z) + N * C1;
126
127  if (ten)
128    z *= C3;
129
130  return (z);
131}
132
133#endif /* _DOUBLE_IS_32BITS */
Note: See TracBrowser for help on using the repository browser.