boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-26 16:52:27 +00:00
Code Issues Projects Releases Wiki Activity
Files
f5cb11eb803aa3a2ff1d73bee18a38c505f3166e
math/test/test_jacobi_theta.hpp
Evan Miller d7141cd353 Jacobi Theta functions (#394) 
* 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.
2020-08-15 18:51:47 -04:00

749 lines
32 KiB
C++
Raw Blame History

// Copyright John Maddock 2006.
// Copyright Evan Miller 2020
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/special_functions/ellint_rf.hpp>
#include <boost/math/special_functions/jacobi_elliptic.hpp>
#include <boost/math/special_functions/jacobi_theta.hpp>
#include <boost/math/special_functions/zeta.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/quadrature/exp_sinh.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/array.hpp>
#include "functor.hpp"
#include "handle_test_result.hpp"
#include "table_type.hpp"
#ifndef SC_
#define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
#endif
template <class Real, class T>
void do_test_jacobi_theta1(const T& data, const char* type_name, const char* test_name) {
typedef Real value_type;
typedef value_type (*pg)(value_type, value_type);
std::cout << "Testing: " << test_name << std::endl;
#ifdef JACOBI_THETA1_FUNCTION_TO_TEST
pg fp2 = JACOBI_THETA1_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
pg fp2 = boost::math::jacobi_theta1<value_type, value_type>;
#else
pg fp2 = boost::math::jacobi_theta1;
#endif
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp2, 0, 1),
extract_result<Real>(2));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta1", test_name);
std::cout << std::endl;
}
template <class Real, class T>
void do_test_jacobi_theta2(const T& data, const char* type_name, const char* test_name) {
typedef Real value_type;
typedef value_type (*pg)(value_type, value_type);
std::cout << "Testing: " << test_name << std::endl;
#ifdef JACOBI_THETA2_FUNCTION_TO_TEST
pg fp2 = JACOBI_THETA2_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
pg fp2 = boost::math::jacobi_theta2<value_type, value_type>;
#else
pg fp2 = boost::math::jacobi_theta2;
#endif
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp2, 0, 1),
extract_result<Real>(2));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta2", test_name);
std::cout << std::endl;
}
template <class Real, class T>
void do_test_jacobi_theta3(const T& data, const char* type_name, const char* test_name) {
typedef Real value_type;
typedef value_type (*pg)(value_type, value_type);
std::cout << "Testing: " << test_name << std::endl;
#ifdef JACOBI_THETA3_FUNCTION_TO_TEST
pg fp2 = JACOBI_THETA3_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
pg fp2 = boost::math::jacobi_theta3<value_type, value_type>;
#else
pg fp2 = boost::math::jacobi_theta3;
#endif
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp2, 0, 1),
extract_result<Real>(2));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta3", test_name);
std::cout << std::endl;
}
template <class Real, class T>
void do_test_jacobi_theta4(const T& data, const char* type_name, const char* test_name) {
typedef Real value_type;
typedef value_type (*pg)(value_type, value_type);
std::cout << "Testing: " << test_name << std::endl;
#ifdef JACOBI_THETA4_FUNCTION_TO_TEST
pg fp2 = JACOBI_THETA4_FUNCTION_TO_TEST;
#elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
pg fp2 = boost::math::jacobi_theta4<value_type, value_type>;
#else
pg fp2 = boost::math::jacobi_theta4;
#endif
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp2, 0, 1),
extract_result<Real>(2));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta4", test_name);
std::cout << std::endl;
}
template <class Real, class T>
void do_test_jacobi_theta_tau(const T& data, const char* type_name, const char* test_name) {
typedef Real value_type;
typedef value_type (*pg)(value_type, value_type);
std::cout << "Testing: " << test_name << std::endl;
#if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
pg fp1 = boost::math::jacobi_theta1tau<value_type, value_type>;
pg fp2 = boost::math::jacobi_theta2tau<value_type, value_type>;
pg fp3 = boost::math::jacobi_theta3tau<value_type, value_type>;
pg fp4 = boost::math::jacobi_theta4tau<value_type, value_type>;
#else
pg fp1 = boost::math::jacobi_theta1tau;
pg fp2 = boost::math::jacobi_theta2tau;
pg fp3 = boost::math::jacobi_theta3tau;
pg fp4 = boost::math::jacobi_theta4tau;
#endif
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp1, 0, 1),
extract_result<Real>(2));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta1tau", test_name);
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp2, 0, 1),
extract_result<Real>(3));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta2tau", test_name);
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp3, 0, 1),
extract_result<Real>(4));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta3tau", test_name);
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp4, 0, 1),
extract_result<Real>(5));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "jacobi_theta4tau", test_name);
std::cout << std::endl;
}
template <typename T>
void test_spots(T, const char* type_name)
{
// Function values calculated on https://wolframalpha.com/
// EllipticTheta[1, z, q]
static const boost::array<boost::array<typename table_type<T>::type, 3>, 22> data1 = {{
{{ SC_(0.25), SC_(0.5), SC_(0.1540230688155610552510349122197994458164480291364308) }},
{{ SC_(0.5), SC_(0.5), SC_(0.402768575854814314826394321410682828786027207014725) }},
{{ SC_(1.0), SC_(0.5), SC_(1.330378498179274650272750199052730280058943456725878763411) }},
{{ SC_(2.0), SC_(0.5), SC_(1.632025902952598833772353216268208997235004608799433589257) }},
{{ SC_(4.0), SC_(0.5), SC_(-1.02330632494166272025903454492708687979388431668700575889) }},
{{ SC_(10.0), SC_(0.5), SC_(-0.506725689219604598643739369857898454980617737596340) }},
{{ SC_(0.25), SC_(0.0078125), SC_(0.147082536469061213804178649159394420990352754783257117514) }},
{{ SC_(0.5), SC_(0.0078125), SC_(0.285031930001354337576834900191853014429316815453397057542) }},
{{ SC_(1.0), SC_(0.0078125), SC_(0.500336519612853406200885943502694674511080381314798446615) }},
{{ SC_(2.0), SC_(0.0078125), SC_(0.540681625286624428872671041984657483294251553554422250931) }},
{{ SC_(4.0), SC_(0.0078125), SC_(-0.44997798288122032252205899314355127447793127666845934004) }},
{{ SC_(10.0), SC_(0.0078125), SC_(-0.32344103052261772036606444750254248788211904932304982822) }},
{{ SC_(1.5), SC_(0.9375), SC_(6.455616598043074010374387709748431776849714249419980311428) }},
{{ SC_(1.5), SC_(0.96875), SC_(8.494748959742732967152146642136570398472098411842769875811) }},
{{ SC_(1.5), SC_(0.984375), SC_(10.27394960641515008799474668007745853274050682503839614462) }},
{{ SC_(1.5), SC_(0.9921875), SC_(10.56322418602486655878408856907764433121023675569880599668) }},
{{ SC_(0.0), SC_(0.0078125), SC_(0.0) }},
{{ SC_(0.0), SC_(0.5), SC_(0.0) }},
{{ SC_(0.0), SC_(0.9375), SC_(0.0) }},
{{ SC_(0.0), SC_(0.96875), SC_(0.0) }},
{{ SC_(0.0), SC_(0.984375), SC_(0.0) }},
{{ SC_(0.0), SC_(0.9921875), SC_(0.0) }},
}};
// EllipticTheta[2, z, q]
static const boost::array<boost::array<typename table_type<T>::type, 3>, 22> data2 = {{
{{ SC_(0.25), SC_(0.5), SC_(1.945359087094512287818938605108992884433591043123906291186) }},
{{ SC_(0.5), SC_(0.5), SC_(1.484216087659583107499509464625356597654932790316228596683) }},
{{ SC_(1.0), SC_(0.5), SC_(0.500198138514456200618643558666164246520575297293771869190) }},
{{ SC_(2.0), SC_(0.5), SC_(-0.31816282165462356641101267196568721591143313305914629995) }},
{{ SC_(4.0), SC_(0.5), SC_(-0.73416812893190753892245332974270105112878210782122749389) }},
{{ SC_(10.0), SC_(0.5), SC_(-1.32067302962326803616213092008760707610192812421263609239) }},
{{ SC_(0.25), SC_(0.0078125), SC_(0.576145327104766930654951565363812166642791552142749996541) }},
{{ SC_(0.5), SC_(0.0078125), SC_(0.521816280475855206768414007009226207248962173223429288460) }},
{{ SC_(1.0), SC_(0.0078125), SC_(0.321229744663905222607893889592775291250157432381753805004) }},
{{ SC_(2.0), SC_(0.0078125), SC_(-0.24740754322178854446297728714315695332313383557595139996) }},
{{ SC_(4.0), SC_(0.0078125), SC_(-0.38862819739105110692686845441417046539892089911988041243) }},
{{ SC_(10.0), SC_(0.0078125), SC_(-0.49890931813624527541778224812869556711737103111567319268) }},
{{ SC_(3.0), SC_(0.9375), SC_(-5.11392816786538016035153925334241975210394670067130953352) }},
{{ SC_(3.0), SC_(0.96875), SC_(-5.29012912680048642016398403857878652305705819627897707515) }},
{{ SC_(3.0), SC_(0.984375), SC_(-3.95437444890235969862463149250591376362876865922981295660) }},
{{ SC_(3.0), SC_(0.9921875), SC_(-1.55309936234390798246955842243578030972976727248025068776) }},
{{ SC_(0.0), SC_(0.0078125), SC_(0.594639849222534631954791856184512118943710280851563329851) }},
{{ SC_(0.0), SC_(0.5), SC_(2.128931250513027558591613402575350180853805396958448940968) }},
{{ SC_(0.0), SC_(0.9375), SC_(6.976947123071698246084428957201676843908839940030780606850) }},
{{ SC_(0.0), SC_(0.96875), SC_(9.947454796978382607130245293535173220560623896730670936855) }},
{{ SC_(0.0), SC_(0.984375), SC_(14.12398706491126638681068088410435889335374416476169364474) }},
{{ SC_(0.0), SC_(0.9921875), SC_(20.01377050922054864039693105446212500742825758308336203414) }},
}};
// EllipticTheta[3, z, q]
static const boost::array<boost::array<typename table_type<T>::type, 3>, 22> data3 = {{
{{ SC_(0.25), SC_(0.5), SC_(1.945383919743612326705943032930976804537995814958244156964) }},
{{ SC_(0.5), SC_(0.5), SC_(1.484396862425166928164115914328477415075581759665236164625) }},
{{ SC_(1.0), SC_(0.5), SC_(0.505893885730484607919474452677852065978820023168006719298) }},
{{ SC_(2.0), SC_(0.5), SC_(0.331435978324530423856445870208355989399154547338364678855) }},
{{ SC_(4.0), SC_(0.5), SC_(0.736474899103717622792604836948158645655031914452730542597) }},
{{ SC_(10.0), SC_(0.5), SC_(1.320991123572679837556511698539830878932277973655257733968) }},
{{ SC_(0.25), SC_(0.0078125), SC_(1.013712231555102950279020764520600278509143561359073286704) }},
{{ SC_(0.5), SC_(0.0078125), SC_(1.008442220428654137020348371834353416731512133091516925412) }},
{{ SC_(1.0), SC_(0.0078125), SC_(0.993497700808926421501904684196885527471205379480110676883) }},
{{ SC_(2.0), SC_(0.0078125), SC_(0.989786817339946335270527719280367827578464526769825119042) }},
{{ SC_(4.0), SC_(0.0078125), SC_(0.997726554836621271192712333586918701234666824689471509074) }},
{{ SC_(10.0), SC_(0.0078125), SC_(1.006376277246758468079371354566456368238591100751823541324) }},
{{ SC_(3.0), SC_(0.9375), SC_(5.113928167865380160351539253342419752103946700671309533529) }},
{{ SC_(3.0), SC_(0.96875), SC_(5.290129126800486420163984038578786523057058196278977075158) }},
{{ SC_(3.0), SC_(0.984375), SC_(3.954374448902359698624631492505913763628768659229812956604) }},
{{ SC_(3.0), SC_(0.9921875), SC_(1.553099362343907982469558422435780309729767272480250687761) }},
{{ SC_(0.0), SC_(0.0078125), SC_(1.015625007450580597140668559497101271987479437621158995242) }},
{{ SC_(0.0), SC_(0.5), SC_(2.128936827211877158669458548544951324612516539940878092889) }},
{{ SC_(0.0), SC_(0.9375), SC_(6.976947123071698246084428957201676843908839940030780606850) }},
{{ SC_(0.0), SC_(0.96875), SC_(9.947454796978382607130245293535173220560623896730670936855) }},
{{ SC_(0.0), SC_(0.984375), SC_(14.12398706491126638681068088410435889335374416476169364474) }},
{{ SC_(0.0), SC_(0.9921875), SC_(20.01377050922054864039693105446212500742825758308336203414) }},
}};
// EllipticTheta[4, z, q]
static const boost::array<boost::array<typename table_type<T>::type, 3>, 20> data4 = {{
{{ SC_(0.25), SC_(0.5), SC_(0.189666257078605856907477593562312286776776156459895303534) }},
{{ SC_(0.5), SC_(0.5), SC_(0.411526533253405515206323680892825857445581901774756902114) }},
{{ SC_(1.0), SC_(0.5), SC_(1.330686328485433289294314954726283002076056588770122570003) }},
{{ SC_(2.0), SC_(0.5), SC_(1.632130562351990831100773069064726698266264072056233877584) }},
{{ SC_(4.0), SC_(0.5), SC_(1.024161147731330827103564503229671566751499815308607281771) }},
{{ SC_(10.0), SC_(0.5), SC_(0.512267623558970547872956225451763774220817155272291299008) }},
{{ SC_(0.25), SC_(0.0078125), SC_(0.986287776496028802869709593974763104913490530471778951141) }},
{{ SC_(0.5), SC_(0.0078125), SC_(0.991557773370274771280909909075180928934876164736866576523) }},
{{ SC_(1.0), SC_(0.0078125), SC_(1.006502289451024620679171207071841795644771897492095505061) }},
{{ SC_(2.0), SC_(0.0078125), SC_(1.010213180491934207237038406874868068814754041886231359135) }},
{{ SC_(4.0), SC_(0.0078125), SC_(1.002273430893140443692155306258136796608858362734966941965) }},
{{ SC_(10.0), SC_(0.0078125), SC_(0.993623712815089968927968772900270119031604138870447126142) }},
{{ SC_(1.5), SC_(0.9375), SC_(6.455616598043074010374387709748431776849714249419980311428) }},
{{ SC_(1.5), SC_(0.96875), SC_(8.494748959742732967152146642136570398472098411842769875811) }},
{{ SC_(1.5), SC_(0.984375), SC_(10.27394960641515008799474668007745853274050682503839614462) }},
{{ SC_(1.5), SC_(0.9921875), SC_(10.56322418602486655878408856907764433121023675569880599668) }},
{{ SC_(0.0), SC_(0.0078125), SC_(0.984375007450580596706987690502899498384498317273182227148) }},
{{ SC_(0.0), SC_(0.5), SC_(0.121124208002580502460849293181867505809858246820960597233) }},
{{ SC_(0.0), SC_(0.9375), SC_(3.4752705802238602772154431927173524635861732774491547e-16) }},
{{ SC_(0.0), SC_(0.96875), SC_(3.5224799214778391114074447929704379790977217720041291e-33) }}
// {{ SC_(0.0), SC_(0.984375), SC_(0.0) }},
// {{ SC_(0.0), SC_(0.9921875), SC_(0.0) }},
}};
do_test_jacobi_theta1<T>(data1, type_name, "Jacobi Theta 1: Wolfrom Alpha Data");
do_test_jacobi_theta2<T>(data2, type_name, "Jacobi Theta 2: Wolfram Alpha Data");
do_test_jacobi_theta3<T>(data3, type_name, "Jacobi Theta 3: Wolfram Alpha Data");
do_test_jacobi_theta4<T>(data4, type_name, "Jacobi Theta 4: Wolfram Alpha Data");
#include "jacobi_theta_data.ipp"
do_test_jacobi_theta_tau<T>(jacobi_theta_data, type_name, "Jacobi Theta: Random Data");
#include "jacobi_theta_small_tau.ipp"
do_test_jacobi_theta_tau<T>(jacobi_theta_small_tau_data, type_name, "Jacobi Theta: Random Data (Small Tau)");
}
#define _check_close(a, b, eps) \
if (abs(a) < 2 * eps * eps || abs(b) < 2 * eps * eps) { \
BOOST_CHECK_SMALL((a) - (b), eps); \
} else { \
BOOST_CHECK_CLOSE_FRACTION(a, b, eps); \
}
template <typename RealType>
inline void test_periodicity(RealType z, RealType q, RealType eps) {
using namespace boost::math;
_check_close(
jacobi_theta1(z, q),
jacobi_theta1(z + constants::two_pi<RealType>(), q),
eps);
_check_close(
jacobi_theta2(z, q),
jacobi_theta2(z + constants::two_pi<RealType>(), q),
eps);
_check_close(
jacobi_theta3(z, q),
jacobi_theta3(z + constants::pi<RealType>(), q),
eps);
_check_close(
jacobi_theta4(z, q),
jacobi_theta4(z + constants::pi<RealType>(), q),
eps);
}
// DLMF 20.2(iii) Translation of the Argument by Half-Periods
template <typename RealType>
inline void test_argument_translation(RealType z, RealType q, RealType eps) {
using namespace boost::math;
_check_close( // DLMF 20.2.11
jacobi_theta1(z, q),
-jacobi_theta2(z + constants::half_pi<RealType>(), q),
eps);
_check_close( // DLMF 20.2.12
jacobi_theta2(z, q),
jacobi_theta1(z + constants::half_pi<RealType>(), q),
eps);
_check_close( // DLMF 20.2.13
jacobi_theta3(z, q),
jacobi_theta4(z + constants::half_pi<RealType>(), q),
eps);
_check_close( // DLMF 20.2.14
jacobi_theta4(z, q),
jacobi_theta3(z + constants::half_pi<RealType>(), q),
eps);
}
// DLMF 20.7(i) Sums of Squares
template <typename RealType>
inline void test_sums_of_squares(RealType z, RealType q, RealType eps) {
using namespace boost::math;
_check_close( // DLMF 20.7.1
jacobi_theta3(RealType(0), q) * jacobi_theta3(z, q),
hypot(
jacobi_theta4(RealType(0), q) * jacobi_theta4(z, q),
jacobi_theta2(RealType(0), q) * jacobi_theta2(z, q)),
eps);
_check_close( // DLMF 20.7.2
jacobi_theta3(RealType(0), q) * jacobi_theta4(z, q),
hypot(
jacobi_theta2(RealType(0), q) * jacobi_theta1(z, q),
jacobi_theta4(RealType(0), q) * jacobi_theta3(z, q)),
eps);
_check_close( // DLMF 20.7.3
jacobi_theta2(RealType(0), q) * jacobi_theta4(z, q),
hypot(
jacobi_theta3(RealType(0), q) * jacobi_theta1(z, q),
jacobi_theta4(RealType(0), q) * jacobi_theta2(z, q)),
eps);
_check_close( // DLMF 20.7.4
jacobi_theta2(RealType(0), q) * jacobi_theta3(z, q),
hypot(
jacobi_theta4(RealType(0), q) * jacobi_theta1(z, q),
jacobi_theta3(RealType(0), q) * jacobi_theta2(z, q)),
eps);
_check_close( // DLMF 20.7.5 (no z)
jacobi_theta3(RealType(0), q) * jacobi_theta3(RealType(0), q),
hypot(
jacobi_theta2(RealType(0), q) * jacobi_theta2(RealType(0), q),
jacobi_theta4(RealType(0), q) * jacobi_theta4(RealType(0), q)),
eps);
}
// DLMF 20.7(ii) Addition Formulas
template <typename RealType>
inline void test_addition_formulas(RealType z, RealType w, RealType q, RealType eps) {
using namespace boost::math;
_check_close( // DLMF 20.7.6
jacobi_theta4(RealType(0), q) * jacobi_theta4(RealType(0), q) *
jacobi_theta1(w + z, q) * jacobi_theta1(w - z, q),
jacobi_theta3(w, q) * jacobi_theta3(w, q) * jacobi_theta2(z, q) * jacobi_theta2(z, q) -
jacobi_theta2(w, q) * jacobi_theta2(w, q) * jacobi_theta3(z, q) * jacobi_theta3(z, q),
eps);
_check_close( // DLMF 20.7.7
jacobi_theta4(RealType(0), q) * jacobi_theta4(RealType(0), q) *
jacobi_theta2(w + z, q) * jacobi_theta2(w - z, q),
jacobi_theta4(w, q) * jacobi_theta4(w, q) * jacobi_theta2(z, q) * jacobi_theta2(z, q) -
jacobi_theta1(w, q) * jacobi_theta1(w, q) * jacobi_theta3(z, q) * jacobi_theta3(z, q),
eps);
_check_close( // DLMF 20.7.8
jacobi_theta4(RealType(0), q) * jacobi_theta4(RealType(0), q) *
jacobi_theta3(w + z, q) * jacobi_theta3(w - z, q),
jacobi_theta4(w, q) * jacobi_theta4(w, q) * jacobi_theta3(z, q) * jacobi_theta3(z, q) -
jacobi_theta1(w, q) * jacobi_theta1(w, q) * jacobi_theta2(z, q) * jacobi_theta2(z, q),
eps);
_check_close( // DLMF 20.7.9
jacobi_theta4(RealType(0), q) * jacobi_theta4(RealType(0), q) *
jacobi_theta4(w + z, q) * jacobi_theta4(w - z, q),
jacobi_theta3(w, q) * jacobi_theta3(w, q) * jacobi_theta3(z, q) * jacobi_theta3(z, q) -
jacobi_theta2(w, q) * jacobi_theta2(w, q) * jacobi_theta2(z, q) * jacobi_theta2(z, q),
eps);
}
// DLMF 20.7(iii) Duplication Formula
template <typename RealType>
inline void test_duplication_formula(RealType z, RealType q, RealType eps) {
using namespace boost::math;
_check_close( // DLMF 20.7.10
jacobi_theta1(z + z, q) * jacobi_theta2(RealType(0), q) * jacobi_theta3(RealType(0), q) * jacobi_theta4(RealType(0), q),
RealType(2) * jacobi_theta1(z, q) * jacobi_theta2(z, q) * jacobi_theta3(z, q) * jacobi_theta4(z, q),
eps);
}
// DLMF 20.7(iv) Transformations of Nome
template <typename RealType>
inline void test_transformations_of_nome(RealType z, RealType q, RealType eps) {
using namespace boost::math;
_check_close( // DLMF 20.7.11
jacobi_theta1(z, q) * jacobi_theta2(z, q) * jacobi_theta4(z + z, q * q),
jacobi_theta3(z, q) * jacobi_theta4(z, q) * jacobi_theta1(z + z, q * q),
eps);
_check_close( // DLMF 20.7.12
jacobi_theta1(z, q * q) * jacobi_theta4(z, q * q) * jacobi_theta2(z, q),
jacobi_theta2(z, q * q) * jacobi_theta3(z, q * q) * jacobi_theta1(z, q),
eps);
}
// DLMF 20.7(v) Watson's Identities
template <typename RealType>
inline void test_watsons_identities(RealType z, RealType w, RealType q, RealType eps) {
using namespace boost::math;
// Rearrange DLMF equations to get q*q on each side of the equality
_check_close( // DLMF 20.7.13
jacobi_theta1(z, q) * jacobi_theta1(w, q)
+ jacobi_theta2(z + w, q * q) * jacobi_theta3(z - w, q * q),
jacobi_theta3(z + w, q * q) * jacobi_theta2(z - w, q * q),
eps);
_check_close( // DLMF 20.7.14
jacobi_theta3(z, q) * jacobi_theta3(w, q)
- jacobi_theta2(z + w, q * q) * jacobi_theta2(z - w, q * q),
jacobi_theta3(z + w, q * q) * jacobi_theta3(z - w, q * q),
eps);
_check_close( // MathWorld Eqn. 48
jacobi_theta3(z, q) - jacobi_theta2(z + z, q*q*q*q),
jacobi_theta3(z + z, q*q*q*q),
eps);
_check_close( // MathWorld Eqn. 49
jacobi_theta4(z, q) + jacobi_theta2(z + z, q*q*q*q),
jacobi_theta3(z + z, q*q*q*q),
eps);
}
template <typename RealType>
inline void test_landen_transformations(RealType z, RealType tau, RealType eps) {
// A and B below are the reciprocals of their DLMF definitions
using namespace boost::math;
// DLMF 20.7.15 (reciprocal)
RealType A = jacobi_theta4tau(RealType(0), tau + tau);
_check_close( // DLMF 20.7.16
jacobi_theta1tau(z + z, tau + tau) * A,
jacobi_theta1tau(z, tau) * jacobi_theta2tau(z, tau),
eps);
_check_close( // DLMF 20.7.17
jacobi_theta2tau(z + z, tau + tau) * A,
jacobi_theta1tau(constants::quarter_pi<RealType>() - z, tau) * jacobi_theta1tau(constants::quarter_pi<RealType>() + z, tau),
eps);
_check_close( // DLMF 20.7.18
jacobi_theta3tau(z + z, tau + tau) * A,
jacobi_theta3tau(constants::quarter_pi<RealType>() - z, tau) * jacobi_theta3tau(constants::quarter_pi<RealType>() + z, tau),
eps);
_check_close( // DLMF 20.7.19
jacobi_theta4tau(z + z, tau + tau) * A,
jacobi_theta3tau(z, tau) * jacobi_theta4tau(z, tau),
eps);
// DLMF 20.7.20 (reciprocal)
RealType B = jacobi_theta3tau(RealType(0), tau) * jacobi_theta4tau(RealType(0), tau) * jacobi_theta3tau(constants::quarter_pi<RealType>(), tau);
_check_close( // DLMF 20.7.21
jacobi_theta1tau(4*z, 4*tau) * B,
jacobi_theta1tau(z, tau)
* jacobi_theta1tau(constants::quarter_pi<RealType>() - z, tau)
* jacobi_theta1tau(constants::quarter_pi<RealType>() + z, tau)
* jacobi_theta2tau(z, tau),
eps);
_check_close( // DLMF 20.7.22
jacobi_theta2tau(4*z, 4*tau) * B,
jacobi_theta2tau(constants::quarter_pi<RealType>()/2 - z, tau)
* jacobi_theta2tau(constants::quarter_pi<RealType>()/2 + z, tau)
* jacobi_theta2tau(constants::three_quarters_pi<RealType>()/2 - z, tau)
* jacobi_theta2tau(constants::three_quarters_pi<RealType>()/2 + z, tau),
eps);
_check_close( // DLMF 20.7.23
jacobi_theta3tau(4*z, 4*tau) * B,
jacobi_theta3tau(constants::quarter_pi<RealType>()/2 - z, tau)
* jacobi_theta3tau(constants::quarter_pi<RealType>()/2 + z, tau)
* jacobi_theta3tau(constants::three_quarters_pi<RealType>()/2 - z, tau)
* jacobi_theta3tau(constants::three_quarters_pi<RealType>()/2 + z, tau),
eps);
_check_close( // DLMF 20.7.24
jacobi_theta4tau(4*z, 4*tau) * B,
jacobi_theta4tau(z, tau)
* jacobi_theta4tau(constants::quarter_pi<RealType>() - z, tau)
* jacobi_theta4tau(constants::quarter_pi<RealType>() + z, tau)
* jacobi_theta3tau(z, tau),
eps);
}
template <typename RealType>
inline void test_special_values(RealType eps) {
// https://mathworld.wolfram.com/JacobiThetaFunctions.html (Eq. 45)
using namespace boost::math;
BOOST_MATH_STD_USING
_check_close(
tgamma(RealType(0.75)) * jacobi_theta3tau(RealType(0), RealType(1)),
sqrt(constants::root_pi<RealType>()),
eps);
_check_close(
tgamma(RealType(1.25))
* constants::root_pi<RealType>()
* sqrt(sqrt(constants::root_two<RealType>()))
* jacobi_theta3tau(RealType(0), constants::root_two<RealType>()),
tgamma(RealType(1.125))
* sqrt(tgamma(RealType(0.25))),
eps);
_check_close(
tgamma(RealType(0.75))
* sqrt(constants::root_two<RealType>())
* jacobi_theta4tau(RealType(0), RealType(1)),
sqrt(constants::root_pi<RealType>()),
eps);
}
template <typename RealType>
inline void test_mellin_transforms(RealType s, RealType integration_eps, RealType result_eps) {
using namespace boost::math;
BOOST_MATH_STD_USING
boost::math::quadrature::exp_sinh<RealType> integrator;
auto f2 = [&, s](RealType t)
{
if (t*t == 0.f)
return RealType(0);
if (t > sqrt(sqrt(std::numeric_limits<RealType>::max())))
return RealType(0);
return pow(t, s-1) * jacobi_theta2tau(RealType(0), t*t);
};
auto f3 = [&, s](RealType t)
{
if (t*t == 0.f)
return RealType(0);
if (t > sqrt(sqrt(std::numeric_limits<RealType>::max())))
return RealType(0);
return pow(t, s-1) * jacobi_theta3m1tau(RealType(0), t*t);
};
auto f4 = [&, s](RealType t)
{
if (t*t == 0.f)
return RealType(0);
if (t > sqrt(sqrt(std::numeric_limits<RealType>::max())))
return RealType(0);
return -pow(t, s-1) * jacobi_theta4m1tau(RealType(0), t*t);
};
_check_close( // DLMF 20.10.1
integrator.integrate(f2, integration_eps),
(pow(RealType(2), s) - 1) * pow(constants::pi<RealType>(), -0.5*s) * tgamma(0.5*s) * zeta(s),
result_eps);
_check_close( // DLMF 20.10.2
integrator.integrate(f3, integration_eps),
pow(constants::pi<RealType>(), -0.5*s) * tgamma(0.5*s) * zeta(s),
result_eps);
_check_close( // DLMF 20.10.3
integrator.integrate(f4, integration_eps),
(1 - pow(RealType(2), 1 - s)) * pow(constants::pi<RealType>(), -0.5*s) * tgamma(0.5*s) * zeta(s),
result_eps);
}
template <typename RealType>
inline void test_laplace_transforms(RealType s, RealType integration_eps, RealType result_eps) {
using namespace boost::math;
BOOST_MATH_STD_USING
RealType beta = -0.5;
RealType l = sinh(abs(beta)) + 1.0;
boost::math::quadrature::exp_sinh<RealType> integrator;
auto f1 = [&, s, l, beta](RealType t)
{
return exp(-s*t) * jacobi_theta1tau(0.5 * beta * constants::pi<RealType>() / l,
constants::pi<RealType>() * t / l / l);
};
auto f4 = [&, s, l, beta](RealType t)
{
return exp(-s*t) * jacobi_theta4tau(0.5 * beta * constants::pi<RealType>() / l,
constants::pi<RealType>() * t / l / l);
};
_check_close( // DLMF 20.10.4 says the RHS should be negative?
integrator.integrate(f1, integration_eps),
l/sqrt(s)*sinh(beta*sqrt(s))/cosh(l*sqrt(s)),
result_eps);
_check_close( // DLMF 20.10.5
integrator.integrate(f4, integration_eps),
l/sqrt(s)*cosh(beta*sqrt(s))/sinh(l*sqrt(s)),
result_eps);
// TODO - DLMF defines two additional relations for theta2 and theta3, but
// these do not match the computed values at all.
}
template <typename RealType>
inline void test_elliptic_functions(RealType z, RealType q, RealType eps) {
using namespace boost::math;
BOOST_MATH_STD_USING
RealType t2 = jacobi_theta2(RealType(0), q);
RealType t3 = jacobi_theta3(RealType(0), q);
RealType t4 = jacobi_theta4(RealType(0), q);
RealType k = t2 * t2 / (t3 * t3);
RealType xi = z / (t3 * t3);
_check_close( // DLMF 22.2.4
jacobi_sn(k, z) *
t2 * jacobi_theta4(xi, q),
t3 * jacobi_theta1(xi, q),
eps);
_check_close( // DLMF 22.2.5
jacobi_cn(k, z) *
t2 * jacobi_theta4(xi, q),
t4 * jacobi_theta2(xi, q),
eps);
_check_close( // DLMF 22.2.6
jacobi_dn(k, z) *
t3 * jacobi_theta4(xi, q),
t4 * jacobi_theta3(xi, q),
eps);
_check_close( // DLMF 22.2.7
jacobi_sd(k, z) *
t2 * t4 * jacobi_theta3(xi, q),
t3 * t3 * jacobi_theta1(xi, q),
eps);
_check_close( // DLMF 22.2.8
jacobi_cd(k, z) *
t2 * jacobi_theta3(xi, q),
t3 * jacobi_theta2(xi, q),
eps);
_check_close( // DLMF 22.2.9
jacobi_cd(k, z) *
t2 * jacobi_theta3(xi, q),
t3 * jacobi_theta2(xi, q),
eps);
_check_close( // DLMF 22.2.9
jacobi_sc(k, z) *
t4 * jacobi_theta2(xi, q),
t3 * jacobi_theta1(xi, q),
eps);
}
template <typename RealType>
inline void test_elliptic_integrals(RealType q, RealType eps) {
using namespace boost::math;
BOOST_MATH_STD_USING
RealType t2 = jacobi_theta2(RealType(0), q);
RealType t3 = jacobi_theta3(RealType(0), q);
RealType t4 = jacobi_theta4(RealType(0), q);
RealType k = t2*t2 / (t3*t3);
RealType k1= t4*t4 / (t3*t3);
if (t3*t3*t3*t3 != 0.f && t4*t4*t4*t4 != 0.f) {
_check_close( // DLMF 20.9.4
ellint_rf(RealType(0), t3*t3*t3*t3, t4*t4*t4*t4),
constants::half_pi<RealType>(),
eps);
}
if (k*k != 0.f && k1*k1 != 0.f) {
_check_close( // DLMF 20.9.5
ellint_rf(RealType(0), k1*k1, RealType(1))
* log(q) / constants::pi<RealType>(),
-ellint_rf(RealType(0), k*k, RealType(1)),
eps);
}
}
Reference in New Issue View Git Blame Copy Permalink