src/libm/e_rem_pio2.c
author Sam Lantinga <slouken@libsdl.org>
Mon, 22 Jun 2015 23:36:06 -0700
changeset 9776 952ff8a5076f
parent 6044 35448a5ea044
child 10650 b6ec7005ca15
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_rem_pio2.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@2756
    15
    "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
slouken@2756
    16
#endif
slouken@2756
    17
slouken@2756
    18
/* __ieee754_rem_pio2(x,y)
slouken@2756
    19
 *
slouken@2756
    20
 * return the remainder of x rem pi/2 in y[0]+y[1]
slouken@2756
    21
 * use __kernel_rem_pio2()
slouken@2756
    22
 */
slouken@2756
    23
slouken@6044
    24
#include "math_libm.h"
slouken@2756
    25
#include "math_private.h"
slouken@2756
    26
slouken@2756
    27
libm_hidden_proto(fabs)
slouken@2756
    28
slouken@2756
    29
/*
slouken@2756
    30
 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
slouken@2756
    31
 */
slouken@2756
    32
#ifdef __STDC__
slouken@2756
    33
     static const int32_t two_over_pi[] = {
slouken@2756
    34
#else
slouken@2756
    35
     static int32_t two_over_pi[] = {
slouken@2756
    36
#endif
slouken@2756
    37
         0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
slouken@2756
    38
         0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
slouken@2756
    39
         0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
slouken@2756
    40
         0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
slouken@2756
    41
         0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
slouken@2756
    42
         0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
slouken@2756
    43
         0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
slouken@2756
    44
         0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
slouken@2756
    45
         0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
slouken@2756
    46
         0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
slouken@2756
    47
         0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
slouken@2756
    48
     };
slouken@2756
    49
slouken@2756
    50
#ifdef __STDC__
slouken@2756
    51
static const int32_t npio2_hw[] = {
slouken@2756
    52
#else
slouken@2756
    53
static int32_t npio2_hw[] = {
slouken@2756
    54
#endif
slouken@2756
    55
    0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
slouken@2756
    56
    0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
slouken@2756
    57
    0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
slouken@2756
    58
    0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
slouken@2756
    59
    0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
slouken@2756
    60
    0x404858EB, 0x404921FB,
slouken@2756
    61
};
slouken@2756
    62
slouken@2756
    63
/*
slouken@2756
    64
 * invpio2:  53 bits of 2/pi
slouken@2756
    65
 * pio2_1:   first  33 bit of pi/2
slouken@2756
    66
 * pio2_1t:  pi/2 - pio2_1
slouken@2756
    67
 * pio2_2:   second 33 bit of pi/2
slouken@2756
    68
 * pio2_2t:  pi/2 - (pio2_1+pio2_2)
slouken@2756
    69
 * pio2_3:   third  33 bit of pi/2
slouken@2756
    70
 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
slouken@2756
    71
 */
slouken@2756
    72
slouken@2756
    73
#ifdef __STDC__
slouken@2756
    74
static const double
slouken@2756
    75
#else
slouken@2756
    76
static double
slouken@2756
    77
#endif
slouken@2756
    78
  zero = 0.00000000000000000000e+00,    /* 0x00000000, 0x00000000 */
slouken@2756
    79
    half = 5.00000000000000000000e-01,  /* 0x3FE00000, 0x00000000 */
slouken@2756
    80
    two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
slouken@2756
    81
    invpio2 = 6.36619772367581382433e-01,       /* 0x3FE45F30, 0x6DC9C883 */
slouken@2756
    82
    pio2_1 = 1.57079632673412561417e+00,        /* 0x3FF921FB, 0x54400000 */
slouken@2756
    83
    pio2_1t = 6.07710050650619224932e-11,       /* 0x3DD0B461, 0x1A626331 */
slouken@2756
    84
    pio2_2 = 6.07710050630396597660e-11,        /* 0x3DD0B461, 0x1A600000 */
slouken@2756
    85
    pio2_2t = 2.02226624879595063154e-21,       /* 0x3BA3198A, 0x2E037073 */
slouken@2756
    86
    pio2_3 = 2.02226624871116645580e-21,        /* 0x3BA3198A, 0x2E000000 */
slouken@2756
    87
    pio2_3t = 8.47842766036889956997e-32;       /* 0x397B839A, 0x252049C1 */
slouken@2756
    88
slouken@2756
    89
#ifdef __STDC__
slouken@2756
    90
int32_t attribute_hidden
slouken@2756
    91
__ieee754_rem_pio2(double x, double *y)
slouken@2756
    92
#else
slouken@2756
    93
int32_t attribute_hidden
slouken@2756
    94
__ieee754_rem_pio2(x, y)
slouken@2756
    95
     double x, y[];
slouken@2756
    96
#endif
slouken@2756
    97
{
slouken@2756
    98
    double z = 0.0, w, t, r, fn;
slouken@2756
    99
    double tx[3];
slouken@2756
   100
    int32_t e0, i, j, nx, n, ix, hx;
slouken@2756
   101
    u_int32_t low;
slouken@2756
   102
slouken@2756
   103
    GET_HIGH_WORD(hx, x);       /* high word of x */
slouken@2756
   104
    ix = hx & 0x7fffffff;
slouken@2756
   105
    if (ix <= 0x3fe921fb) {     /* |x| ~<= pi/4 , no need for reduction */
slouken@2756
   106
        y[0] = x;
slouken@2756
   107
        y[1] = 0;
slouken@2756
   108
        return 0;
slouken@2756
   109
    }
slouken@2756
   110
    if (ix < 0x4002d97c) {      /* |x| < 3pi/4, special case with n=+-1 */
slouken@2756
   111
        if (hx > 0) {
slouken@2756
   112
            z = x - pio2_1;
slouken@2756
   113
            if (ix != 0x3ff921fb) {     /* 33+53 bit pi is good enough */
slouken@2756
   114
                y[0] = z - pio2_1t;
slouken@2756
   115
                y[1] = (z - y[0]) - pio2_1t;
slouken@2756
   116
            } else {            /* near pi/2, use 33+33+53 bit pi */
slouken@2756
   117
                z -= pio2_2;
slouken@2756
   118
                y[0] = z - pio2_2t;
slouken@2756
   119
                y[1] = (z - y[0]) - pio2_2t;
slouken@2756
   120
            }
slouken@2756
   121
            return 1;
slouken@2756
   122
        } else {                /* negative x */
slouken@2756
   123
            z = x + pio2_1;
slouken@2756
   124
            if (ix != 0x3ff921fb) {     /* 33+53 bit pi is good enough */
slouken@2756
   125
                y[0] = z + pio2_1t;
slouken@2756
   126
                y[1] = (z - y[0]) + pio2_1t;
slouken@2756
   127
            } else {            /* near pi/2, use 33+33+53 bit pi */
slouken@2756
   128
                z += pio2_2;
slouken@2756
   129
                y[0] = z + pio2_2t;
slouken@2756
   130
                y[1] = (z - y[0]) + pio2_2t;
slouken@2756
   131
            }
slouken@2756
   132
            return -1;
slouken@2756
   133
        }
slouken@2756
   134
    }
slouken@2756
   135
    if (ix <= 0x413921fb) {     /* |x| ~<= 2^19*(pi/2), medium size */
slouken@2756
   136
        t = fabs(x);
slouken@2756
   137
        n = (int32_t) (t * invpio2 + half);
slouken@2756
   138
        fn = (double) n;
slouken@2756
   139
        r = t - fn * pio2_1;
slouken@2756
   140
        w = fn * pio2_1t;       /* 1st round good to 85 bit */
slouken@2756
   141
        if (n < 32 && ix != npio2_hw[n - 1]) {
slouken@2756
   142
            y[0] = r - w;       /* quick check no cancellation */
slouken@2756
   143
        } else {
slouken@2756
   144
            u_int32_t high;
slouken@2756
   145
            j = ix >> 20;
slouken@2756
   146
            y[0] = r - w;
slouken@2756
   147
            GET_HIGH_WORD(high, y[0]);
slouken@2756
   148
            i = j - ((high >> 20) & 0x7ff);
slouken@2756
   149
            if (i > 16) {       /* 2nd iteration needed, good to 118 */
slouken@2756
   150
                t = r;
slouken@2756
   151
                w = fn * pio2_2;
slouken@2756
   152
                r = t - w;
slouken@2756
   153
                w = fn * pio2_2t - ((t - r) - w);
slouken@2756
   154
                y[0] = r - w;
slouken@2756
   155
                GET_HIGH_WORD(high, y[0]);
slouken@2756
   156
                i = j - ((high >> 20) & 0x7ff);
slouken@2756
   157
                if (i > 49) {   /* 3rd iteration need, 151 bits acc */
slouken@2756
   158
                    t = r;      /* will cover all possible cases */
slouken@2756
   159
                    w = fn * pio2_3;
slouken@2756
   160
                    r = t - w;
slouken@2756
   161
                    w = fn * pio2_3t - ((t - r) - w);
slouken@2756
   162
                    y[0] = r - w;
slouken@2756
   163
                }
slouken@2756
   164
            }
slouken@2756
   165
        }
slouken@2756
   166
        y[1] = (r - y[0]) - w;
slouken@2756
   167
        if (hx < 0) {
slouken@2756
   168
            y[0] = -y[0];
slouken@2756
   169
            y[1] = -y[1];
slouken@2756
   170
            return -n;
slouken@2756
   171
        } else
slouken@2756
   172
            return n;
slouken@2756
   173
    }
slouken@2756
   174
    /*
slouken@2756
   175
     * all other (large) arguments
slouken@2756
   176
     */
slouken@2756
   177
    if (ix >= 0x7ff00000) {     /* x is inf or NaN */
slouken@2756
   178
        y[0] = y[1] = x - x;
slouken@2756
   179
        return 0;
slouken@2756
   180
    }
slouken@2756
   181
    /* set z = scalbn(|x|,ilogb(x)-23) */
slouken@2756
   182
    GET_LOW_WORD(low, x);
slouken@2756
   183
    SET_LOW_WORD(z, low);
slouken@2756
   184
    e0 = (ix >> 20) - 1046;     /* e0 = ilogb(z)-23; */
slouken@2756
   185
    SET_HIGH_WORD(z, ix - ((int32_t) (e0 << 20)));
slouken@2756
   186
    for (i = 0; i < 2; i++) {
slouken@2756
   187
        tx[i] = (double) ((int32_t) (z));
slouken@2756
   188
        z = (z - tx[i]) * two24;
slouken@2756
   189
    }
slouken@2756
   190
    tx[2] = z;
slouken@2756
   191
    nx = 3;
slouken@2756
   192
    while (tx[nx - 1] == zero)
slouken@2756
   193
        nx--;                   /* skip zero term */
slouken@2756
   194
    n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
slouken@2756
   195
    if (hx < 0) {
slouken@2756
   196
        y[0] = -y[0];
slouken@2756
   197
        y[1] = -y[1];
slouken@2756
   198
        return -n;
slouken@2756
   199
    }
slouken@2756
   200
    return n;
slouken@2756
   201
}