|
|
1.1 root 1: /*
2: * Copyright (c) 1985 Regents of the University of California.
3: * All rights reserved.
4: *
5: * Redistribution and use in source and binary forms are permitted
6: * provided that: (1) source distributions retain this entire copyright
7: * notice and comment, and (2) distributions including binaries display
8: * the following acknowledgement: ``This product includes software
9: * developed by the University of California, Berkeley and its contributors''
10: * in the documentation or other materials provided with the distribution
11: * and in all advertising materials mentioning features or use of this
12: * software. Neither the name of the University nor the names of its
13: * contributors may be used to endorse or promote products derived
14: * from this software without specific prior written permission.
15: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
16: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
17: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
18: *
19: * All recipients should regard themselves as participants in an ongoing
20: * research project and hence should feel obligated to report their
21: * experiences (good or bad) with these elementary function codes, using
22: * the sendbug(8) program, to the authors.
23: */
24:
25: #ifndef lint
26: static char sccsid[] = "@(#)acosh.c 5.5 (Berkeley) 6/1/90";
27: #endif /* not lint */
28:
29: /* ACOSH(X)
30: * RETURN THE INVERSE HYPERBOLIC COSINE OF X
31: * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS)
32: * CODED IN C BY K.C. NG, 2/16/85;
33: * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85.
34: *
35: * Required system supported functions :
36: * sqrt(x)
37: *
38: * Required kernel function:
39: * log1p(x) ...return log(1+x)
40: *
41: * Method :
42: * Based on
43: * acosh(x) = log [ x + sqrt(x*x-1) ]
44: * we have
45: * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else
46: * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) .
47: * These formulae avoid the over/underflow complication.
48: *
49: * Special cases:
50: * acosh(x) is NaN with signal if x<1.
51: * acosh(NaN) is NaN without signal.
52: *
53: * Accuracy:
54: * acosh(x) returns the exact inverse hyperbolic cosine of x nearly
55: * rounded. In a test run with 512,000 random arguments on a VAX, the
56: * maximum observed error was 3.30 ulps (units of the last place) at
57: * x=1.0070493753568216 .
58: *
59: * Constants:
60: * The hexadecimal values are the intended ones for the following constants.
61: * The decimal values may be used, provided that the compiler will convert
62: * from decimal to binary accurately enough to produce the hexadecimal values
63: * shown.
64: */
65:
66: #include "mathimpl.h"
67:
68: vc(ln2hi, 6.9314718055829871446E-1 ,7217,4031,0000,f7d0, 0, .B17217F7D00000)
69: vc(ln2lo, 1.6465949582897081279E-12 ,bcd5,2ce7,d9cc,e4f1, -39, .E7BCD5E4F1D9CC)
70:
71: ic(ln2hi, 6.9314718036912381649E-1, -1, 1.62E42FEE00000)
72: ic(ln2lo, 1.9082149292705877000E-10,-33, 1.A39EF35793C76)
73:
74: #ifdef vccast
75: #define ln2hi vccast(ln2hi)
76: #define ln2lo vccast(ln2lo)
77: #endif
78:
79: double acosh(x)
80: double x;
81: {
82: double t,big=1.E20; /* big+1==big */
83:
84: #if !defined(vax)&&!defined(tahoe)
85: if(x!=x) return(x); /* x is NaN */
86: #endif /* !defined(vax)&&!defined(tahoe) */
87:
88: /* return log1p(x) + log(2) if x is large */
89: if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);}
90:
91: t=sqrt(x-1.0);
92: return(log1p(t*(t+sqrt(x+1.0))));
93: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.