[Bug libstdc++/51133] New: Incorrect implementation of std::tr1::hermite()
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51133 Bug #: 51133 Summary: Incorrect implementation of std::tr1::hermite() Classification: Unclassified Product: gcc Version: 4.6.1 Status: UNCONFIRMED Severity: normal Priority: P3 Component: libstdc++ AssignedTo: unassig...@gcc.gnu.org ReportedBy: dick...@gmail.com Created attachment 25820 --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25820 Patch fixing error in poly_hermite.tcc recursion relation The hermite() function in std::tr1, defined in libstdc++-v3/include/tr1/poly_hermite.tcc uses a recursion relation which has a sign error. The following line: __H_n = 2 * (__x * __H_nm1 + (__i - 1) * __H_nm2); should be changed to: __H_n = 2 * (__x * __H_nm1 - (__i - 1) * __H_nm2); I have attached a patch (generated inside poly_hermite.tcc's folder) which makes this change. This patch has been tested with gcc 4.6.1 by checking via numerical integration that the resultant polynomial has the correct orthogonality. As implemented, the result of this numerical integration was extremely accurate to at least 14 decimal places at up to n=84 for double precision.
[Bug libstdc++/51133] Incorrect implementation of std::tr1::hermite()
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51133 --- Comment #3 from Jason Dick 2011-11-15 05:48:37 UTC --- Created attachment 25823 --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25823 Test case for hermite(), requires GSL This test verifies the Hermite polynomials using numerical integration (using the GSL software package for the integration). An error is returned on the first integral which does not return the expected result.
[Bug libstdc++/51133] Incorrect implementation of std::tr1::hermite()
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51133 --- Comment #4 from Jason Dick 2011-11-15 05:59:27 UTC --- I have attached a tiny program which I used to test these polynomials. I believe the test is entirely robust provided the first two Hermite polynomials are correct. On my machine, this test fails at n=85 because the normalization factor inside the square root becomes too large. Computing the normalization factor in logarithmic space instead allows us to push up to n=87, but then GSL can no longer integrate the function (likely to to a limit of floating point accuracy).
[Bug libstdc++/51133] Incorrect implementation of std::tr1::hermite()
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51133 Jason Dick changed: What|Removed |Added Attachment #25823|0 |1 is obsolete|| --- Comment #10 from Jason Dick 2011-11-15 10:57:49 UTC --- Created attachment 25825 --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=25825 New test case for hermite() using only standard library Okay, I've coded up a really quick and dirty recursive integration routine and used that instead of the GSL routine. This routine is much slower and less accurate, but it works more or less. The only difficulty is that this particular integration routine isn't able to guarantee the total error very well. But as long as you use a tolerance on the integration that is a few orders of magnitude better than you check, the result works well enough. Here I've used a tolerance of 1e-10, and checked to ensure that the results agree with the expected value at 1e-7. With these parameters, and an integration range of -16.5 to 16.6, I get the expected result up to n=104 (after which an inf is encountered in the integration). As I understand it, you need to be much more clever to produce an integration routine that is able to guarantee total error even on very complex functions. This routine is also susceptible to the possibility of choosing integration values that just happen to be zero, even though the rest of the integration interval is far from zero (this is why my integration interval is not even). If you would like, I would happily support the use of this very simple integration function. If this is not good enough as it is, I wouldn't mind putting a little bit of work into implementing a more complicated routine. I am a bit unfamiliar with the practices and conventions at use within GNU development, though. So please let me know if anything in particular would be required on my part to do this. P.S. I've also noticed that I (foolishly) used both the names hermite_test and test_hermite, so I've fixed that to only use the name hermite_test.
[Bug libstdc++/51133] Incorrect implementation of std::tr1::hermite()
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=51133 --- Comment #12 from Jason Dick 2011-11-15 11:08:16 UTC --- (In reply to comment #11) > Thanks, but please, instead of doing this here, for this specific testcase, > please collaborate with Ed at integrating this code in the testsuite. Sounds good!
