sysdeps/ieee754/flt-32/e_rem_pio2f.c

   1 /* e_rem_pio2f.c -- float version of e_rem_pio2.c
   2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   3  */
   4
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #if defined(LIBM_SCCS) && !defined(lint)
  17 static char rcsid[] = "$NetBSD: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
  18 #endif
  19
  20 /* __ieee754_rem_pio2f(x,y)
  21  *
  22  * return the remainder of x rem pi/2 in y[0]+y[1]
  23  * use __kernel_rem_pio2f()
  24  */
  25
  26 #include "math.h"
  27 #include "math_private.h"
  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 0xA2, 0xF9, 0x83, 0x6E, 0x4E, 0x44, 0x15, 0x29, 0xFC,
  38 0x27, 0x57, 0xD1, 0xF5, 0x34, 0xDD, 0xC0, 0xDB, 0x62,
  39 0x95, 0x99, 0x3C, 0x43, 0x90, 0x41, 0xFE, 0x51, 0x63,
  40 0xAB, 0xDE, 0xBB, 0xC5, 0x61, 0xB7, 0x24, 0x6E, 0x3A,
  41 0x42, 0x4D, 0xD2, 0xE0, 0x06, 0x49, 0x2E, 0xEA, 0x09,
  42 0xD1, 0x92, 0x1C, 0xFE, 0x1D, 0xEB, 0x1C, 0xB1, 0x29,
  43 0xA7, 0x3E, 0xE8, 0x82, 0x35, 0xF5, 0x2E, 0xBB, 0x44,
  44 0x84, 0xE9, 0x9C, 0x70, 0x26, 0xB4, 0x5F, 0x7E, 0x41,
  45 0x39, 0x91, 0xD6, 0x39, 0x83, 0x53, 0x39, 0xF4, 0x9C,
  46 0x84, 0x5F, 0x8B, 0xBD, 0xF9, 0x28, 0x3B, 0x1F, 0xF8,
  47 0x97, 0xFF, 0xDE, 0x05, 0x98, 0x0F, 0xEF, 0x2F, 0x11,
  48 0x8B, 0x5A, 0x0A, 0x6D, 0x1F, 0x6D, 0x36, 0x7E, 0xCF,
  49 0x27, 0xCB, 0x09, 0xB7, 0x4F, 0x46, 0x3F, 0x66, 0x9E,
  50 0x5F, 0xEA, 0x2D, 0x75, 0x27, 0xBA, 0xC7, 0xEB, 0xE5,
  51 0xF1, 0x7B, 0x3D, 0x07, 0x39, 0xF7, 0x8A, 0x52, 0x92,
  52 0xEA, 0x6B, 0xFB, 0x5F, 0xB1, 0x1F, 0x8D, 0x5D, 0x08,
  53 0x56, 0x03, 0x30, 0x46, 0xFC, 0x7B, 0x6B, 0xAB, 0xF0,
  54 0xCF, 0xBC, 0x20, 0x9A, 0xF4, 0x36, 0x1D, 0xA9, 0xE3,
  55 0x91, 0x61, 0x5E, 0xE6, 0x1B, 0x08, 0x65, 0x99, 0x85,
  56 0x5F, 0x14, 0xA0, 0x68, 0x40, 0x8D, 0xFF, 0xD8, 0x80,
  57 0x4D, 0x73, 0x27, 0x31, 0x06, 0x06, 0x15, 0x56, 0xCA,
  58 0x73, 0xA8, 0xC9, 0x60, 0xE2, 0x7B, 0xC0, 0x8C, 0x6B,
  59 };
  60
  61 /* This array is like the one in e_rem_pio2.c, but the numbers are
  62    single precision and the last 8 bits are forced to 0.  */
  63 #ifdef __STDC__
  64 static const int32_t npio2_hw[] = {
  65 #else
  66 static int32_t npio2_hw[] = {
  67 #endif
  68 0x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
  69 0x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
  70 0x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
  71 0x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
  72 0x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
  73 0x4242c700, 0x42490f00
  74 };
  75
  76 /*
  77  * invpio2:  24 bits of 2/pi
  78  * pio2_1:   first  17 bit of pi/2
  79  * pio2_1t:  pi/2 - pio2_1
  80  * pio2_2:   second 17 bit of pi/2
  81  * pio2_2t:  pi/2 - (pio2_1+pio2_2)
  82  * pio2_3:   third  17 bit of pi/2
  83  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  84  */
  85
  86 #ifdef __STDC__
  87 static const float
  88 #else
  89 static float
  90 #endif
  91 zero =  0.0000000000e+00, /* 0x00000000 */
  92 half =  5.0000000000e-01, /* 0x3f000000 */
  93 two8 =  2.5600000000e+02, /* 0x43800000 */
  94 invpio2 =  6.3661980629e-01, /* 0x3f22f984 */
  95 pio2_1  =  1.5707855225e+00, /* 0x3fc90f80 */
  96 pio2_1t =  1.0804334124e-05, /* 0x37354443 */
  97 pio2_2  =  1.0804273188e-05, /* 0x37354400 */
  98 pio2_2t =  6.0770999344e-11, /* 0x2e85a308 */
  99 pio2_3  =  6.0770943833e-11, /* 0x2e85a300 */
 100 pio2_3t =  6.1232342629e-17; /* 0x248d3132 */
 101
 102 #ifdef __STDC__
 103         int32_t __ieee754_rem_pio2f(float x, float *y)
 104 #else
 105         int32_t __ieee754_rem_pio2f(x,y)
 106         float x,y[];
 107 #endif
 108 {
 109         float z,w,t,r,fn;
 110         float tx[3];
 111         int32_t e0,i,j,nx,n,ix,hx;
 112
 113         GET_FLOAT_WORD(hx,x);
 114         ix = hx&0x7fffffff;
 115         if(ix<=0x3f490fd8)   /* |x| ~<= pi/4 , no need for reduction */
 116             {y[0] = x; y[1] = 0; return 0;}
 117         if(ix<0x4016cbe4) {  /* |x| < 3pi/4, special case with n=+-1 */
 118             if(hx>0) {
 119                 z = x - pio2_1;
 120                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
 121                     y[0] = z - pio2_1t;
 122                     y[1] = (z-y[0])-pio2_1t;
 123                 } else {                /* near pi/2, use 24+24+24 bit pi */
 124                     z -= pio2_2;
 125                     y[0] = z - pio2_2t;
 126                     y[1] = (z-y[0])-pio2_2t;
 127                 }
 128                 return 1;
 129             } else {    /* negative x */
 130                 z = x + pio2_1;
 131                 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
 132                     y[0] = z + pio2_1t;
 133                     y[1] = (z-y[0])+pio2_1t;
 134                 } else {                /* near pi/2, use 24+24+24 bit pi */
 135                     z += pio2_2;
 136                     y[0] = z + pio2_2t;
 137                     y[1] = (z-y[0])+pio2_2t;
 138                 }
 139                 return -1;
 140             }
 141         }
 142         if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
 143             t  = fabsf(x);
 144             n  = (int32_t) (t*invpio2+half);
 145             fn = (float)n;
 146             r  = t-fn*pio2_1;
 147             w  = fn*pio2_1t;    /* 1st round good to 40 bit */
 148             if(n<32&&(int32_t)(ix&0xffffff00)!=npio2_hw[n-1]) {
 149                 y[0] = r-w;     /* quick check no cancellation */
 150             } else {
 151                 u_int32_t high;
 152                 j  = ix>>23;
 153                 y[0] = r-w;
 154                 GET_FLOAT_WORD(high,y[0]);
 155                 i = j-((high>>23)&0xff);
 156                 if(i>8) {  /* 2nd iteration needed, good to 57 */
 157                     t  = r;
 158                     w  = fn*pio2_2;
 159                     r  = t-w;
 160                     w  = fn*pio2_2t-((t-r)-w);
 161                     y[0] = r-w;
 162                     GET_FLOAT_WORD(high,y[0]);
 163                     i = j-((high>>23)&0xff);
 164                     if(i>25)  { /* 3rd iteration need, 74 bits acc */
 165                         t  = r; /* will cover all possible cases */
 166                         w  = fn*pio2_3;
 167                         r  = t-w;
 168                         w  = fn*pio2_3t-((t-r)-w);
 169                         y[0] = r-w;
 170                     }
 171                 }
 172             }
 173             y[1] = (r-y[0])-w;
 174             if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 175             else         return n;
 176         }
 177     /*
 178      * all other (large) arguments
 179      */
 180         if(ix>=0x7f800000) {            /* x is inf or NaN */
 181             y[0]=y[1]=x-x; return 0;
 182         }
 183     /* set z = scalbn(|x|,ilogb(x)-7) */
 184         e0      = (ix>>23)-134;         /* e0 = ilogb(z)-7; */
 185         SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
 186         for(i=0;i<2;i++) {
 187                 tx[i] = (float)((int32_t)(z));
 188                 z     = (z-tx[i])*two8;
 189         }
 190         tx[2] = z;
 191         nx = 3;
 192         while(tx[nx-1]==zero) nx--;     /* skip zero term */
 193         n  =  __kernel_rem_pio2f(tx,y,e0,nx,2,two_over_pi);
 194         if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 195         return n;
 196 }