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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [log10q.c] - Blame information for rev 740

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn
/*                                                      log10l.c
2
 *
3
 *      Common logarithm, 128-bit long double precision
4
 *
5
 *
6
 *
7
 * SYNOPSIS:
8
 *
9
 * long double x, y, log10l();
10
 *
11
 * y = log10l( x );
12
 *
13
 *
14
 *
15
 * DESCRIPTION:
16
 *
17
 * Returns the base 10 logarithm of x.
18
 *
19
 * The argument is separated into its exponent and fractional
20
 * parts.  If the exponent is between -1 and +1, the logarithm
21
 * of the fraction is approximated by
22
 *
23
 *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
24
 *
25
 * Otherwise, setting  z = 2(x-1)/x+1),
26
 *
27
 *     log(x) = z + z^3 P(z)/Q(z).
28
 *
29
 *
30
 *
31
 * ACCURACY:
32
 *
33
 *                      Relative error:
34
 * arithmetic   domain     # trials      peak         rms
35
 *    IEEE      0.5, 2.0     30000      2.3e-34     4.9e-35
36
 *    IEEE     exp(+-10000)  30000      1.0e-34     4.1e-35
37
 *
38
 * In the tests over the interval exp(+-10000), the logarithms
39
 * of the random arguments were uniformly distributed over
40
 * [-10000, +10000].
41
 *
42
 */
43
 
44
/*
45
   Cephes Math Library Release 2.2:  January, 1991
46
   Copyright 1984, 1991 by Stephen L. Moshier
47
   Adapted for glibc November, 2001
48
 
49
    This library is free software; you can redistribute it and/or
50
    modify it under the terms of the GNU Lesser General Public
51
    License as published by the Free Software Foundation; either
52
    version 2.1 of the License, or (at your option) any later version.
53
 
54
    This library is distributed in the hope that it will be useful,
55
    but WITHOUT ANY WARRANTY; without even the implied warranty of
56
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
57
    Lesser General Public License for more details.
58
 
59
    You should have received a copy of the GNU Lesser General Public
60
    License along with this library; if not, write to the Free Software
61
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
62
 
63
 */
64
 
65
#include "quadmath-imp.h"
66
 
67
/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
68
 * 1/sqrt(2) <= x < sqrt(2)
69
 * Theoretical peak relative error = 5.3e-37,
70
 * relative peak error spread = 2.3e-14
71
 */
72
static const __float128 P[13] =
73
{
74
  1.313572404063446165910279910527789794488E4Q,
75
  7.771154681358524243729929227226708890930E4Q,
76
  2.014652742082537582487669938141683759923E5Q,
77
  3.007007295140399532324943111654767187848E5Q,
78
  2.854829159639697837788887080758954924001E5Q,
79
  1.797628303815655343403735250238293741397E5Q,
80
  7.594356839258970405033155585486712125861E4Q,
81
  2.128857716871515081352991964243375186031E4Q,
82
  3.824952356185897735160588078446136783779E3Q,
83
  4.114517881637811823002128927449878962058E2Q,
84
  2.321125933898420063925789532045674660756E1Q,
85
  4.998469661968096229986658302195402690910E-1Q,
86
  1.538612243596254322971797716843006400388E-6Q
87
};
88
static const __float128 Q[12] =
89
{
90
  3.940717212190338497730839731583397586124E4Q,
91
  2.626900195321832660448791748036714883242E5Q,
92
  7.777690340007566932935753241556479363645E5Q,
93
  1.347518538384329112529391120390701166528E6Q,
94
  1.514882452993549494932585972882995548426E6Q,
95
  1.158019977462989115839826904108208787040E6Q,
96
  6.132189329546557743179177159925690841200E5Q,
97
  2.248234257620569139969141618556349415120E5Q,
98
  5.605842085972455027590989944010492125825E4Q,
99
  9.147150349299596453976674231612674085381E3Q,
100
  9.104928120962988414618126155557301584078E2Q,
101
  4.839208193348159620282142911143429644326E1Q
102
/* 1.000000000000000000000000000000000000000E0Q, */
103
};
104
 
105
/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
106
 * where z = 2(x-1)/(x+1)
107
 * 1/sqrt(2) <= x < sqrt(2)
108
 * Theoretical peak relative error = 1.1e-35,
109
 * relative peak error spread 1.1e-9
110
 */
111
static const __float128 R[6] =
112
{
113
  1.418134209872192732479751274970992665513E5Q,
114
 -8.977257995689735303686582344659576526998E4Q,
115
  2.048819892795278657810231591630928516206E4Q,
116
 -2.024301798136027039250415126250455056397E3Q,
117
  8.057002716646055371965756206836056074715E1Q,
118
 -8.828896441624934385266096344596648080902E-1Q
119
};
120
static const __float128 S[6] =
121
{
122
  1.701761051846631278975701529965589676574E6Q,
123
 -1.332535117259762928288745111081235577029E6Q,
124
  4.001557694070773974936904547424676279307E5Q,
125
 -5.748542087379434595104154610899551484314E4Q,
126
  3.998526750980007367835804959888064681098E3Q,
127
 -1.186359407982897997337150403816839480438E2Q
128
/* 1.000000000000000000000000000000000000000E0Q, */
129
};
130
 
131
static const __float128
132
/* log10(2) */
133
L102A = 0.3125Q,
134
L102B = -1.14700043360188047862611052755069732318101185E-2Q,
135
/* log10(e) */
136
L10EA = 0.5Q,
137
L10EB = -6.570551809674817234887108108339491770560299E-2Q,
138
/* sqrt(2)/2 */
139
SQRTH = 7.071067811865475244008443621048490392848359E-1Q;
140
 
141
 
142
 
143
/* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
144
 
145
static __float128
146
neval (__float128 x, const __float128 *p, int n)
147
{
148
  __float128 y;
149
 
150
  p += n;
151
  y = *p--;
152
  do
153
    {
154
      y = y * x + *p--;
155
    }
156
  while (--n > 0);
157
  return y;
158
}
159
 
160
 
161
/* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
162
 
163
static __float128
164
deval (__float128 x, const __float128 *p, int n)
165
{
166
  __float128 y;
167
 
168
  p += n;
169
  y = x + *p--;
170
  do
171
    {
172
      y = y * x + *p--;
173
    }
174
  while (--n > 0);
175
  return y;
176
}
177
 
178
 
179
 
180
__float128
181
log10q (__float128 x)
182
{
183
  __float128 z;
184
  __float128 y;
185
  int e;
186
  int64_t hx, lx;
187
 
188
/* Test for domain */
189
  GET_FLT128_WORDS64 (hx, lx, x);
190
  if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
191
    return (-1.0Q / (x - x));
192
  if (hx < 0)
193
    return (x - x) / (x - x);
194
  if (hx >= 0x7fff000000000000LL)
195
    return (x + x);
196
 
197
/* separate mantissa from exponent */
198
 
199
/* Note, frexp is used so that denormal numbers
200
 * will be handled properly.
201
 */
202
  x = frexpq (x, &e);
203
 
204
 
205
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
206
 * where z = 2(x-1)/x+1)
207
 */
208
  if ((e > 2) || (e < -2))
209
    {
210
      if (x < SQRTH)
211
        {                       /* 2( 2x-1 )/( 2x+1 ) */
212
          e -= 1;
213
          z = x - 0.5Q;
214
          y = 0.5Q * z + 0.5Q;
215
        }
216
      else
217
        {                       /*  2 (x-1)/(x+1)   */
218
          z = x - 0.5Q;
219
          z -= 0.5Q;
220
          y = 0.5Q * x + 0.5Q;
221
        }
222
      x = z / y;
223
      z = x * x;
224
      y = x * (z * neval (z, R, 5) / deval (z, S, 5));
225
      goto done;
226
    }
227
 
228
 
229
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
230
 
231
  if (x < SQRTH)
232
    {
233
      e -= 1;
234
      x = 2.0 * x - 1.0Q;       /*  2x - 1  */
235
    }
236
  else
237
    {
238
      x = x - 1.0Q;
239
    }
240
  z = x * x;
241
  y = x * (z * neval (x, P, 12) / deval (x, Q, 11));
242
  y = y - 0.5 * z;
243
 
244
done:
245
 
246
  /* Multiply log of fraction by log10(e)
247
   * and base 2 exponent by log10(2).
248
   */
249
  z = y * L10EB;
250
  z += x * L10EB;
251
  z += e * L102B;
252
  z += y * L10EA;
253
  z += x * L10EA;
254
  z += e * L102A;
255
  return (z);
256
}

powered by: WebSVN 2.1.0

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