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