OpenCores
URL https://opencores.org/ocsvn/openrisc_me/openrisc_me/trunk

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [newlib-1.18.0/] [newlib/] [libm/] [mathfp/] [s_logarithm.c] - Blame information for rev 258

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 207 jeremybenn
 
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
/*
11
FUNCTION
12
       <<log>>, <<logf>>, <<log10>>, <<log10f>>, <<logarithm>>, <<logarithmf>>---natural or base 10 logarithms
13
 
14
INDEX
15
    log
16
INDEX
17
    logf
18
INDEX
19
    log10
20
INDEX
21
    log10f
22
 
23
ANSI_SYNOPSIS
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
 
30
TRAD_SYNOPSIS
31
       #include <math.h>
32
       double log(<[x]>);
33
       double <[x]>;
34
 
35
       float logf(<[x]>);
36
       float <[x]>;
37
 
38
       double log10(<[x]>);
39
       double <[x]>;
40
 
41
       float log10f(<[x]>);
42
       float <[x]>;
43
 
44
DESCRIPTION
45
Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e
46
(where e is the base of the natural system of logarithms, 2.71828@dots{}) or
47
base 10.
48
<<log>> and <<logf>> are identical save for the return and argument types.
49
<<log10>> and <<log10f>> are identical save for the return and argument types.
50
 
51
RETURNS
52
Normally, returns the calculated value.  When <[x]> is zero, the
53
returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>.
54
When <[x]> is negative, the returned value is <<-HUGE_VAL>> and
55
<<errno>> is set to <<EDOM>>.  You can control the error behavior via
56
<<matherr>>.
57
 
58
PORTABILITY
59
<<log>> is ANSI. <<logf>> is an extension.
60
 
61
<<log10>> is ANSI. <<log10f>> is an extension.
62
*/
63
 
64
 
65
/******************************************************************
66
 * Logarithm
67
 *
68
 * Input:
69
 *   x - floating point value
70
 *   ten - indicates base ten numbers
71
 *
72
 * Output:
73
 *   logarithm of x
74
 *
75
 * Description:
76
 *   This routine calculates logarithms.
77
 *
78
 *****************************************************************/
79
 
80
#include "fdlibm.h"
81
#include "zmath.h"
82
 
83
#ifndef _DOUBLE_IS_32BITS
84
 
85
static const double a[] = { -0.64124943423745581147e+02,
86
                             0.16383943563021534222e+02,
87
                            -0.78956112887481257267 };
88
static const double b[] = { -0.76949932108494879777e+03,
89
                             0.31203222091924532844e+03,
90
                            -0.35667977739034646171e+02 };
91
static const double C1 =  22713.0 / 32768.0;
92
static const double C2 =  1.428606820309417232e-06;
93
static const double C3 =  0.43429448190325182765;
94
 
95
double
96
_DEFUN (logarithm, (double, int),
97
        double x _AND
98
        int ten)
99
{
100
  int N;
101
  double f, w, z;
102
 
103
  /* Check for range and domain errors here. */
104
  if (x == 0.0)
105
    {
106
      errno = ERANGE;
107
      return (-z_infinity.d);
108
    }
109
  else if (x < 0.0)
110
    {
111
      errno = EDOM;
112
      return (z_notanum.d);
113
    }
114
  else if (!isfinite(x))
115
    {
116
      if (isnan(x))
117
        return (z_notanum.d);
118
      else
119
        return (z_infinity.d);
120
    }
121
 
122
  /* Get the exponent and mantissa where x = f * 2^N. */
123
  f = frexp (x, &N);
124
 
125
  z = f - 0.5;
126
 
127
  if (f > __SQRT_HALF)
128
    z = (z - 0.5) / (f * 0.5 + 0.5);
129
  else
130
    {
131
      N--;
132
      z /= (z * 0.5 + 0.5);
133
    }
134
  w = z * z;
135
 
136
  /* Use Newton's method with 4 terms. */
137
  z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]);
138
 
139
  if (N != 0)
140
    z = (N * C2 + z) + N * C1;
141
 
142
  if (ten)
143
    z *= C3;
144
 
145
  return (z);
146
}
147
 
148
#endif /* _DOUBLE_IS_32BITS */

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.