Skip to content

Commit

Permalink
Update C files from GSL 2.7.1 to 2.8
Browse files Browse the repository at this point in the history
  • Loading branch information
@watanabe-j
watanabe-j committed Jul 19, 2024
1 parent 6ed02f0 commit 4dfda92
Show file tree
Hide file tree
Showing 6 changed files with 45 additions and 23 deletions.
19 changes: 19 additions & 0 deletions src/gsl/integration/gsl_integration.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -390,6 +390,25 @@ int gsl_integration_qagiu (gsl_function * f,
// int gsl_integration_fixed(const gsl_function * func, double * result, // edited for qfratio
// const gsl_integration_fixed_workspace * w); // edited for qfratio

// /* Lebedev quadrature */ // edited for qfratio

// typedef struct // edited for qfratio
// { // edited for qfratio
// size_t n; /* number of nodes/weights */ // edited for qfratio
// double *weights; /* quadrature weights */ // edited for qfratio
// double *x; /* x quadrature nodes */ // edited for qfratio
// double *y; /* y quadrature nodes */ // edited for qfratio
// double *z; /* z quadrature nodes */ // edited for qfratio
// double *theta; /* theta quadrature nodes */ // edited for qfratio
// double *phi; /* phi quadrature nodes */ // edited for qfratio
// } gsl_integration_lebedev_workspace; // edited for qfratio

// gsl_integration_lebedev_workspace * gsl_integration_lebedev_alloc(const size_t n); // edited for qfratio

// void gsl_integration_lebedev_free(gsl_integration_lebedev_workspace * w); // edited for qfratio

// size_t gsl_integration_lebedev_n(const gsl_integration_lebedev_workspace * w); // edited for qfratio

__END_DECLS

#endif /* __GSL_INTEGRATION_H__ */
3 changes: 2 additions & 1 deletion src/gsl/specfunc/exp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
/* Author: G. Jungman */

#include <config.h>
#include <stdlib.h>
#include "../gsl_math.h" // edited for qfratio
#include "../err/gsl_errno.h" // edited for qfratio
#include "gsl_sf_gamma.h" // edited for qfratio
Expand Down Expand Up @@ -218,7 +219,7 @@ int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
// const double sy = GSL_SIGN(y); // edited for qfratio
// const int N = (int) floor(l10_val); // edited for qfratio
// const double arg_val = (l10_val - N) * M_LN10; // edited for qfratio
// const double arg_err = 2.0 * GSL_DBL_EPSILON * (fabs(x) + fabs(ly) + M_LN10*fabs(N)); // edited for qfratio
// const double arg_err = 2.0 * GSL_DBL_EPSILON * (fabs(x) + fabs(ly) + M_LN10*abs(N)); // edited for qfratio

// result->val = sy * exp(arg_val); // edited for qfratio
// result->err = arg_err * fabs(result->val); // edited for qfratio
Expand Down
2 changes: 1 addition & 1 deletion src/gsl/specfunc/gamma.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1616,7 +1616,7 @@ int gsl_sf_lnchoose_e(unsigned int n, unsigned int m, gsl_sf_result * result)
// prod *= tk; // edited for qfratio
// } // edited for qfratio
// result->val = prod; // edited for qfratio
// result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs(n-m); // edited for qfratio
// result->err = 2.0 * GSL_DBL_EPSILON * prod * fabs((double) (n-m)); // edited for qfratio
// return GSL_SUCCESS; // edited for qfratio
// } // edited for qfratio
// else // edited for qfratio
Expand Down
16 changes: 8 additions & 8 deletions src/gsl/specfunc/hyperg_1F1.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1006,7 +1006,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r
Mn = Mnm1;
}
result->val = Ma/Mn;
result->err = 2.0 * GSL_DBL_EPSILON * (fabs((double)(a)) + 1.0) * fabs(Ma/Mn); // edited for qfratio
result->err = 2.0 * GSL_DBL_EPSILON * (abs(a) + 1.0) * fabs(Ma/Mn);
return stat_CF1;
}
else if(b > a && b < 2*a + x && b > x) {
Expand Down Expand Up @@ -1037,7 +1037,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r
stat_ex = gsl_sf_exp_e(x, &ex); /* 1F1(b,b,x) */
result->val = ex.val * Ma/Mn;
result->err = ex.err * fabs(Ma/Mn);
result->err += 4.0 * GSL_DBL_EPSILON * (fabs((double)(b-a))+1.0) * fabs(result->val); // edited for qfratio
result->err += 4.0 * GSL_DBL_EPSILON * (abs(b-a)+1.0) * fabs(result->val);
return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
}
else if(x >= 0.0) {
Expand All @@ -1060,7 +1060,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r
}
result->val = Mn;
result->err = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
result->err *= fabs((double)(a-b))+1.0; // edited for qfratio
result->err *= abs(a-b)+1.0;
return GSL_SUCCESS;
}
else {
Expand Down Expand Up @@ -1089,7 +1089,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r
Mn = Mnp1;
}
result->val = Mn;
result->err = fabs(Mn) * (1.0 + fabs((double)(a))) * fabs(r_Mn.err / r_Mn.val); // edited for qfratio
result->err = fabs(Mn) * (1.0 + abs(a)) * fabs(r_Mn.err / r_Mn.val);
result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
return GSL_SUCCESS;
}
Expand All @@ -1116,7 +1116,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r
}
result->val = Man;
result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
result->err *= fabs((double)(b-a))+1.0; // edited for qfratio
result->err *= abs(b-a)+1.0;
return GSL_SUCCESS;
}
else {
Expand Down Expand Up @@ -1172,7 +1172,7 @@ hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * r

result->val = Mn;
result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
result->err *= fabs((double)(b-a))+1.0; // edited for qfratio
result->err *= abs(b-a)+1.0;
return GSL_SUCCESS;
}
}
Expand Down Expand Up @@ -1816,11 +1816,11 @@ gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result
/* Standard domain error due to singularity. */
DOMAIN_ERROR(result);
}
else if(x > 100.0 && GSL_MAX_DBL(1.0,fabs((double)(b-a)))*GSL_MAX_DBL(1.0,fabs((double)(1-a))) < 0.5 * x) { // edited for qfratio
else if(x > 100.0 && GSL_MAX_DBL(1.0,abs(b-a))*GSL_MAX_DBL(1.0,abs(1-a)) < 0.5 * x) {
/* x -> +Inf asymptotic */
return hyperg_1F1_asymp_posx(a, b, x, result);
}
else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs((double)(a)))*GSL_MAX_DBL(1.0,fabs((double)(1+a-b))) < 0.5 * fabs(x)) { // edited for qfratio
else if(x < -100.0 && GSL_MAX_DBL(1.0,abs(a))*GSL_MAX_DBL(1.0,abs(1+a-b)) < 0.5 * fabs(x)) {
/* x -> -Inf asymptotic */
return hyperg_1F1_asymp_negx(a, b, x, result);
}
Expand Down
23 changes: 12 additions & 11 deletions src/gsl/specfunc/hyperg_U.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
/* Author: G. Jungman */

#include <config.h>
#include <stdlib.h>
#include "../gsl_math.h" // edited for qfratio
#include "../err/gsl_errno.h" // edited for qfratio
#include "gsl_sf_exp.h" // edited for qfratio
Expand All @@ -40,9 +41,9 @@

#define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON)

#define SERIES_EVAL_OK(a,b,x) ((fabs((double)(a)) < 5 && b < 5 && x < 2.0) || (fabs((double)(a)) < 10 && b < 10 && x < 1.0)) // edited for qfratio
#define SERIES_EVAL_OK(a,b,x) ((fabs((double) a) < 5 && b < 5 && x < 2.0) || (fabs((double) a) < 10 && b < 10 && x < 1.0))

#define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs((double)(a)),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)) // edited for qfratio
#define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs((double) a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x))

/* Log[U(a,2a,x)]
* [Abramowitz+stegun, 13.6.21]
Expand Down Expand Up @@ -287,7 +288,7 @@ hyperg_U_finite_sum(int N, double a, double b, double x, double xeps,

result->val = sum_val * poch.val;
result->err = fabs(sum_val) * poch.err + sum_err * fabs(poch.val);
result->err += fabs(poch.val) * (fabs((double)(N)) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val); // edited for qfratio
result->err += fabs(poch.val) * (abs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val);
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
return stat_poch;
Expand Down Expand Up @@ -873,7 +874,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
}
else if(a == -1) {
result->val = -b + x;
result->err = 2.0 * GSL_DBL_EPSILON * (fabs((double)(b)) + fabs(x)); // edited for qfratio
result->err = 2.0 * GSL_DBL_EPSILON * (abs(b) + fabs(x));
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
result->e10 = 0;
return GSL_SUCCESS;
Expand Down Expand Up @@ -925,7 +926,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
lnm.val = scale_count*lnscale;
lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val);
y.val = Ua;
y.err = 4.0 * GSL_DBL_EPSILON * (fabs((double)(a))+1.0) * fabs(Ua); // edited for qfratio
y.err = 4.0 * GSL_DBL_EPSILON * (abs(a)+1.0) * fabs(Ua);
return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
}
else if(b >= 2.0*a + x) {
Expand Down Expand Up @@ -959,7 +960,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale));
y.val = Ua;
y.err = fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
y.err += 2.0 * GSL_DBL_EPSILON * (fabs((double)(a)) + 1.0) * fabs(Ua); // edited for qfratio
y.err += 2.0 * GSL_DBL_EPSILON * (abs(a) + 1.0) * fabs(Ua);
stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
return GSL_ERROR_SELECT_2(stat_e, stat_1);
}
Expand Down Expand Up @@ -1010,7 +1011,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
double lnscl = -scale_count*log(scale_factor);
double lnpre_val = lnU_target + lnscl;
double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl));
double oUa_err = 2.0 * (fabs((double)(a_target-a)) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua); // edited for qfratio
double oUa_err = 2.0 * (abs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua);
int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err,
1.0/Ua, oUa_err,
result);
Expand Down Expand Up @@ -1054,7 +1055,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
}
Ua1_bck_val = Ua;
Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs((double)(a1-a))+CF1_count+1.0) * fabs(Ua); // edited for qfratio
Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (abs(a1-a)+CF1_count+1.0) * fabs(Ua);
stat_bck = stat_CF1;
}

Expand Down Expand Up @@ -1120,7 +1121,7 @@ hyperg_U_int_bge1(const int a, const int b, const double x,
}
Ua1_for_val = Ua;
Ua1_for_err = fabs(Ua) * fabs(r_Ua.err/r_Ua.val);
Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs((double)(a1-a0))+1.0) * fabs(Ua1_for_val); // edited for qfratio
Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (abs(a1-a0)+1.0) * fabs(Ua1_for_val);
}

/* Now do the matching to produce the final result.
Expand Down Expand Up @@ -1451,7 +1452,7 @@ hyperg_U_bge1(const double a, const double b, const double x,

lnscale = log(scale_factor);
lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale;
lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs((double)(scale_count_for - scale_count_bck))*fabs(lnscale)); // edited for qfratio
lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + abs(scale_count_for - scale_count_bck)*fabs(lnscale));
y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck;
y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val);
stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
Expand Down Expand Up @@ -1680,7 +1681,7 @@ gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x,
int stat_e;
int stat_U = hyperg_U_int_bge1(ap, bp, x, &U);
double ln_pre_val = (1.0-b)*ln_x;
double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs((double)(b))+1.0) * fabs(ln_x); // edited for qfratio
double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (abs(b)+1.0) * fabs(ln_x);
ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */
stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
U.val, U.err,
Expand Down
5 changes: 3 additions & 2 deletions src/gsl/specfunc/poch.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
/* Author: G. Jungman */

#include <config.h>
#include <stdlib.h>
#include "../gsl_math.h" // edited for qfratio
#include "../err/gsl_errno.h" // edited for qfratio
#include "gsl_sf_exp.h" // edited for qfratio
Expand Down Expand Up @@ -166,7 +167,7 @@ pochrel_smallx(const double a, const double x, gsl_sf_result * result)

if(bp == a) {
result->val = dpoch1;
result->err = 2.0 * GSL_DBL_EPSILON * (fabs((double)(incr)) + 1.0) * fabs(result->val); // edited for qfratio
result->err = 2.0 * GSL_DBL_EPSILON * (abs(incr) + 1.0) * fabs(result->val);
return GSL_SUCCESS;
}
else {
Expand All @@ -181,7 +182,7 @@ pochrel_smallx(const double a, const double x, gsl_sf_result * result)
double trig = t1 - t2;
result->val = dpoch1 * (1.0 + x*trig) + trig;
result->err = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));
result->err += 2.0 * GSL_DBL_EPSILON * (fabs((double)(incr)) + 1.0) * fabs(result->val); // edited for qfratio
result->err += 2.0 * GSL_DBL_EPSILON * (abs(incr) + 1.0) * fabs(result->val);
return GSL_SUCCESS;
}
}
Expand Down

0 comments on commit 4dfda92

Please sign in to comment.