xref: /illumos-gate/usr/src/lib/libm/common/C/acosh.c (revision fd6d41c5025e9fb45a115fc82d86e9983d1e9fd6)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24  */
25 /*
26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27  * Use is subject to license terms.
28  */
29 
30 #pragma weak __acosh = acosh
31 
32 /* INDENT OFF */
33 /*
34  * acosh(x)
35  * Method :
36  *	Based on
37  *		acosh(x) = log [ x + sqrt(x*x-1) ]
38  *	we have
39  *		acosh(x) := log(x)+ln2,	if x is large; else
40  *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x > 2; else
41  *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t = x-1.
42  *
43  * Special cases:
44  *	acosh(x) is NaN with signal if x < 1.
45  *	acosh(NaN) is NaN without signal.
46  */
47 /* INDENT ON */
48 
49 #include "libm_protos.h"	/* _SVID_libm_error */
50 #include "libm_macros.h"
51 #include <math.h>
52 
53 static const double
54 	one = 1.0,
55 	ln2 = 6.93147180559945286227e-01;	/* 3FE62E42, FEFA39EF */
56 
57 double
58 acosh(double x) {
59 	double t;
60 	int hx;
61 
62 	hx = ((int *) &x)[HIWORD];
63 	if (hx < 0x3ff00000) {	/* x < 1 */
64 		if (isnan(x))
65 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
66 			return (hx >= 0xfff80000 ? x : (x - x) / (x - x));
67 			/* assumes sparc-like QNaN */
68 #else
69 			return (x - x) / (x - x);
70 #endif
71 		else
72 			return (_SVID_libm_err(x, x, 29));
73 	} else if (hx >= 0x41b00000) {
74 		/* x > 2**28 */
75 		if (hx >= 0x7ff00000) {	/* x is inf of NaN */
76 #if defined(FPADD_TRAPS_INCOMPLETE_ON_NAN)
77 			return (hx >= 0x7ff80000 ? x : x + x);
78 			/* assumes sparc-like QNaN */
79 #else
80 			return (x + x);
81 #endif
82 		} else	/* acosh(huge)=log(2x) */
83 			return (log(x) + ln2);
84 	} else if (((hx - 0x3ff00000) | ((int *) &x)[LOWORD]) == 0) {
85 		return (0.0);	/* acosh(1) = 0 */
86 	} else if (hx > 0x40000000) {
87 		/* 2**28 > x > 2 */
88 		t = x * x;
89 		return (log(2.0 * x - one / (x + sqrt(t - one))));
90 	} else {
91 		/* 1 < x < 2 */
92 		t = x - one;
93 		return (log1p(t + sqrt(2.0 * t + t * t)));
94 	}
95 }
96