Date: Fri, 7 Jun 2013 19:16:38 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: numerics@freebsd.org Subject: acoshl, asinhl, atanhl almost finished Message-ID: <20130607185037.F2756@besplex.bde.org>
next in thread | raw e-mail | index | archive | help
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 <sys/cdefs.h>
% -__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 <float.h>
% +#ifdef __i386__
% +#include <ieeefp.h>
% +#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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130607185037.F2756>