source: trunk/sys/libm/e_acosh.c @ 186

Last change on this file since 186 was 1, checked in by alain, 7 years ago

First import

File size: 1.6 KB
Line 
1
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
15/* __ieee754_acosh(x)
16 * Method :
17 *      Based on
18 *              acosh(x) = log [ x + sqrt(x*x-1) ]
19 *      we have
20 *              acosh(x) := log(x)+ln2, if x is large; else
21 *              acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
22 *              acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
23 *
24 * Special cases:
25 *      acosh(x) is NaN with signal if x<1.
26 *      acosh(NaN) is NaN without signal.
27 */
28
29#include <libm/fdlibm.h>
30
31#ifdef __STDC__
32static const double 
33#else
34static double 
35#endif
36one     = 1.0,
37ln2     = 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
38
39#ifdef __STDC__
40        double __ieee754_acosh(double x)
41#else
42        double __ieee754_acosh(x)
43        double x;
44#endif
45{       
46        double t;
47        int n0,hx;
48        n0 = ((*(int*)&one)>>29)^1;
49        hx = *(n0+(int*)&x);
50        if(hx<0x3ff00000) {             /* x < 1 */
51            return (x-x)/(x-x);
52        } else if(hx >=0x41b00000) {    /* x > 2**28 */
53            if(hx >=0x7ff00000) {       /* x is inf of NaN */
54                return x+x;
55            } else 
56                return __ieee754_log(x)+ln2;    /* acosh(huge)=log(2x) */
57        } else if(((hx-0x3ff00000)|*(1-n0+(int*)&x))==0) {
58            return 0.0;                 /* acosh(1) = 0 */
59        } else if (hx > 0x40000000) {   /* 2**28 > x > 2 */
60            t=x*x;
61            return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
62        } else {                        /* 1<x<2 */
63            t = x-one;
64            return log1p(t+sqrt(2.0*t+t*t));
65        }
66}
Note: See TracBrowser for help on using the repository browser.