help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: sisyphus1
Subject: Re: [Help-gsl] test release for GSL 2.4
Date: Thu, 22 Jun 2017 12:30:49 +1000

-----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






reply via email to

[Prev in Thread] Current Thread [Next in Thread]