path: root/src/math/acoshl.c

/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.c */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */
/* acoshl(x)
 * Method :
 *      Based on
 *              acoshl(x) = logl [ x + sqrtl(x*x-1) ]
 *      we have
 *              acoshl(x) := logl(x)+ln2,       if x is large; else
 *              acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
 *              acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
 *
 * Special cases:
 *      acoshl(x) is NaN with signal if x<1.
 *      acoshl(NaN) is NaN without signal.
 */

#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double acoshl(long double x)
{
	return acosh(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
static const long double
one = 1.0,
ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */

long double acoshl(long double x)
{
	long double t;
	uint32_t se,i0,i1;

	GET_LDOUBLE_WORDS(se, i0, i1, x);
	if (se < 0x3fff || se & 0x8000) {  /* x < 1 */
		return (x-x)/(x-x);
	} else if (se >= 0x401d) {  /* x > 2**30 */
		if (se >= 0x7fff)  /* x is inf or NaN */
			return x+x;
		return logl(x) + ln2;  /* acoshl(huge) = logl(2x) */
	} else if (((se-0x3fff)|i0|i1) == 0) {
		return 0.0;            /* acosh(1) = 0 */
	} else if (se > 0x4000) {  /* x > 2 */
		t = x*x;
		return logl(2.0*x - one/(x + sqrtl(t - one)));
	}
	/* 1 < x <= 2 */
	t = x - one;
	return log1pl(t + sqrtl(2.0*t + t*t));
}
#endif