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