There is a bug in error estimation for gsl_sf_exp_mult_e10_e:
gsl_sf_exp_e10_e(10000) gives error estimate 3.911404048e-11
but
gsl_sf_exp_mult_e10_e(10000, 1) gives 3.911012947e-15
so it looks like just multiplying by 1 gives us 4 more digits of
precision. True error is about 1.6884e-11 (using MPFR
http://ex-cs.sist.ac.jp/~tkouya/try_mpfr.html), so the error is
clearly underestimated. Changing line:
const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(ly);
to
const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(x);
seems to fix the problem, but I'm not sure this is correct.
Also, there is a discrepancy in result format between exp_e10 and
exp_mult_e10 for x = 100 and y = 1.
The former returns 2.688 with e10=43, the latter 2.688e43 with e10=0.
It would be nicer to return the first result in both cases, but since
both are correct it might not be worth fixing.
Attached is a file reproducing both cases.
Szymon Jaroszewicz
#include <stdio.h>
#include <gsl/gsl_sf_exp.h>
void test_x_y_e10(double x, double y)
{
gsl_sf_result_e10 result;
int status;
printf("\nx=%f, y=%f\n", x, y);
status = gsl_sf_exp_e10_e(x, &result);
printf("gsl_exp_e10: (%.15g +/- %.15g) * 10^%d\n", result.val, result.err, result.e10);
status = gsl_sf_exp_mult_e10_e(x, y, &result);
printf("gsl_exp_mult_e10: (%.15g +/- %.15g) * 10^%d\n", result.val, result.err, result.e10);
}
int
main(void)
{
test_x_y_e10(10000, 1);
test_x_y_e10(100, 1);
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
