Skip to content

Commit e5e12d6

Browse files
steppisteppi
committed
BUG: Make fixes to get tests to pass
1 parent 48d659f commit e5e12d6

File tree

12 files changed

+168
-156
lines changed

12 files changed

+168
-156
lines changed

‎include/xsf/cephes/cbrt.h‎

Lines changed: 70 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -52,80 +52,81 @@ namespace cephes {
5252

5353
namespace detail {
5454

55-
constexpr double CBRT2 = 1.2599210498948731647672;
56-
constexpr double CBRT4 = 1.5874010519681994747517;
57-
constexpr double CBRT2I = 0.79370052598409973737585;
58-
constexpr double CBRT4I = 0.62996052494743658238361;
55+
inline constexpr double CBRT2 = 1.2599210498948731647672;
56+
inline constexpr double CBRT4 = 1.5874010519681994747517;
57+
inline constexpr double CBRT2I = 0.79370052598409973737585;
58+
inline constexpr double CBRT4I = 0.62996052494743658238361;
5959

60-
XSF_HOST_DEVICE inline double cbrt(double x) {
61-
int e, rem, sign;
62-
double z;
60+
}
6361

64-
if (!std::isfinite(x)) {
65-
return x;
66-
}
67-
if (x == 0) {
68-
return (x);
69-
}
70-
if (x > 0) {
71-
sign = 1;
72-
} else {
73-
sign = -1;
74-
x = -x;
75-
}
62+
XSF_HOST_DEVICE inline double cbrt(double x) {
63+
int e, rem, sign;
64+
double z;
7665

77-
z = x;
78-
/* extract power of 2, leaving
79-
* mantissa between 0.5 and 1
80-
*/
81-
x = std::frexp(x, &e);
66+
if (!std::isfinite(x)) {
67+
return x;
68+
}
69+
if (x == 0) {
70+
return (x);
71+
}
72+
if (x > 0) {
73+
sign = 1;
74+
} else {
75+
sign = -1;
76+
x = -x;
77+
}
78+
79+
z = x;
80+
/* extract power of 2, leaving
81+
* mantissa between 0.5 and 1
82+
*/
83+
x = std::frexp(x, &e);
84+
85+
/* Approximate cube root of number between .5 and 1,
86+
* peak relative error = 9.2e-6
87+
*/
88+
x = (((-1.3466110473359520655053e-1 * x + 5.4664601366395524503440e-1) * x - 9.5438224771509446525043e-1) *
89+
x +
90+
1.1399983354717293273738e0) *
91+
x +
92+
4.0238979564544752126924e-1;
93+
94+
/* exponent divided by 3 */
95+
if (e >= 0) {
96+
rem = e;
97+
e /= 3;
98+
rem -= 3 * e;
99+
if (rem == 1) {
100+
x *= detail::CBRT2;
101+
} else if (rem == 2) {
102+
x *= detail::CBRT4;
103+
}
104+
}
105+
/* argument less than 1 */
106+
else {
107+
e = -e;
108+
rem = e;
109+
e /= 3;
110+
rem -= 3 * e;
111+
if (rem == 1) {
112+
x *= detail::CBRT2I;
113+
} else if (rem == 2) {
114+
x *= detail::CBRT4I;
115+
}
116+
e = -e;
117+
}
82118

83-
/* Approximate cube root of number between .5 and 1,
84-
* peak relative error = 9.2e-6
85-
*/
86-
x = (((-1.3466110473359520655053e-1 * x + 5.4664601366395524503440e-1) * x - 9.5438224771509446525043e-1) *
87-
x +
88-
1.1399983354717293273738e0) *
89-
x +
90-
4.0238979564544752126924e-1;
119+
/* multiply by power of 2 */
120+
x = std::ldexp(x, e);
121+
122+
/* Newton iteration */
123+
x -= (x - (z / (x * x))) * 0.33333333333333333333;
124+
x -= (x - (z / (x * x))) * 0.33333333333333333333;
91125

92-
/* exponent divided by 3 */
93-
if (e >= 0) {
94-
rem = e;
95-
e /= 3;
96-
rem -= 3 * e;
97-
if (rem == 1) {
98-
x *= CBRT2;
99-
} else if (rem == 2) {
100-
x *= CBRT4;
101-
}
102-
}
103-
/* argument less than 1 */
104-
else {
105-
e = -e;
106-
rem = e;
107-
e /= 3;
108-
rem -= 3 * e;
109-
if (rem == 1) {
110-
x *= CBRT2I;
111-
} else if (rem == 2) {
112-
x *= CBRT4I;
113-
}
114-
e = -e;
115-
}
116-
117-
/* multiply by power of 2 */
118-
x = std::ldexp(x, e);
119-
120-
/* Newton iteration */
121-
x -= (x - (z / (x * x))) * 0.33333333333333333333;
122-
x -= (x - (z / (x * x))) * 0.33333333333333333333;
123-
124-
if (sign < 0)
125-
x = -x;
126-
return (x);
127-
}
128-
} // namespace detail
126+
if (sign < 0)
127+
x = -x;
128+
return (x);
129+
}
129130

130131
} // namespace cephes
131132
} // namespace xsf

‎include/xsf/cephes/exp10.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,8 @@
5757
#include "../config.h"
5858
#include "../error.h"
5959

60+
#include "polevl.h"
61+
6062
namespace xsf {
6163
namespace cephes {
6264

‎include/xsf/cephes/gdtr.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,7 @@ namespace cephes {
111111
XSF_HOST_DEVICE inline double gdtr(double a, double b, double x) {
112112

113113
if (x < 0.0) {
114-
sf_error("gdtr", SF_ERROR_DOMAIN, NULL);
114+
set_error("gdtr", SF_ERROR_DOMAIN, NULL);
115115
return (std::numeric_limits<double>::quiet_NaN());
116116
}
117117
return (igam(b, a * x));
@@ -129,7 +129,7 @@ namespace cephes {
129129
XSF_HOST_DEVICE inline double gdtri(double a, double b, double y) {
130130

131131
if ((y < 0.0) || (y > 1.0) || (a <= 0.0) || (b < 0.0)) {
132-
sf_error("gdtri", SF_ERROR_DOMAIN, NULL);
132+
set_error("gdtri", SF_ERROR_DOMAIN, NULL);
133133
return (std::numeric_limits<double>::quiet_NaN());
134134
}
135135

‎include/xsf/cephes/ndtr.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -135,6 +135,7 @@
135135
#pragma once
136136

137137
#include "../config.h"
138+
#include "../error.h"
138139

139140
#include "const.h"
140141
#include "polevl.h"

‎include/xsf/cephes/shichi.h‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -151,7 +151,7 @@ namespace cephes {
151151
} // namespace detail
152152

153153
/* Sine and cosine integrals */
154-
XSF_HOST_DEVICE inline int shichi(double x, double *si, double *ci) {
154+
XSF_HOST_DEVICE inline int shichi(double x, double &si, double &ci) {
155155
double k, z, c, s, a, b;
156156
short sign;
157157

@@ -163,8 +163,8 @@ namespace cephes {
163163
}
164164

165165
if (x == 0.0) {
166-
*si = 0.0;
167-
*ci = -std::numeric_limits<double>::infinity();
166+
si = 0.0;
167+
ci = -std::numeric_limits<double>::infinity();
168168
return (0);
169169
}
170170

@@ -216,20 +216,20 @@ namespace cephes {
216216

217217
asymp:
218218
if (x > 1000) {
219-
*si = std::numeric_limits<double>::infinity();
220-
*ci = std::numeric_limits<double>::infinity();
219+
si = std::numeric_limits<double>::infinity();
220+
ci = std::numeric_limits<double>::infinity();
221221
} else {
222222
/* Asymptotic expansions
223223
* http://functions.wolfram.com/GammaBetaErf/CoshIntegral/06/02/
224224
* http://functions.wolfram.com/GammaBetaErf/SinhIntegral/06/02/0001/
225225
*/
226226
a = detail::hyp3f0(0.5, 1, 1, 4.0 / (x * x));
227227
b = detail::hyp3f0(1, 1, 1.5, 4.0 / (x * x));
228-
*si = std::cosh(x) / x * a + std::sinh(x) / (x * x) * b;
229-
*ci = std::sinh(x) / x * a + std::cosh(x) / (x * x) * b;
228+
si = std::cosh(x) / x * a + std::sinh(x) / (x * x) * b;
229+
ci = std::sinh(x) / x * a + std::cosh(x) / (x * x) * b;
230230
}
231231
if (sign) {
232-
*si = -*si;
232+
si = -si;
233233
}
234234
return 0;
235235

@@ -238,9 +238,9 @@ namespace cephes {
238238
s = -s;
239239
}
240240

241-
*si = s;
241+
si = s;
242242

243-
*ci = detail::SCIPY_EULER + std::log(x) + c;
243+
ci = detail::SCIPY_EULER + std::log(x) + c;
244244
return (0);
245245
}
246246

‎include/xsf/cephes/sici.h‎

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -138,7 +138,7 @@ namespace cephes {
138138

139139
} // namespace detail
140140

141-
XSF_HOST_DEVICE inline int sici(double x, double *si, double *ci) {
141+
XSF_HOST_DEVICE inline int sici(double x, double &si, double &ci) {
142142
double z, c, s, f, g;
143143
short sign;
144144

@@ -150,24 +150,24 @@ namespace cephes {
150150
}
151151

152152
if (x == 0.0) {
153-
*si = 0.0;
154-
*ci = -std::numeric_limits<double>::infinity();
153+
si = 0.0;
154+
ci = -std::numeric_limits<double>::infinity();
155155
return (0);
156156
}
157157

158158
if (x > 1.0e9) {
159159
if (std::isinf(x)) {
160160
if (sign == -1) {
161-
*si = -M_PI_2;
162-
*ci = std::numeric_limits<double>::quiet_NaN();
161+
si = -M_PI_2;
162+
ci = std::numeric_limits<double>::quiet_NaN();
163163
} else {
164-
*si = M_PI_2;
165-
*ci = 0;
164+
si = M_PI_2;
165+
ci = 0;
166166
}
167167
return 0;
168168
}
169-
*si = M_PI_2 - std::cos(x) / x;
170-
*ci = std::sin(x) / x;
169+
si = M_PI_2 - std::cos(x) / x;
170+
ci = std::sin(x) / x;
171171
}
172172

173173
if (x > 4.0) {
@@ -181,8 +181,8 @@ namespace cephes {
181181
if (sign) {
182182
s = -s;
183183
}
184-
*si = s;
185-
*ci = detail::SCIPY_EULER + std::log(x) + c; /* real part if x < 0 */
184+
si = s;
185+
ci = detail::SCIPY_EULER + std::log(x) + c; /* real part if x < 0 */
186186
return (0);
187187

188188
/* The auxiliary functions are:
@@ -211,11 +211,11 @@ namespace cephes {
211211
f = polevl(z, detail::sici_FN8, 8) / (x * p1evl(z, detail::sici_FD8, 8));
212212
g = z * polevl(z, detail::sici_GN8, 8) / p1evl(z, detail::sici_GD8, 9);
213213
}
214-
*si = M_PI_2 - f * c - g * s;
214+
si = M_PI_2 - f * c - g * s;
215215
if (sign) {
216-
*si = -(*si);
216+
si = -(si);
217217
}
218-
*ci = f * s - g * c;
218+
ci = f * s - g * c;
219219

220220
return (0);
221221
}

‎include/xsf/exp.h‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,11 @@
11
#pragma once
22

3-
#include "xsf/cephes/exp10.h"
4-
#include "xsf/cephes/exp2.h"
3+
#include "config.h"
4+
#include "trig.h"
5+
6+
#include "cephes/exp10.h"
7+
#include "cephes/exp2.h"
8+
#include "cephes/unity.h"
59

610
namespace xsf {
711

@@ -45,12 +49,12 @@ inline std::complex<float> expm1(std::complex<float> z) {
4549
return static_cast<std::complex<float>>(expm1(static_cast<std::complex<double>>(z)));
4650
}
4751

48-
double exp2(double x) { return cephes::exp2(x); }
52+
inline double exp2(double x) { return cephes::exp2(x); }
4953

50-
float exp2(float x) { return exp2(static_cast<double>(x)); }
54+
inline float exp2(float x) { return exp2(static_cast<double>(x)); }
5155

52-
double exp10(double x) { return cephes::exp10(x); }
56+
inline double exp10(double x) { return cephes::exp10(x); }
5357

54-
float exp10(float x) { return exp10(static_cast<double>(x)); }
58+
inline float exp10(float x) { return exp10(static_cast<double>(x)); }
5559

5660
} // namespace xsf

‎include/xsf/mathieu.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22

3+
#include "error.h"
34
#include "specfun/specfun.h"
45

56
namespace xsf {

‎include/xsf/par_cyl.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22

3+
#include "error.h"
34
#include "specfun/specfun.h"
45

56
namespace xsf {

0 commit comments

Comments
 (0)