Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
4973 right-hear 1
/* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */
2
/* @(#)e_acosh.c 5.1 93/09/24 */
3
/*
4
 * ====================================================
5
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6
 *
7
 * Developed at SunPro, a Sun Microsystems, Inc. business.
8
 * Permission to use, copy, modify, and distribute this
9
 * software is freely granted, provided that this notice
10
 * is preserved.
11
 * ====================================================
12
 */
13
 
14
#if defined(LIBM_SCCS) && !defined(lint)
15
static char rcsid[] = "$Id: e_acosh.c,v 1.6 1994/08/18 23:04:54 jtc Exp $";
16
#endif
17
 
18
/* __ieee754_acosh(x)
19
 * Method :
20
 *	Based on
21
 *		acosh(x) = log [ x + sqrt(x*x-1) ]
22
 *	we have
23
 *		acosh(x) := log(x)+ln2,	if x is large; else
24
 *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
25
 *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
26
 *
27
 * Special cases:
28
 *	acosh(x) is NaN with signal if x<1.
29
 *	acosh(NaN) is NaN without signal.
30
 */
31
 
32
#include "math.h"
33
#include "math_private.h"
34
 
35
#ifdef __STDC__
36
static const double
37
#else
38
static double
39
#endif
40
one	= 1.0,
41
ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
42
 
43
#ifdef __STDC__
44
	double __ieee754_acosh(double x)
45
#else
46
	double __ieee754_acosh(x)
47
	double x;
48
#endif
49
{
50
	double t;
51
	int32_t hx;
52
	u_int32_t lx;
53
	EXTRACT_WORDS(hx,lx,x);
54
	if(hx<0x3ff00000) {		/* x < 1 */
55
	    return (x-x)/(x-x);
56
	} else if(hx >=0x41b00000) {	/* x > 2**28 */
57
	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
58
	        return x+x;
59
	    } else
60
		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
61
	} else if(((hx-0x3ff00000)|lx)==0) {
62
	    return 0.0;			/* acosh(1) = 0 */
63
	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
64
	    t=x*x;
65
	    return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
66
	} else {			/* 1
67
	    t = x-one;
68
	    return log1p(t+sqrt(2.0*t+t*t));
69
	}
70
}