-----Original Message-----
From: Max Gaspa
Sent: Wednesday, June 21, 2017 11:26 PM
To: [email protected]
Subject: Re: [Help-gsl] test release for GSL 2.4
[snip]
- char filename[] = "test.dat";
+ FILE *f = tmpfile();
in few words the change imply using tmpfile(). Unfortunately in Windows
(XP is fine! 7 fails) that call is creating a temporary file in C:\ that
is usually not writable for security reason. The call fails and the
pointer f is NULL.
Interestingly, while the call fails with my 64-bit gcc-6.3.0, it succeeds
with the 32-bit gcc-6.3.0. At least, f is not NULL - though it seems that
the file is un-writeable.
Presumably this accounts for the difference in the behaviour between my
32-bit and 64-bit builds.
[snip]
I'd like to change the compilation flags (using SSE2 and or -mfloat-store)
just to see what will happen...and then using MPFR I'd like to understand
better the reason of the failure.
The manual says that gsl_sf_bessel_j2_e() computes ((3/x^2 - 1) * sin(x) - 3
* cos(x)/x)/x
Using MPFR (with a sufficiently large precision) shows that the *exact*
value of that computation (for x = 1048576.0), rounded to 17 decimal digits
is -3.1518539455252413e-007.
When I perform that calculation (with double precision) using 64-bit
gcc-6.3.0, I get a result of -3.1518539455252539e-007.
This is the value "obtained" in the test-suite.log - and it's simply wrong.
Using MPFR, I find that performing that calculation with double (53-bit)
precision, gives -3.1518539455252417e-7.
The culprit for the error is, of course, the gcc-6.3.0 trig functions. For
1048576.0, according to gcc-6.3.0:
sin: 3.3049314002173596e-001
cos: 9.4380839390131155e-001
According to MPFR:
sin: 3.3049314002173469e-1
cos: 9.4380839390131199e-1
And the specfunc failure that we see is because gcc-6.3.0 gives us the wrong
trig values.
Switching to gcc-7.1.0 fixes this problem for me - as gcc-7.1.0 agrees with
MPFR.
mpfr_jn() returns a completely different value for the value 1048576.0 (n =
2).
I'm getting -7.02097920391228138906e-4.
I guess (hope) that it's a different function to gsl_sf_bessel_j2_e().
Cheers,
Rob
