1251599Sdas/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */ 2251599Sdas 32116Sjkh/* @(#)s_asinh.c 5.1 93/09/24 */ 42116Sjkh/* 52116Sjkh * ==================================================== 62116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 72116Sjkh * 82116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 92116Sjkh * Permission to use, copy, modify, and distribute this 108870Srgrimes * software is freely granted, provided that this notice 112116Sjkh * is preserved. 122116Sjkh * ==================================================== 132116Sjkh */ 142116Sjkh 15176451Sdas#include <sys/cdefs.h> 16176451Sdas__FBSDID("$FreeBSD$"); 172116Sjkh 18251599Sdas/* 19251599Sdas * See s_asinh.c for complete comments. 20251599Sdas * 21251599Sdas * Converted to long double by David Schultz <das@FreeBSD.ORG> and 22251599Sdas * Bruce D. Evans. 232116Sjkh */ 242116Sjkh 25251599Sdas#include <float.h> 26251599Sdas#ifdef __i386__ 27251599Sdas#include <ieeefp.h> 28251599Sdas#endif 29251599Sdas 30251599Sdas#include "fpmath.h" 312116Sjkh#include "math.h" 322116Sjkh#include "math_private.h" 332116Sjkh 34251599Sdas/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */ 35251599Sdas/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */ 36251599Sdas#if LDBL_MANT_DIG == 64 37251599Sdas#define EXP_LARGE 34 38251599Sdas#define EXP_TINY -34 39251599Sdas#elif LDBL_MANT_DIG == 113 40251599Sdas#define EXP_LARGE 58 41251599Sdas#define EXP_TINY -58 42251599Sdas#else 43251599Sdas#error "Unsupported long double format" 44251599Sdas#endif 45251599Sdas 46251599Sdas#if LDBL_MAX_EXP != 0x4000 47251599Sdas/* We also require the usual expsign encoding. */ 48251599Sdas#error "Unsupported long double format" 49251599Sdas#endif 50251599Sdas 51251599Sdas#define BIAS (LDBL_MAX_EXP - 1) 52251599Sdas 538870Srgrimesstatic const double 542116Sjkhone = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 558870Srgrimeshuge= 1.00000000000000000000e+300; 562116Sjkh 57251599Sdas#if LDBL_MANT_DIG == 64 58251599Sdasstatic const union IEEEl2bits 59251599Sdasu_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); 60251599Sdas#define ln2 u_ln2.e 61251599Sdas#elif LDBL_MANT_DIG == 113 62251599Sdasstatic const long double 63251599Sdasln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ 64251599Sdas#else 65251599Sdas#error "Unsupported long double format" 66251599Sdas#endif 67251599Sdas 68251599Sdaslong double 69251599Sdasasinhl(long double x) 708870Srgrimes{ 71251599Sdas long double t, w; 72251599Sdas uint16_t hx, ix; 73251599Sdas 74251599Sdas ENTERI(); 75251599Sdas GET_LDBL_EXPSIGN(hx, x); 76251599Sdas ix = hx & 0x7fff; 77251599Sdas if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */ 78251599Sdas if (ix < BIAS + EXP_TINY) { /* |x| < TINY, or misnormal */ 79251599Sdas if (huge + x > one) RETURNI(x); /* return x inexact except 0 */ 808870Srgrimes } 81251599Sdas if (ix >= BIAS + EXP_LARGE) { /* |x| >= LARGE, or misnormal */ 82251599Sdas w = logl(fabsl(x))+ln2; 83251599Sdas } else if (ix >= 0x4000) { /* LARGE > |x| >= 2.0, or misnormal */ 84251599Sdas t = fabsl(x); 85251599Sdas w = logl(2.0*t+one/(sqrtl(x*x+one)+t)); 86251599Sdas } else { /* 2.0 > |x| >= TINY, or misnormal */ 872116Sjkh t = x*x; 88251599Sdas w =log1pl(fabsl(x)+t/(one+sqrtl(one+t))); 892116Sjkh } 90251599Sdas RETURNI((hx & 0x8000) == 0 ? w : -w); 912116Sjkh} 92