src/libm/e_pow.c
author Sam Lantinga <slouken@libsdl.org>
Mon, 22 Jun 2015 23:36:06 -0700
changeset 9776 952ff8a5076f
parent 7678 286c42d7c5ed
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_pow.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@2756
    14
static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $";
slouken@2756
    15
#endif
slouken@2756
    16
slouken@2756
    17
/* __ieee754_pow(x,y) return x**y
slouken@2756
    18
 *
slouken@2756
    19
 *		      n
slouken@2756
    20
 * Method:  Let x =  2   * (1+f)
slouken@2756
    21
 *	1. Compute and return log2(x) in two pieces:
slouken@2756
    22
 *		log2(x) = w1 + w2,
slouken@2756
    23
 *	   where w1 has 53-24 = 29 bit trailing zeros.
slouken@2756
    24
 *	2. Perform y*log2(x) = n+y' by simulating muti-precision
slouken@2756
    25
 *	   arithmetic, where |y'|<=0.5.
slouken@2756
    26
 *	3. Return x**y = 2**n*exp(y'*log2)
slouken@2756
    27
 *
slouken@2756
    28
 * Special cases:
slouken@2756
    29
 *	1.  (anything) ** 0  is 1
slouken@2756
    30
 *	2.  (anything) ** 1  is itself
slouken@2756
    31
 *	3.  (anything) ** NAN is NAN
slouken@2756
    32
 *	4.  NAN ** (anything except 0) is NAN
slouken@2756
    33
 *	5.  +-(|x| > 1) **  +INF is +INF
slouken@2756
    34
 *	6.  +-(|x| > 1) **  -INF is +0
slouken@2756
    35
 *	7.  +-(|x| < 1) **  +INF is +0
slouken@2756
    36
 *	8.  +-(|x| < 1) **  -INF is +INF
slouken@2756
    37
 *	9.  +-1         ** +-INF is NAN
slouken@2756
    38
 *	10. +0 ** (+anything except 0, NAN)               is +0
slouken@2756
    39
 *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
slouken@2756
    40
 *	12. +0 ** (-anything except 0, NAN)               is +INF
slouken@2756
    41
 *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
slouken@2756
    42
 *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
slouken@2756
    43
 *	15. +INF ** (+anything except 0,NAN) is +INF
slouken@2756
    44
 *	16. +INF ** (-anything except 0,NAN) is +0
slouken@2756
    45
 *	17. -INF ** (anything)  = -0 ** (-anything)
slouken@2756
    46
 *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
slouken@2756
    47
 *	19. (-anything except 0 and inf) ** (non-integer) is NAN
slouken@2756
    48
 *
slouken@2756
    49
 * Accuracy:
slouken@2756
    50
 *	pow(x,y) returns x**y nearly rounded. In particular
slouken@2756
    51
 *			pow(integer,integer)
slouken@2756
    52
 *	always returns the correct integer provided it is
slouken@2756
    53
 *	representable.
slouken@2756
    54
 *
slouken@2756
    55
 * Constants :
slouken@2756
    56
 * The hexadecimal values are the intended ones for the following
slouken@2756
    57
 * constants. The decimal values may be used, provided that the
slouken@2756
    58
 * compiler will convert from decimal to binary accurately enough
slouken@2756
    59
 * to produce the hexadecimal values shown.
slouken@2756
    60
 */
slouken@2756
    61
slouken@6044
    62
#include "math_libm.h"
slouken@2756
    63
#include "math_private.h"
slouken@2756
    64
slouken@2756
    65
libm_hidden_proto(scalbn)
slouken@2756
    66
    libm_hidden_proto(fabs)
slouken@2756
    67
#ifdef __STDC__
slouken@2756
    68
     static const double
slouken@2756
    69
#else
slouken@2756
    70
     static double
slouken@2756
    71
#endif
slouken@2756
    72
       bp[] = { 1.0, 1.5, }, dp_h[] = {
slouken@2756
    73
     0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
slouken@2756
    74
slouken@2756
    75
         dp_l[] = {
slouken@2756
    76
     0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
slouken@2756
    77
slouken@2756
    78
         zero = 0.0, one = 1.0, two = 2.0, two53 = 9007199254740992.0,  /* 0x43400000, 0x00000000 */
slouken@3337
    79
         huge_val = 1.0e300, tiny = 1.0e-300,
slouken@2756
    80
         /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
slouken@2756
    81
         L1 = 5.99999999999994648725e-01,       /* 0x3FE33333, 0x33333303 */
slouken@2756
    82
         L2 = 4.28571428578550184252e-01,       /* 0x3FDB6DB6, 0xDB6FABFF */
slouken@2756
    83
         L3 = 3.33333329818377432918e-01,       /* 0x3FD55555, 0x518F264D */
slouken@2756
    84
         L4 = 2.72728123808534006489e-01,       /* 0x3FD17460, 0xA91D4101 */
slouken@2756
    85
         L5 = 2.30660745775561754067e-01,       /* 0x3FCD864A, 0x93C9DB65 */
slouken@2756
    86
         L6 = 2.06975017800338417784e-01,       /* 0x3FCA7E28, 0x4A454EEF */
slouken@2756
    87
         P1 = 1.66666666666666019037e-01,       /* 0x3FC55555, 0x5555553E */
slouken@2756
    88
         P2 = -2.77777777770155933842e-03,      /* 0xBF66C16C, 0x16BEBD93 */
slouken@2756
    89
         P3 = 6.61375632143793436117e-05,       /* 0x3F11566A, 0xAF25DE2C */
slouken@2756
    90
         P4 = -1.65339022054652515390e-06,      /* 0xBEBBBD41, 0xC5D26BF1 */
slouken@2756
    91
         P5 = 4.13813679705723846039e-08,       /* 0x3E663769, 0x72BEA4D0 */
slouken@2756
    92
         lg2 = 6.93147180559945286227e-01,      /* 0x3FE62E42, 0xFEFA39EF */
slouken@2756
    93
         lg2_h = 6.93147182464599609375e-01,    /* 0x3FE62E43, 0x00000000 */
slouken@2756
    94
         lg2_l = -1.90465429995776804525e-09,   /* 0xBE205C61, 0x0CA86C39 */
slouken@2756
    95
         ovt = 8.0085662595372944372e-0017,     /* -(1024-log2(ovfl+.5ulp)) */
slouken@2756
    96
         cp = 9.61796693925975554329e-01,       /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
slouken@2756
    97
         cp_h = 9.61796700954437255859e-01,     /* 0x3FEEC709, 0xE0000000 =(float)cp */
slouken@2756
    98
         cp_l = -7.02846165095275826516e-09,    /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */
slouken@2756
    99
         ivln2 = 1.44269504088896338700e+00,    /* 0x3FF71547, 0x652B82FE =1/ln2 */
slouken@2756
   100
         ivln2_h = 1.44269502162933349609e+00,  /* 0x3FF71547, 0x60000000 =24b 1/ln2 */
slouken@2756
   101
         ivln2_l = 1.92596299112661746887e-08;  /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */
slouken@2756
   102
slouken@2756
   103
#ifdef __STDC__
slouken@2756
   104
     double attribute_hidden __ieee754_pow(double x, double y)
slouken@2756
   105
#else
slouken@2756
   106
     double attribute_hidden __ieee754_pow(x, y)
slouken@2756
   107
     double x, y;
slouken@2756
   108
#endif
slouken@2756
   109
     {
slouken@2756
   110
         double z, ax, z_h, z_l, p_h, p_l;
slouken@2756
   111
         double y1, t1, t2, r, s, t, u, v, w;
slouken@2756
   112
         int32_t i, j, k, yisint, n;
slouken@2756
   113
         int32_t hx, hy, ix, iy;
slouken@2756
   114
         u_int32_t lx, ly;
slouken@2756
   115
slouken@2756
   116
         EXTRACT_WORDS(hx, lx, x);
slouken@2756
   117
         EXTRACT_WORDS(hy, ly, y);
slouken@2756
   118
         ix = hx & 0x7fffffff;
slouken@2756
   119
         iy = hy & 0x7fffffff;
slouken@2756
   120
slouken@2756
   121
         /* y==zero: x**0 = 1 */
slouken@2756
   122
         if ((iy | ly) == 0)
slouken@2756
   123
             return one;
slouken@2756
   124
slouken@2756
   125
         /* +-NaN return x+y */
slouken@2756
   126
         if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) ||
slouken@2756
   127
             iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0)))
slouken@2756
   128
             return x + y;
slouken@2756
   129
slouken@2756
   130
         /* determine if y is an odd int when x < 0
slouken@2756
   131
          * yisint = 0       ... y is not an integer
slouken@2756
   132
          * yisint = 1       ... y is an odd int
slouken@2756
   133
          * yisint = 2       ... y is an even int
slouken@2756
   134
          */
slouken@2756
   135
         yisint = 0;
slouken@2756
   136
         if (hx < 0) {
slouken@2756
   137
             if (iy >= 0x43400000)
slouken@2756
   138
                 yisint = 2;    /* even integer y */
slouken@2756
   139
             else if (iy >= 0x3ff00000) {
slouken@2756
   140
                 k = (iy >> 20) - 0x3ff;        /* exponent */
slouken@2756
   141
                 if (k > 20) {
slouken@2756
   142
                     j = ly >> (52 - k);
slouken@2756
   143
                     if ((j << (52 - k)) == ly)
slouken@2756
   144
                         yisint = 2 - (j & 1);
slouken@2756
   145
                 } else if (ly == 0) {
slouken@2756
   146
                     j = iy >> (20 - k);
slouken@2756
   147
                     if ((j << (20 - k)) == iy)
slouken@2756
   148
                         yisint = 2 - (j & 1);
slouken@2756
   149
                 }
slouken@2756
   150
             }
slouken@2756
   151
         }
slouken@2756
   152
slouken@2756
   153
         /* special value of y */
slouken@2756
   154
         if (ly == 0) {
slouken@2756
   155
             if (iy == 0x7ff00000) {    /* y is +-inf */
slouken@2756
   156
                 if (((ix - 0x3ff00000) | lx) == 0)
slouken@2756
   157
                     return y - y;      /* inf**+-1 is NaN */
slouken@2756
   158
                 else if (ix >= 0x3ff00000)     /* (|x|>1)**+-inf = inf,0 */
slouken@2756
   159
                     return (hy >= 0) ? y : zero;
slouken@2756
   160
                 else           /* (|x|<1)**-,+inf = inf,0 */
slouken@2756
   161
                     return (hy < 0) ? -y : zero;
slouken@2756
   162
             }
slouken@2756
   163
             if (iy == 0x3ff00000) {    /* y is  +-1 */
slouken@2756
   164
                 if (hy < 0)
slouken@2756
   165
                     return one / x;
slouken@2756
   166
                 else
slouken@2756
   167
                     return x;
slouken@2756
   168
             }
slouken@2756
   169
             if (hy == 0x40000000)
slouken@2756
   170
                 return x * x;  /* y is  2 */
slouken@2756
   171
             if (hy == 0x3fe00000) {    /* y is  0.5 */
slouken@2756
   172
                 if (hx >= 0)   /* x >= +0 */
slouken@2756
   173
                     return __ieee754_sqrt(x);
slouken@2756
   174
             }
slouken@2756
   175
         }
slouken@2756
   176
slouken@2756
   177
         ax = fabs(x);
slouken@2756
   178
         /* special value of x */
slouken@2756
   179
         if (lx == 0) {
slouken@2756
   180
             if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) {
gabomdq@7678
   181
                 z = ax;        /* x is +-0,+-inf,+-1 */
slouken@2756
   182
                 if (hy < 0)
slouken@2756
   183
                     z = one / z;       /* z = (1/|x|) */
slouken@2756
   184
                 if (hx < 0) {
slouken@2756
   185
                     if (((ix - 0x3ff00000) | yisint) == 0) {
slouken@2756
   186
                         z = (z - z) / (z - z); /* (-1)**non-int is NaN */
slouken@2756
   187
                     } else if (yisint == 1)
slouken@2756
   188
                         z = -z;        /* (x<0)**odd = -(|x|**odd) */
slouken@2756
   189
                 }
slouken@2756
   190
                 return z;
slouken@2756
   191
             }
slouken@2756
   192
         }
slouken@2756
   193
slouken@2756
   194
         /* (x<0)**(non-int) is NaN */
slouken@2756
   195
         if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
slouken@2756
   196
             return (x - x) / (x - x);
slouken@2756
   197
slouken@2756
   198
         /* |y| is huge */
slouken@2756
   199
         if (iy > 0x41e00000) { /* if |y| > 2**31 */
slouken@2756
   200
             if (iy > 0x43f00000) {     /* if |y| > 2**64, must o/uflow */
slouken@2756
   201
                 if (ix <= 0x3fefffff)
slouken@3337
   202
                     return (hy < 0) ? huge_val * huge_val : tiny * tiny;
slouken@2756
   203
                 if (ix >= 0x3ff00000)
slouken@3337
   204
                     return (hy > 0) ? huge_val * huge_val : tiny * tiny;
slouken@2756
   205
             }
slouken@2756
   206
             /* over/underflow if x is not close to one */
slouken@2756
   207
             if (ix < 0x3fefffff)
slouken@3337
   208
                 return (hy < 0) ? huge_val * huge_val : tiny * tiny;
slouken@2756
   209
             if (ix > 0x3ff00000)
slouken@3337
   210
                 return (hy > 0) ? huge_val * huge_val : tiny * tiny;
slouken@2756
   211
             /* now |1-x| is tiny <= 2**-20, suffice to compute
slouken@2756
   212
                log(x) by x-x^2/2+x^3/3-x^4/4 */
slouken@2756
   213
             t = x - 1;         /* t has 20 trailing zeros */
slouken@2756
   214
             w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25));
slouken@2756
   215
             u = ivln2_h * t;   /* ivln2_h has 21 sig. bits */
slouken@2756
   216
             v = t * ivln2_l - w * ivln2;
slouken@2756
   217
             t1 = u + v;
slouken@2756
   218
             SET_LOW_WORD(t1, 0);
slouken@2756
   219
             t2 = v - (t1 - u);
slouken@2756
   220
         } else {
slouken@2756
   221
             double s2, s_h, s_l, t_h, t_l;
slouken@2756
   222
             n = 0;
slouken@2756
   223
             /* take care subnormal number */
slouken@2756
   224
             if (ix < 0x00100000) {
slouken@2756
   225
                 ax *= two53;
slouken@2756
   226
                 n -= 53;
slouken@2756
   227
                 GET_HIGH_WORD(ix, ax);
slouken@2756
   228
             }
slouken@2756
   229
             n += ((ix) >> 20) - 0x3ff;
slouken@2756
   230
             j = ix & 0x000fffff;
slouken@2756
   231
             /* determine interval */
slouken@2756
   232
             ix = j | 0x3ff00000;       /* normalize ix */
slouken@2756
   233
             if (j <= 0x3988E)
slouken@2756
   234
                 k = 0;         /* |x|<sqrt(3/2) */
slouken@2756
   235
             else if (j < 0xBB67A)
slouken@2756
   236
                 k = 1;         /* |x|<sqrt(3)   */
slouken@2756
   237
             else {
slouken@2756
   238
                 k = 0;
slouken@2756
   239
                 n += 1;
slouken@2756
   240
                 ix -= 0x00100000;
slouken@2756
   241
             }
slouken@2756
   242
             SET_HIGH_WORD(ax, ix);
slouken@2756
   243
slouken@2756
   244
             /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
slouken@2756
   245
             u = ax - bp[k];    /* bp[0]=1.0, bp[1]=1.5 */
slouken@2756
   246
             v = one / (ax + bp[k]);
slouken@2756
   247
             s = u * v;
slouken@2756
   248
             s_h = s;
slouken@2756
   249
             SET_LOW_WORD(s_h, 0);
slouken@2756
   250
             /* t_h=ax+bp[k] High */
slouken@2756
   251
             t_h = zero;
slouken@2756
   252
             SET_HIGH_WORD(t_h,
slouken@2756
   253
                           ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
slouken@2756
   254
             t_l = ax - (t_h - bp[k]);
slouken@2756
   255
             s_l = v * ((u - s_h * t_h) - s_h * t_l);
slouken@2756
   256
             /* compute log(ax) */
slouken@2756
   257
             s2 = s * s;
slouken@2756
   258
             r = s2 * s2 * (L1 +
slouken@2756
   259
                            s2 * (L2 +
slouken@2756
   260
                                  s2 * (L3 +
slouken@2756
   261
                                        s2 * (L4 + s2 * (L5 + s2 * L6)))));
slouken@2756
   262
             r += s_l * (s_h + s);
slouken@2756
   263
             s2 = s_h * s_h;
slouken@2756
   264
             t_h = 3.0 + s2 + r;
slouken@2756
   265
             SET_LOW_WORD(t_h, 0);
slouken@2756
   266
             t_l = r - ((t_h - 3.0) - s2);
slouken@2756
   267
             /* u+v = s*(1+...) */
slouken@2756
   268
             u = s_h * t_h;
slouken@2756
   269
             v = s_l * t_h + t_l * s;
slouken@2756
   270
             /* 2/(3log2)*(s+...) */
slouken@2756
   271
             p_h = u + v;
slouken@2756
   272
             SET_LOW_WORD(p_h, 0);
slouken@2756
   273
             p_l = v - (p_h - u);
slouken@2756
   274
             z_h = cp_h * p_h;  /* cp_h+cp_l = 2/(3*log2) */
slouken@2756
   275
             z_l = cp_l * p_h + p_l * cp + dp_l[k];
slouken@2756
   276
             /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
slouken@2756
   277
             t = (double) n;
slouken@2756
   278
             t1 = (((z_h + z_l) + dp_h[k]) + t);
slouken@2756
   279
             SET_LOW_WORD(t1, 0);
slouken@2756
   280
             t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
slouken@2756
   281
         }
slouken@2756
   282
slouken@2756
   283
         s = one;               /* s (sign of result -ve**odd) = -1 else = 1 */
slouken@2756
   284
         if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
slouken@2756
   285
             s = -one;          /* (-ve)**(odd int) */
slouken@2756
   286
slouken@2756
   287
         /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
slouken@2756
   288
         y1 = y;
slouken@2756
   289
         SET_LOW_WORD(y1, 0);
slouken@2756
   290
         p_l = (y - y1) * t1 + y * t2;
slouken@2756
   291
         p_h = y1 * t1;
slouken@2756
   292
         z = p_l + p_h;
slouken@2756
   293
         EXTRACT_WORDS(j, i, z);
slouken@2756
   294
         if (j >= 0x40900000) { /* z >= 1024 */
slouken@2756
   295
             if (((j - 0x40900000) | i) != 0)   /* if z > 1024 */
slouken@3337
   296
                 return s * huge_val * huge_val;        /* overflow */
slouken@2756
   297
             else {
slouken@2756
   298
                 if (p_l + ovt > z - p_h)
slouken@3337
   299
                     return s * huge_val * huge_val;    /* overflow */
slouken@2756
   300
             }
slouken@2756
   301
         } else if ((j & 0x7fffffff) >= 0x4090cc00) {   /* z <= -1075 */
slouken@2756
   302
             if (((j - 0xc090cc00) | i) != 0)   /* z < -1075 */
slouken@2756
   303
                 return s * tiny * tiny;        /* underflow */
slouken@2756
   304
             else {
slouken@2756
   305
                 if (p_l <= z - p_h)
slouken@2756
   306
                     return s * tiny * tiny;    /* underflow */
slouken@2756
   307
             }
slouken@2756
   308
         }
slouken@2756
   309
         /*
slouken@2756
   310
          * compute 2**(p_h+p_l)
slouken@2756
   311
          */
slouken@2756
   312
         i = j & 0x7fffffff;
slouken@2756
   313
         k = (i >> 20) - 0x3ff;
slouken@2756
   314
         n = 0;
slouken@2756
   315
         if (i > 0x3fe00000) {  /* if |z| > 0.5, set n = [z+0.5] */
slouken@2756
   316
             n = j + (0x00100000 >> (k + 1));
slouken@2756
   317
             k = ((n & 0x7fffffff) >> 20) - 0x3ff;      /* new k for n */
slouken@2756
   318
             t = zero;
slouken@2756
   319
             SET_HIGH_WORD(t, n & ~(0x000fffff >> k));
slouken@2756
   320
             n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
slouken@2756
   321
             if (j < 0)
slouken@2756
   322
                 n = -n;
slouken@2756
   323
             p_h -= t;
slouken@2756
   324
         }
slouken@2756
   325
         t = p_l + p_h;
slouken@2756
   326
         SET_LOW_WORD(t, 0);
slouken@2756
   327
         u = t * lg2_h;
slouken@2756
   328
         v = (p_l - (t - p_h)) * lg2 + t * lg2_l;
slouken@2756
   329
         z = u + v;
slouken@2756
   330
         w = v - (z - u);
slouken@2756
   331
         t = z * z;
slouken@2756
   332
         t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
slouken@2756
   333
         r = (z * t1) / (t1 - two) - (w + z * w);
slouken@2756
   334
         z = one - (r - z);
slouken@2756
   335
         GET_HIGH_WORD(j, z);
slouken@2756
   336
         j += (n << 20);
slouken@2756
   337
         if ((j >> 20) <= 0)
slouken@2756
   338
             z = scalbn(z, n);  /* subnormal output */
slouken@2756
   339
         else
slouken@2756
   340
             SET_HIGH_WORD(z, j);
slouken@2756
   341
         return s * z;
slouken@2756
   342
     }