Skip to content

Jacobi Theta functions #394

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
NAThompson merged 33 commits into from
Aug 15, 2020
Merged

Jacobi Theta functions #394

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
33 commits
Select commit Hold shift + click to select a range
898c3ad
Jacobi Theta functions
evanmiller Jul 7, 2020
c78dcd7
[CI SKIP] Jacobi theta: Add special-value tests and more
evanmiller Jul 7, 2020
10fc561
Jacobi theta: Test two more of Watson's identities [CI SKIP]
evanmiller Jul 7, 2020
4b6701d
Improve precision of Jacobi theta functions [CI SKIP]
evanmiller Jul 8, 2020
9c68c11
Jacobi theta: Make changes suggested in #394 [CI SKIP]
evanmiller Jul 8, 2020
aa382f5
Add quadrature tests to Jacobi theta functions [CI SKIP]
evanmiller Jul 8, 2020
45e8ab9
Test Jacobi thetas against elliptic functions and elliptic integrals …
evanmiller Jul 9, 2020
0d3bc43
Test Jacobi Thetas against their Laplace transforms [CI SKIP]
evanmiller Jul 9, 2020
c293d48
Add a note on using log1p with Jacobi theta functions [CI SKIP]
evanmiller Jul 9, 2020
6b8dd5a
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 9, 2020
08391cd
Add random data tests to Jacobi Theta functions [CI SKIP]
evanmiller Jul 9, 2020
d73472e
Add small-tau tests and simplify Jacobi Theta code [CI SKIP]
evanmiller Jul 10, 2020
78fc9e2
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 13, 2020
98bc16c
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 17, 2020
2ecdf32
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 23, 2020
32a4d73
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 27, 2020
a2499bb
Add user documentation for Jacobi Theta functions [CI SKIP]
evanmiller Jul 29, 2020
041ca09
Add function graphs to Jacobi Theta docs [CI SKIP]
evanmiller Jul 30, 2020
3a0e3f3
Define Jacobi Theta test tolerances [CI SKIP]
evanmiller Jul 31, 2020
4630855
Merge branch 'develop' into jacobi-theta [CI SKIP]
evanmiller Jul 31, 2020
cf7ff40
Add implementation note on Jacobi theta functions [CI SKIP]
evanmiller Jul 31, 2020
b116c75
Consolidate Jacobi Theta ULPs plotting programs [CI SKIP]
evanmiller Aug 1, 2020
9fab61b
Fix q domain checking of jacobi_theta4 [CI SKIP]
evanmiller Aug 1, 2020
fb8d3f3
Add ULPs plots to Jacobi Theta docs [CI SKIP]
evanmiller Aug 1, 2020
dd2fcd1
Add missing Jacobi Theta ULPs plots [CI SKIP]
evanmiller Aug 1, 2020
d629f80
Add LaTeX source for Jacobi Theta equations [CI SKIP]
evanmiller Aug 1, 2020
ef39393
Remove unused Jacobi Theta PNG equations [CI SKIP]
evanmiller Aug 1, 2020
7ecda0f
Add Jacobi Theta performance script [CI SKIP]
evanmiller Aug 2, 2020
8cea3e9
Remove vestigial eps*eps check from jacobi_theta3 [CI SKIP]
evanmiller Aug 2, 2020
2902061
Update Jacobi Theta docs per code review comments [CI SKIP]
evanmiller Aug 2, 2020
5bea1e5
Enable arg promotion for Jacobi Theta functions [CI SKIP]
evanmiller Aug 3, 2020
46c59f8
Fix Jacobi Theta plotting script [CI SKIP]
evanmiller Aug 3, 2020
c9ac3a6
Change Jacobi Theta convergence criterion [CI SKIP]
evanmiller Aug 4, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
[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
  • Loading branch information
@evanmiller
evanmiller committed Jul 7, 2020
commit c78dcd72aa0fb9e6f0f484c10f68c8b262d9ace6
37 changes: 30 additions & 7 deletions include/boost/math/special_functions/jacobi_theta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,10 +166,25 @@ template <class RealType, class Policy>
inline RealType jacobi_theta4m1tau(RealType z, RealType tau, const Policy& pol);

#ifdef __JACOBI_THETA_USE_IMAGINARY
// The following _IMAGINARY theta functions assume imaginary z and are for
// internal use only. The z argument is scaled by tau, and sines and cosines
// are replaced with hyperbolic sines and cosines to accommodate the imaginary
// argument. (Recall sinh=(exp(x)-exp(-x))/2 and cosh=(exp(x)+exp(-x))/2)
//
// The return values are scaled by exp(-tau*z²/π)/sqrt(tau).
//
// Using the convention τ'=-1/τ and noting that both τ and τ' are always
// imaginary (at least in our applications), the _IMAGINARY functions are used
// to implement the following four relations:
//
// [20.7.30] sqrt(-iτ)θ₁(z|τ) = -i*exp(iτ'z²/π)*θ₁(zτ'|τ')
// [20.7.31] sqrt(-iτ)θ₂(z|τ) = exp(iτ'z²/π)*θ₄(zτ'|τ')
// [20.7.32] sqrt(-iτ)θ₃(z|τ) = exp(iτ'z²/π)*θ₃(zτ'|τ')
// [20.7.33] sqrt(-iτ)θ₄(z|τ) = exp(iτ'z²/π)*θ₂(zτ'|τ')
template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = RealType(0);
Expand Down Expand Up @@ -200,7 +215,7 @@ _IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy& pol) {
template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = RealType(0);
Expand All @@ -227,7 +242,7 @@ _IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy& pol) {
template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 1;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = exp(-z*z*tau/constants::pi<RealType>());
Expand All @@ -254,7 +269,7 @@ _IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy& pol) {
template <class RealType, class Policy>
inline RealType
_IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy& pol) {
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 1;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType /* q_n, */ delta, delta1, delta2, result = exp(-z*z*tau/constants::pi<RealType>());
Expand Down Expand Up @@ -289,7 +304,7 @@ template <class RealType, class Policy>
inline RealType
jacobi_theta1tau(RealType z, RealType tau, const Policy& pol)
Copy link
Collaborator

@jzmaddock jzmaddock Aug 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Point of pedantry: I'd prefer it if the actual entry points followed the logic in https://www.boost.org/doc/libs/1_73_0/libs/math/doc/html/math_toolkit/special_tut/special_tut_impl.html, what we gain is the ability to handle arguments of different types (including integers), plus expression templates when used with multiprecision types.

On a related note, I'd like to see the entry points added to this test: https://github.com/boostorg/math/blob/develop/test/compile_test/instantiate.hpp,

Also note that each special function appears in multiple sections in that file - you will very likely get test failures when doing this, but it does catch a lot of issues - including now, not forwarding policies correctly ;)

Copy link
Contributor Author

@evanmiller evanmiller Aug 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be clear - should I be using promote_args and friends in the public API? (I was advised against this previously.)

I've added the entry points (privately) to instantiate.hpp - lots of failures - will keep me busy tomorrow.

Copy link
Contributor Author

@evanmiller evanmiller Aug 3, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How's it look now?

Copy link
Collaborator

@jzmaddock jzmaddock Aug 3, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be clear - should I be using promote_args and friends in the public API? (I was advised against this previously.)

Yeah, my bad, you had a complex number interface before and promote_args doesn't work (easily) with that. Now that it's all real-valued, might as well be consistent with everything else.

aside 1: promote_args is mis-named, should really be "make_common_floating_point_type" or something.
aside 2: I guess it would help if I updated it to function with complex number types as well!

I've added the entry points (privately) to instantiate.hpp - lots of failures - will keep me busy tomorrow.

Likely most will disappear once the argument type handling is in place.

How's it look now?

You're way ahead of me, will catch up with this later!

{
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType q_n, delta, result = RealType(0);
Expand Down Expand Up @@ -341,6 +356,7 @@ inline T jacobi_theta1tau(T z, T tau) {
template <class RealType, class Policy>
inline RealType
jacobi_theta1(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (q <= 0.0 || q >= 1.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta1<%1%>(%1%)",
"q must be greater than 0 and less than 1 but got %1%.", q, pol);
Expand All @@ -359,7 +375,7 @@ template <class RealType, class Policy>
inline RealType
jacobi_theta2tau(RealType z, RealType tau, const Policy& pol)
{
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING
unsigned n = 0;
RealType eps = policies::get_epsilon<RealType, Policy>();
RealType q_n, delta, result = RealType(0);
Expand Down Expand Up @@ -407,6 +423,7 @@ inline T jacobi_theta2tau(T z, T tau) {
template <class RealType, class Policy>
inline RealType
jacobi_theta2(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (q <= 0.0 || q >= 1.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta2<%1%>(%1%)",
"q must be greater than 0 and less than 1 but got %1%.", q, pol);
Expand Down Expand Up @@ -458,6 +475,7 @@ template <class RealType, class Policy>
inline RealType
jacobi_theta3tau(RealType z, RealType tau, const Policy& pol)
{
BOOST_MATH_STD_USING
if (tau <= 0.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta3tau<%1%>(%1%)",
"tau must be greater than 0 but got %1%.", tau, pol);
Expand Down Expand Up @@ -488,6 +506,7 @@ inline T jacobi_theta3tau(T z, T tau) {
template <class RealType, class Policy>
inline RealType
jacobi_theta3m1(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (q <= 0.0 || q >= 1.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta3m1<%1%>(%1%)",
"q must be greater than 0 and less than 1 but got %1%.", q, pol);
Expand All @@ -505,6 +524,7 @@ inline T jacobi_theta3m1(T z, T q) {
template <class RealType, class Policy>
inline RealType
jacobi_theta3(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (q <= 0.0 || q >= 1.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta3<%1%>(%1%)",
"q must be greater than 0 and less than 1 but got %1%.", q, pol);
Expand All @@ -524,7 +544,7 @@ template <class RealType, class Policy>
inline RealType
jacobi_theta4m1tau(RealType z, RealType tau, const Policy& pol)
{
BOOST_MATH_STD_USING
BOOST_MATH_STD_USING

RealType eps = policies::get_epsilon<RealType, Policy>();
RealType q_n, delta, result = RealType(0);
Expand Down Expand Up @@ -557,6 +577,7 @@ template <class RealType, class Policy>
inline RealType
jacobi_theta4tau(RealType z, RealType tau, const Policy& pol)
{
BOOST_MATH_STD_USING
if (tau <= 0.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta4tau<%1%>(%1%)",
"tau must be greater than 0 but got %1%.", tau, pol);
Expand Down Expand Up @@ -589,6 +610,7 @@ inline T jacobi_theta4tau(T z, T tau) {
template <class RealType, class Policy>
inline RealType
jacobi_theta4m1(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (q <= 0.0 || q >= 1.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta4m1<%1%>(%1%)",
"q must be greater than 0 and less than 1 but got %1%.", q, pol);
Expand All @@ -606,6 +628,7 @@ inline T jacobi_theta4m1(T z, T q) {
template <class RealType, class Policy>
inline RealType
jacobi_theta4(RealType z, RealType q, const Policy& pol) {
BOOST_MATH_STD_USING
if (abs(q) >= 1.0 || abs(q) == 0.0) {
return policies::raise_domain_error<RealType>("boost::math::jacobi_theta4<%1%>(%1%)",
"|q| must be greater than zero and less than 1, but got %1%.", q, pol);
Expand Down
12 changes: 8 additions & 4 deletions test/test_jacobi_theta.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@

#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
Expand Down Expand Up @@ -34,6 +35,7 @@ BOOST_AUTO_TEST_CASE( test_main )
{
expected_results();
BOOST_MATH_CONTROL_FP;
BOOST_MATH_STD_USING

using namespace boost::math;

Expand Down Expand Up @@ -61,27 +63,29 @@ BOOST_AUTO_TEST_CASE( test_main )
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) {
double eps = std::numeric_limits<double>::epsilon();
test_periodicity(z, q, sqrt(eps));
test_argument_translation(z, q, 1.1 * sqrt(eps));
test_sums_of_squares(z, q, sqrt(eps));
// The addition formula is complicated, cut it some extra slack
test_addition_formulas(z, M_LN2, q, sqrt(sqrt(eps)));
test_addition_formulas(z, constants::ln_two<double>(), q, sqrt(sqrt(eps)));
test_duplication_formula(z, q, sqrt(eps));
test_transformations_of_nome(z, q, sqrt(eps));
test_watsons_identities(z, 0.5, q, sqrt(eps));
test_landen_transformations(z, -log(q)/M_PI, sqrt(eps));
test_landen_transformations(z, -log(q)/constants::pi<double>(), sqrt(eps));
}
}

test_special_values(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(boost::math::concepts::real_concept(0), "real_concept");
test_spots(concepts::real_concept(0), "real_concept");
#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
Expand Down
29 changes: 28 additions & 1 deletion test/test_jacobi_theta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
// Copyright John Maddock 2006.
// Copyright Evan Miller 2020
// #include <boost/math/concepts/real_concept.hpp>
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
Expand Down Expand Up @@ -478,3 +477,31 @@ inline void test_landen_transformations(RealType z, RealType tau, RealType eps)
* 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);
}