SDL  2.0
e_log.c
Go to the documentation of this file.
1 /* @(#)e_log.c 5.1 93/09/24 */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static const char rcsid[] =
15  "$NetBSD: e_log.c,v 1.8 1995/05/10 20:45:49 jtc Exp $";
16 #endif
17 
18 /* __ieee754_log(x)
19  * Return the logrithm of x
20  *
21  * Method :
22  * 1. Argument Reduction: find k and f such that
23  * x = 2^k * (1+f),
24  * where sqrt(2)/2 < 1+f < sqrt(2) .
25  *
26  * 2. Approximation of log(1+f).
27  * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
28  * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
29  * = 2s + s*R
30  * We use a special Reme algorithm on [0,0.1716] to generate
31  * a polynomial of degree 14 to approximate R The maximum error
32  * of this polynomial approximation is bounded by 2**-58.45. In
33  * other words,
34  * 2 4 6 8 10 12 14
35  * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
36  * (the values of Lg1 to Lg7 are listed in the program)
37  * and
38  * | 2 14 | -58.45
39  * | Lg1*s +...+Lg7*s - R(z) | <= 2
40  * | |
41  * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
42  * In order to guarantee error in log below 1ulp, we compute log
43  * by
44  * log(1+f) = f - s*(f - R) (if f is not too large)
45  * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
46  *
47  * 3. Finally, log(x) = k*ln2 + log(1+f).
48  * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
49  * Here ln2 is split into two floating point number:
50  * ln2_hi + ln2_lo,
51  * where n*ln2_hi is always exact for |n| < 2000.
52  *
53  * Special cases:
54  * log(x) is NaN with signal if x < 0 (including -INF) ;
55  * log(+INF) is +INF; log(0) is -INF with signal;
56  * log(NaN) is that NaN with no signal.
57  *
58  * Accuracy:
59  * according to an error analysis, the error is always less than
60  * 1 ulp (unit in the last place).
61  *
62  * Constants:
63  * The hexadecimal values are the intended ones for the following
64  * constants. The decimal values may be used, provided that the
65  * compiler will convert from decimal to binary accurately enough
66  * to produce the hexadecimal values shown.
67  */
68 
69 #include "math_libm.h"
70 #include "math_private.h"
71 
72 #ifdef __STDC__
73 static const double
74 #else
75 static double
76 #endif
77  ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
78  ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
79  two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
80  Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
81  Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
82  Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
83  Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
84  Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
85  Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
86  Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
87 
88 #ifdef __STDC__
89 static const double zero = 0.0;
90 #else
91 static double zero = 0.0;
92 #endif
93 
94 #ifdef __STDC__
95 double attribute_hidden
96 __ieee754_log(double x)
97 #else
98 double attribute_hidden
100  double x;
101 #endif
102 {
103  double hfsq, f, s, z, R, w, t1, t2, dk;
104  int32_t k, hx, i, j;
105  u_int32_t lx;
106 
107  EXTRACT_WORDS(hx, lx, x);
108 
109  k = 0;
110  if (hx < 0x00100000) { /* x < 2**-1022 */
111  if (((hx & 0x7fffffff) | lx) == 0)
112  return -two54 / zero; /* log(+-0)=-inf */
113  if (hx < 0)
114  return (x - x) / zero; /* log(-#) = NaN */
115  k -= 54;
116  x *= two54; /* subnormal number, scale up x */
117  GET_HIGH_WORD(hx, x);
118  }
119  if (hx >= 0x7ff00000)
120  return x + x;
121  k += (hx >> 20) - 1023;
122  hx &= 0x000fffff;
123  i = (hx + 0x95f64) & 0x100000;
124  SET_HIGH_WORD(x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
125  k += (i >> 20);
126  f = x - 1.0;
127  if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */
128  if (f == zero) {
129  if (k == 0)
130  return zero;
131  else {
132  dk = (double) k;
133  return dk * ln2_hi + dk * ln2_lo;
134  }
135  }
136  R = f * f * (0.5 - 0.33333333333333333 * f);
137  if (k == 0)
138  return f - R;
139  else {
140  dk = (double) k;
141  return dk * ln2_hi - ((R - dk * ln2_lo) - f);
142  }
143  }
144  s = f / (2.0 + f);
145  dk = (double) k;
146  z = s * s;
147  i = hx - 0x6147a;
148  w = z * z;
149  j = 0x6b851 - hx;
150  t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
151  t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
152  i |= j;
153  R = t2 + t1;
154  if (i > 0) {
155  hfsq = 0.5 * f * f;
156  if (k == 0)
157  return f - (hfsq - s * (hfsq + R));
158  else
159  return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) -
160  f);
161  } else {
162  if (k == 0)
163  return f - s * (f - R);
164  else
165  return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f);
166  }
167 }
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:103
static double zero
Definition: e_log.c:91
GLdouble GLdouble z
GLdouble s
Definition: SDL_opengl.h:2063
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int int int return Display Window Cursor return Display Window return Display Drawable GC int int unsigned int unsigned int return Display Drawable GC int int _Xconst char int return Display Drawable GC int int unsigned int unsigned int return Display return Display Cursor return Display GC return XModifierKeymap return char Display Window int return Display return Display Atom return Display Window XWindowAttributes return Display Window return Display XEvent Bool(*) XPointer return Display Window Bool unsigned int int int Window Cursor Time return Display Window int return KeySym return Display _Xconst char Bool return Display _Xconst char return XKeyEvent char int KeySym XComposeStatus return Display int int int XVisualInfo return Display Window int int return _Xconst char return Display XEvent return Display Drawable GC XImage int int int int unsigned int unsigned int return Display Window Window Window int int int int unsigned int return Display Window Window int int return Display Window unsigned int unsigned int return Display Window Bool long XEvent return Display GC unsigned long return Display Window int Time return Display Window Window return Display Window unsigned long return Display Window XSizeHints Display Colormap XColor int return char int XTextProperty return XFontStruct _Xconst char int int int int XCharStruct return Display Window return Display Time return Display Colormap return Display Window Window int int unsigned int unsigned int int int return Display Window int return XExtensionInfo Display char XExtensionHooks int XPointer return XExtensionInfo XExtensionInfo Display return Display return Display unsigned long Display GC Display char long Display xReply int Bool return Display Bool return Display int SDL_X11_XESetEventToWireRetType return Display Window Window Window Window unsigned int return Display XShmSegmentInfo return Display Drawable GC XImage int int int int unsigned int unsigned int Boo k)
Definition: SDL_x11sym.h:213
static double Lg1
Definition: e_log.c:80
static double ln2_lo
Definition: e_log.c:78
GLfloat f
static double Lg6
Definition: e_log.c:85
#define attribute_hidden
Definition: math_private.h:24
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:131
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int in i)
Definition: SDL_x11sym.h:50
unsigned int u_int32_t
Definition: math_private.h:29
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:93
static double Lg5
Definition: e_log.c:84
GLubyte GLubyte GLubyte GLubyte w
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int int in j)
Definition: SDL_x11sym.h:50
static double Lg4
Definition: e_log.c:83
static double Lg3
Definition: e_log.c:82
double attribute_hidden __ieee754_log(double x)
Definition: e_log.c:99
static double Lg7
Definition: e_log.c:86
static double two54
Definition: e_log.c:79
static double Lg2
Definition: e_log.c:81
static double ln2_hi
Definition: e_log.c:77