Re: [Help-gsl] test release for GSL 2.4

[sisyphus1]

-----Original Message----- 
From: Max Gaspa
Sent: Wednesday, June 21, 2017 11:26 PM
To: address@hidden
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