From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from gtd-gmbh.de (mail.gtd.eu [46.24.46.35]) by sourceware.org (Postfix) with ESMTPS id CED9F3871029 for ; Tue, 17 Mar 2020 14:57:33 +0000 (GMT) X-MDAV-Result: clean X-MDAV-Processed: gtd-gmbh.de, Tue, 17 Mar 2020 15:56:24 +0100 Received: from localhost.localdomain [(78.42.121.247)] by gtd-gmbh.de (MDaemon PRO v19.5.4) with ESMTPA id md50013616789.msg; Tue, 17 Mar 2020 15:56:23 +0100 X-Spam-Processed: gtd-gmbh.de, Tue, 17 Mar 2020 15:56:23 +0100 (not processed: message from trusted or authenticated source) X-MDRemoteIP: 78.42.121.247 X-MDHelo: localhost.localdomain X-MDArrival-Date: Tue, 17 Mar 2020 15:56:23 +0100 X-Authenticated-Sender: fabian.schriever@gtd-gmbh.de X-Return-Path: prvs=1345bad6eb=fabian.schriever@gtd-gmbh.de X-Envelope-From: fabian.schriever@gtd-gmbh.de X-MDaemon-Deliver-To: newlib@sourceware.org From: Fabian Schriever To: newlib@sourceware.org Cc: Fabian Schriever Subject: [PATCH] Fix for k_tan.c specific inputs Date: Tue, 17 Mar 2020 15:48:44 +0100 Message-Id: <20200317144844.1075-1-fabian.schriever@gtd-gmbh.de> X-Mailer: git-send-email 2.24.1.windows.2 MIME-Version: 1.0 Content-Transfer-Encoding: 8bit X-Spam-Status: No, score=-23.2 required=5.0 tests=GIT_PATCH_0, GIT_PATCH_1, GIT_PATCH_2, GIT_PATCH_3, KAM_DMARC_STATUS, KAM_LAZY_DOMAIN_SECURITY, RCVD_IN_DNSWL_NONE, SPF_HELO_NONE, SPF_NONE, SUBJ_OBFU_PUNCT_FEW, UNPARSEABLE_RELAY autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: newlib@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Newlib mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 17 Mar 2020 14:57:37 -0000 This fix for k_tan.c is a copy from fdlibm version 5.3 (see also http://www.netlib.org/fdlibm/readme), adjusted to use the macros available in newlib (SET_LOW_WORD). This fix reduces the ULP error of the value shown in the fdlibm readme (tan(1.7765241907548024E+269)) to 0.45 (thereby reducing the error by 1). This issue only happens for large numbers that get reduced by the range reduction to a value smaller in magnitude than 2^-28, that is also reduced an uneven number of times. This seems rather unlikely given that one ULP is (much) larger than 2^-28 for the values that may cause an issue. Although given the sheer number of values a double can represent, it is still possible that there are more affected values, finding them however will be quite hard, if not impossible. We also took a look at how another library (libm in FreeBSD) handles the issue: In FreeBSD the complete if branch which checks for values smaller than 2^-28 (or rather 2^-27, another change done by FreeBSD) is moved out of the kernel function and into the external function. This means that the value that gets checked for this condition is the unreduced value. Therefore the input value which caused a problem in the fdlibm/newlib kernel tan will run through the full polynomial, including the careful calculation of -1/(x+r). So the difference is really whether r or y is used. r = y + p with p being the result of the polynomial with 1/3*x^3 being the largest (and magnitude defining) value. With x being <2^-27 we therefore know that p is smaller than y (y has to be at least the size of the value of x last mantissa bit divided by 2, which is at least x*2^-51 for doubles) by enough to warrant saying that r ~ y. So we can conclude that the general implementation of this special case is the same, FreeBSD simply has a different philosophy on when to handle especially small numbers. --- newlib/libm/math/k_tan.c | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/newlib/libm/math/k_tan.c b/newlib/libm/math/k_tan.c index 9f5b30760..4be82d599 100644 --- a/newlib/libm/math/k_tan.c +++ b/newlib/libm/math/k_tan.c @@ -84,14 +84,27 @@ T[] = { __int32_t ix,hx; GET_HIGH_WORD(hx,x); ix = hx&0x7fffffff; /* high word of |x| */ - if(ix<0x3e300000) /* x < 2**-28 */ - {if((int)x==0) { /* generate inexact */ - __uint32_t low; - GET_LOW_WORD(low,x); - if(((ix|low)|(iy+1))==0) return one/fabs(x); - else return (iy==1)? x: -one/x; - } - } + if(ix<0x3e300000) { /* x < 2**-28 */ + if((int)x==0) { /* generate inexact */ + __uint32_t low; + GET_LOW_WORD(low,x); + if(((ix|low)|(iy+1))==0) return one/fabs(x); + else { + if(iy==1) + return x; + else { + double a, t; + z = w = x + y; + SET_LOW_WORD(z,0); + v = y - (z - x); + t = a = -one / w; + SET_LOW_WORD(t,0); + s = one + t * z; + return t + a * (s + t * v); + } + } + } + } if(ix>=0x3FE59428) { /* |x|>=0.6744 */ if(hx<0) {x = -x; y = -y;} z = pio4-x; -- 2.24.1.windows.2