|
1 #include "Python.h" |
|
2 |
|
3 #ifndef HAVE_HYPOT |
|
4 double hypot(double x, double y) |
|
5 { |
|
6 double yx; |
|
7 |
|
8 x = fabs(x); |
|
9 y = fabs(y); |
|
10 if (x < y) { |
|
11 double temp = x; |
|
12 x = y; |
|
13 y = temp; |
|
14 } |
|
15 if (x == 0.) |
|
16 return 0.; |
|
17 else { |
|
18 yx = y/x; |
|
19 return x*sqrt(1.+yx*yx); |
|
20 } |
|
21 } |
|
22 #endif /* HAVE_HYPOT */ |
|
23 |
|
24 #ifndef HAVE_COPYSIGN |
|
25 static double |
|
26 copysign(double x, double y) |
|
27 { |
|
28 /* use atan2 to distinguish -0. from 0. */ |
|
29 if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) { |
|
30 return fabs(x); |
|
31 } else { |
|
32 return -fabs(x); |
|
33 } |
|
34 } |
|
35 #endif /* HAVE_COPYSIGN */ |
|
36 |
|
37 #ifndef HAVE_LOG1P |
|
38 #include <float.h> |
|
39 |
|
40 double |
|
41 log1p(double x) |
|
42 { |
|
43 /* For x small, we use the following approach. Let y be the nearest |
|
44 float to 1+x, then |
|
45 |
|
46 1+x = y * (1 - (y-1-x)/y) |
|
47 |
|
48 so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny, |
|
49 the second term is well approximated by (y-1-x)/y. If abs(x) >= |
|
50 DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest |
|
51 then y-1-x will be exactly representable, and is computed exactly |
|
52 by (y-1)-x. |
|
53 |
|
54 If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be |
|
55 round-to-nearest then this method is slightly dangerous: 1+x could |
|
56 be rounded up to 1+DBL_EPSILON instead of down to 1, and in that |
|
57 case y-1-x will not be exactly representable any more and the |
|
58 result can be off by many ulps. But this is easily fixed: for a |
|
59 floating-point number |x| < DBL_EPSILON/2., the closest |
|
60 floating-point number to log(1+x) is exactly x. |
|
61 */ |
|
62 |
|
63 double y; |
|
64 if (fabs(x) < DBL_EPSILON/2.) { |
|
65 return x; |
|
66 } else if (-0.5 <= x && x <= 1.) { |
|
67 /* WARNING: it's possible than an overeager compiler |
|
68 will incorrectly optimize the following two lines |
|
69 to the equivalent of "return log(1.+x)". If this |
|
70 happens, then results from log1p will be inaccurate |
|
71 for small x. */ |
|
72 y = 1.+x; |
|
73 return log(y)-((y-1.)-x)/y; |
|
74 } else { |
|
75 /* NaNs and infinities should end up here */ |
|
76 return log(1.+x); |
|
77 } |
|
78 } |
|
79 #endif /* HAVE_LOG1P */ |
|
80 |
|
81 /* |
|
82 * ==================================================== |
|
83 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
|
84 * |
|
85 * Developed at SunPro, a Sun Microsystems, Inc. business. |
|
86 * Permission to use, copy, modify, and distribute this |
|
87 * software is freely granted, provided that this notice |
|
88 * is preserved. |
|
89 * ==================================================== |
|
90 */ |
|
91 |
|
92 static const double ln2 = 6.93147180559945286227E-01; |
|
93 static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */ |
|
94 static const double two_pow_p28 = 268435456.0; /* 2**28 */ |
|
95 static const double zero = 0.0; |
|
96 |
|
97 /* asinh(x) |
|
98 * Method : |
|
99 * Based on |
|
100 * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] |
|
101 * we have |
|
102 * asinh(x) := x if 1+x*x=1, |
|
103 * := sign(x)*(log(x)+ln2)) for large |x|, else |
|
104 * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else |
|
105 * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) |
|
106 */ |
|
107 |
|
108 #ifndef HAVE_ASINH |
|
109 double |
|
110 asinh(double x) |
|
111 { |
|
112 double w; |
|
113 double absx = fabs(x); |
|
114 |
|
115 if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) { |
|
116 return x+x; |
|
117 } |
|
118 if (absx < two_pow_m28) { /* |x| < 2**-28 */ |
|
119 return x; /* return x inexact except 0 */ |
|
120 } |
|
121 if (absx > two_pow_p28) { /* |x| > 2**28 */ |
|
122 w = log(absx)+ln2; |
|
123 } |
|
124 else if (absx > 2.0) { /* 2 < |x| < 2**28 */ |
|
125 w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx)); |
|
126 } |
|
127 else { /* 2**-28 <= |x| < 2= */ |
|
128 double t = x*x; |
|
129 w = log1p(absx + t / (1.0 + sqrt(1.0 + t))); |
|
130 } |
|
131 return copysign(w, x); |
|
132 |
|
133 } |
|
134 #endif /* HAVE_ASINH */ |
|
135 |
|
136 /* acosh(x) |
|
137 * Method : |
|
138 * Based on |
|
139 * acosh(x) = log [ x + sqrt(x*x-1) ] |
|
140 * we have |
|
141 * acosh(x) := log(x)+ln2, if x is large; else |
|
142 * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else |
|
143 * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. |
|
144 * |
|
145 * Special cases: |
|
146 * acosh(x) is NaN with signal if x<1. |
|
147 * acosh(NaN) is NaN without signal. |
|
148 */ |
|
149 |
|
150 #ifndef HAVE_ACOSH |
|
151 double |
|
152 acosh(double x) |
|
153 { |
|
154 if (Py_IS_NAN(x)) { |
|
155 return x+x; |
|
156 } |
|
157 if (x < 1.) { /* x < 1; return a signaling NaN */ |
|
158 errno = EDOM; |
|
159 #ifdef Py_NAN |
|
160 return Py_NAN; |
|
161 #else |
|
162 return (x-x)/(x-x); |
|
163 #endif |
|
164 } |
|
165 else if (x >= two_pow_p28) { /* x > 2**28 */ |
|
166 if (Py_IS_INFINITY(x)) { |
|
167 return x+x; |
|
168 } else { |
|
169 return log(x)+ln2; /* acosh(huge)=log(2x) */ |
|
170 } |
|
171 } |
|
172 else if (x == 1.) { |
|
173 return 0.0; /* acosh(1) = 0 */ |
|
174 } |
|
175 else if (x > 2.) { /* 2 < x < 2**28 */ |
|
176 double t = x*x; |
|
177 return log(2.0*x - 1.0 / (x + sqrt(t - 1.0))); |
|
178 } |
|
179 else { /* 1 < x <= 2 */ |
|
180 double t = x - 1.0; |
|
181 return log1p(t + sqrt(2.0*t + t*t)); |
|
182 } |
|
183 } |
|
184 #endif /* HAVE_ACOSH */ |
|
185 |
|
186 /* atanh(x) |
|
187 * Method : |
|
188 * 1.Reduced x to positive by atanh(-x) = -atanh(x) |
|
189 * 2.For x>=0.5 |
|
190 * 1 2x x |
|
191 * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) |
|
192 * 2 1 - x 1 - x |
|
193 * |
|
194 * For x<0.5 |
|
195 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) |
|
196 * |
|
197 * Special cases: |
|
198 * atanh(x) is NaN if |x| >= 1 with signal; |
|
199 * atanh(NaN) is that NaN with no signal; |
|
200 * |
|
201 */ |
|
202 |
|
203 #ifndef HAVE_ATANH |
|
204 double |
|
205 atanh(double x) |
|
206 { |
|
207 double absx; |
|
208 double t; |
|
209 |
|
210 if (Py_IS_NAN(x)) { |
|
211 return x+x; |
|
212 } |
|
213 absx = fabs(x); |
|
214 if (absx >= 1.) { /* |x| >= 1 */ |
|
215 errno = EDOM; |
|
216 #ifdef Py_NAN |
|
217 return Py_NAN; |
|
218 #else |
|
219 return x/zero; |
|
220 #endif |
|
221 } |
|
222 if (absx < two_pow_m28) { /* |x| < 2**-28 */ |
|
223 return x; |
|
224 } |
|
225 if (absx < 0.5) { /* |x| < 0.5 */ |
|
226 t = absx+absx; |
|
227 t = 0.5 * log1p(t + t*absx / (1.0 - absx)); |
|
228 } |
|
229 else { /* 0.5 <= |x| <= 1.0 */ |
|
230 t = 0.5 * log1p((absx + absx) / (1.0 - absx)); |
|
231 } |
|
232 return copysign(t, x); |
|
233 } |
|
234 #endif /* HAVE_ATANH */ |