GNU Scientific Library - Bugs: bug #25413, [patch] monte test failure under...
You are not allowed to post comments on this tracker with your current authentication level.
bug #25413: [patch] monte test failure under mingw
Submitter: | Guido De Rosa <gderosa> | ||
Submitted: | Mon 26 Jan 2009 07:52:59 PM UTC | ||
Category: | Accuracy problem | Severity: | 3 - Normal |
Operating System: | Уеб-?траница | Status: | Fixed |
Assigned to: | None | Open/Closed: | Closed |
Release: | 1.12 |
( Jump to the original submission )
Tue 10 Feb 2009 05:18:16 PM UTC, comment #9: |
-Deleted Account- <bjg> |
Mon 02 Feb 2009 02:24:43 AM UTC, comment #8: Another patch (vegas-types.diff).
As far as I tested it fixes the problem, and it's not a workaround.
Looking at the original monte/vegas.c :
I agree that FPU issues may lead only to tiny, platform-dependent inaccuracies but... when floating-point variables are mixed/implicitly converted to integers, errors get sometimes amplified and sometimes hidden.
A trivial example: convert 1.499999999999 and 1.500000000001 to the nearest integer, and what you get? 1 and 2, which is a huge discrepancy!
Even better than my patch, a kind of gsl_int_pow_uint() should be implemented to calculate an integer elevated to a an unsigned int power, returning another integer:
So, the above lines of vegas.c would look like:
which is clean code, with no type messing at all.
|
Guido De Rosa <gderosa> |
Sun 01 Feb 2009 05:32:33 PM UTC, comment #7: Thank you for the logs. I do think that is a compiler or system library problem.
Whether they are in double or extended registers, the results of the intermediate operations would not be different enough to change the final result assuming they are computed correctly.
If the pow function pow(x,1) doesn't return an exact result of x it could cause the final incorrect rounding.
|
-Deleted Account- <bjg> |
Sun 01 Feb 2009 06:48:28 AM UTC, comment #6: Two gdb sessions uploaded:
gdb.txt on optimized code (gcc -O2) which make the test fail
and, for comparison
gdb-no-optim.txt on non-optimized code (gcc -O0) which correctly sets 'calls' to 300
|
Guido De Rosa <gderosa> |
Sat 31 Jan 2009 09:04:31 PM UTC, comment #5: The first lines of your test.log show
The limits of integation are:
xl[0]=0.000000 xu[0]=1.000000
num_dim=1, calls=299, it_num=1, max_it_num=5 veb=1, alph=1.50,
The correct value of calls should be 300.
There are only a few steps in the calculation of calls, if you step through them you should find where it is wrong. Starting from line 166 of vegas.c I get
boxes=166
|
-Deleted Account- <bjg> |
Sat 31 Jan 2009 04:35:10 PM UTC, comment #4: Excess precision problems are not not MinGW specific, they are common to the x86 architecture. This sort of "bug" is regularly reported to the gcc tracker, and regularly closed or suspended in a "won't fix" state.
As far as I understand, relying on extended precision on x86 is not safe; serious practical consequences may be rare, yet unpredictable.
http://arxiv.org/PS_cache/cs/pdf/0701/0701192v5.pdf (sec 3.1)
There are several (tested) way to fix/workaround this in GSL without the performance loss of -ffloat-store.
1) -msse[{2|3}] -mfpmath=sse (on Pentium3, Athlon XP and newer)
When SSE is not available:
2) force the FPU to IEEE standard precision
2.1) for example:
--- vegas.c 2008-11-29 17:42:43.000000000 +0100
/* gsl headers */
+ /* 53-bit mantissa */
where $(MINGW_prefix)/include/fenv.h reads:
/* The floating point environment set by MSVCRT _fpreset (53-bit mantissa) */
2.2) equivalently, link libglsmonte to $(MINGW_prefix)/lib/CRT_fp8.o
2.3) better, write ieee-utils/fp-win32.c to implement gsl_ieee_set_mode(), then set GSL_IEEE_MODE as needed.
Naturally, there should be a way to integrate all (or part of) this stuff into the autoconf infrastructure.
Comments?
|
Guido De Rosa <gderosa> |
Thu 29 Jan 2009 06:48:24 PM UTC, comment #3:
|
Guido De Rosa <gderosa> |
Thu 29 Jan 2009 01:59:26 PM UTC, comment #2: Thanks for the link. I don't believe the test failure should be ignored or worked around with -ffloat-store as the relative size of the error is huge O(1e-3).
This suggests that there is a more serious underlying problem with extended precision registers that is specific to Mingw as it doesn't with the same versions of Gcc on GNU/Linux or other operating systems.
I'd suggest that you rerun the test with tracing for the vegas method by adding s->verbose=1 before the call to vegas and see where the error is occurring.
#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { s->verbose=1; gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; } |
-Deleted Account- <bjg> |
Tue 27 Jan 2009 05:55:36 AM UTC, comment #1: The patch is just a workaround "that works" -- but for something "deeper", here are some suggestions recently received from MinGW developers:
|
Guido De Rosa <gderosa> |
Mon 26 Jan 2009 07:52:59 PM UTC, original submission:
This small patch fixes a long-standing problem when building on Windows with MinGW port of gcc:
http://lists.gnu.org/archive/html/bug-gsl/2006-09/msg00004.html
Thanks to Danny Smith for suggestions (my previous work-around was much more convoluted):
https://sourceforge.net/tracker2/?func=detail&aid=2535618&group_id=2435&atid=102435
Apply this patch and run ./autogen.sh to rebuild configure script and other files.
|
Guido De Rosa <gderosa> |
Depends on the following items: None found
Items that depend on this one: None found
Follow 8 latest changes.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2009-02-10 | bjg | Status | Need Info | ![]() |
Fixed |
Open/Closed | Open | ![]() |
Closed | ||
2009-02-02 | gderosa | Attached File | - | ![]() |
Added vegas-types.diff, #17380 |
2009-02-01 | gderosa | Attached File | - | ![]() |
Added gdb.txt, #17370 |
Attached File | - | ![]() |
Added gdb-no-opim.txt, #17371 | ||
2009-01-29 | gderosa | Attached File | - | ![]() |
Added test.log, #17351 |
2009-01-29 | bjg | Status | None | ![]() |
Need Info |
2009-01-26 | gderosa | Attached File | - | ![]() |
Added gsl-mingw.diff, #17333 |
There is some discussion in the mingw math.h header related to this -- there is a bug in the MSCVRT pow function where it returns incorrect results for some integer arguments, as in this case. They suggest some purely Mingw based workarounds there. However since the gsl_pow_int function fixes the problem and we have it available I will switch to that, thanks for the suggestion.