xref: /dflybsd-src/contrib/openbsd_libm/src/e_acosh.c (revision 4382f29d99a100bd77a81697c2f699c11f6a472a)
1*05a0b428SJohn Marino /* @(#)e_acosh.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino  * ====================================================
4*05a0b428SJohn Marino  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino  *
6*05a0b428SJohn Marino  * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino  * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino  * software is freely granted, provided that this notice
9*05a0b428SJohn Marino  * is preserved.
10*05a0b428SJohn Marino  * ====================================================
11*05a0b428SJohn Marino  */
12*05a0b428SJohn Marino 
13*05a0b428SJohn Marino /* acosh(x)
14*05a0b428SJohn Marino  * Method :
15*05a0b428SJohn Marino  *	Based on
16*05a0b428SJohn Marino  *		acosh(x) = log [ x + sqrt(x*x-1) ]
17*05a0b428SJohn Marino  *	we have
18*05a0b428SJohn Marino  *		acosh(x) := log(x)+ln2,	if x is large; else
19*05a0b428SJohn Marino  *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
20*05a0b428SJohn Marino  *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
21*05a0b428SJohn Marino  *
22*05a0b428SJohn Marino  * Special cases:
23*05a0b428SJohn Marino  *	acosh(x) is NaN with signal if x<1.
24*05a0b428SJohn Marino  *	acosh(NaN) is NaN without signal.
25*05a0b428SJohn Marino  */
26*05a0b428SJohn Marino 
27*05a0b428SJohn Marino #include <float.h>
28*05a0b428SJohn Marino #include <math.h>
29*05a0b428SJohn Marino 
30*05a0b428SJohn Marino #include "math_private.h"
31*05a0b428SJohn Marino 
32*05a0b428SJohn Marino static const double
33*05a0b428SJohn Marino one	= 1.0,
34*05a0b428SJohn Marino ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
35*05a0b428SJohn Marino 
36*05a0b428SJohn Marino double
acosh(double x)37*05a0b428SJohn Marino acosh(double x)
38*05a0b428SJohn Marino {
39*05a0b428SJohn Marino 	double t;
40*05a0b428SJohn Marino 	int32_t hx;
41*05a0b428SJohn Marino 	u_int32_t lx;
42*05a0b428SJohn Marino 	EXTRACT_WORDS(hx,lx,x);
43*05a0b428SJohn Marino 	if(hx<0x3ff00000) {		/* x < 1 */
44*05a0b428SJohn Marino 	    return (x-x)/(x-x);
45*05a0b428SJohn Marino 	} else if(hx >=0x41b00000) {	/* x > 2**28 */
46*05a0b428SJohn Marino 	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
47*05a0b428SJohn Marino 	        return x+x;
48*05a0b428SJohn Marino 	    } else
49*05a0b428SJohn Marino 		return log(x)+ln2;	/* acosh(huge)=log(2x) */
50*05a0b428SJohn Marino 	} else if(((hx-0x3ff00000)|lx)==0) {
51*05a0b428SJohn Marino 	    return 0.0;			/* acosh(1) = 0 */
52*05a0b428SJohn Marino 	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
53*05a0b428SJohn Marino 	    t=x*x;
54*05a0b428SJohn Marino 	    return log(2.0*x-one/(x+sqrt(t-one)));
55*05a0b428SJohn Marino 	} else {			/* 1<x<2 */
56*05a0b428SJohn Marino 	    t = x-one;
57*05a0b428SJohn Marino 	    return log1p(t+sqrt(2.0*t+t*t));
58*05a0b428SJohn Marino 	}
59*05a0b428SJohn Marino }
60*05a0b428SJohn Marino 
61*05a0b428SJohn Marino #if	LDBL_MANT_DIG == DBL_MANT_DIG
62*05a0b428SJohn Marino __strong_alias(acoshl, acosh);
63*05a0b428SJohn Marino #endif	/* LDBL_MANT_DIG == DBL_MANT_DIG */
64