From owner-freebsd-numerics@FreeBSD.ORG Fri Jun 7 09:16:42 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id B48B15BD for ; Fri, 7 Jun 2013 09:16:42 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail108.syd.optusnet.com.au (mail108.syd.optusnet.com.au [211.29.132.59]) by mx1.freebsd.org (Postfix) with ESMTP id 09AAF129B for ; Fri, 7 Jun 2013 09:16:42 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail108.syd.optusnet.com.au (Postfix) with ESMTPS id 2B78F1A1AAC for ; Fri, 7 Jun 2013 19:16:40 +1000 (EST) Date: Fri, 7 Jun 2013 19:16:38 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: numerics@freebsd.org Subject: acoshl, asinhl, atanhl almost finished Message-ID: <20130607185037.F2756@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=Q6eKePKa c=1 sm=1 a=tYOxqLTWYZkA:10 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=pVaANWmWtz0A:10 a=i0XaXLn4f-YvG9u7eWMA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 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, 07 Jun 2013 09:16:42 -0000 I converted the fdlibm double precision inverse hyperbolic functions to long double precision and am now testing them. atanhl has passed a lot of test on ld80: % --- e_atanh.c Thu May 30 18:14:16 2013 % +++ e_atanhl.c Fri Jun 7 18:07:32 2013 % @@ -1,2 +1,4 @@ % +/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */ % +/* Conversion to float by Bruce D. Evans. */ % % /* @(#)e_atanh.c 1.3 95/01/18 */ % @@ -14,49 +16,44 @@ % % #include % -__FBSDID("$FreeBSD: src/lib/msun/src/e_atanh.c,v 1.9 2012/11/17 01:50:07 svnexp Exp $"); % +__FBSDID("$FreeBSD$"); % % -/* __ieee754_atanh(x) % - * Method : % - * 1.Reduced x to positive by atanh(-x) = -atanh(x) % - * 2.For x>=0.5 % - * 1 2x x % - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) % - * 2 1 - x 1 - x % - * % - * For x<0.5 % - * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) % - * % - * Special cases: % - * atanh(x) is NaN if |x| > 1 with signal; % - * atanh(NaN) is that NaN with no signal; % - * atanh(+-1) is +-INF with signal. % - * % - */ % +#include % +#ifdef __i386__ % +#include % +#endif % % #include "math.h" % #include "math_private.h" % % -static const double one = 1.0, huge = 1e300; % +#if LDBL_MAX_EXP != 0x4000 % +/* We also require the usual expsign encoding. */ % +#error "Unsupported long double format" % +#endif % + % +#define BIAS (LDBL_MAX_EXP - 1) % + % +static const double one = 1.0; % +static volatile long double huge = 0x1p10000L; % static const double zero = 0.0; % % -double % -__ieee754_atanh(double x) % +long double % +atanhl(long double x) % { This was supposed to look as much like the double version as the float version does, including style bugs, but almost every line from here on needed changing so I fixed most of the horizontal whitespace. % - double t; % - int32_t hx,ix; % - u_int32_t lx; % - EXTRACT_WORDS(hx,lx,x); % - ix = hx&0x7fffffff; % - if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ % - return (x-x)/(x-x); Extracting all the bits for this test is too painful, and this test is inefficient anyway, so use an FP comparison for classifying +-1. Only classificatins based on the exponent alone are easy to use and efficient, especially in long double prcision. % - if(ix==0x3ff00000) % - return x/zero; % - if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ This threshold is badly chosen (and became worse when it was blindly copied for float precision). Change to a higher threshold corresponding to 2**-26 here. x is not a very good approximation to atanh(x) when x is in the changed range, but it is much better than the general approximation used for the rest of the range <= 0.5 (the latter has roundoff errors of 1.0-1.5 ulps over the whole range). % - SET_HIGH_WORD(x,ix); % - if(ix<0x3fe00000) { /* x < 0.5 */ % + long double t; % + uint16_t hx, ix; % + % + ENTERI(); % + GET_LDBL_EXPSIGN(hx, x); % + ix = hx & 0x7fff; % + if (ix >= 0x3fff) /* |x| >= 1 or x is NaN or misnormal */ % + RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x)); The RETURNI()s give larger diffs. Everything else changed literally on these 2 lines (previously 4 more complicated lines), so change the logic a little too. % + if (ix < BIAS - (LDBL_MANT_DIG + 1) / 2 && (huge + x) > zero) % + RETURNI(x); /* |x| < 2**-(MANT_DIG+1)/2 */ Misnormals are handled in all cases by arranging to do arithmetic on them so that the get converted to a qNaN for return. For example, if x is unnormal here, (huge + x) > zero is false here (x acts like an sNaN), so we don't just return it here, but fall through to cases where it gets converted. % + SET_LDBL_EXPSIGN(x, ix); % + if (ix < 0x3ffe) { /* |x| < 0.5 or x is misnormal */ % t = x+x; % - t = 0.5*log1p(t+t*x/(one-x)); % + t = 0.5*log1pl(t+t*x/(one-x)); % } else Preserving the style bugs includes not adding a comment for this case. % - t = 0.5*log1p((x+x)/(one-x)); % - if(hx>=0) return t; else return -t; % + t = 0.5*log1pl((x+x)/(one-x)); % + RETURNI((hx & 0x8000) == 0 ? t : -t); % } As expected since the algorithm wasn't changed and there few possibilities for numerical accidents, the errors in ulps seem to be about the same as in lower precisions (actually almost half an ulp better because log1pl() is better). Bruce