mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
d7141cd353
* Jacobi Theta functions Implementations, tests, and ULP plotting programs are provided for the four Jacobi Theta functions per #373. Twenty-four public C++ functions are provided in all, covering various precision-preserving scenarios. Documentation for collaborators is provided in the code comments. Proper documentation for end users will be provided when the implementation and APIs are finalized. Some tests are failing; this implementation is meant to start a conversation. The core dilemma faced by the author was that large values of |q| resulted in slow convergence, and sometimes wildly inaccurate results. Following the implementation note in DLMF 20.14, I added code to switch over to the imaginary versions of the theta functions when |q| > 0.85. This restored accuracy such that all of the identity tests passed for a loose-enough epsilon, but then lost precision to the point that the Wolfram Alpha spot checks failed. It is the author's hope that someone with floating-point experience can tame the exponential dragons and squeeze the ULPs back down to a reasonable range when |q| is large. When #392 is merged I will add more thorough value tests, although I fully expect them to fail until the underlying precision issues are resolved. As a final note, the precision issues do not affect the z=0 case - the ULP plots indicate these return values within 2 ULP across all valid |q|. So that's a start. * [CI SKIP] Jacobi theta: Add special-value tests and more * Add tests covering z=0 special values from MathWorld * Add missing real_concept header * Replace M_PI and friends with constants::pi etc * Use BOOST_MATH_STD_USING in more places * Jacobi theta: Test two more of Watson's identities [CI SKIP] See https://mathworld.wolfram.com/JacobiThetaFunctions.html (Equations 48 and 49) * Improve precision of Jacobi theta functions [CI SKIP] Rewrite the private imaginary versions to use double-sided summations following DLMF 20.13.4 and 20.13.5. This cuts down the worst of the precision issues by a factor of 10, and gets more of the tests to pass. I am confident enough in the code path to eliminate the compile-time __JACOBI_THETA_USE_IMAGINARY flag. In fact the imaginary-z code paths are now enabled for all |q| > 0.04, i.e. most legal values of q. More extensive tests will be needed to illuminate any remaining precision issues. * Jacobi theta: Make changes suggested in #394 [CI SKIP] * Add LICENSE notice to main file * Document convergence criteria * Eliminate eps*eps = 0 logic. This causes some disagreement with the zero returned by Wolfram Alpha for z=0, q > 0.99 in the fourth function. Mathematically, the fourth function is never exactly zero, so I don't trust Wolfram here. * Per code-review comments, remove multiplications by floating-point 2. * Tweak the plotting programs to display their titles, and to uniformly use `float` as their CoarseType and `long double` as their `PreciseType`. * Add quadrature tests to Jacobi theta functions [CI SKIP] The quadrature tests revealed a problem in the m1 functions: they too should switch to the _IMAGINARY logic for q > exp(-pi), or will suffer from slow convergence. Fix them. Also tighten tolerances for many tests from sqrt(eps) to 100 * eps. * Test Jacobi thetas against elliptic functions and elliptic integrals [CI SKIP] See: * https://dlmf.nist.gov/22.2 * https://dlmf.nist.gov/20.9#i * Test Jacobi Thetas against their Laplace transforms [CI SKIP] See: * https://dlmf.nist.gov/20.10#ii I did find some disagreement, and dropped the negative sign from the theta1 equation. DLMF's theta2 and theta3 Laplace transform equations do not agree at all with the computed values - will need to investigate. In the meantime, the two implemented equations agree to 4 EPS so I am keeping them. * Add a note on using log1p with Jacobi theta functions [CI SKIP] See discussion: * https://github.com/boostorg/math/pull/394#issuecomment-655871762 * Add random data tests to Jacobi Theta functions [CI SKIP] Add a test data generator program for the Jacobi theta functions. This program will produce data for the tau parameterization, so that precision isn't lost during the log-transformation. This distinguishes it from the Wolfram Alpha data, which is parameterized by q. A few of these new random-data tests are failing, but not by obscene margins (< 100 EPS). These failures will be addressed when the test tolerances are finalized. * Add small-tau tests and simplify Jacobi Theta code [CI SKIP] Add tests for small tau (i.e. large q). The tests are failing with mean ~ 200 EPS and max ~ 800 EPS. These look like worst-case input, and should be the focus of future accuracy improvements. This commit also simplifies the _IMAGINARY code by abstracting all of the loops into a single svelte function. * Add user documentation for Jacobi Theta functions [CI SKIP] * Add function graphs to Jacobi Theta docs [CI SKIP] * Define Jacobi Theta test tolerances [CI SKIP] * Add implementation note on Jacobi theta functions [CI SKIP] * Consolidate Jacobi Theta ULPs plotting programs [CI SKIP] * Fix q domain checking of jacobi_theta4 [CI SKIP] * Add ULPs plots to Jacobi Theta docs [CI SKIP] Also add the built HTML files for easy evaluation. A full rebuild is needed for the new docs to appear in the indexes. * Add missing Jacobi Theta ULPs plots [CI SKIP] * Add LaTeX source for Jacobi Theta equations [CI SKIP] * Remove unused Jacobi Theta PNG equations [CI SKIP] * Add Jacobi Theta performance script [CI SKIP] Provided by @NAThompson. * Remove vestigial eps*eps check from jacobi_theta3 [CI SKIP] * Update Jacobi Theta docs per code review comments [CI SKIP] * Enable arg promotion for Jacobi Theta functions [CI SKIP] Add Jacobi theta functions to the instantiation tests and fix up everything needed to make them pass. This changes the function signatures to use promote_args. * Fix Jacobi Theta plotting script [CI SKIP] This script broke when the promote_args API was added. * Change Jacobi Theta convergence criterion [CI SKIP] Compare the non-oscillating part of the delta to the previous one. This avoids some headaches comparing the delta to the partial sum, because the partial sum can be a small number due to the oscillating component alternating signs. Because successive terms involve either q^n^2 or exp(-(pi*n)^2), convergence should still happen pretty quickly. Graphs have been updated and tests still passs with no noticeable difference.
120 lines
4.4 KiB
C++
120 lines
4.4 KiB
C++
|
|
#include <pch_light.hpp>
|
|
#include <boost/math/concepts/real_concept.hpp>
|
|
#include "test_jacobi_theta.hpp"
|
|
|
|
// Test file for the Jacobi Theta functions, a.k.a the four horsemen of the
|
|
// Jacobi elliptic integrals. At the moment only Wolfrma Alpha spot checks are
|
|
// used. We should generate extra-precise numbers with NTL::RR or some such.
|
|
|
|
void expected_results()
|
|
{
|
|
//
|
|
// Define the max and mean errors expected for
|
|
// various compilers and platforms.
|
|
//
|
|
//
|
|
add_expected_result(
|
|
".*", // compiler
|
|
".*", // stdlib
|
|
".*", // platform
|
|
".*", // test type(s)
|
|
".*Small Tau.*", // test data group
|
|
".*", 1000, 200); // test function
|
|
|
|
add_expected_result(
|
|
".*", // compiler
|
|
".*", // stdlib
|
|
".*", // platform
|
|
".*", // test type(s)
|
|
".*Wolfram Alpha.*", // test data group
|
|
".*", 60, 15); // test function
|
|
|
|
// Catch all cases come last:
|
|
//
|
|
add_expected_result(
|
|
".*", // compiler
|
|
".*", // stdlib
|
|
".*", // platform
|
|
".*", // test type(s)
|
|
".*", // test data group
|
|
".*", 20, 5); // test function
|
|
//
|
|
// Finish off by printing out the compiler/stdlib/platform names,
|
|
// we do this to make it easier to mark up expected error rates.
|
|
//
|
|
std::cout << "Tests run with " << BOOST_COMPILER << ", "
|
|
<< BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
|
|
}
|
|
|
|
BOOST_AUTO_TEST_CASE( test_main )
|
|
{
|
|
expected_results();
|
|
BOOST_MATH_CONTROL_FP;
|
|
BOOST_MATH_STD_USING
|
|
|
|
using namespace boost::math;
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta1(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta1(0.0, 1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta2(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta2(0.0, 1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta3(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta3(0.0, 1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta4(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta4(0.0, 1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta1tau(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta1tau(0.0, -1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta2tau(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta2tau(0.0, -1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta3tau(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta3tau(0.0, -1.0), std::domain_error);
|
|
|
|
BOOST_CHECK_THROW(jacobi_theta4tau(0.0, 0.0), std::domain_error);
|
|
BOOST_CHECK_THROW(jacobi_theta4tau(0.0, -1.0), std::domain_error);
|
|
|
|
double eps = std::numeric_limits<double>::epsilon();
|
|
for (double q=0.0078125; q<1.0; q += 0.0078125) { // = 1/128
|
|
for (double z=-8.0; z<=8.0; z += 0.125) {
|
|
test_periodicity(z, q, 100 * eps);
|
|
test_argument_translation(z, q, 100 * eps);
|
|
test_sums_of_squares(z, q, 100 * eps);
|
|
// The addition formula is complicated, cut it some extra slack
|
|
test_addition_formulas(z, constants::ln_two<double>(), q, sqrt(sqrt(eps)));
|
|
test_duplication_formula(z, q, 100 * eps);
|
|
test_transformations_of_nome(z, q, 100 * eps);
|
|
test_watsons_identities(z, 0.5, q, 101 * eps);
|
|
test_landen_transformations(z, -log(q)/constants::pi<double>(), sqrt(eps));
|
|
test_elliptic_functions(z, q, 5 * sqrt(eps));
|
|
}
|
|
test_elliptic_integrals(q, 10 * eps);
|
|
}
|
|
|
|
test_special_values(eps);
|
|
|
|
for (double s=0.125; s<3.0; s+=0.125) {
|
|
test_mellin_transforms(2.0 + s, eps, 3 * eps);
|
|
test_laplace_transforms(s, eps, 4 * eps);
|
|
}
|
|
|
|
test_spots(0.0F, "float");
|
|
test_spots(0.0, "double");
|
|
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
|
|
test_spots(0.0L, "long double");
|
|
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
|
|
test_spots(concepts::real_concept(0), "real_concept");
|
|
#endif
|
|
#else
|
|
std::cout << "<note>The long double tests have been disabled on this platform "
|
|
"either because the long double overloads of the usual math functions are "
|
|
"not available at all, or because they are too inaccurate for these tests "
|
|
"to pass.</note>" << std::endl;
|
|
#endif
|
|
}
|