e_rem_pio2f.c revision 176356
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 * Debugged and optimized by Bruce D. Evans. 4 */ 5 6/* 7 * ==================================================== 8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 9 * 10 * Developed at SunPro, a Sun Microsystems, Inc. business. 11 * Permission to use, copy, modify, and distribute this 12 * software is freely granted, provided that this notice 13 * is preserved. 14 * ==================================================== 15 */ 16 17#ifndef lint 18static char rcsid[] = "$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 176356 2008-02-17 07:31:59Z das $"; 19#endif 20 21/* __ieee754_rem_pio2f(x,y) 22 * 23 * return the remainder of x rem pi/2 in y[0]+y[1] 24 * use double precision internally 25 * use __kernel_rem_pio2() for large x 26 */ 27 28#include "math.h" 29#include "math_private.h" 30 31/* 32 * invpio2: 53 bits of 2/pi 33 * pio2_1: first 33 bit of pi/2 34 * pio2_1t: pi/2 - pio2_1 35 */ 36 37static const double 38zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 39half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 40two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 41invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 42pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 43pio2_1t = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */ 44 45 int32_t __ieee754_rem_pio2f(float x, float *y) 46{ 47 double w,t,r,fn; 48 double tx[1],ty[2]; 49 float z; 50 int32_t e0,n,ix,hx; 51 52 GET_FLOAT_WORD(hx,x); 53 ix = hx&0x7fffffff; 54 /* 33+53 bit pi is good enough for medium size */ 55 if(ix<=0x49490f80) { /* |x| ~<= 2^19*(pi/2), medium size */ 56 t = fabsf(x); 57 n = (int32_t) (t*invpio2+half); 58 fn = (double)n; 59 r = t-fn*pio2_1; 60 w = fn*pio2_1t; 61 y[0] = r-w; 62 y[1] = (r-y[0])-w; 63 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 64 else return n; 65 } 66 /* 67 * all other (large) arguments 68 */ 69 if(ix>=0x7f800000) { /* x is inf or NaN */ 70 y[0]=y[1]=x-x; return 0; 71 } 72 /* set z = scalbn(|x|,ilogb(|x|)-23) */ 73 e0 = (ix>>23)-150; /* e0 = ilogb(|x|)-23; */ 74 SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23))); 75 tx[0] = z; 76 n = __kernel_rem_pio2(tx,ty,e0,1,1); 77 y[0] = ty[0]; 78 y[1] = ty[0] - y[0]; 79 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 80 return n; 81} 82