Skip to content

Size the T_l recurrence passes of the analytic add-back by their growth rate #836

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

Open
CharlesCNorton wants to merge 2 commits into
base: main
Choose a base branch
from
+94 −137
Open

Size the T_l recurrence passes of the analytic add-back by their growth rate #836

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4b8e4b2
Size the T_l recurrence passes of the analytic add-back by their grow…
CharlesCNorton Sep 8, 2026
4827da7
Merge branch 'main' into nestor-extended-precision
CharlesCNorton Sep 9, 2026
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
168 changes: 58 additions & 110 deletions src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,66 @@
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <numbers>
#include <vector>

namespace vmecpp {

namespace {

// Logarithm of the growth of the homogeneous solutions of the T_l recurrence
// over a forward pass, sqrt(B/A)^kL, above which the recurrence is run
// backward instead; a forward pass amplifies the rounding of its inputs by at
// most this factor, ten.
constexpr double kMaxForwardLogGrowth = std::numbers::ln10;

// The zero seed of a backward pass contaminates T_l by (A/B)^{(top - l)/2} of
// T_top; the pass starts far enough above kL to bring that below 1e-17 at kL.
constexpr double kMinSeedLogDecay = 17.0 * std::numbers::ln10;

// T_l = int_{-1}^{1} t^l / sqrt(A t^2 + 2 d t + B) dt for l = 0, ..., kL at
// tangential grid point kl from the three-term recurrence
// (l + 1) A T_{l+1} + (2 l + 1) d T_l + l B T_{l-1} = sqrtc2 - (-1)^l sqrta2
// with T_0 given. The characteristic roots of the homogeneous recurrence are a
// complex pair of modulus sqrt(B/A), so a forward pass amplifies the rounding
// of T_0 and of the right-hand sides by sqrt(B/A)^l and is used while that
// stays below exp(kMaxForwardLogGrowth); otherwise the recurrence runs
// backward from a zero seed, whose contamination decays by sqrt(A/B) per step.
void ComputeTl(double A, double B, double d, double sqrtc2, double sqrta2,
double T0, int kL, int kl, std::vector<Eigen::VectorXd>& m_T) {
const auto rhs = [sqrtc2, sqrta2](int l) {
return sqrtc2 + (l % 2 == 0 ? -sqrta2 : sqrta2);
};
m_T[0][kl] = T0;
const double log_ratio = std::log(B / A);
if (kL * log_ratio <= 2.0 * kMaxForwardLogGrowth) {
double T_prev = 0.0; // T_{-1}
for (int l = 0; l < kL; ++l) {
const double T_next =
(rhs(l) - (2 * l + 1) * d * m_T[l][kl] - l * B * T_prev) /
((l + 1) * A);
T_prev = m_T[l][kl];
m_T[l + 1][kl] = T_next;
}
return;
}
const int tail =
static_cast<int>(std::ceil(2.0 * kMinSeedLogDecay / log_ratio));
double T_hi = 0.0; // T_{l+1}
double T_cur = 0.0; // T_l
for (int l = kL + tail; l >= 1; --l) {
const double T_lo =
(rhs(l) - (2 * l + 1) * d * T_cur - (l + 1) * A * T_hi) / (l * B);
T_hi = T_cur;
T_cur = T_lo;
if (1 <= l - 1 && l - 1 <= kL) {
m_T[l - 1][kl] = T_lo;
}
}
}

} // namespace

SingularIntegrals::SingularIntegrals(const Sizes* s,
const FourierBasisFastToroidal* fb,
const TangentialPartitioning* tp,
Expand Down Expand Up @@ -220,117 +276,9 @@ void SingularIntegrals::prepareUpdate(
(sqrtam * sqrta2[kl] - am[kl] + d[kl])) /
sqrtam;

// Fill all Tlp[0..L] and Tlm[0..L] by picking the numerically stable
// direction of the three-term recurrence on a per-(+/-), per-kl basis.
//
// The characteristic roots of the homogeneous recurrence satisfy
// A*r^2 + 2*d*r + B = 0 -> |r1 r2| = B/A.
// For T^+: (A, B) = (ap, am), so |r1 r2| = am/ap.
// For T^-: (A, B) = (am, ap), so |r1 r2| = ap/am.
// If B > A (at least one |r| > 1), forward iteration is unstable and
// backward (Miller's algorithm) is used instead; otherwise forward is fine.
//
// T^{\pm}_0 is analytic (above); T^{\pm}_{-1} = 0. Forward produces
// T_{l+1} from T_l and T_{l-1}; backward produces T_{l-1} from T_l and
// T_{l+1} via the same recurrence solved in reverse. For backward,
// iteration starts from a zero seed far above the required L; the result
// is then normalized to match the analytic T^{\pm}_0.
//
// rhs(l+1) = sqrtc2 + (-1)^{l+1}*sqrta2 (same for T^+ and T^-).
const int kL = mf + nf;
// The spurious solution is damped by (A/B)^kTailExtra per pass.
// For the worst realistic ratio (A/B ~ 0.5) suppression is ~0.5^50 ~ 1e-16.
const int kTailExtra = 50;
const int kLtail = kL + kTailExtra;

// Only switch to backward when the forward spurious-mode growth
// (|r1 r2| = B/A) would actually exceed double precision over kL steps.
// Threshold: forward is considered stable as long as (B/A)^kL < 1e10,
// i.e. spurious amplitude stays within ~1e10 of the particular solution.
// Near-degenerate kl (|r1|~|r2|~1) fall in the forward branch, where
// zero-seed Miller is known to misconverge (spurious modes never damp).
// Formula: kL * ln(B/A) < ln(1e10) -> B/A < exp(ln(1e10)/kL).
constexpr double kLogGrowthThreshold = 10.0 * 2.30258509299; // ln(1e10)
const double logRatioP =
(am[kl] > ap[kl] && ap[kl] > 0.0) ? std::log(am[kl] / ap[kl]) : 0.0;
const bool useBackwardP =
static_cast<double>(kL) * logRatioP > kLogGrowthThreshold;
const double logRatioM =
(ap[kl] > am[kl] && am[kl] > 0.0) ? std::log(ap[kl] / am[kl]) : 0.0;
const bool useBackwardM =
static_cast<double>(kL) * logRatioM > kLogGrowthThreshold;

// --- T^+: A = ap, B = am ---
Tlp[0][kl] = T0p;
if (useBackwardP) {
// forward unstable -> use backward recurrence.
double T_hi = 0.0;
double T_cur = 1.0e-300;
for (int l = kLtail; l >= 1; --l) {
const double rhs = sqrtc2[kl] + (l % 2 == 0 ? -1.0 : 1.0) * sqrta2[kl];
const double T_lo =
(rhs - (2 * l + 1) * d[kl] * T_cur - (l + 1) * ap[kl] * T_hi) /
(l * am[kl]);
T_hi = T_cur;
T_cur = T_lo;
if (l - 1 <= kL) {
Tlp[l - 1][kl] = T_lo;
}
}
const double scaleP = T0p / Tlp[0][kl];
for (int l = 0; l <= kL; ++l) {
Tlp[l][kl] *= scaleP;
}
} else {
// forward stable.
double T_prev = 0.0; // T^+_{-1}
int sgn = 1;
for (int fl = 0; fl < kL; ++fl) {
sgn = -sgn;
const double rhs = sqrtc2[kl] + sgn * sqrta2[kl];
const double T_next =
(rhs - (2 * fl + 1) * d[kl] * Tlp[fl][kl] - fl * am[kl] * T_prev) /
(ap[kl] * (fl + 1));
T_prev = Tlp[fl][kl];
Tlp[fl + 1][kl] = T_next;
}
}

// --- T^-: A = am, B = ap ---
Tlm[0][kl] = T0m;
if (useBackwardM) {
// forward unstable -> use backward recurrence.
double T_hi = 0.0;
double T_cur = 1.0e-300;
for (int l = kLtail; l >= 1; --l) {
const double rhs = sqrtc2[kl] + (l % 2 == 0 ? -1.0 : 1.0) * sqrta2[kl];
const double T_lo =
(rhs - (2 * l + 1) * d[kl] * T_cur - (l + 1) * am[kl] * T_hi) /
(l * ap[kl]);
T_hi = T_cur;
T_cur = T_lo;
if (l - 1 <= kL) {
Tlm[l - 1][kl] = T_lo;
}
}
const double scaleM = T0m / Tlm[0][kl];
for (int l = 0; l <= kL; ++l) {
Tlm[l][kl] *= scaleM;
}
} else {
// forward stable.
double T_prev = 0.0; // T^-_{-1}
int sgn = 1;
for (int fl = 0; fl < kL; ++fl) {
sgn = -sgn;
const double rhs = sqrtc2[kl] + sgn * sqrta2[kl];
const double T_next =
(rhs - (2 * fl + 1) * d[kl] * Tlm[fl][kl] - fl * ap[kl] * T_prev) /
(am[kl] * (fl + 1));
T_prev = Tlm[fl][kl];
Tlm[fl + 1][kl] = T_next;
}
}
ComputeTl(ap[kl], am[kl], d[kl], sqrtc2[kl], sqrta2[kl], T0p, kL, kl, Tlp);
ComputeTl(am[kl], ap[kl], d[kl], sqrtc2[kl], sqrta2[kl], T0m, kL, kl, Tlm);
} // kl
} // prepareUpdate

Expand Down
63 changes: 36 additions & 27 deletions src/vmecpp/cpp/vmecpp/free_boundary/singular_integrals/singular_integrals_test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -164,24 +164,33 @@ TEST(TestSingularIntegrals, CheckConstants) {
} // CheckConstants

// Verify that prepareUpdate computes T^{+/-}_l accurately for ALL l in [0, L]
// at a range of Fourier resolutions.
// at a range of Fourier resolutions and metric coefficients.
//
// Resolutions span the low-kL regime (forward recurrence is accurate, so the
// forward branch of prepareUpdate is exercised) all the way up to high-kL
// (Miller backward recurrence fires; forward would lose all digits).
// The cases span the low-kL regime, where the forward recurrence is accurate,
// up to kL = 45, where one direction of each recurrence loses all digits and
// the other has to run backward from a seed far above kL.
//
// Reference: 64-point Gauss-Legendre quadrature of the defining integral.
// This is independent of both recurrence directions and reaches ~1e-13
// relative accuracy at the chosen coefficients up to l = 45.
class TlpTlmAccuracyTest
: public ::testing::TestWithParam<std::pair<int, int>> {};
struct TlCase {
int mpol;
int ntor;
// metric coefficients (a, b2, c) = (guu, guv, gvv) of the tangent plane
double a;
double b2;
double c;
};

class TlpTlmAccuracyTest : public ::testing::TestWithParam<TlCase> {};

TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) {
// GL-64 reference noise at l near kL is ~1e-13; 1e-11 leaves a safe margin
// while staying far below the forward-recurrence error that Miller corrects.
// while staying far below the error of a recurrence run in the wrong
// direction or started too close above kL.
static constexpr double kTolerance = 1.0e-11;

const auto [mpol, ntor] = GetParam();
const auto [mpol, ntor, a_val, b2_val, c_val] = GetParam();

const bool lasym = false;
const int nfp = 5;
Expand All @@ -200,21 +209,14 @@ TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) {
const int kL = mf + nf;
SingularIntegrals si(&s, &fb, &tp, &sg, nf, mf);

// Geometry coefficients chosen so am/ap ~ 4.7 and the discriminant
// d^2 - ap*am = 0.16 - 2.31 < 0 (so the integrand is smooth on [-1, 1]).
// ap = a + b2 + c = 0.7
// am = a - b2 + c = 3.3
// d = c - a = 0.4
// At these coefficients the forward recurrence loses ~kL * log10(am/ap)
// ~ 0.67 * kL significant digits; for kL >= 15 it loses > 10 digits and the
// Miller activation threshold (growth > 1e10) triggers.
const double a_val = 0.8;
const double b2_val = -1.3;
const double c_val = 1.2;
const double ap = a_val + b2_val + c_val;
const double am = a_val - b2_val + c_val;
const double d = c_val - a_val;
ASSERT_GT(am, ap) << "test setup: need am > ap to exercise Miller path";
// The homogeneous solutions of the T^+ (T^-) recurrence grow by
// sqrt(am/ap) (sqrt(ap/am)) per forward step, so one of the two has to run
// backward once kL is large enough.
ASSERT_GT(std::max(am / ap, ap / am), 2.0)
<< "test setup: need one unstable forward direction";
ASSERT_LT(d * d, ap * am) << "test setup: need smooth integrand on [-1,1]";

const int numLocal = tp.ztMax - tp.ztMin;
Expand Down Expand Up @@ -253,15 +255,22 @@ TEST_P(TlpTlmAccuracyTest, MatchesQuadrature) {
}
}

// (mpol, ntor) pairs spanning the forward-stable regime (kL=15, Miller just
// barely fires) up to high-kL where forward would lose >30 digits.
// The first three cases share the coefficients ap = 0.7, am = 3.3, d = 0.4
// (am/ap = 4.7) at kL = 15, 27 and 45; the last one takes the metric of the
// cth_like_free_bdy boundary at the point of its largest cross term
// (ap/am = 2.2) at kL = 33, where a backward pass seeded 50 steps above kL
// leaves 2e-9 of the seed in T^-_33.
INSTANTIATE_TEST_SUITE_P(
ResolutionSweep, TlpTlmAccuracyTest,
::testing::Values(std::pair<int, int>{6, 8}, std::pair<int, int>{12, 14},
std::pair<int, int>{20, 24}),
[](const ::testing::TestParamInfo<std::pair<int, int>>& info) {
return "mpol" + std::to_string(info.param.first) + "_ntor" +
std::to_string(info.param.second);
::testing::Values(TlCase{6, 8, 0.8, -1.3, 1.2},
TlCase{12, 14, 0.8, -1.3, 1.2},
TlCase{20, 24, 0.8, -1.3, 1.2},
TlCase{16, 16, 1.360e-02, 2.241e-02, 4.626e-02}),
[](const ::testing::TestParamInfo<TlCase>& info) {
return "mpol" + std::to_string(info.param.mpol) + "_ntor" +
std::to_string(info.param.ntor) + "_ap" +
std::to_string(static_cast<int>(
1000.0 * (info.param.a + info.param.b2 + info.param.c)));
});

} // namespace vmecpp
Loading