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