Jacobi Theta functions

[evanmiller]

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.

Any feedback appreciated.

[jzmaddock]

That looks really good!

You can save me some time by providing a couple of concrete examples which currently have low accuracy?

[evanmiller]

@jzmaddock I can provide specific examples later today but the failing spot tests (../../../b2 test_jacobi_theta) should point you in the right direction. Also see the ULP plotting programs added to reporting/accuracy. Try large q / small tau (say q = 0.99 or tau = 0.001) with non-zero z (say 1.5) to trigger the _IMAGINARY code paths on any of the public functions.

To be clear: the real-valued code paths are inaccurate for large q. The _IMAGINARY code paths are accurate but imprecise.

Just so you know: the _IMAGINARY functions are rescaled and rearranged a bit - they assume imaginary z, so the sines and cosines became hyperbolic. I then tried to be clever by consolidating the exponentials. Because the identity tests pass I believe the expressions are mathematically correct, but they're losing precision somehow. I suspect either the exp(-z*z) or sqrt(tau) but I am honestly out of my depth here.

Will provide actual examples tonight if I haven't heard from you.

[jzmaddock]

OK I have a test case, thanks.

Minor nits while I think of them:

• test_jacobi_theta.cpp needs to include boost/math/concepts/real_concept.hpp to build.
• The tests use some constants like M_PI which are Unix'isms: please use boost::math::constants::pi() for portability.
• There are a bunch of places where you're accidentally calling ::log(double) (for example) rather than std::log(T) which will quite adversely effect the long double precisions. Adding the using declarations or BOOST_MATH_STD_USING will fix these. Instantiating the functions with boost::math::concepts::std_real_concept will help find all the locations.

After that the error rates don't look bad (200eps or so), but.... as an example jacobi_theta_2(3, 0.9921875) is so totally trivial to calculate (done in one iteration, and I don't see any obvious things that would loose precision) that I would like to investigate some more.

One possible improvement is in:

template <class RealType, class Policy>
inline RealType
jacobi_theta2(RealType z, RealType q, const Policy& pol) {
    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);
    }
    return jacobi_theta2tau(z, -log(q)/constants::pi<RealType>(), pol);
}

where log(q) might be better expressed as log1p(q-1) when q is ultra close to 1, but makes no difference in the case above. And I'm not completely convinced it ever would make a difference, I need to think about that ;)

[jzmaddock]

One more nit: if you instantiate with for example number<debug_adaptor<cpp_bin_float_50::backend_type>> ie a type that uses expression templates, then there are quite a few build errors where expressions are passed as function arguments.

The error in the case above come entirely from the line

delta1 = exp(tau*(z*(RealType(2*n) - z/constants::pi<RealType>())
                -constants::pi<RealType>() * RealType(n)*RealType(n)));

The issue being that z(2n-z/pi) happens to evaluate to very close to pi in this case. I don't see any algebraic simplification though?

[evanmiller]

• test_jacobi_theta.cpp needs to include boost/math/concepts/real_concept.hpp to build.

Thanks

• The tests use some constants like M_PI which are Unix'isms: please use boost::math::constants::pi() for portability.

Will do

• There are a bunch of places where you're accidentally calling ::log(double) (for example) rather than std::log(T) which will quite adversely effect the long double precisions. Adding the using declarations or BOOST_MATH_STD_USING will fix these. Instantiating the functions with boost::math::concepts::std_real_concept will help find all the locations.

Aha! I never would have thought of that. It likely affects the use of fmod as well - this would basically reduce long doubles to doubles (for the z argument) on the _IMAGINARY path. Since my ULP plots use long doubles for the high-precision case, this would basically wreck all the plots.

I'll make these changes and report back.

[NAThompson]

@evanmiller : It'll probably a few iterations on this, so make sure to tag your commits with [CI SKIP] so we don't DDOS the CI system.

Also, could you drop a performance benchmark into reporting/performance? This one is a good template, but I'd also template on float, double and long double.

[evanmiller]

I've made the changes suggested and added a few more tests covering special values (did you know the third theta function is related to the gamma function?). I don't see any perceptible change to the test results though.

@NAThompson I will look into adding a performance benchmark later tonight.

[evanmiller]

Good news! I think I've just about cracked the precision issues. I completed the square on the exponentials, and rewrote them as double-sided summations following DLMF 20.13.4 and 20.13.5.

The mean errors on the MathWorld data are now down to <10 (eps? ULPs? whatever test_hetero reports). They were 50-200 before. The ULP plots show the precision getting steadily worse as z grows - so there's still room for improvement. I should have a better picture on the overall accuracy once the 1000-bit code gets plugged in.

BTW What are we benchmarking against exactly? I don't see Jacobi theta functions in either GSL or RMath.

Commit incoming.

[NAThompson]

BTW What are we benchmarking against exactly?

Ourselves! Run the benchmark under perf record and examine the asm; you'll enjoy it.

[NAThompson]

If you want, you can use my "math_unit_test.hpp" and use CHECK_ULP_CLOSE if you think it'll be a bit more ergonomic. It also compiles much faster that Boost.Test.

[NAThompson]

It looks like your implementation is now correct; all values are basically in the condition number envelope.

[NAThompson]

It's jacobi3 that confuses me. The error is pretty good, but it's uncorrelated with the condition number envelope, which feels very strange.

[NAThompson]

I'm not sure this comment is very informative; I was hoping for something like "expected error is cond(f; x)eps/2, so when delta <..." and then some perturbative analysis.

[evanmiller]

The comment represents the extent of my thinking and, well, education. I'm happy to incorporate a more formal convergence analysis but it's simply outside my ken.

[NAThompson]

lol, alright np. I'll see if I can get a minute to work out some error analysis for these functions.

[jzmaddock]

IMO there's no need to require delta < eps, normally we would just use abs(result * eps) > abs(delta) as if that's true, then adding delta can have no effect on result. The only exception would be if delta can grow as well as shrink.

[evanmiller]

The sign of delta can alternate, which could theoretically put result at 0 after two iterations. So I wanted to handle that somehow, and distinguish the situation from an underflow.

[NAThompson]

Wait then aren't all these calls to abs superfluous?

[evanmiller]

Sorry, I should have said: the numerator is q_n (always positive) but the denominator is result (which may oscillate). It would probably be better to change the denominator to previous_q_n.

[NAThompson]

I personally might track \sum q^{n^2} since it provides a nice L1 bound on the sum.

[evanmiller]

I'm not sure the sum will work because of the alternating signs. In addition to the oscillating component there is the (-1)^n in two of the formulas. With q close to 1, the straight sum could be large relative to the final result, and so an error that is small relative to the sum may be large relative to the result.

[evanmiller]

I have updated the convergence criterion to compare q_n to the previous one. Everything including ULP plots look pretty much the same. Let me know what you think.

[NAThompson]

Looking at this graph, it appears that the x->0 accuracy is not limited by the condition number of function evaluation and hence it could be implemented better. Is this a case to recommend use of jacobi_theta1tau?

[evanmiller]

tau can help with precision issues on q but not x.

[NAThompson]

I'm confused. You say you enforce that tau is purely imaginary, yet your function signature says tau is the same type as x. . .

[evanmiller]

tau the C++ argument is real-typed. The provided number is implicitly multiplied by i. This seemed simpler than trying to provide a complex-valued API.

[NAThompson]

I don't understand what this file is doing. It appears you are spitting out values of the function from your own code; but isn't that tautologous?

[evanmiller]

This is what the other tools/*_data.cpp files do – use a high-precision implementation to compare against the low-precision implementations.

[jzmaddock]

Sort of, here's how I approach this:

Let's say we have a principal method, plus a bunch of special case handling (asymptotic expansions at one end of the range say as an example). The main things that can go wrong:

1. The principal method has an outright bug (eg terms in series not correctly calculated).
2. The principal method looses precision due to either cancellation or an incorrect termination condition.
3. There are bugs in the special cases, or the selection logic between them.

We can deal with (1) by comparison to a "known good" - so job done on that I think.

For (2) and (3) what I try and do is calculate spot values at very high precision using the principal method alone, the important thing is to remove all the special case handling as otherwise as @NAThompson says you're not really testing much. In an ideal world you would use a completely different method, but those are really hard to find. When I remember, what I should have done (but didn't always) is archive a "naive" version of the principal method in the test value generator file and call that rather than the actual method - I suspect I may have been sloppy and just quickly commented stuff out in the header occasionally though :(

However, with regard to different methods: do the defining series for these functions satisfy that? They don't have to be fast for test value generation as they're only run the once. In extreme cases (1F1) I've resorted to interval arithmetic using mpfi_float to handle series with lots of cancellation when generating spot tests.

[evanmiller]

That makes sense. We have essentially two implementations of each function, the IMAGINARY path (for tau < 1) and the regular path. The code chooses the one that is expected to converge more quickly. So one possibility would be to invert the branch logic so that the test data always chooses the slower path when generating test data.

[NAThompson]

@evanmiller : Could you drop this file into reporting/performance/jacobi_theta_performance.cpp?

//  (C) Copyright Nick Thompson 2020.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <random>
#include <benchmark/benchmark.h>
#include <boost/math/special_functions/jacobi_theta.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

using boost::multiprecision::number;
using boost::multiprecision::mpfr_float_backend;
using boost::multiprecision::float128;
using boost::multiprecision::cpp_bin_float_50;
using boost::multiprecision::cpp_bin_float_100;
using boost::math::jacobi_theta1;
using boost::math::jacobi_theta1tau;

template<class Real>
void JacobiTheta1(benchmark::State& state)
{
    std::random_device rd;
    std::mt19937_64 mt(rd());
    std::uniform_real_distribution<long double> unif(0,0.01);

    Real x = static_cast<Real>(unif(mt));
    Real q = static_cast<Real>(unif(mt));
    for (auto _ : state)
    {
        benchmark::DoNotOptimize(jacobi_theta1(x, q));
        x += std::numeric_limits<Real>::epsilon();
    }
}

BENCHMARK_TEMPLATE(JacobiTheta1, float);
BENCHMARK_TEMPLATE(JacobiTheta1, double);
BENCHMARK_TEMPLATE(JacobiTheta1, long double);
BENCHMARK_TEMPLATE(JacobiTheta1, float128);
BENCHMARK_TEMPLATE(JacobiTheta1, number<mpfr_float_backend<100>>);
BENCHMARK_TEMPLATE(JacobiTheta1, number<mpfr_float_backend<200>>);
BENCHMARK_TEMPLATE(JacobiTheta1, number<mpfr_float_backend<300>>);
BENCHMARK_TEMPLATE(JacobiTheta1, number<mpfr_float_backend<400>>);
BENCHMARK_TEMPLATE(JacobiTheta1, number<mpfr_float_backend<1000>>);
BENCHMARK_TEMPLATE(JacobiTheta1, cpp_bin_float_50);
BENCHMARK_TEMPLATE(JacobiTheta1, cpp_bin_float_100);

template<class Real>
void JacobiTheta1Tau(benchmark::State& state)
{
    std::random_device rd;
    std::mt19937_64 mt(rd());
    std::uniform_real_distribution<long double> unif(0,0.01);

    Real x = static_cast<Real>(unif(mt));
    Real q = static_cast<Real>(unif(mt));
    for (auto _ : state)
    {
        benchmark::DoNotOptimize(jacobi_theta1tau(x, q));
        x += std::numeric_limits<Real>::epsilon();
    }
}

BENCHMARK_TEMPLATE(JacobiTheta1Tau, float);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, double);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, long double);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, float128);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, number<mpfr_float_backend<100>>);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, number<mpfr_float_backend<200>>);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, number<mpfr_float_backend<300>>);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, number<mpfr_float_backend<400>>);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, number<mpfr_float_backend<1000>>);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, cpp_bin_float_50);
BENCHMARK_TEMPLATE(JacobiTheta1Tau, cpp_bin_float_100);

BENCHMARK_MAIN();

Also, running this file, I find that something has gone wrong with the ADL for mpfr types, there are some warnings from unused policies from g++-10.

[NAThompson]

Is there equivalent Mathematica syntax that has the same nome/tau/scaling conventions? It might be nice to put that in here . . .

I see EllipticTheta[a,q,z]

[evanmiller]

I don't think Mathematica provides a tau parameterization. I'll add a link in the docs to EllipticTheta, though.

[NAThompson]

I'm getting more confused about the tau parametrization. It now looks to me like you are trading a multiplication on each iteration for a logarithm and a call to std::exp on each iteration. . .

[evanmiller]

Implementing purely in q, without a log, would require a pow in place of the exp. The private IMAGINARY code paths have to do the exp anyway since they are implementing hyperbolic sine and cosine.

Exposing tau has the added benefit of letting the user provide a pre-logged argument for the sake of precision, e.g. with q close to 1 (so the user can use log1p) or if the argument provided would have been an exponential such as exp(-a).

[NAThompson]

You don't need a call to pow, since q^(n+1)^2 = q^n^2qnqn*q. Here's some pseudocode:

template<typename Real>
Real new_theta1(Real x, Real q)
{
    using std::sqrt;
    using std::sin;
    using std::abs;
    Real scale = 2*sqrt(sqrt(q));
    Real term = sin(x);
    Real sum = term;
    size_t n = 1;
    Real qn = q;
    Real qnsq = q;
    Real prefactor = qnsq*qn;
    while (prefactor > abs(sum)*std::numeric_limits<Real>::epsilon())
    {
        term = prefactor*sin((2*n+1)*x);
        if (n & 1) {
            sum -= term;
        } else {
            sum += term;
        }
        qn *= q;
        qnsq = qnsq*q*qn*qn;
        prefactor = qnsq*qn;
        ++n;
    }
    return scale*sum;
}

This executes in 13ns on my machine. More improvements are available from a Clenshaw recurrence.

[evanmiller]

Very cool! Happy to explore performance improvements after everything else is settled.

[evanmiller]

@evanmiller : Could you drop this file into reporting/performance/jacobi_theta_performance.cpp?

Will do!

[NAThompson]

Maybe precompute tau/constants::pi<RealType>() so that there isn't a division on each iteration.

People don't realize that elementary special function evaluations are in the same complexity class as division!

[evanmiller]

If we're going to squeeze out divisions like that, it might make sense to divide tau by pi further up the call stack.

[NAThompson]

Calling sin on each iteration seems expensive; can't it be done via a Clenshaw recurrence?

[evanmiller]

The iteration count is typically between 2 and 5 so I'm not sure a recurrence formula will buy much (also because there's already an exp in each iteration).

[NAThompson]

Ok, so most Clenshaw recurrence discussions descend rather than ascend, so we instead can use

sin((2(n+1)+1)x) = sin((2n+1)x)*cos(2x) + cos((2n+1)x)*sin(2x)

[NAThompson]

Combine this with the recurrence for q^n^2 and you can do the whole thing with a single call to std::cos and std::sin and then only muls and adds after that.

[NAThompson]

Clenshaw recurrence instead of cos on each iteration?

[jzmaddock]

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 ;)

[evanmiller]

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.

[evanmiller]

How's it look now?

[jzmaddock]

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!

[jzmaddock]

This looks like it's ready to go? @NAThompson ?

[NAThompson]

I'm happy with this; the only potential point of improvement is using the trigonometric recurrence and the q^{n^2} recurrence to increase the speed. But perfect is the enemy of the good.

[evanmiller]

@NAThompson's ideas may also yield accuracy improvements – but those can always take the form of a future PR. So you have my blessing to merge as well.

[pulver]

Are the build errors on https://travis-ci.org/github/boostorg/math/builds/718417184 related to this PR?

[NAThompson]

My bad; I told @evanmiller to [CI SKIP] until he was ready to merge, and then I just merged it. So looks like this one is my fault.

[jzmaddock]

Testing the "obvious" fixes now....

[jzmaddock]

Tentative fixes pushed: just looks like the usual concept-failures. Glad someone was paying attention to develop results as I missed that!