1251599Sdas/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */ 2141296Sdas 3141296Sdas/* @(#)e_atanh.c 1.3 95/01/18 */ 42116Sjkh/* 52116Sjkh * ==================================================== 62116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 72116Sjkh * 8141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 92116Sjkh * Permission to use, copy, modify, and distribute this 10141296Sdas * software is freely granted, provided that this notice 112116Sjkh * is preserved. 122116Sjkh * ==================================================== 13141296Sdas * 142116Sjkh */ 152116Sjkh 16176451Sdas#include <sys/cdefs.h> 17176451Sdas__FBSDID("$FreeBSD$"); 182116Sjkh 19251599Sdas/* 20251599Sdas * See e_atanh.c for complete comments. 212116Sjkh * 22251599Sdas * Converted to long double by David Schultz <das@FreeBSD.ORG> and 23251599Sdas * Bruce D. Evans. 242116Sjkh */ 252116Sjkh 26251599Sdas#include <float.h> 27251599Sdas#ifdef __i386__ 28251599Sdas#include <ieeefp.h> 29251599Sdas#endif 30251599Sdas 31251599Sdas#include "fpmath.h" 322116Sjkh#include "math.h" 332116Sjkh#include "math_private.h" 342116Sjkh 35251599Sdas/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */ 36251599Sdas#if LDBL_MANT_DIG == 64 37251599Sdas#define EXP_TINY -34 38251599Sdas#elif LDBL_MANT_DIG == 113 39251599Sdas#define EXP_TINY -58 40251599Sdas#else 41251599Sdas#error "Unsupported long double format" 42251599Sdas#endif 43251599Sdas 44251599Sdas#if LDBL_MAX_EXP != 0x4000 45251599Sdas/* We also require the usual expsign encoding. */ 46251599Sdas#error "Unsupported long double format" 47251599Sdas#endif 48251599Sdas 49251599Sdas#define BIAS (LDBL_MAX_EXP - 1) 50251599Sdas 512116Sjkhstatic const double one = 1.0, huge = 1e300; 522116Sjkhstatic const double zero = 0.0; 532116Sjkh 54251599Sdaslong double 55251599Sdasatanhl(long double x) 562116Sjkh{ 57251599Sdas long double t; 58251599Sdas uint16_t hx, ix; 59251599Sdas 60251599Sdas ENTERI(); 61251599Sdas GET_LDBL_EXPSIGN(hx, x); 62251599Sdas ix = hx & 0x7fff; 63251599Sdas if (ix >= 0x3fff) /* |x| >= 1, or NaN or misnormal */ 64251599Sdas RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x)); 65251599Sdas if (ix < BIAS + EXP_TINY && (huge + x) > zero) 66251599Sdas RETURNI(x); /* x is tiny */ 67251599Sdas SET_LDBL_EXPSIGN(x, ix); 68251599Sdas if (ix < 0x3ffe) { /* |x| < 0.5, or misnormal */ 692116Sjkh t = x+x; 70251599Sdas t = 0.5*log1pl(t+t*x/(one-x)); 71141296Sdas } else 72251599Sdas t = 0.5*log1pl((x+x)/(one-x)); 73251599Sdas RETURNI((hx & 0x8000) == 0 ? t : -t); 742116Sjkh} 75