From owner-svn-src-head@FreeBSD.ORG Fri Apr 29 23:13:43 2011 Return-Path: Delivered-To: svn-src-head@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id BD2061065677; Fri, 29 Apr 2011 23:13:43 +0000 (UTC) (envelope-from kargl@FreeBSD.org) Received: from svn.freebsd.org (svn.freebsd.org [IPv6:2001:4f8:fff6::2c]) by mx1.freebsd.org (Postfix) with ESMTP id A8EA08FC17; Fri, 29 Apr 2011 23:13:43 +0000 (UTC) Received: from svn.freebsd.org (localhost [127.0.0.1]) by svn.freebsd.org (8.14.4/8.14.4) with ESMTP id p3TNDhVf019726; Fri, 29 Apr 2011 23:13:43 GMT (envelope-from kargl@svn.freebsd.org) Received: (from kargl@localhost) by svn.freebsd.org (8.14.4/8.14.4/Submit) id p3TNDhD6019720; Fri, 29 Apr 2011 23:13:43 GMT (envelope-from kargl@svn.freebsd.org) Message-Id: <201104292313.p3TNDhD6019720@svn.freebsd.org> From: Steve Kargl Date: Fri, 29 Apr 2011 23:13:43 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org X-SVN-Group: head MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cc: Subject: svn commit: r221234 - in head/lib/msun: ld128 ld80 src X-BeenThere: svn-src-head@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: SVN commit messages for the src tree for head/-current List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 29 Apr 2011 23:13:43 -0000 Author: kargl Date: Fri Apr 29 23:13:43 2011 New Revision: 221234 URL: http://svn.freebsd.org/changeset/base/221234 Log: Improve the accuracy from a max ULP of ~2000 to max ULP < 0.79 on i386-class hardware for sinl and cosl. The hand-rolled argument reduction have been replaced by e_rem_pio2l() implementations. To preserve history the following commands have been executed: svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h The ld80 version has been tested by bde, das, and kargl over the last few years (bde, das) and few months (kargl). An older ld128 version was tested by das. The committed version has only been compiled tested via 'make universe'. Approved by: das (mentor) Obtained from: bde Added: head/lib/msun/ld128/e_rem_pio2l.h - copied, changed from r221219, head/lib/msun/src/e_rem_pio2.c head/lib/msun/ld80/e_rem_pio2l.h - copied, changed from r221219, head/lib/msun/src/e_rem_pio2.c Modified: head/lib/msun/src/s_cosl.c head/lib/msun/src/s_sinl.c head/lib/msun/src/s_tanl.c Copied and modified: head/lib/msun/ld128/e_rem_pio2l.h (from r221219, head/lib/msun/src/e_rem_pio2.c) ============================================================================== --- head/lib/msun/src/e_rem_pio2.c Fri Apr 29 18:49:36 2011 (r221219, copy source) +++ head/lib/msun/ld128/e_rem_pio2l.h Fri Apr 29 23:13:43 2011 (r221234) @@ -1,8 +1,8 @@ - -/* @(#)e_rem_pio2.c 1.4 95/01/18 */ +/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this @@ -16,7 +16,7 @@ #include __FBSDID("$FreeBSD$"); -/* __ieee754_rem_pio2(x,y) +/* ld128 version of __ieee754_rem_pio2l(x,y) * * return the remainder of x rem pi/2 in y[0]+y[1] * use __kernel_rem_pio2() @@ -26,132 +26,85 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +#include "fpmath.h" + +#define BIAS (LDBL_MAX_EXP - 1) /* - * invpio2: 53 bits of 2/pi - * pio2_1: first 33 bit of pi/2 + * XXX need to verify that nonzero integer multiples of pi/2 within the + * range get no closer to a long double than 2**-140, or that + * ilogb(x) + ilogb(min_delta) < 45 - -140. + */ +/* + * invpio2: 113 bits of 2/pi + * pio2_1: first 68 bits of pi/2 * pio2_1t: pi/2 - pio2_1 - * pio2_2: second 33 bit of pi/2 + * pio2_2: second 68 bits of pi/2 * pio2_2t: pi/2 - (pio2_1+pio2_2) - * pio2_3: third 33 bit of pi/2 + * pio2_3: third 68 bits of pi/2 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const double zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ -two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ -invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ -pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ -pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ -pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ -pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ -pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ -pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ +two24 = 1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */ -#ifndef INLINE_REM_PIO2 -extern +static const long double +invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */ +pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */ +pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */ +pio2_2 = 2.0222662487959507323994779168837751e-21L, /* 0x13198a2e03707344a400000000000.0p-181 */ +pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df531d89cd91.0p-254 */ +pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */ +pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */ + +#ifdef INLINE_REM_PIO2L +static inline __always_inline #endif -__inline int -__ieee754_rem_pio2(double x, double *y) +int +__ieee754_rem_pio2l(long double x, long double *y) { - double z,w,t,r,fn; - double tx[3],ty[2]; - int32_t e0,i,j,nx,n,ix,hx; - u_int32_t low; - - GET_HIGH_WORD(hx,x); /* high word of x */ - ix = hx&0x7fffffff; -#if 0 /* Must be handled in caller. */ - if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ - {y[0] = x; y[1] = 0; return 0;} -#endif - if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ - if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ - goto medium; /* cancellation -- use medium case */ - if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ - if (hx > 0) { - z = x - pio2_1; /* one round good to 85 bits */ - y[0] = z - pio2_1t; - y[1] = (z-y[0])-pio2_1t; - return 1; - } else { - z = x + pio2_1; - y[0] = z + pio2_1t; - y[1] = (z-y[0])+pio2_1t; - return -1; - } - } else { - if (hx > 0) { - z = x - 2*pio2_1; - y[0] = z - 2*pio2_1t; - y[1] = (z-y[0])-2*pio2_1t; - return 2; - } else { - z = x + 2*pio2_1; - y[0] = z + 2*pio2_1t; - y[1] = (z-y[0])+2*pio2_1t; - return -2; - } - } - } - if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ - if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ - if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ - goto medium; - if (hx > 0) { - z = x - 3*pio2_1; - y[0] = z - 3*pio2_1t; - y[1] = (z-y[0])-3*pio2_1t; - return 3; - } else { - z = x + 3*pio2_1; - y[0] = z + 3*pio2_1t; - y[1] = (z-y[0])+3*pio2_1t; - return -3; - } - } else { - if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ - goto medium; - if (hx > 0) { - z = x - 4*pio2_1; - y[0] = z - 4*pio2_1t; - y[1] = (z-y[0])-4*pio2_1t; - return 4; - } else { - z = x + 4*pio2_1; - y[0] = z + 4*pio2_1t; - y[1] = (z-y[0])+4*pio2_1t; - return -4; - } - } - } - if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ -medium: + union IEEEl2bits u,u1; + long double z,w,t,r,fn; + double tx[5],ty[3]; + int64_t n; + int e0,ex,i,j,nx; + int16_t expsign; + + u.e = x; + expsign = u.xbits.expsign; + ex = expsign & 0x7fff; + if (ex < BIAS + 45 || ex == BIAS + 45 && + u.bits.manh < 0x921fb54442d1LL) { + /* |x| ~< 2^45*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); - fn = fn-0x1.8p52; -#ifdef HAVE_EFFICIENT_IRINT - n = irint(fn); + fn = x*invpio2+0x1.8p112; + fn = fn-0x1.8p112; +#ifdef HAVE_EFFICIENT_I64RINT + n = i64rint(fn); #else - n = (int32_t)fn; + n = fn; #endif r = x-fn*pio2_1; - w = fn*pio2_1t; /* 1st round good to 85 bit */ + w = fn*pio2_1t; /* 1st round good to 180 bit */ { - u_int32_t high; - j = ix>>20; + union IEEEl2bits u2; + int ex1; + j = ex; y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>16) { /* 2nd iteration needed, good to 118 */ + u2.e = y[0]; + ex1 = u2.xbits.expsign & 0x7fff; + i = j-ex1; + if(i>51) { /* 2nd iteration needed, good to 248 */ t = r; w = fn*pio2_2; r = t-w; w = fn*pio2_2t-((t-r)-w); y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>49) { /* 3rd iteration need, 151 bits acc */ + u2.e = y[0]; + ex1 = u2.xbits.expsign & 0x7fff; + i = j-ex1; + if(i>119) { /* 3rd iteration need, 316 bits acc */ t = r; /* will cover all possible cases */ w = fn*pio2_3; r = t-w; @@ -166,22 +119,25 @@ medium: /* * all other (large) arguments */ - if(ix>=0x7ff00000) { /* x is inf or NaN */ + if(ex==0x7fff) { /* x is inf or NaN */ y[0]=y[1]=x-x; return 0; } /* set z = scalbn(|x|,ilogb(x)-23) */ - GET_LOW_WORD(low,x); - SET_LOW_WORD(z,low); - e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ - SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); - for(i=0;i<2;i++) { + u1.e = x; + e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */ + u1.xbits.expsign = ex - e0; + z = u1.e; + for(i=0;i<4;i++) { tx[i] = (double)((int32_t)(z)); z = (z-tx[i])*two24; } - tx[2] = z; - nx = 3; + tx[4] = z; + nx = 5; while(tx[nx-1]==zero) nx--; /* skip zero term */ - n = __kernel_rem_pio2(tx,ty,e0,nx,1); - if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} - y[0] = ty[0]; y[1] = ty[1]; return n; + n = __kernel_rem_pio2(tx,ty,e0,nx,3); + t = (long double)ty[2] + ty[1]; + r = t + ty[0]; + w = ty[0] - (r - t); + if(expsign<0) {y[0] = -r; y[1] = -w; return -n;} + y[0] = r; y[1] = w; return n; } Copied and modified: head/lib/msun/ld80/e_rem_pio2l.h (from r221219, head/lib/msun/src/e_rem_pio2.c) ============================================================================== --- head/lib/msun/src/e_rem_pio2.c Fri Apr 29 18:49:36 2011 (r221219, copy source) +++ head/lib/msun/ld80/e_rem_pio2l.h Fri Apr 29 23:13:43 2011 (r221234) @@ -1,8 +1,8 @@ - -/* @(#)e_rem_pio2.c 1.4 95/01/18 */ +/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this @@ -16,7 +16,7 @@ #include __FBSDID("$FreeBSD$"); -/* __ieee754_rem_pio2(x,y) +/* ld80 version of __ieee754_rem_pio2l(x,y) * * return the remainder of x rem pi/2 in y[0]+y[1] * use __kernel_rem_pio2() @@ -26,132 +26,95 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +#include "fpmath.h" + +#define BIAS (LDBL_MAX_EXP - 1) /* - * invpio2: 53 bits of 2/pi - * pio2_1: first 33 bit of pi/2 + * invpio2: 64 bits of 2/pi + * pio2_1: first 39 bits of pi/2 * pio2_1t: pi/2 - pio2_1 - * pio2_2: second 33 bit of pi/2 + * pio2_2: second 39 bits of pi/2 * pio2_2t: pi/2 - (pio2_1+pio2_2) - * pio2_3: third 33 bit of pi/2 + * pio2_3: third 39 bits of pi/2 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ static const double zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ -invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ -pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ -pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ -pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ -pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ -pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ -pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ +pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */ +pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */ +pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */ + +#if defined(__amd64__) || defined(__i386__) +/* Long double constants are slow on these arches, and broken on i386. */ +static const volatile double +invpio2hi = 6.3661977236758138e-01, /* 0x145f306dc9c883.0p-53 */ +invpio2lo = -3.9356538861223811e-17, /* -0x16b00000000000.0p-107 */ +pio2_1thi = -1.0746346554971943e-12, /* -0x12e7b9676733af.0p-92 */ +pio2_1tlo = 8.8451028997905949e-29, /* 0x1c080000000000.0p-146 */ +pio2_2thi = 6.3683171635109499e-25, /* 0x18a2e03707344a.0p-133 */ +pio2_2tlo = 2.3183081793789774e-41, /* 0x10280000000000.0p-187 */ +pio2_3thi = -2.7529965190440717e-37, /* -0x176b7ed8fbbacc.0p-174 */ +pio2_3tlo = -4.2006647512740502e-54; /* -0x19c00000000000.0p-230 */ +#define invpio2 ((long double)invpio2hi + invpio2lo) +#define pio2_1t ((long double)pio2_1thi + pio2_1tlo) +#define pio2_2t ((long double)pio2_2thi + pio2_2tlo) +#define pio2_3t ((long double)pio2_3thi + pio2_3tlo) +#else +static const long double +invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */ +pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */ +pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */ +pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */ +#endif -#ifndef INLINE_REM_PIO2 -extern +#ifdef INLINE_REM_PIO2L +static inline __always_inline #endif -__inline int -__ieee754_rem_pio2(double x, double *y) +int +__ieee754_rem_pio2l(long double x, long double *y) { - double z,w,t,r,fn; + union IEEEl2bits u,u1; + long double z,w,t,r,fn; double tx[3],ty[2]; - int32_t e0,i,j,nx,n,ix,hx; - u_int32_t low; + int e0,ex,i,j,nx,n; + int16_t expsign; - GET_HIGH_WORD(hx,x); /* high word of x */ - ix = hx&0x7fffffff; -#if 0 /* Must be handled in caller. */ - if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ - {y[0] = x; y[1] = 0; return 0;} -#endif - if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */ - if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */ - goto medium; /* cancellation -- use medium case */ - if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */ - if (hx > 0) { - z = x - pio2_1; /* one round good to 85 bits */ - y[0] = z - pio2_1t; - y[1] = (z-y[0])-pio2_1t; - return 1; - } else { - z = x + pio2_1; - y[0] = z + pio2_1t; - y[1] = (z-y[0])+pio2_1t; - return -1; - } - } else { - if (hx > 0) { - z = x - 2*pio2_1; - y[0] = z - 2*pio2_1t; - y[1] = (z-y[0])-2*pio2_1t; - return 2; - } else { - z = x + 2*pio2_1; - y[0] = z + 2*pio2_1t; - y[1] = (z-y[0])+2*pio2_1t; - return -2; - } - } - } - if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */ - if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */ - if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */ - goto medium; - if (hx > 0) { - z = x - 3*pio2_1; - y[0] = z - 3*pio2_1t; - y[1] = (z-y[0])-3*pio2_1t; - return 3; - } else { - z = x + 3*pio2_1; - y[0] = z + 3*pio2_1t; - y[1] = (z-y[0])+3*pio2_1t; - return -3; - } - } else { - if (ix == 0x401921fb) /* |x| ~= 4pi/2 */ - goto medium; - if (hx > 0) { - z = x - 4*pio2_1; - y[0] = z - 4*pio2_1t; - y[1] = (z-y[0])-4*pio2_1t; - return 4; - } else { - z = x + 4*pio2_1; - y[0] = z + 4*pio2_1t; - y[1] = (z-y[0])+4*pio2_1t; - return -4; - } - } - } - if(ix<0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ -medium: + u.e = x; + expsign = u.xbits.expsign; + ex = expsign & 0x7fff; + if (ex < BIAS + 25 || ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2) { + /* |x| ~< 2^25*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52); - fn = fn-0x1.8p52; + fn = x*invpio2+0x1.8p63; + fn = fn-0x1.8p63; #ifdef HAVE_EFFICIENT_IRINT n = irint(fn); #else - n = (int32_t)fn; + n = fn; #endif r = x-fn*pio2_1; - w = fn*pio2_1t; /* 1st round good to 85 bit */ + w = fn*pio2_1t; /* 1st round good to 102 bit */ { - u_int32_t high; - j = ix>>20; + union IEEEl2bits u2; + int ex1; + j = ex; y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>16) { /* 2nd iteration needed, good to 118 */ + u2.e = y[0]; + ex1 = u2.xbits.expsign & 0x7fff; + i = j-ex1; + if(i>22) { /* 2nd iteration needed, good to 141 */ t = r; w = fn*pio2_2; r = t-w; w = fn*pio2_2t-((t-r)-w); y[0] = r-w; - GET_HIGH_WORD(high,y[0]); - i = j-((high>>20)&0x7ff); - if(i>49) { /* 3rd iteration need, 151 bits acc */ + u2.e = y[0]; + ex1 = u2.xbits.expsign & 0x7fff; + i = j-ex1; + if(i>61) { /* 3rd iteration need, 180 bits acc */ t = r; /* will cover all possible cases */ w = fn*pio2_3; r = t-w; @@ -166,14 +129,14 @@ medium: /* * all other (large) arguments */ - if(ix>=0x7ff00000) { /* x is inf or NaN */ + if(ex==0x7fff) { /* x is inf or NaN */ y[0]=y[1]=x-x; return 0; } /* set z = scalbn(|x|,ilogb(x)-23) */ - GET_LOW_WORD(low,x); - SET_LOW_WORD(z,low); - e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ - SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); + u1.e = x; + e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */ + u1.xbits.expsign = ex - e0; + z = u1.e; for(i=0;i<2;i++) { tx[i] = (double)((int32_t)(z)); z = (z-tx[i])*two24; @@ -181,7 +144,9 @@ medium: tx[2] = z; nx = 3; while(tx[nx-1]==zero) nx--; /* skip zero term */ - n = __kernel_rem_pio2(tx,ty,e0,nx,1); - if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;} - y[0] = ty[0]; y[1] = ty[1]; return n; + n = __kernel_rem_pio2(tx,ty,e0,nx,2); + r = (long double)ty[0] + ty[1]; + w = ty[1] - (r - ty[0]); + if(expsign<0) {y[0] = -r; y[1] = -w; return -n;} + y[0] = r; y[1] = w; return n; } Modified: head/lib/msun/src/s_cosl.c ============================================================================== --- head/lib/msun/src/s_cosl.c Fri Apr 29 23:00:23 2011 (r221233) +++ head/lib/msun/src/s_cosl.c Fri Apr 29 23:13:43 2011 (r221234) @@ -28,7 +28,6 @@ __FBSDID("$FreeBSD$"); /* - * Compute cos(x) for x where x is reduced to y = x - k * pi / 2. * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows * an accuracy of <= 0.7412 ULP. */ @@ -36,27 +35,22 @@ __FBSDID("$FreeBSD$"); #include #include "math.h" +#define INLINE_REM_PIO2L #include "math_private.h" -#include "fpmath.h" - #if LDBL_MANT_DIG == 64 -#define NX 3 -#define PREC 2 +#include "../ld80/e_rem_pio2l.h" #elif LDBL_MANT_DIG == 113 -#define NX 5 -#define PREC 3 +#include "../ld128/e_rem_pio2l.h" #else #error "Unsupported long double format" #endif -static const long double two24 = 1.67772160000000000000e+07L; - long double cosl(long double x) { union IEEEl2bits z; - int i, e0; - double xd[NX], yd[PREC]; + int e0; + long double y[2]; long double hi, lo; z.e = x; @@ -74,26 +68,9 @@ cosl(long double x) if (z.e < M_PI_4) return (__kernel_cosl(z.e, 0)); - /* Split z.e into a 24-bit representation. */ - e0 = ilogbl(z.e) - 23; - z.e = scalbnl(z.e, -e0); - for (i = 0; i < NX; i++) { - xd[i] = (double)((int32_t)z.e); - z.e = (z.e - xd[i]) * two24; - } - - /* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */ - e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC); - -#if PREC == 2 - hi = (long double)yd[0] + yd[1]; - lo = yd[1] - (hi - yd[0]); -#else /* PREC == 3 */ - long double t; - t = (long double)yd[2] + yd[1]; - hi = t + yd[0]; - lo = yd[0] - (hi - t); -#endif + e0 = __ieee754_rem_pio2l(x, y); + hi = y[0]; + lo = y[1]; switch (e0 & 3) { case 0: Modified: head/lib/msun/src/s_sinl.c ============================================================================== --- head/lib/msun/src/s_sinl.c Fri Apr 29 23:00:23 2011 (r221233) +++ head/lib/msun/src/s_sinl.c Fri Apr 29 23:13:43 2011 (r221234) @@ -27,34 +27,25 @@ #include __FBSDID("$FreeBSD$"); -/* - * Compute sin(x) for x where x is reduced to y = x - k * pi / 2. - */ - #include #include "math.h" +#define INLINE_REM_PIO2L #include "math_private.h" -#include "fpmath.h" - #if LDBL_MANT_DIG == 64 -#define NX 3 -#define PREC 2 +#include "../ld80/e_rem_pio2l.h" #elif LDBL_MANT_DIG == 113 -#define NX 5 -#define PREC 3 +#include "../ld128/e_rem_pio2l.h" #else #error "Unsupported long double format" #endif -static const long double two24 = 1.67772160000000000000e+07L; - long double sinl(long double x) { union IEEEl2bits z; - int i, e0, s; - double xd[NX], yd[PREC]; + int e0, s; + long double y[2]; long double hi, lo; z.e = x; @@ -75,26 +66,9 @@ sinl(long double x) return (s ? -hi : hi); } - /* Split z.e into a 24-bit representation. */ - e0 = ilogbl(z.e) - 23; - z.e = scalbnl(z.e, -e0); - for (i = 0; i < NX; i++) { - xd[i] = (double)((int32_t)z.e); - z.e = (z.e - xd[i]) * two24; - } - - /* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */ - e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC); - -#if PREC == 2 - hi = (long double)yd[0] + yd[1]; - lo = yd[1] - (hi - yd[0]); -#else /* PREC == 3 */ - long double t; - t = (long double)yd[2] + yd[1]; - hi = t + yd[0]; - lo = yd[0] - (hi - t); -#endif + e0 = __ieee754_rem_pio2l(x, y); + hi = y[0]; + lo = y[1]; switch (e0 & 3) { case 0: @@ -111,5 +85,5 @@ sinl(long double x) break; } - return (s ? -hi : hi); + return (hi); } Modified: head/lib/msun/src/s_tanl.c ============================================================================== --- head/lib/msun/src/s_tanl.c Fri Apr 29 23:00:23 2011 (r221233) +++ head/lib/msun/src/s_tanl.c Fri Apr 29 23:13:43 2011 (r221234) @@ -28,7 +28,6 @@ __FBSDID("$FreeBSD$"); /* - * Compute tan(x) for x where x is reduced to y = x - k * pi / 2. * Limited testing on pseudorandom numbers drawn within [0:4e8] shows * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%). @@ -37,27 +36,22 @@ __FBSDID("$FreeBSD$"); #include #include "math.h" +#define INLINE_REM_PIO2L #include "math_private.h" -#include "fpmath.h" - #if LDBL_MANT_DIG == 64 -#define NX 3 -#define PREC 2 +#include "../ld80/e_rem_pio2l.h" #elif LDBL_MANT_DIG == 113 -#define NX 5 -#define PREC 3 +#include "../ld128/e_rem_pio2l.h" #else #error "Unsupported long double format" #endif -static const long double two24 = 1.67772160000000000000e+07L; - long double tanl(long double x) { union IEEEl2bits z; - int i, e0, s; - double xd[NX], yd[PREC]; + int e0, s; + long double y[2]; long double hi, lo; z.e = x; @@ -78,26 +72,9 @@ tanl(long double x) return (s ? -hi : hi); } - /* Split z.e into a 24-bit representation. */ - e0 = ilogbl(z.e) - 23; - z.e = scalbnl(z.e, -e0); - for (i = 0; i < NX; i++) { - xd[i] = (double)((int32_t)z.e); - z.e = (z.e - xd[i]) * two24; - } - - /* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */ - e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC); - -#if PREC == 2 - hi = (long double)yd[0] + yd[1]; - lo = yd[1] - (hi - yd[0]); -#else /* PREC == 3 */ - long double t; - t = (long double)yd[2] + yd[1]; - hi = t + yd[0]; - lo = yd[0] - (hi - t); -#endif + e0 = __ieee754_rem_pio2l(x, y); + hi = y[0]; + lo = y[1]; switch (e0 & 3) { case 0: @@ -110,5 +87,5 @@ tanl(long double x) break; } - return (s ? -hi : hi); + return (hi); }