elliptic_theta.hpp

Go to the documentation of this file.

#pragma once
 
#include "math.hpp"
#include <cstddef>
#include <limits>
#include <stdexcept>
 
namespace uni20
{
template <class R> struct theta_constant_result
{
    R scaled_value;
    R log_scale;
};
 
namespace detail
{
template <class R>
concept ThetaReal = native_math::NativeReal<R>
#if UNI20_HAS_MPREAL
                    || std::same_as<R, mpreal>
#endif
    ;
 
template <ThetaReal R> theta_constant_result<R> theta_constant_series(unsigned kind, R const& t, std::size_t max_terms)
{
  auto const p = precision_of(t);
  auto number = [p](auto n) { return scalar_at_precision<R>(n, p); };
  R const zero = number(0), one = number(1), two = number(2);
  R const pi_value = uni20::pi<R> * one;
  bool const dual = t < one;
  // S-transform: theta2 <-> theta4, theta3 -> theta3, with factor 1/sqrt(t).
  unsigned const series_kind = dual && kind != 3 ? 6 - kind : kind;
  R scale = dual ? -math::log(t) / two : zero;
  if (series_kind == 2)
  {
    // Divide pi first: pi*t can overflow while pi*t/4 is representable.
    R const leading = dual ? (pi_value / number(4)) / t : (pi_value / number(4)) * t;
    scale -= leading;
  }
  if (!uni20::isfinite(scale)) throw std::overflow_error("theta constant logarithmic scale exceeds scalar range");
 
  // An infinite decay is harmless: every nonconstant series term then rounds
  // to zero. Never form exp(+decay) or subtract nearly equal theta constants.
  R const decay = dual ? pi_value / t : pi_value * t;
  R const target = numeric_limits<R>::epsilon(t) / number(16);
  if (!(target > zero)) throw std::underflow_error("theta constant series tolerance exceeds scalar range");
  R sum = series_kind == 2 ? two : one, correction = zero, n = one;
  for (std::size_t index = 0; index < max_terms; ++index)
  {
    R const power = series_kind == 2 ? n * (n + one) : n * n;
    R const term = two * math::exp(-decay * power);
    // In the reduced domain t>=1 successive absolute terms have ratio at
    // most exp(-3*pi)<1/2. The entire omitted tail is less than twice term.
    // This bounds truncation in exact arithmetic, not floating-point rounding.
    if (term <= target) return {sum + correction, scale};
    R const signed_term = series_kind == 4 && index % 2 == 0 ? -term : term;
    R const next = sum + signed_term;
    correction += (sum - next) + signed_term; // |sum| > |signed_term|.
    sum = next;
    n += one;
  }
  throw std::runtime_error("theta constant series exhausted its term budget");
}
} // namespace detail
 
template <detail::ThetaReal R>
theta_constant_result<R> elliptic_theta_constant(unsigned kind, R const& t, std::size_t max_terms = 4096)
{
  if (kind < 2 || kind > 4 || max_terms == 0)
    throw std::invalid_argument("theta constant requires kind 2..4 and a positive term budget");
  auto const p = precision_of(t);
  if (p.is_exact()) throw std::invalid_argument("theta constant requires explicit finite precision");
  if (!uni20::isfinite(t) || !(t > 0)) throw std::invalid_argument("theta constant requires finite t>0");
#if UNI20_HAS_MPREAL
  if constexpr (std::same_as<R, mpreal>)
  {
    if (p.bit_count() > std::numeric_limits<Precision::bit_count_type>::max() - 32)
      throw std::overflow_error("theta constant guard precision exceeds supported range");
    auto const work = Precision::bits(p.bit_count() + 32);
    auto result = detail::theta_constant_series(kind, t.at(work), max_terms);
    theta_constant_result<R> rounded{result.scaled_value.at(p), result.log_scale.at(p)};
    if (!uni20::isfinite(rounded.log_scale))
      throw std::overflow_error("theta constant logarithmic scale exceeds result range");
    return rounded;
  }
  else
#endif
    return detail::theta_constant_series(kind, t, max_terms);
}
 
template <detail::ThetaReal R>
theta_constant_result<R> elliptic_theta_constant(unsigned kind, R const& t, Precision p, std::size_t max_terms = 4096)
{
  if (p.is_exact()) throw std::invalid_argument("theta constant requires explicit finite precision");
  return elliptic_theta_constant(kind, scalar_at_precision<R>(t, p), max_terms);
}
} // namespace uni20