@@ -11,9 +11,12 @@ using std::isinf;
11
11
#include < stdio.h>
12
12
#include < stdint.h>
13
13
#include < stdlib.h>
14
+ #include < memory>
15
+ #include < limits>
14
16
15
17
#include " xsf/binom.h"
16
18
#include " xsf/lambertw.h"
19
+ #include " sf_error.h"
17
20
18
21
19
22
/* Stirling numbers of the second kind
@@ -32,53 +35,56 @@ using std::isinf;
32
35
33
36
// Dynamic programming
34
37
35
- double _stirling2_dp (double n, double k){
38
+ double _stirling2_dp (double n, double k) {
36
39
if ((n == 0 && k == 0 ) || (n==1 && k==1 )) {
37
40
return 1 .;
38
41
}
39
- if (k <= 0 || k > n || n < 0 ){
42
+ if (k <= 0 || k > n || n < 0 ) {
40
43
return 0 .;
41
44
}
42
45
int arraySize = k <= n - k + 1 ? k : n - k + 1 ;
43
- double *curr = (double *) malloc (arraySize * sizeof (double ));
44
- for (int i = 0 ; i < arraySize; i++){
46
+ auto curr = std::unique_ptr<double []>{new (std::nothrow) double [arraySize]};
47
+ if (curr == nullptr ) {
48
+ sf_error (" stirling2" , SF_ERROR_MEMORY, NULL );
49
+ return std::numeric_limits<double >::quiet_NaN ();
50
+ }
51
+ for (int i = 0 ; i < arraySize; i++) {
45
52
curr[i] = 1 .;
46
53
}
47
54
if (k <= n - k + 1 ) {
48
- for (int i = 1 ; i < n - k + 1 ; i++){
49
- for (int j = 1 ; j < k; j++){
55
+ for (int i = 1 ; i < n - k + 1 ; i++) {
56
+ for (int j = 1 ; j < k; j++) {
50
57
curr[j] = (j + 1 ) * curr[j] + curr[j - 1 ];
51
- if (isinf (curr[j])){
52
- free (curr );
58
+ if (isinf (curr[j])) {
59
+ sf_error ( " stirling2 " , SF_ERROR_OVERFLOW, NULL );
53
60
return INFINITY; // numeric overflow
54
61
}
55
62
}
56
63
}
57
64
} else {
58
- for (int i = 1 ; i < k; i++){
59
- for (int j = 1 ; j < n - k + 1 ; j++){
65
+ for (int i = 1 ; i < k; i++) {
66
+ for (int j = 1 ; j < n - k + 1 ; j++) {
60
67
curr[j] = (i + 1 ) * curr[j - 1 ] + curr[j];
61
- if (isinf (curr[j])){
62
- free (curr );
68
+ if (isinf (curr[j])) {
69
+ sf_error ( " stirling2 " , SF_ERROR_OVERFLOW, NULL );
63
70
return INFINITY; // numeric overflow
64
71
}
65
72
}
66
73
}
67
74
}
68
75
double output = curr[arraySize - 1 ];
69
- free (curr);
70
76
return output;
71
77
}
72
78
73
79
74
80
75
81
// second order Temme approximation
76
82
77
- double _stirling2_temme (double n, double k){
78
- if ((n == k && n >= 0 ) || (n > 0 && k==1 )){
83
+ double _stirling2_temme (double n, double k) {
84
+ if ((n == k && n >= 0 ) || (n > 0 && k==1 )) {
79
85
return 1 .;
80
86
}
81
- if (k <= 0 || k > n || n < 0 ){
87
+ if (k <= 0 || k > n || n < 0 ) {
82
88
return 0 .;
83
89
}
84
90
double mu = (double )k / (double )n;
0 commit comments