src/libm/e_rem_pio2.c

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