From owner-freebsd-numerics@FreeBSD.ORG Fri Mar 13 05:34:41 2015 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id A6D574A2 for ; Fri, 13 Mar 2015 05:34:41 +0000 (UTC) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id 88DF1C0D for ; Fri, 13 Mar 2015 05:34:41 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id t2D5Ho8A063320 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Thu, 12 Mar 2015 22:17:50 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id t2D5HocW063319; Thu, 12 Mar 2015 22:17:50 -0700 (PDT) (envelope-from sgk) Date: Thu, 12 Mar 2015 22:17:50 -0700 From: Steve Kargl To: enh Subject: Re: libm functions sensitive to current rounding mode Message-ID: <20150313051749.GA63174@troutmask.apl.washington.edu> References: MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.18-1 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 13 Mar 2015 05:34:41 -0000 On Thu, Mar 12, 2015 at 08:58:06PM -0700, enh via freebsd-numerics wrote: > various SoC vendors have been sending us (Android) assembler > implementations of various libm functions, and we noticed when looking > at one change that it was slower on certain inputs than the FreeBSD C > implementation we were already using. the response was that the cost > was coping with other rounding modes, and they gave this example: > > #include > #include > #include > > int main(){ > double x = -0xb.6365e22ee46fp4; > fesetround(FE_UPWARD); > printf("sin(%f) = %.16e\n", x, sin(x)); > fesetround(FE_DOWNWARD); > printf("sin(%f) = %.16e\n", x, sin(x)); > > x = 0x1.921fb54442d16p0; > fesetround(FE_UPWARD); > printf("sin(%f) = %.16e\n", x, sin(x)); > fesetround(FE_DOWNWARD); > printf("sin(%f) = %.16e\n", x, sin(x)); > return 0; > } > > see https://android-review.googlesource.com/112700 for the original > change and related conversation. > > glibc 2.19: > > sin(-182.212373) = -2.4759225463534308e-18 > sin(-182.212374) = -2.4759225463534309e-18 > sin(1.570797) = 1.0000000000000000e+00 > sin(1.570796) = 1.0000000000000000e+00 > > (glibc's fixed bugs here in the past. see > https://sourceware.org/bugzilla/show_bug.cgi?id=3976 for example.) > > the FreeBSD libm: > > sin(-182.212374) = -2.4759225463534304e-18 > sin(-182.212374) = 2.2575413760606011e-11 > sin(1.570796) = 1.0000000000000000e+00 > sin(1.570796) = 1.0000000002522276e+00 > > it looks like e_sqrtl.c, s_fmal.c, s_fma.c, and s_fmaf.c save/restore > the rounding mode but nothing else does. > > let me know if you'd like me to send trivial patches for any of the > affected functions. (do all the libm functions need this?) > http://wiki.freebsd.org/Numerics documents the progress on implementing missing libm functions. It's somewhat out of date as only powl and tgammal as required by C99/11 are missing for the long double functions. I don't know the status of complex.h functions. I don't do wiki, so haven't even tried to update the page. At the bottom of the page, you'll see a list of future directions. Dealing with rounding modes is on the list. As of now, FreeBSD libm should work with round-to-nearest; other modes probably need to be tested. e_sqrtl.c was implemented to handle the rounding mode for 2 reasons: (1) IEEE-754 mandates that square root is correctly rounded in all rounding modes; and (2), fdlibm's src/e_sqrt.c (might be wrong filename) has a long comment that explains an algorithm that can be used to implement software sqrt() in all rounding modes. That's what I implemented and then Bruce Evans and David Schultz fixed what I did. I suspect all other functions, which are implemented in software, need to be tested and patches developed. If you have trivial patches, the best thing to do is submit a bug report via FreeBSD's bugzilla so that the patch doesn't get lost. Someday after I get powl and tgammal written, and I become independently wealthy, I may have time to work on libm. BTW, n1785.pdf is adding a number new functions to the C library (e.g., sinpi(x), etc). I suspect I would work on those before dealing with rounding modes; as all of my personal codes are expecting round-to-nearest. -- Steve