src/libm/e_log.c
author Sam Lantinga <slouken@libsdl.org>
Mon, 22 Jun 2015 23:36:06 -0700
changeset 9776 952ff8a5076f
parent 6044 35448a5ea044
child 11683 48bcba563d9c
permissions -rw-r--r--
Fixed bug 3030 - SDL_RecreateWindow fails to restore title, icon, etc.

Adam M.

It loses the title and icon when window recreation fails. For instance, this may happen when trying to create an OpenGL ES window on a system that doesn't support it. But at that point, the title and icon have already been lost.
     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
    99 __ieee754_log(x)
   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 }