Bug 66689 - comp_ellint_3 and ellint_3 return garbage values
Summary: comp_ellint_3 and ellint_3 return garbage values
Status: RESOLVED FIXED
Alias: None
Product: gcc
Classification: Unclassified
Component: libstdc++ (show other bugs)
Version: 5.1.0
: P3 normal
Target Milestone: 8.0
Assignee: Ed Smith-Rowland
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2015-06-27 08:49 UTC by John Maddock
Modified: 2018-04-23 09:01 UTC (History)
3 users (show)

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed: 2015-06-27 00:00:00


Attachments
Patch with regenerated testcases for all. tr1 and std. (236.02 KB, application/x-bzip2)
2017-11-17 18:14 UTC, Ed Smith-Rowland
Details

Note You need to log in before you can comment on or make changes to this bug.
Description John Maddock 2015-06-27 08:49:42 UTC
The std::tr1 functions comp_ellint_3 and ellint_3 appear to return garbage values whenever the nu parameter is non zero.

For example:

comp_ellint_3(0.75, 0) returns 1.91099 which is correct, but
comp_ellint_3(0.75, 0.5) returns 1.52852 when it should be 2.80011 (see for example http://www.wolframalpha.com/input/?i=EllipticPi[0.5%2C+0.75^2] but note different parameterization on the wolfram site)

Similar errors occur with ellint_3 rendering these functions essentially useless as far as I can see?
Comment 1 John Maddock 2015-06-28 12:32:38 UTC
I think I've found the problem - make the nu parameter negative and you get the correct answers.

These functions appear to use the same definition for ellint_3 as GSL uses (same code?), but TR1 uses a different one which matches A&S: note the change of sign in the nu parameter in the integral:

GSL has Pi(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))

TR1 has Pi(\phi,k,n) = \int_0^\phi dt 1/((1 - n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))
Comment 2 Jonathan Wakely 2015-06-28 13:33:46 UTC
Thanks for the analysis.

I'd like to re-use the TR1 code for an implementation of ISO/IEC 29124:2010 (whether or not it also goes into C++17) so we should fix this first.
Comment 3 Ed Smith-Rowland 2017-03-29 06:02:00 UTC
Well, for whatever reason, TR29124 and C++17 chose the - in front of \nu.

As long as we document and warn I think this isn't a defect.

I'll dig through the docs and get back.
Comment 4 Ed Smith-Rowland 2017-11-17 06:27:19 UTC
My last comment is nuts.  I was thrown by the fact that GSL, against which I've been testing, and the Carlson papers that form the basis if the implementation use the +nu convention.  I must change this.  This is a huge bug.
Comment 5 Ed Smith-Rowland 2017-11-17 06:29:18 UTC
in other news I've switched to boost to test this.
Comment 6 Ed Smith-Rowland 2017-11-17 18:14:23 UTC
Created attachment 42635 [details]
Patch with regenerated testcases for all. tr1 and std.

2017-11-17  Edward Smith-Rowland  <3dw4rd@verizon.net>

	PR libstdc++/pr66689 - comp_ellint_3 and ellint_3 return garbage values
	* include/tr1/ell_integral.tcc: Correct the nu sign convention
	in ellint_3 and comp_ellint_3.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	06_comp_ellint_3/check_value.cc: Regen with correct values.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	14_ellint_3/check_value.cc: Ditto.
	* testsuite/special_functions/06_comp_ellint_3/check_value.cc: Ditto.
	* testsuite/special_functions/13_ellint_3/check_value.cc: Ditto.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	01_assoc_laguerre/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	02_assoc_legendre/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	03_beta/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	04_comp_ellint_1/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	05_comp_ellint_2/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	07_conf_hyperg/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	08_cyl_bessel_i/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	09_cyl_bessel_j/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	10_cyl_bessel_k/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	11_cyl_neumann/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	12_ellint_1/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	13_ellint_2/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	15_expint/check_value_neg.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	16_hermite/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	17_hyperg/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	18_laguerre/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	19_legendre/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	20_riemann_zeta/check_value_neg.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	21_sph_bessel/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	22_sph_legendre/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	23_sph_neumann/check_value.cc: Regen.
	* testsuite/ext/special_functions/conf_hyperg/check_value.cc: Regen.
	* testsuite/ext/special_functions/hyperg/check_value.cc: Regen.
	* testsuite/special_functions/01_assoc_laguerre/check_value.cc: Regen.
	* testsuite/special_functions/02_assoc_legendre/check_value.cc: Regen.
	* testsuite/special_functions/03_beta/check_value.cc: Regen.
	* testsuite/special_functions/04_comp_ellint_1/check_value.cc: Regen.
	* testsuite/special_functions/05_comp_ellint_2/check_value.cc: Regen.
	* testsuite/special_functions/07_cyl_bessel_i/check_value.cc: Regen.
	* testsuite/special_functions/08_cyl_bessel_j/check_value.cc: Regen.
	* testsuite/special_functions/09_cyl_bessel_k/check_value.cc: Regen.
 	* testsuite/special_functions/10_cyl_neumann/check_value.cc: Regen.
	* testsuite/special_functions/11_ellint_1/check_value.cc: Regen.
	* testsuite/special_functions/12_ellint_2/check_value.cc: Regen.
	* testsuite/special_functions/14_expint/check_value.cc: Regen.
	* testsuite/special_functions/15_hermite/check_value.cc: Regen.
	* testsuite/special_functions/16_laguerre/check_value.cc: Regen.
	* testsuite/special_functions/17_legendre/check_value.cc: Regen.
	* testsuite/special_functions/18_riemann_zeta/check_value.cc: Regen.
	* testsuite/special_functions/19_sph_bessel/check_value.cc: Regen.
	* testsuite/special_functions/20_sph_legendre/check_value.cc: Regen.
	* testsuite/special_functions/21_sph_neumann/check_value.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	06_comp_ellint_3/pr66689.cc: New.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	14_ellint_3/pr66689.cc: New.
	* testsuite/special_functions/06_comp_ellint_3/pr66689.cc: New.
	* testsuite/special_functions/13_ellint_3/pr66689.cc: New.
Comment 7 Ed Smith-Rowland 2018-04-20 19:35:41 UTC
We fixed this back in 1017-11-18.