From owner-svn-src-head@FreeBSD.ORG Sun Aug 31 21:38:04 2014 Return-Path: Delivered-To: svn-src-head@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id A9821322; Sun, 31 Aug 2014 21:38:04 +0000 (UTC) Received: from svn.freebsd.org (svn.freebsd.org [IPv6:2001:1900:2254:2068::e6a:0]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client did not present a certificate) by mx1.freebsd.org (Postfix) with ESMTPS id 88EDA1CCE; Sun, 31 Aug 2014 21:38:04 +0000 (UTC) Received: from svn.freebsd.org ([127.0.1.70]) by svn.freebsd.org (8.14.9/8.14.9) with ESMTP id s7VLc4GB082704; Sun, 31 Aug 2014 21:38:04 GMT (envelope-from kargl@FreeBSD.org) Received: (from kargl@localhost) by svn.freebsd.org (8.14.9/8.14.9/Submit) id s7VLc4MA082700; Sun, 31 Aug 2014 21:38:04 GMT (envelope-from kargl@FreeBSD.org) Message-Id: <201408312138.s7VLc4MA082700@svn.freebsd.org> X-Authentication-Warning: svn.freebsd.org: kargl set sender to kargl@FreeBSD.org using -f From: Steve Kargl Date: Sun, 31 Aug 2014 21:38:04 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org Subject: svn commit: r270893 - head/lib/msun/src X-SVN-Group: head MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit X-BeenThere: svn-src-head@freebsd.org X-Mailman-Version: 2.1.18-1 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: Sun, 31 Aug 2014 21:38:04 -0000 Author: kargl Date: Sun Aug 31 21:38:03 2014 New Revision: 270893 URL: http://svnweb.freebsd.org/changeset/base/270893 Log: Compute sin(pi*x) without actually doing the pi*x multiplication. sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is the precision of x. The new argument reduction is an optimization compared to the old code, and it removes a chunk of dead code. Accuracy tests in the intervals (-21,-20), (-20,-19), ... (-1,0) show no differences between the old and new code. Obtained from: bde Modified: head/lib/msun/src/e_lgamma_r.c head/lib/msun/src/e_lgammaf_r.c Modified: head/lib/msun/src/e_lgamma_r.c ============================================================================== --- head/lib/msun/src/e_lgamma_r.c Sun Aug 31 21:18:23 2014 (r270892) +++ head/lib/msun/src/e_lgamma_r.c Sun Aug 31 21:38:03 2014 (r270893) @@ -156,37 +156,35 @@ w6 = -1.63092934096575273989e-03; /* 0x static const double zero= 0.00000000000000000000e+00; - static double sin_pi(double x) +/* + * Compute sin(pi*x) without actually doing the pi*x multiplication. + * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is + * the precision of x. + */ +static double +sin_pi(double x) { + volatile double vz; double y,z; - int n,ix; + int n; - GET_HIGH_WORD(ix,x); - ix &= 0x7fffffff; + y = -x; - if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0); - y = -x; /* x is assume negative */ + vz = y+0x1p52; /* depend on 0 <= y < 0x1p52 */ + z = vz-0x1p52; /* rint(y) for the above range */ + if (z == y) + return (zero); + + vz = y+0x1p50; + GET_LOW_WORD(n,vz); /* bits for rounded y (units 0.25) */ + z = vz-0x1p50; /* y rounded to a multiple of 0.25 */ + if (z > y) { + z -= 0.25; /* adjust to round down */ + n--; + } + n &= 7; /* octant of y mod 2 */ + y = y - z + n * 0.25; /* y mod 2 */ - /* - * argument reduction, make sure inexact flag not raised if input - * is an integer - */ - z = floor(y); - if(z!=y) { /* inexact anyway */ - y *= 0.5; - y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */ - n = (int) (y*4.0); - } else { - if(ix>=0x43400000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x43300000) z = y+two52; /* exact */ - GET_LOW_WORD(n,z); - n &= 1; - y = n; - n<<= 2; - } - } switch (n) { case 0: y = __kernel_sin(pi*y,zero,0); break; case 1: Modified: head/lib/msun/src/e_lgammaf_r.c ============================================================================== --- head/lib/msun/src/e_lgammaf_r.c Sun Aug 31 21:18:23 2014 (r270892) +++ head/lib/msun/src/e_lgammaf_r.c Sun Aug 31 21:38:03 2014 (r270893) @@ -89,37 +89,35 @@ w6 = -1.6309292987e-03; /* 0xbad5c4e8 * static const float zero= 0.0000000000e+00; - static float sin_pif(float x) +/* + * Compute sin(pi*x) without actually doing the pi*x multiplication. + * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is + * the precision of x. + */ +static float +sin_pi(float x) { + volatile float vz; float y,z; - int n,ix; + int n; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; + y = -x; - if(ix<0x3e800000) return __kernel_sindf(pi*x); - y = -x; /* x is assume negative */ + vz = y+0x1p23F; /* depend on 0 <= y < 0x1p23 */ + z = vz-0x1p23F; /* rintf(y) for the above range */ + if (z == y) + return (zero); + + vz = y+0x1p21F; + GET_FLOAT_WORD(n,vz); /* bits for rounded y (units 0.25) */ + z = vz-0x1p21F; /* y rounded to a multiple of 0.25 */ + if (z > y) { + z -= 0.25F; /* adjust to round down */ + n--; + } + n &= 7; /* octant of y mod 2 */ + y = y - z + n * 0.25F; /* y mod 2 */ - /* - * argument reduction, make sure inexact flag not raised if input - * is an integer - */ - z = floorf(y); - if(z!=y) { /* inexact anyway */ - y *= (float)0.5; - y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ - n = (int) (y*(float)4.0); - } else { - if(ix>=0x4b800000) { - y = zero; n = 0; /* y must be even */ - } else { - if(ix<0x4b000000) z = y+two23; /* exact */ - GET_FLOAT_WORD(n,z); - n &= 1; - y = n; - n<<= 2; - } - } switch (n) { case 0: y = __kernel_sindf(pi*y); break; case 1: @@ -157,7 +155,7 @@ __ieee754_lgammaf_r(float x, int *signga if(hx<0) { if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ return one/zero; - t = sin_pif(x); + t = sin_pi(x); if(t==zero) return one/zero; /* -integer */ nadj = __ieee754_logf(pi/fabsf(t*x)); if(t