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

Subversion Repositories openrisc

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

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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