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

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: 

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


181:  {
182:    double tot_boxes = pow ((double) boxes, (double) dim);
183:    state->calls_per_box = GSL_MAX (calls / tot_boxes, 2);
184:    calls = state->calls_per_box * tot_boxes;
185:  }


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


int gsl_int_pow_uint(int x, unsigned int n);
unsigned int gsl_uint_pow_uint(unsigned int x, unsigned int n);
/* with long int variants, etc. */


So, the above lines of vegas.c would look like:


      {
        unsigned int tot_boxes = gsl_uint_pow_uint (boxes, dim);
        state->calls_per_box = GSL_MAX (calls / tot_boxes, 2);
        calls = state->calls_per_box * tot_boxes;
      }


which is clean code, with no type messing at all.



(file #17380)

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



(file #17370, file #17371)

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
box_per_bin=max(166/50,1)=3
bins=min(166/3,50)=50
boxes=50*3=150
tot_boxes=pow(150,1)=150
calls_per_box=max(333/150,2)=2
calls=2*150=300

-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)
http://gcc.gnu.org/ml/gcc/2003-08/msg01195.html
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323
http://www.arl.wustl.edu/~lockwood/class/cs306/books/artofasm/Chapter_14/CH14-3.html

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);
+
   return GSL_SUCCESS;
 }


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.

Comments?

Guido De Rosa <gderosa>
Thu 29 Jan 2009 06:48:24 PM UTC, 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 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:

https://sourceforge.net/tracker2/?func=detail&aid=2535618&group_id=2435&atid=102435#artifact_comment_3418073

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>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attached Files
file #17380:  vegas-types.diff added by gderosa (517B - application/octet-stream)
file #17370:  gdb.txt added by gderosa (6KiB - text/plain)
file #17371:  gdb-no-opim.txt added by gderosa (5KiB - text/plain)
file #17351:  test.log added by gderosa (339KiB - 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
  • -email is unavailable- added by gderosa (Submitted the item)
  •  

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2009-02-10 bjg StatusNeed Info Fixed
        Open/ClosedOpen 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 StatusNone Need Info
    2009-01-26 gderosa Attached File- Added gsl-mingw.diff, #17333

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code