e_rem_pio2.c (176467) | e_rem_pio2.c (176476) |
---|---|
1 2/* @(#)e_rem_pio2.c 1.4 95/01/18 */ 3/* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 * | 1 2/* @(#)e_rem_pio2.c 1.4 95/01/18 */ 3/* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 * |
13 * Optimized by Bruce D. Evans. |
|
13 */ 14 15#include <sys/cdefs.h> | 14 */ 15 16#include <sys/cdefs.h> |
16__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2.c 176467 2008-02-22 18:43:23Z bde $"); | 17__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2.c 176476 2008-02-23 12:53:21Z bde $"); |
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 <float.h> --- 97 unchanged lines hidden (view full) --- 122 y[0] = z + 4*pio2_1t; 123 y[1] = (z-y[0])+4*pio2_1t; 124 return -4; 125 } 126 } 127 } 128 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 129medium: | 18 19/* __ieee754_rem_pio2(x,y) 20 * 21 * return the remainder of x rem pi/2 in y[0]+y[1] 22 * use __kernel_rem_pio2() 23 */ 24 25#include <float.h> --- 97 unchanged lines hidden (view full) --- 123 y[0] = z + 4*pio2_1t; 124 y[1] = (z-y[0])+4*pio2_1t; 125 return -4; 126 } 127 } 128 } 129 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 130medium: |
130 t = fabs(x); | |
131 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ | 131 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ |
132 STRICT_ASSIGN(double,fn,t*invpio2+0x1.8p52); | 132 STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); |
133 fn = fn-0x1.8p52; 134#ifdef HAVE_EFFICIENT_IRINT 135 n = irint(fn); 136#else 137 n = (int32_t)fn; 138#endif | 133 fn = fn-0x1.8p52; 134#ifdef HAVE_EFFICIENT_IRINT 135 n = irint(fn); 136#else 137 n = (int32_t)fn; 138#endif |
139 r = t-fn*pio2_1; | 139 r = x-fn*pio2_1; |
140 w = fn*pio2_1t; /* 1st round good to 85 bit */ 141 { 142 u_int32_t high; 143 j = ix>>20; 144 y[0] = r-w; 145 GET_HIGH_WORD(high,y[0]); 146 i = j-((high>>20)&0x7ff); 147 if(i>16) { /* 2nd iteration needed, good to 118 */ --- 9 unchanged lines hidden (view full) --- 157 w = fn*pio2_3; 158 r = t-w; 159 w = fn*pio2_3t-((t-r)-w); 160 y[0] = r-w; 161 } 162 } 163 } 164 y[1] = (r-y[0])-w; | 140 w = fn*pio2_1t; /* 1st round good to 85 bit */ 141 { 142 u_int32_t high; 143 j = ix>>20; 144 y[0] = r-w; 145 GET_HIGH_WORD(high,y[0]); 146 i = j-((high>>20)&0x7ff); 147 if(i>16) { /* 2nd iteration needed, good to 118 */ --- 9 unchanged lines hidden (view full) --- 157 w = fn*pio2_3; 158 r = t-w; 159 w = fn*pio2_3t-((t-r)-w); 160 y[0] = r-w; 161 } 162 } 163 } 164 y[1] = (r-y[0])-w; |
165 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 166 else return n; | 165 return n; |
167 } 168 /* 169 * all other (large) arguments 170 */ 171 if(ix>=0x7ff00000) { /* x is inf or NaN */ 172 y[0]=y[1]=x-x; return 0; 173 } 174 /* set z = scalbn(|x|,ilogb(x)-23) */ --- 15 unchanged lines hidden --- | 166 } 167 /* 168 * all other (large) arguments 169 */ 170 if(ix>=0x7ff00000) { /* x is inf or NaN */ 171 y[0]=y[1]=x-x; return 0; 172 } 173 /* set z = scalbn(|x|,ilogb(x)-23) */ --- 15 unchanged lines hidden --- |