e_rem_pio2f.c (176410) | e_rem_pio2f.c (176465) |
---|---|
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#include <sys/cdefs.h> | 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#include <sys/cdefs.h> |
18__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 176410 2008-02-19 15:42:46Z bde $"); | 18__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 176465 2008-02-22 15:55:14Z bde $"); |
19 20/* __ieee754_rem_pio2f(x,y) 21 * 22 * return the remainder of x rem pi/2 in y[0]+y[1] 23 * use double precision internally 24 * use __kernel_rem_pio2() for large x 25 */ 26 | 19 20/* __ieee754_rem_pio2f(x,y) 21 * 22 * return the remainder of x rem pi/2 in y[0]+y[1] 23 * use double precision internally 24 * use __kernel_rem_pio2() for large x 25 */ 26 |
27#include <float.h> 28 |
|
27#include "math.h" 28#include "math_private.h" 29 30/* 31 * invpio2: 53 bits of 2/pi 32 * pio2_1: first 33 bit of pi/2 33 * pio2_1t: pi/2 - pio2_1 34 */ --- 12 unchanged lines hidden (view full) --- 47 float z; 48 int32_t e0,n,ix,hx; 49 50 GET_FLOAT_WORD(hx,x); 51 ix = hx&0x7fffffff; 52 /* 33+53 bit pi is good enough for medium size */ 53 if(ix<=0x49490f80) { /* |x| ~<= 2^19*(pi/2), medium size */ 54 t = fabsf(x); | 29#include "math.h" 30#include "math_private.h" 31 32/* 33 * invpio2: 53 bits of 2/pi 34 * pio2_1: first 33 bit of pi/2 35 * pio2_1t: pi/2 - pio2_1 36 */ --- 12 unchanged lines hidden (view full) --- 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#ifdef HAVE_EFFICIENT_IRINT 58 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 59 STRICT_ASSIGN(double,fn,t*invpio2+0x1.8p52); 60 fn = fn-0x1.8p52; 61 n = irint(fn); 62#else |
|
55 n = (int32_t) (t*invpio2+half); 56 fn = (double)n; | 63 n = (int32_t) (t*invpio2+half); 64 fn = (double)n; |
65#endif |
|
57 r = t-fn*pio2_1; 58 w = fn*pio2_1t; 59 y[0] = r-w; 60 y[1] = (r-y[0])-w; 61 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 62 else return n; 63 } 64 /* --- 15 unchanged lines hidden --- | 66 r = t-fn*pio2_1; 67 w = fn*pio2_1t; 68 y[0] = r-w; 69 y[1] = (r-y[0])-w; 70 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 71 else return n; 72 } 73 /* --- 15 unchanged lines hidden --- |