overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// ArcHyperbolicDD.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// long double acoshl(long double x);
29// long double asinhl(long double x);
30// long double atanhl(long double x);
31//
32
33#include "math.h"
34#include "fp_private.h"
35#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
36#include "DD.h"
37
38
39static const DblDbl kLn2DD = {{0x3fe62e42, 0xfefa39ef,
40 0x3c7abc9e, 0x3b39803f}};
41static const Float kInfF = {0x7f800000};
42
43long double fabsl(long double x);
44long double logl(long double x);
45long double sqrtl(long double x);
46long double log1pl(long double x);
47long double copysignl(long double x, long double y);
48
49
50long double acoshl(long double x)
51{
52 DblDbl xDD;
53 long double y;
54 uint32_t xhi;
55 double fpenv;
56
57 xDD.f = x;
58 xhi = xDD.i[0];
59
60 // return x if x is a NaN or x is +INF
61 if ((xDD.d[0] != xDD.d[0]) || (xDD.d[0] == kInfF.f))
62 return x;
63
64 FEGETENVD(fpenv);
65 FESETENVD(0.0);
66
67 // calculate ArcCosh for x >= 1.0
68 if (x >= 1.0L) {
69 if (xhi > 0x5fefffffu) { // huge case (x >= 2.0^512)
70 y = logl(x) + kLn2DD.f;
71 FESETENVD(fpenv);
72 return y;
73 }
74 else if (xhi > 0x3ff40000u) { // intermediate case (x > 1.25)
75 y = logl(2.0L*x - 1.0L/(x + sqrtl(x*x -1.0L)));
76 FESETENVD(fpenv);
77 return y;
78 }
79 else { // small case (x <= 1.25)
80 y = x - 1.0L;
81 y = y + sqrtl(y * (2.0L + y));
82 y = log1pl(y);
83 FESETENVD(fpenv);
84 return y;
85 }
86 }
87
88 xDD.i[0] = 0x7ff80000u; // return value is NaN
89 xDD.i[1] = 0u;
90 xDD.i[2] = 0x7ff80000u;
91 xDD.i[3] = 0u;
92 FESETENVD(fpenv);
93 return (xDD.f); // return NaN
94}
95
96
97long double asinhl(long double x)
98{
99 DblDbl xDD;
100 long double absx, y;
101 uint32_t xhi;
102 double fpenv;
103
104 xDD.f = x;
105 absx = fabsl(x);
106 xhi = xDD.i[0] & 0x7fffffffu;
107
108 if (xhi >= 0x7ff00000u) // NaN or INF is returned
109 return x;
110
111 if (xhi < 0x29900000u) // absx < 2.0^(-358)
112 return x;
113
114 FEGETENVD(fpenv);
115 FESETENVD(0.0);
116
117 // x is normal; calculate ArcSinh(x) based on magnitude of |x|
118 if (xhi <= 0x3ff55555u) { // absx <= 4.0/3.0
119 y = 1.0L / absx;
120 y = log1pl(absx + absx/(y + sqrtl(1.0L + y*y)));
121 }
122 else if (xhi < 0x5ff00000u) // absx < 2.0^512
123 y = logl(2.0L*absx + 1.0/(absx + sqrtl(1.0L + absx*absx)));
124 else // absx >= 2.0^512
125 y = logl(absx) + kLn2DD.f;
126
127 y = copysignl(y,x);
128 FESETENVD(fpenv);
129 return y; // fix sign of result
130}
131
132
133long double atanhl(long double x)
134{
135 DblDbl xDD;
136 long double absx, y;
137 uint32_t xhi;
138 double fpenv;
139
140 xDD.f = x;
141 absx = fabsl(x); // absx = |x|
142 xhi = xDD.i[0] & 0x7fffffffu;
143
144 if (xDD.d[0] != xDD.d[0]) // NaN is returned
145 return x;
146
147 if (xhi < 0x29900000u) { // absx < 2.0^(-358)
148 return x;
149 }
150
151 FEGETENVD(fpenv);
152 FESETENVD(0.0);
153
154 if (absx <= 1.0L) { // valid argument
155 y = 0.5L * log1pl(2.0L*absx / (1.0L - absx));
156 y = copysignl(y,x);
157 FESETENVD(fpenv);
158 return y;
159 }
160 else { // invalid argument
161 y = log1pl(-absx);
162 FESETENVD(fpenv);
163 return y;
164 }
165}
166#endif