1141296Sdas
2141296Sdas/* @(#)e_acosh.c 1.3 95/01/18 */
32116Sjkh/*
42116Sjkh * ====================================================
52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
62116Sjkh *
7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business.
82116Sjkh * Permission to use, copy, modify, and distribute this
9141296Sdas * software is freely granted, provided that this notice
102116Sjkh * is preserved.
112116Sjkh * ====================================================
12141296Sdas *
132116Sjkh */
142116Sjkh
15176451Sdas#include <sys/cdefs.h>
16176451Sdas__FBSDID("$FreeBSD$");
172116Sjkh
182116Sjkh/* __ieee754_acosh(x)
192116Sjkh * Method :
20141296Sdas *	Based on
212116Sjkh *		acosh(x) = log [ x + sqrt(x*x-1) ]
222116Sjkh *	we have
232116Sjkh *		acosh(x) := log(x)+ln2,	if x is large; else
242116Sjkh *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
252116Sjkh *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
262116Sjkh *
272116Sjkh * Special cases:
282116Sjkh *	acosh(x) is NaN with signal if x<1.
292116Sjkh *	acosh(NaN) is NaN without signal.
302116Sjkh */
312116Sjkh
32251599Sdas#include <float.h>
33251599Sdas
342116Sjkh#include "math.h"
352116Sjkh#include "math_private.h"
362116Sjkh
378870Srgrimesstatic const double
382116Sjkhone	= 1.0,
392116Sjkhln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
402116Sjkh
4197407Salfreddouble
4297407Salfred__ieee754_acosh(double x)
438870Srgrimes{
442116Sjkh	double t;
452116Sjkh	int32_t hx;
462116Sjkh	u_int32_t lx;
472116Sjkh	EXTRACT_WORDS(hx,lx,x);
482116Sjkh	if(hx<0x3ff00000) {		/* x < 1 */
492116Sjkh	    return (x-x)/(x-x);
502116Sjkh	} else if(hx >=0x41b00000) {	/* x > 2**28 */
512116Sjkh	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
522116Sjkh	        return x+x;
53141296Sdas	    } else
542116Sjkh		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
552116Sjkh	} else if(((hx-0x3ff00000)|lx)==0) {
562116Sjkh	    return 0.0;			/* acosh(1) = 0 */
572116Sjkh	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
582116Sjkh	    t=x*x;
59141296Sdas	    return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
602116Sjkh	} else {			/* 1<x<2 */
612116Sjkh	    t = x-one;
62141296Sdas	    return log1p(t+sqrt(2.0*t+t*t));
632116Sjkh	}
642116Sjkh}
65251599Sdas
66251599Sdas#if LDBL_MANT_DIG == 53
67251599Sdas__weak_reference(acosh, acoshl);
68251599Sdas#endif
69