1251599Sdas/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */ 2141296Sdas 3141296Sdas/* @(#)e_acosh.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_acosh.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_LARGE is the threshold above which we use acosh(x) ~= log(2x). */ 36251599Sdas#if LDBL_MANT_DIG == 64 37251599Sdas#define EXP_LARGE 34 38251599Sdas#elif LDBL_MANT_DIG == 113 39251599Sdas#define EXP_LARGE 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 518870Srgrimesstatic const double 52251599Sdasone = 1.0; 532116Sjkh 54251599Sdas#if LDBL_MANT_DIG == 64 55251599Sdasstatic const union IEEEl2bits 56251599Sdasu_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); 57251599Sdas#define ln2 u_ln2.e 58251599Sdas#elif LDBL_MANT_DIG == 113 59251599Sdasstatic const long double 60251599Sdasln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ 61251599Sdas#else 62251599Sdas#error "Unsupported long double format" 63251599Sdas#endif 64251599Sdas 65251599Sdaslong double 66251599Sdasacoshl(long double x) 678870Srgrimes{ 68251599Sdas long double t; 69251599Sdas int16_t hx; 70251599Sdas 71251599Sdas ENTERI(); 72251599Sdas GET_LDBL_EXPSIGN(hx, x); 73251599Sdas if (hx < 0x3fff) { /* x < 1, or misnormal */ 74251599Sdas RETURNI((x-x)/(x-x)); 75251599Sdas } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */ 76251599Sdas if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */ 77251599Sdas RETURNI(x+x); 78141296Sdas } else 79251599Sdas RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */ 80251599Sdas } else if (hx == 0x3fff && x == 1) { 81251599Sdas RETURNI(0.0); /* acosh(1) = 0 */ 82251599Sdas } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */ 832116Sjkh t=x*x; 84251599Sdas RETURNI(logl(2.0*x-one/(x+sqrtl(t-one)))); 852116Sjkh } else { /* 1<x<2 */ 862116Sjkh t = x-one; 87251599Sdas RETURNI(log1pl(t+sqrtl(2.0*t+t*t))); 882116Sjkh } 892116Sjkh} 90