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/test_jacobi_theta.cpp

@@ -75,7 +75,9 @@ BOOST_AUTO_TEST_CASE( test_main )
             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);

test/test_jacobi_theta.hpp

@@ -3,6 +3,8 @@
 #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>
@@ -570,3 +572,84 @@ inline void test_mellin_transforms(RealType s, RealType integration_eps, RealTyp
             (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_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);
+    }
+}