bugGNU 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

Submitted by:  Guido De Rosa <gderosa>
Submitted on:  Mon Jan 26 19:52:59 2009  
Category: Accuracy problemSeverity: 3 - Normal
Operating System: Уеб-?траницаStatus: Fixed
Assigned to: NoneOpen/Closed: Closed
Release: 1.12

(Jump to the original submission Jump to the original submission)

Tue Feb 10 17:18:16 2009, comment #9:

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.

-Deleted Account- <bjg>
Mon Feb 2 02:24:43 2009, 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 :

  • two integers are casted to double, in order to call pow() and assign the result to another double (line 182)
  • an integer is divided by a double, then compared to an integer constant, and the result is assigned to an integer (line 183)
  • the product of an integer and a double is assigned to an integer (line 184)

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.

(file #17380)

Guido De Rosa <gderosa>
Sun Feb 1 17:32:33 2009, 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 Feb 1 06:48:28 2009, 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

(file #17370, file #17371)

Guido De Rosa <gderosa>
Sat Jan 31 21:04:31 2009, 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


-Deleted Account- <bjg>
Sat Jan 31 16:35:10 2009, 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)
SSE registers work better than x87 FPU and give better performance. This is the best choice, when possible (and, for the records, is the default gcc behavior on x86_64).

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
+++ vegas.c.new 2009-01-31 15:31:43.703125000 +0100
@@ -56,6 +56,7 @@
/* standard headers */
#include <math.h>
#include <stdio.h>
+#include <fenv.h>

/* gsl headers */
#include <gsl/gsl_math.h>
@@ -502,6 +503,9 @@
state->bins = state->bins_max;
state->ostream = stdout;

+ /* 53-bit mantissa */
+ feupdateenv(FE_PC53_ENV);

where $(MINGW_prefix)/include/fenv.h reads:

/* The floating point environment set by MSVCRT _fpreset (53-bit mantissa) */
#define FE_PC53_ENV ((const fenv_t *)-2)

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.


Guido De Rosa <gderosa>
Thu Jan 29 18:48:24 2009, comment #3:

> I'd suggest that you rerun the test with tracing for the vegas > method by adding s->verbose=1 before the call to vegas

Here's the result (test.log attached).

(file #17351)

Guido De Rosa <gderosa>
Thu Jan 29 13:59:26 2009, 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 Jan 27 05:55:36 2009, 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 Jan 26 19:52:59 2009, original submission:

This small patch fixes a long-standing problem when building on Windows with MinGW port of gcc:


Thanks to Danny Smith for suggestions (my previous work-around was much more convoluted):


Apply this patch and run ./autogen.sh to rebuild configure script and other files.

Guido De Rosa <gderosa>


Attached Files
file #17380:  vegas-types.diff added by gderosa (517B - application/octet-stream)
file #17370:  gdb.txt added by gderosa (6kB - text/plain)
file #17371:  gdb-no-opim.txt added by gderosa (5kB - text/plain)
file #17351:  test.log added by gderosa (339kB - application/octet-stream - monte test output with s->verbose=1)
file #17333:  gsl-mingw.diff added by gderosa (440B - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by gderosa (Submitted the item)

    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 8 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue Feb 10 17:18:16 2009bjgStatusNeed Info=>Fixed
    Mon Feb 2 02:24:43 2009gderosaAttached File-=>Added vegas-types.diff, #17380
    Sun Feb 1 06:48:28 2009gderosaAttached File-=>Added gdb.txt, #17370
      Attached File-=>Added gdb-no-opim.txt, #17371
    Thu Jan 29 18:48:24 2009gderosaAttached File-=>Added test.log, #17351
    Thu Jan 29 13:59:26 2009bjgStatusNone=>Need Info
    Mon Jan 26 19:52:59 2009gderosaAttached File-=>Added gsl-mingw.diff, #17333

    Back to the top

    Powered by Savane 3.1-cleanup1