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