Skip to content

Commit 7f0c47d

Browse files
steppisteppi
authored
Merge pull request #11 from steppi/scipy-updates
MAINT: Incorporate relevant updates which have been made in SciPy since migration
2 parents cd55a14 + 401d2e0 commit 7f0c47d

File tree

10 files changed

+621
-216
lines changed

10 files changed

+621
-216
lines changed

‎include/xsf/amos.h

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,9 +5,12 @@
55

66
namespace xsf {
77

8+
//
9+
// Return sf_error equivalents for AMOS ierr values.
10+
// 'ierr' refers to the last parameter of the AMOS functions
11+
// airy(), besh(), besi(), besj(), besk(), besy(), and biry().
12+
//
813
inline sf_error_t ierr_to_sferr(int nz, int ierr) {
9-
/* Return sf_error equivalents for amos ierr values */
10-
1114
if (nz != 0) {
1215
return SF_ERROR_UNDERFLOW;
1316
}
@@ -23,6 +26,8 @@ inline sf_error_t ierr_to_sferr(int nz, int ierr) {
2326
return SF_ERROR_NO_RESULT;
2427
case 5: /* Algorithm termination condition not met */
2528
return SF_ERROR_NO_RESULT;
29+
case 6: /* Memory allocation failed */
30+
return SF_ERROR_MEMORY;
2631
}
2732

2833
return SF_ERROR_OK;

‎include/xsf/amos/amos.h

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,7 @@
9696

9797
#include <math.h>
9898
#include <complex.h>
99+
#include <memory> // unique_ptr
99100

100101
namespace xsf {
101102
namespace amos {
@@ -2362,6 +2363,7 @@ inline int besy(
23622363
// CANCE BY ARGUMENT REDUCTION
23632364
// IERR=5, ERROR - NO COMPUTATION,
23642365
// ALGORITHM TERMINATION CONDITION NOT MET
2366+
// IERR=6, Memory allocation failed.
23652367
//
23662368
//***LONG DESCRIPTION
23672369
//
@@ -2450,7 +2452,6 @@ inline int besy(
24502452
std::complex<double> c1, c2, hci, st;
24512453
double elim, exr, exi, ey, tay, xx, yy, ascle, rtol, atol, tol, aa, bb, r1m5;
24522454
int i, k, k1, k2, nz, nz1, nz2;
2453-
std::complex<double>* cwrk = new std::complex<double>[n];
24542455

24552456
xx = std::real(z);
24562457
yy = std::imag(z);
@@ -2467,7 +2468,14 @@ inline int besy(
24672468
nz1 = besh(z, fnu, kode, 1, n, cy, ierr);
24682469
if ((*ierr != 0) && (*ierr != 3)) { return 0; }
24692470

2470-
nz2 = besh(z, fnu, kode, 2, n, cwrk, ierr);
2471+
auto cwrk = std::unique_ptr<std::complex<double>[]>
2472+
{new (std::nothrow) std::complex<double>[n]};
2473+
if (cwrk == nullptr) {
2474+
*ierr = 6; // Memory allocation failed.
2475+
return 0;
2476+
}
2477+
2478+
nz2 = besh(z, fnu, kode, 2, n, cwrk.get(), ierr);
24712479
if ((*ierr != 0) && (*ierr != 3)) { return 0; }
24722480

24732481
nz = (nz1 > nz2 ? nz2 : nz1);
@@ -2531,7 +2539,6 @@ inline int besy(
25312539
if ((st == 0.0) && (ey == 0.0)) { nz += 1; }
25322540
}
25332541

2534-
delete [] cwrk;
25352542
return nz;
25362543
}
25372544

‎include/xsf/cephes/chdtr.h

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -162,9 +162,10 @@ namespace xsf {
162162
namespace cephes {
163163

164164
XSF_HOST_DEVICE inline double chdtrc(double df, double x) {
165-
166-
if (x < 0.0)
167-
return 1.0; /* modified by T. Oliphant */
165+
if (x < 0.0) {
166+
set_error("chdtr", SF_ERROR_DOMAIN, NULL);
167+
return (std::numeric_limits<double>::quiet_NaN());
168+
}
168169
return (igamc(df / 2.0, x / 2.0));
169170
}
170171

‎include/xsf/cephes/igam.h

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -328,6 +328,10 @@ namespace cephes {
328328
XSF_HOST_DEVICE inline double igam(double a, double x) {
329329
double absxma_a;
330330

331+
if (std::isnan(a) || std::isnan(x)) {
332+
return std::numeric_limits<double>::quiet_NaN();
333+
}
334+
331335
if (x < 0 || a < 0) {
332336
set_error("gammainc", SF_ERROR_DOMAIN, NULL);
333337
return std::numeric_limits<double>::quiet_NaN();
@@ -367,6 +371,10 @@ namespace cephes {
367371
XSF_HOST_DEVICE double igamc(double a, double x) {
368372
double absxma_a;
369373

374+
if (std::isnan(a) || std::isnan(x)) {
375+
return std::numeric_limits<double>::quiet_NaN();
376+
}
377+
370378
if (x < 0 || a < 0) {
371379
set_error("gammaincc", SF_ERROR_DOMAIN, NULL);
372380
return std::numeric_limits<double>::quiet_NaN();

‎include/xsf/error.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ typedef enum {
1111
SF_ERROR_DOMAIN, /* out of domain */
1212
SF_ERROR_ARG, /* invalid input parameter */
1313
SF_ERROR_OTHER, /* unclassified error */
14+
SF_ERROR_MEMORY, /* memory allocation failed */
1415
SF_ERROR__LAST
1516
} sf_error_t;
1617

‎include/xsf/log_exp.h

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include <cmath>
44

55
#include "config.h"
6+
#include "error.h"
67

78
namespace xsf {
89

@@ -65,4 +66,22 @@ T log_expit(T x) {
6566
return -std::log1p(std::exp(-x));
6667
};
6768

69+
70+
/* Compute log(1 - exp(x)). */
71+
template <typename T>
72+
T log1mexp(T x) {
73+
if (x > 0) {
74+
set_error("_log1mexp", SF_ERROR_DOMAIN, NULL);
75+
return std::numeric_limits<T>::quiet_NaN();
76+
}
77+
if (x == 0) {
78+
set_error("_log1mexp", SF_ERROR_SINGULAR, NULL);
79+
return -std::numeric_limits<T>::infinity();
80+
}
81+
if (x < -1) {
82+
return std::log1p(-std::exp(x));
83+
}
84+
return std::log(-std::expm1(x));
85+
}
86+
6887
} // namespace xsf

‎include/xsf/mathieu.h

Lines changed: 62 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ T cem_cva(T m, T q) {
1717
int int_m, kd = 1;
1818

1919
if ((m < 0) || (m != floor(m))) {
20-
set_error("cem_cva", SF_ERROR_DOMAIN, NULL);
20+
set_error("mathieu_a", SF_ERROR_DOMAIN, NULL);
2121
return std::numeric_limits<T>::quiet_NaN();
2222
}
2323
int_m = (int) m;
@@ -41,7 +41,7 @@ T sem_cva(T m, T q) {
4141
int int_m, kd = 4;
4242

4343
if ((m <= 0) || (m != floor(m))) {
44-
set_error("cem_cva", SF_ERROR_DOMAIN, NULL);
44+
set_error("mathieu_b", SF_ERROR_DOMAIN, NULL);
4545
return std::numeric_limits<T>::quiet_NaN();
4646
}
4747
int_m = (int) m;
@@ -67,7 +67,7 @@ void cem(T m, T q, T x, T &csf, T &csd) {
6767
if ((m < 0) || (m != floor(m))) {
6868
csf = std::numeric_limits<T>::quiet_NaN();
6969
csd = std::numeric_limits<T>::quiet_NaN();
70-
set_error("cem", SF_ERROR_DOMAIN, NULL);
70+
set_error("mathieu_cem", SF_ERROR_DOMAIN, NULL);
7171
} else {
7272
int_m = (int) m;
7373
if (q < 0) {
@@ -85,7 +85,15 @@ void cem(T m, T q, T x, T &csf, T &csd) {
8585
csd = -sgn * d;
8686
}
8787
} else {
88-
specfun::mtu0(kf, int_m, q, x, &csf, &csd);
88+
using specfun::Status;
89+
Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd);
90+
if (status != Status::OK) {
91+
csf = std::numeric_limits<T>::quiet_NaN();
92+
csd = std::numeric_limits<T>::quiet_NaN();
93+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
94+
: SF_ERROR_OTHER;
95+
set_error("mathieu_cem", sf_error, NULL);
96+
}
8997
}
9098
}
9199
}
@@ -97,7 +105,7 @@ void sem(T m, T q, T x, T &csf, T &csd) {
97105
if ((m < 0) || (m != floor(m))) {
98106
csf = std::numeric_limits<T>::quiet_NaN();
99107
csd = std::numeric_limits<T>::quiet_NaN();
100-
set_error("sem", SF_ERROR_DOMAIN, NULL);
108+
set_error("mathieu_sem", SF_ERROR_DOMAIN, NULL);
101109
} else {
102110
int_m = (int) m;
103111
if (int_m == 0) {
@@ -117,7 +125,15 @@ void sem(T m, T q, T x, T &csf, T &csd) {
117125
csd = -sgn * d;
118126
}
119127
} else {
120-
specfun::mtu0(kf, int_m, q, x, &csf, &csd);
128+
using specfun::Status;
129+
Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd);
130+
if (status != Status::OK) {
131+
csf = std::numeric_limits<T>::quiet_NaN();
132+
csd = std::numeric_limits<T>::quiet_NaN();
133+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
134+
: SF_ERROR_OTHER;
135+
set_error("mathieu_sem", sf_error, NULL);
136+
}
121137
}
122138
}
123139
}
@@ -130,10 +146,18 @@ void mcm1(T m, T q, T x, T &f1r, T &d1r) {
130146
if ((m < 0) || (m != floor(m)) || (q < 0)) {
131147
f1r = std::numeric_limits<T>::quiet_NaN();
132148
d1r = std::numeric_limits<T>::quiet_NaN();
133-
set_error("mcm1", SF_ERROR_DOMAIN, NULL);
149+
set_error("mathieu_modcem1", SF_ERROR_DOMAIN, NULL);
134150
} else {
151+
using specfun::Status;
135152
int_m = (int) m;
136-
specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
153+
Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
154+
if (status != Status::OK) {
155+
f1r = std::numeric_limits<T>::quiet_NaN();
156+
d1r = std::numeric_limits<T>::quiet_NaN();
157+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
158+
: SF_ERROR_OTHER;
159+
set_error("mathieu_modcem1", sf_error, NULL);
160+
}
137161
}
138162
}
139163

@@ -145,10 +169,18 @@ void msm1(T m, T q, T x, T &f1r, T &d1r) {
145169
if ((m < 1) || (m != floor(m)) || (q < 0)) {
146170
f1r = std::numeric_limits<T>::quiet_NaN();
147171
d1r = std::numeric_limits<T>::quiet_NaN();
148-
set_error("msm1", SF_ERROR_DOMAIN, NULL);
172+
set_error("mathieu_modsem1", SF_ERROR_DOMAIN, NULL);
149173
} else {
174+
using specfun::Status;
150175
int_m = (int) m;
151-
specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
176+
Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
177+
if (status != Status::OK) {
178+
f1r = std::numeric_limits<T>::quiet_NaN();
179+
d1r = std::numeric_limits<T>::quiet_NaN();
180+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
181+
: SF_ERROR_OTHER;
182+
set_error("mathieu_modsem1", sf_error, NULL);
183+
}
152184
}
153185
}
154186

@@ -160,10 +192,18 @@ void mcm2(T m, T q, T x, T &f2r, T &d2r) {
160192
if ((m < 0) || (m != floor(m)) || (q < 0)) {
161193
f2r = std::numeric_limits<T>::quiet_NaN();
162194
d2r = std::numeric_limits<T>::quiet_NaN();
163-
set_error("mcm2", SF_ERROR_DOMAIN, NULL);
195+
set_error("mathieu_modcem2", SF_ERROR_DOMAIN, NULL);
164196
} else {
197+
using specfun::Status;
165198
int_m = (int) m;
166-
specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
199+
Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
200+
if (status != Status::OK) {
201+
f2r = std::numeric_limits<T>::quiet_NaN();
202+
d2r = std::numeric_limits<T>::quiet_NaN();
203+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
204+
: SF_ERROR_OTHER;
205+
set_error("mathieu_modcem2", sf_error, NULL);
206+
}
167207
}
168208
}
169209

@@ -175,10 +215,18 @@ void msm2(T m, T q, T x, T &f2r, T &d2r) {
175215
if ((m < 1) || (m != floor(m)) || (q < 0)) {
176216
f2r = std::numeric_limits<T>::quiet_NaN();
177217
d2r = std::numeric_limits<T>::quiet_NaN();
178-
set_error("msm2", SF_ERROR_DOMAIN, NULL);
218+
set_error("mathieu_modsem2", SF_ERROR_DOMAIN, NULL);
179219
} else {
220+
using specfun::Status;
180221
int_m = (int) m;
181-
specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
222+
Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r);
223+
if (status != Status::OK) {
224+
f2r = std::numeric_limits<T>::quiet_NaN();
225+
d2r = std::numeric_limits<T>::quiet_NaN();
226+
sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY
227+
: SF_ERROR_OTHER;
228+
set_error("mathieu_modsem2", sf_error, NULL);
229+
}
182230
}
183231
}
184232

‎include/xsf/par_cyl.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -628,7 +628,7 @@ void pbdv(T v, T x, T &pdf, T &pdd) {
628628
num = std::abs((int) v) + 2;
629629
dv = (T *) malloc(sizeof(T) * 2 * num);
630630
if (dv == NULL) {
631-
set_error("pbdv", SF_ERROR_OTHER, "memory allocation error");
631+
set_error("pbdv", SF_ERROR_MEMORY, "memory allocation error");
632632
pdf = std::numeric_limits<T>::quiet_NaN();
633633
pdd = std::numeric_limits<T>::quiet_NaN();
634634
} else {
@@ -653,7 +653,7 @@ void pbvv(T v, T x, T &pvf, T &pvd) {
653653
num = std::abs((int) v) + 2;
654654
vv = (T *) malloc(sizeof(T) * 2 * num);
655655
if (vv == NULL) {
656-
set_error("pbvv", SF_ERROR_OTHER, "memory allocation error");
656+
set_error("pbvv", SF_ERROR_MEMORY, "memory allocation error");
657657
pvf = std::numeric_limits<T>::quiet_NaN();
658658
pvd = std::numeric_limits<T>::quiet_NaN();
659659
} else {

0 commit comments

Comments
 (0)