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