1*998b640bSDavid Schultz /* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */ 2*998b640bSDavid Schultz 3*998b640bSDavid Schultz /* @(#)e_acosh.c 1.3 95/01/18 */ 4*998b640bSDavid Schultz /* 5*998b640bSDavid Schultz * ==================================================== 6*998b640bSDavid Schultz * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 7*998b640bSDavid Schultz * 8*998b640bSDavid Schultz * Developed at SunSoft, a Sun Microsystems, Inc. business. 9*998b640bSDavid Schultz * Permission to use, copy, modify, and distribute this 10*998b640bSDavid Schultz * software is freely granted, provided that this notice 11*998b640bSDavid Schultz * is preserved. 12*998b640bSDavid Schultz * ==================================================== 13*998b640bSDavid Schultz * 14*998b640bSDavid Schultz */ 15*998b640bSDavid Schultz 16*998b640bSDavid Schultz #include <sys/cdefs.h> 17*998b640bSDavid Schultz __FBSDID("$FreeBSD$"); 18*998b640bSDavid Schultz 19*998b640bSDavid Schultz /* 20*998b640bSDavid Schultz * See e_acosh.c for complete comments. 21*998b640bSDavid Schultz * 22*998b640bSDavid Schultz * Converted to long double by David Schultz <das@FreeBSD.ORG> and 23*998b640bSDavid Schultz * Bruce D. Evans. 24*998b640bSDavid Schultz */ 25*998b640bSDavid Schultz 26*998b640bSDavid Schultz #include <float.h> 27*998b640bSDavid Schultz #ifdef __i386__ 28*998b640bSDavid Schultz #include <ieeefp.h> 29*998b640bSDavid Schultz #endif 30*998b640bSDavid Schultz 31*998b640bSDavid Schultz #include "fpmath.h" 32*998b640bSDavid Schultz #include "math.h" 33*998b640bSDavid Schultz #include "math_private.h" 34*998b640bSDavid Schultz 35*998b640bSDavid Schultz /* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */ 36*998b640bSDavid Schultz #if LDBL_MANT_DIG == 64 37*998b640bSDavid Schultz #define EXP_LARGE 34 38*998b640bSDavid Schultz #elif LDBL_MANT_DIG == 113 39*998b640bSDavid Schultz #define EXP_LARGE 58 40*998b640bSDavid Schultz #else 41*998b640bSDavid Schultz #error "Unsupported long double format" 42*998b640bSDavid Schultz #endif 43*998b640bSDavid Schultz 44*998b640bSDavid Schultz #if LDBL_MAX_EXP != 0x4000 45*998b640bSDavid Schultz /* We also require the usual expsign encoding. */ 46*998b640bSDavid Schultz #error "Unsupported long double format" 47*998b640bSDavid Schultz #endif 48*998b640bSDavid Schultz 49*998b640bSDavid Schultz #define BIAS (LDBL_MAX_EXP - 1) 50*998b640bSDavid Schultz 51*998b640bSDavid Schultz static const double 52*998b640bSDavid Schultz one = 1.0; 53*998b640bSDavid Schultz 54*998b640bSDavid Schultz #if LDBL_MANT_DIG == 64 55*998b640bSDavid Schultz static const union IEEEl2bits 56*998b640bSDavid Schultz u_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); 57*998b640bSDavid Schultz #define ln2 u_ln2.e 58*998b640bSDavid Schultz #elif LDBL_MANT_DIG == 113 59*998b640bSDavid Schultz static const long double 60*998b640bSDavid Schultz ln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ 61*998b640bSDavid Schultz #else 62*998b640bSDavid Schultz #error "Unsupported long double format" 63*998b640bSDavid Schultz #endif 64*998b640bSDavid Schultz 65*998b640bSDavid Schultz long double 66*998b640bSDavid Schultz acoshl(long double x) 67*998b640bSDavid Schultz { 68*998b640bSDavid Schultz long double t; 69*998b640bSDavid Schultz int16_t hx; 70*998b640bSDavid Schultz 71*998b640bSDavid Schultz ENTERI(); 72*998b640bSDavid Schultz GET_LDBL_EXPSIGN(hx, x); 73*998b640bSDavid Schultz if (hx < 0x3fff) { /* x < 1, or misnormal */ 74*998b640bSDavid Schultz RETURNI((x-x)/(x-x)); 75*998b640bSDavid Schultz } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */ 76*998b640bSDavid Schultz if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */ 77*998b640bSDavid Schultz RETURNI(x+x); 78*998b640bSDavid Schultz } else 79*998b640bSDavid Schultz RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */ 80*998b640bSDavid Schultz } else if (hx == 0x3fff && x == 1) { 81*998b640bSDavid Schultz RETURNI(0.0); /* acosh(1) = 0 */ 82*998b640bSDavid Schultz } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */ 83*998b640bSDavid Schultz t=x*x; 84*998b640bSDavid Schultz RETURNI(logl(2.0*x-one/(x+sqrtl(t-one)))); 85*998b640bSDavid Schultz } else { /* 1<x<2 */ 86*998b640bSDavid Schultz t = x-one; 87*998b640bSDavid Schultz RETURNI(log1pl(t+sqrtl(2.0*t+t*t))); 88*998b640bSDavid Schultz } 89*998b640bSDavid Schultz } 90