#!/usr/bin/env python3
"""Audit a Heun connection coefficient against its block expansion.

This script implements the nonresonant numerical slice used in Chapter 7,
page 8 of *Advanced Methods of Linear ODEs*.  It evaluates the
Lisovyy--Naidiuk three-term recurrence and its dual continued fraction for
the general Heun equation with

    theta_0 = 1/7, theta_1 = 2/9, theta_t = 1/5,
    theta_infinity = 3/11, omega = 2/5, lambda = 1/t.

The continued-fraction tail, the first-order infinite sum, and the forward
recurrence limit all converge algebraically with the cutoff.  Values at
K, 2K, 4K, ... are therefore extrapolated with the ordinary Richardson
table for a 1/K expansion.  The cutoff table and its last-step estimate are
printed rather than hidden behind a claimed number of correct digits.

The full optional 2 x 2 array is *source-first*: rows are the signs of
theta_0 and columns are the signs of theta_1 in

    psi^[0]_epsilon = sum_epsilon' C(epsilon theta_0,
                                      epsilon' theta_1) psi^[1]_epsilon'.

Status caveats:

* The recurrence/continued-fraction connection formula is rigorous on its
  stated nonresonant Heun chart; this program is only a numerical audit of
  one parameter slice.
* The classical-block expression is used here only through its analytic
  first coefficient.  Agreement does not by itself prove a global block
  identity or justify continuation beyond |lambda| < 1.
* Richardson extrapolation assumes the expected inverse-cutoff asymptotic
  expansion.  Increase both ``--base-cutoff`` and ``--levels`` before
  interpreting more digits than the displayed convergence estimate allows.

Only mpmath and Python's standard library are required.

Primary reference: O. Lisovyy and A. Naidiuk, *Perturbative connection
formulas for Heun equations*, J. Phys. A 55 (2022) 434005,
doi:10.1088/1751-8121/ac9ba7, arXiv:2208.01604.
"""

from __future__ import annotations

import argparse
from dataclasses import dataclass
import platform
import sys
from time import perf_counter
from typing import Callable, Iterable, Sequence

import mpmath as mp


class AuditError(RuntimeError):
    """A parameter lies outside the chart or a numerical pole was met."""


@dataclass(frozen=True)
class HeunParameters:
    """Signed endpoint parameters and the fixed remaining Heun data."""

    theta0: mp.mpf
    theta1: mp.mpf
    thetat: mp.mpf
    thetainf: mp.mpf
    omega: mp.mpf

    @property
    def x(self) -> mp.mpf:
        return mp.mpf("0.5") - self.theta0 + self.theta1


@dataclass(frozen=True)
class RichardsonResult:
    """Cutoffs, triangular table, extrapolate, and last diagonal change."""

    cutoffs: tuple[int, ...]
    table: tuple[tuple[mp.mpf | mp.mpc, ...], ...]

    @property
    def raw(self) -> mp.mpf | mp.mpc:
        return self.table[-1][0]

    @property
    def value(self) -> mp.mpf | mp.mpc:
        return self.table[-1][-1]

    @property
    def error_estimate(self) -> mp.mpf:
        if len(self.table) < 2:
            return mp.inf
        # Compare consecutive diagonal extrapolants R_{j,j} and
        # R_{j-1,j-1}.  ``table[-2][-2]`` would instead select the
        # penultimate entry of the previous row and overstate the change.
        return abs(self.table[-1][-1] - self.table[-2][-1])


def canonical_parameters(theta0_sign: int = 1, theta1_sign: int = 1) -> HeunParameters:
    """Return the book's exact-rational slice with chosen endpoint signs."""

    return HeunParameters(
        theta0=theta0_sign * mp.mpf(1) / 7,
        theta1=theta1_sign * mp.mpf(2) / 9,
        thetat=mp.mpf(1) / 5,
        thetainf=mp.mpf(3) / 11,
        omega=mp.mpf(2) / 5,
    )


def _is_near_integer(value: mp.mpf | mp.mpc, tolerance: mp.mpf) -> bool:
    if abs(mp.im(value)) > tolerance:
        return False
    real = mp.re(value)
    return abs(real - mp.nint(real)) <= tolerance


def _is_nonpositive_integer(value: mp.mpf | mp.mpc, tolerance: mp.mpf) -> bool:
    return _is_near_integer(value, tolerance) and mp.re(value) <= tolerance


def _validation_tolerance() -> mp.mpf:
    return mp.power(10, -max(18, mp.mp.dps // 2))


def validate_chart(params: HeunParameters) -> None:
    """Reject endpoint resonance and static gamma/recurrence poles."""

    tol = _validation_tolerance()
    for label, theta in (("theta_0", params.theta0), ("theta_1", params.theta1)):
        if _is_near_integer(2 * theta, tol):
            raise AuditError(
                f"endpoint resonance: {label}={mp.nstr(theta, 16)} lies in (1/2) Z"
            )

    gamma_arguments = {
        "1 - 2 theta_0": 1 - 2 * params.theta0,
        "2 theta_1": 2 * params.theta1,
        "x + omega": params.x + params.omega,
        "x - omega": params.x - params.omega,
    }
    for label, argument in gamma_arguments.items():
        if _is_nonpositive_integer(argument, tol):
            raise AuditError(
                "gamma pole in the hypergeometric prefactor: "
                f"{label}={mp.nstr(argument, 16)}"
            )

    # Every recurrence denominator is D_k=(k+x)^2-omega^2, k>=0.
    # It vanishes exactly when either -x+omega or -x-omega is a
    # nonnegative integer.
    for sign in (1, -1):
        root = -params.x + sign * params.omega
        if _is_near_integer(root, tol) and mp.re(root) >= -tol:
            raise AuditError(
                "recurrence pole: (k+x)^2-omega^2 vanishes at "
                f"k={mp.nstr(mp.re(root), 16)}"
            )


def validate_lambda(lam: mp.mpf | mp.mpc) -> None:
    if not (mp.isfinite(mp.re(lam)) and mp.isfinite(mp.im(lam))):
        raise AuditError("lambda must be finite")
    if abs(lam) >= 1:
        raise AuditError(
            f"the implemented weak-coupling chart requires |lambda| < 1; got {mp.nstr(abs(lam), 12)}"
        )
    if abs(1 - lam) <= _validation_tolerance():
        raise AuditError("lambda is too close to the singular point lambda=1")


def denominator(k: int, params: HeunParameters) -> mp.mpf:
    """D_k=(k+x)^2-omega^2 in the rescaled Heun recurrence."""

    return (k + params.x) ** 2 - params.omega**2


def alpha(k: int, params: HeunParameters) -> mp.mpf:
    """The diagonal recurrence coefficient alpha_k, k >= 0."""

    if k < 0:
        raise AuditError(f"alpha_k is only defined here for k >= 0, got {k}")
    den = denominator(k, params)
    if abs(den) <= _validation_tolerance():
        raise AuditError(f"recurrence denominator D_{k} is zero or numerically singular")
    numerator = (
        (k + mp.mpf("0.5") - params.theta0 - params.thetat) ** 2
        - params.theta0**2
        - params.thetainf**2
        + params.omega**2
    )
    return -numerator / den


def beta(k: int, params: HeunParameters) -> mp.mpf:
    """The subdiagonal recurrence coefficient beta_k, k >= 0."""

    if k < 0:
        raise AuditError(f"beta_k is only defined here for k >= 0, got {k}")
    if k == 0:
        # beta_0 multiplies a_{-1}=0 and is never needed analytically.
        return mp.mpf("0")
    den = denominator(k, params) * denominator(k - 1, params)
    if abs(den) <= _validation_tolerance():
        raise AuditError(
            f"recurrence denominator D_{k} D_{k - 1} is zero or numerically singular"
        )
    numerator = (
        k
        * (k - 2 * params.theta0)
        * ((k - params.theta0 + params.theta1 - params.thetat) ** 2 - params.thetainf**2)
    )
    return numerator / den


def hypergeometric_coefficient(k: int, params: HeunParameters) -> mp.mpf:
    """h_k=(x+omega)_k(x-omega)_k/[k!(1-2theta_0)_k]."""

    if k < 0:
        raise AuditError(f"hypergeometric coefficient requires k >= 0, got {k}")
    return (
        mp.rf(params.x + params.omega, k)
        * mp.rf(params.x - params.omega, k)
        / (mp.factorial(k) * mp.rf(1 - 2 * params.theta0, k))
    )


def hypergeometric_connection(params: HeunParameters) -> mp.mpf:
    """The lambda=0 Gauss connection coefficient C_HG."""

    validate_chart(params)
    return (
        mp.gamma(1 - 2 * params.theta0)
        * mp.gamma(2 * params.theta1)
        / (mp.gamma(params.x + params.omega) * mp.gamma(params.x - params.omega))
    )


def block_f1(params: HeunParameters) -> mp.mpf:
    """Analytic classical-block coefficient f_1 for the Heun slice."""

    validate_chart(params)
    quarter = mp.mpf(1) / 4
    gap = quarter - params.omega**2
    if abs(params.omega * gap) <= _validation_tolerance():
        raise AuditError(
            "the displayed accessory-inversion chart excludes "
            "omega=0 and omega^2=1/4"
        )
    a_factor = quarter - params.omega**2 + params.theta0**2 - params.theta1**2
    b_factor = quarter - params.omega**2 + params.thetainf**2 - params.thetat**2
    digamma_jump = mp.digamma(params.x + params.omega) - mp.digamma(
        params.x - params.omega
    )
    return (
        -a_factor * b_factor * digamma_jump / (4 * params.omega * gap)
        - (params.theta0 + params.theta1) * b_factor / (2 * gap)
    )


def f1_partial_sum(cutoff: int, params: HeunParameters) -> mp.mpf:
    """Return 1/2+theta_t-sum_{k=1}^K(alpha_{k-1}+beta_k)."""

    total = mp.fsum(alpha(k - 1, params) + beta(k, params) for k in range(1, cutoff + 1))
    return mp.mpf("0.5") + params.thetat - total


def continued_fraction_log_ratio(
    cutoff: int, lam: mp.mpf | mp.mpc, params: HeunParameters
) -> mp.mpf | mp.mpc:
    """Truncate the dual Riccati recurrence and return log(C/C_HG).

    The boundary eta_{K+1}=1 selects the root approached on |lambda|<1.
    The omitted sum over k>K is restored by cutoff Richardson
    extrapolation in :func:`connection_from_backward_recurrence`.
    """

    validate_lambda(lam)
    eta_next: mp.mpf | mp.mpc = mp.mpf("1")
    etas: list[mp.mpf | mp.mpc] = [mp.mpf("0")] * (cutoff + 1)
    dynamic_tol = _validation_tolerance()
    for k in range(cutoff, 0, -1):
        if abs(eta_next) <= dynamic_tol:
            raise AuditError(
                f"continued-fraction denominator eta_{k + 1} is zero or numerically singular"
            )
        eta_k = 1 - lam * alpha(k - 1, params) - lam * beta(k, params) / eta_next
        if abs(eta_k) <= dynamic_tol:
            raise AuditError(f"continued-fraction factor eta_{k} is zero or numerically singular")
        etas[k] = eta_k
        eta_next = eta_k

    return -(mp.mpf("0.5") + params.thetat) * mp.log(1 - lam) + mp.fsum(
        mp.log(etas[k]) for k in range(1, cutoff + 1)
    )


def forward_connection_partial(
    cutoff: int, lam: mp.mpf | mp.mpc, params: HeunParameters
) -> mp.mpf | mp.mpc:
    """Compute C_N from the forward recurrence with a_-1=0, a_0=1."""

    validate_lambda(lam)
    a_previous: mp.mpf | mp.mpc = mp.mpf("0")
    a_current: mp.mpf | mp.mpc = mp.mpf("1")
    for k in range(cutoff):
        a_next = a_current - lam * (
            alpha(k, params) * a_current + beta(k, params) * a_previous
        )
        a_previous, a_current = a_current, a_next
    return (
        hypergeometric_connection(params)
        * mp.power(1 - lam, mp.mpf("0.5") - params.thetat)
        * a_current
    )


def normalized_recurrence_coefficients(
    highest_index: int, lam: mp.mpf | mp.mpc, params: HeunParameters
) -> tuple[mp.mpf | mp.mpc, ...]:
    """Return a_0,...,a_highest_index from the finite-lambda recurrence."""

    if highest_index < 0:
        raise AuditError("highest recurrence index must be nonnegative")
    validate_lambda(lam)
    values: list[mp.mpf | mp.mpc] = [mp.mpf("1")]
    a_previous: mp.mpf | mp.mpc = mp.mpf("0")
    a_current: mp.mpf | mp.mpc = mp.mpf("1")
    for k in range(highest_index):
        a_next = a_current - lam * (
            alpha(k, params) * a_current + beta(k, params) * a_previous
        )
        values.append(a_next)
        a_previous, a_current = a_current, a_next
    return tuple(values)


def build_richardson(
    evaluator: Callable[[int], mp.mpf | mp.mpc], base_cutoff: int, levels: int
) -> RichardsonResult:
    """Build the ordinary Richardson table for an expansion in 1/K."""

    cutoffs = tuple(base_cutoff * 2**level for level in range(levels))
    table: list[tuple[mp.mpf | mp.mpc, ...]] = []
    for row_index, cutoff in enumerate(cutoffs):
        row: list[mp.mpf | mp.mpc] = [evaluator(cutoff)]
        for order in range(1, row_index + 1):
            scale = mp.mpf(2) ** order
            row.append((scale * row[order - 1] - table[row_index - 1][order - 1]) / (scale - 1))
        table.append(tuple(row))
    return RichardsonResult(cutoffs=cutoffs, table=tuple(table))


def connection_from_backward_recurrence(
    lam: mp.mpf | mp.mpc,
    params: HeunParameters,
    base_cutoff: int,
    levels: int,
) -> tuple[mp.mpf | mp.mpc, RichardsonResult]:
    """Return C and the Richardson table for the dual recurrence."""

    prefactor = hypergeometric_connection(params)
    log_table = build_richardson(
        lambda cutoff: continued_fraction_log_ratio(cutoff, lam, params),
        base_cutoff,
        levels,
    )
    return prefactor * mp.exp(log_table.value), log_table


def connection_from_forward_recurrence(
    lam: mp.mpf | mp.mpc,
    params: HeunParameters,
    base_cutoff: int,
    levels: int,
) -> RichardsonResult:
    """Return the Richardson table for the direct forward limit C_N."""

    return build_richardson(
        lambda cutoff: forward_connection_partial(cutoff, lam, params),
        base_cutoff,
        levels,
    )


def signed_parameter_grid() -> tuple[tuple[HeunParameters, HeunParameters], ...]:
    return (
        (canonical_parameters(+1, +1), canonical_parameters(+1, -1)),
        (canonical_parameters(-1, +1), canonical_parameters(-1, -1)),
    )


def connection_matrix(
    lam: mp.mpf | mp.mpc, base_cutoff: int, levels: int
) -> tuple[list[list[mp.mpf | mp.mpc]], list[list[RichardsonResult]]]:
    """Compute the source-first 2 x 2 connection array."""

    matrix: list[list[mp.mpf | mp.mpc]] = []
    tables: list[list[RichardsonResult]] = []
    for parameter_row in signed_parameter_grid():
        value_row: list[mp.mpf | mp.mpc] = []
        table_row: list[RichardsonResult] = []
        for params in parameter_row:
            value, table = connection_from_backward_recurrence(
                lam, params, base_cutoff, levels
            )
            value_row.append(value)
            table_row.append(table)
        matrix.append(value_row)
        tables.append(table_row)
    return matrix, tables


def determinant_2x2(matrix: Sequence[Sequence[mp.mpf | mp.mpc]]) -> mp.mpf | mp.mpc:
    return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]


def _parse_lambda(text: str) -> mp.mpf | mp.mpc:
    try:
        value = mp.mpmathify(text.strip().replace("i", "j"))
    except (TypeError, ValueError) as exc:
        raise AuditError(
            "lambda must be a real or complex numeric literal, for example 0.04 or 0.04+0.01j"
        ) from exc
    if not isinstance(value, (mp.mpf, mp.mpc)):
        raise AuditError("lambda must evaluate to one numeric value")
    return value


def _format(value: mp.mpf | mp.mpc, digits: int) -> str:
    tolerance = mp.power(10, -max(12, digits - 4))
    if abs(mp.im(value)) <= tolerance:
        return mp.nstr(mp.re(value), digits)
    return mp.nstr(value, digits)


def _print_richardson_table(
    title: str, result: RichardsonResult, digits: int, transform: Callable[[mp.mpf | mp.mpc], mp.mpf | mp.mpc] | None = None
) -> None:
    transform = transform or (lambda value: value)
    print(title)
    print("  cutoff             raw cutoff value             Richardson diagonal")
    for cutoff, row in zip(result.cutoffs, result.table):
        print(
            f"  {cutoff:6d}  {_format(transform(row[0]), digits):>31}"
            f"  {_format(transform(row[-1]), digits):>31}"
        )
    print(f"  last diagonal change = {_format(result.error_estimate, digits)}")


def _coefficient_diagonal_change(
    prefactor: mp.mpf | mp.mpc,
    logarithm_table: RichardsonResult,
) -> mp.mpf:
    """Return the last diagonal change after exponentiating to C units."""
    if len(logarithm_table.table) < 2:
        return mp.inf
    latest = logarithm_table.table[-1][-1]
    previous = logarithm_table.table[-2][-1]
    return abs(prefactor * (mp.exp(latest) - mp.exp(previous)))


def _parameter_line(params: HeunParameters) -> str:
    return (
        f"theta0={_format(params.theta0, 12)}, "
        f"theta1={_format(params.theta1, 12)}, "
        f"thetat={_format(params.thetat, 12)}, "
        f"thetainf={_format(params.thetainf, 12)}, "
        f"omega={_format(params.omega, 12)}"
    )


def make_parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description=(
            "Check a fixed general-Heun connection coefficient using the "
            "dual continued fraction, the forward recurrence, and the "
            "analytic first classical-block coefficient."
        )
    )
    parser.add_argument(
        "--lambda",
        dest="lambda_text",
        default="0.04",
        metavar="VALUE",
        help="weak-coupling value lambda=1/t, real or complex (default: 0.04)",
    )
    parser.add_argument(
        "--dps", type=int, default=60, help="mpmath decimal precision, at least 30 (default: 60)"
    )
    parser.add_argument(
        "--base-cutoff",
        type=int,
        default=128,
        help="first recurrence cutoff K, at least 16 (default: 128)",
    )
    parser.add_argument(
        "--levels",
        type=int,
        default=7,
        help="number of doubled cutoffs/Richardson rows, 2--12 (default: 7)",
    )
    parser.add_argument(
        "--show-matrix",
        action="store_true",
        help="print the full source-first 2 x 2 coefficient array",
    )
    parser.add_argument(
        "--show-scaling",
        action="store_true",
        help="compare C(lambda) with C_HG exp(f1 lambda) at 2lambda, lambda, lambda/2, lambda/4",
    )
    parser.add_argument(
        "--check-forward",
        action="store_true",
        help="independently extrapolate the direct forward recurrence",
    )
    return parser


def validate_numerical_options(parser: argparse.ArgumentParser, args: argparse.Namespace) -> None:
    if args.dps < 30:
        parser.error("--dps must be at least 30")
    if args.base_cutoff < 16:
        parser.error("--base-cutoff must be at least 16")
    if not 2 <= args.levels <= 12:
        parser.error("--levels must lie between 2 and 12")
    largest = args.base_cutoff * 2 ** (args.levels - 1)
    if largest > 2_000_000:
        parser.error(
            "largest cutoff exceeds 2,000,000; reduce --base-cutoff or --levels"
        )


def run(args: argparse.Namespace, parser: argparse.ArgumentParser) -> None:
    validate_numerical_options(parser, args)
    mp.mp.dps = args.dps
    lam = _parse_lambda(args.lambda_text)
    validate_lambda(lam)
    if args.show_scaling and 2 * abs(lam) >= 1:
        parser.error("--show-scaling requires |2 lambda| < 1")
    if args.show_scaling and abs(lam) == 0:
        parser.error("--show-scaling requires nonzero lambda")

    started = perf_counter()
    digits = min(40, max(18, args.dps - 12))
    table_digits = min(24, digits)
    params = canonical_parameters()
    for parameter_row in signed_parameter_grid():
        for signed_params in parameter_row:
            validate_chart(signed_params)

    print("General-Heun recurrence / classical-block audit")
    print(f"parameters: {_parameter_line(params)}")
    print(f"lambda={_format(lam, digits)}, |lambda|={_format(abs(lam), 12)}")
    print(
        f"precision={args.dps} dps, cutoffs={args.base_cutoff}*2^j, "
        f"j=0,...,{args.levels - 1}"
    )
    print()

    c_hg = hypergeometric_connection(params)
    f1_exact = block_f1(params)
    f1_table = build_richardson(
        lambda cutoff: f1_partial_sum(cutoff, params), args.base_cutoff, args.levels
    )
    print(f"C_HG = {_format(c_hg, digits)}")
    print(f"f1 (analytic block formula) = {_format(f1_exact, digits)}")
    print(f"f1 (ODE recurrence sum)     = {_format(f1_table.value, digits)}")
    print(f"absolute f1 discrepancy     = {_format(abs(f1_table.value - f1_exact), digits)}")
    print()
    _print_richardson_table("First-order infinite-sum audit:", f1_table, table_digits)
    print()

    first_a = normalized_recurrence_coefficients(2, lam, params)
    print("Finite-lambda recurrence slice (u_k = h_k a_k):")
    for k in range(3):
        h_k = hypergeometric_coefficient(k, params)
        print(
            f"  k={k}: h_k={_format(h_k, table_digits)}, "
            f"a_k={_format(first_a[k], table_digits)}, "
            f"u_k={_format(h_k * first_a[k], table_digits)}"
        )
    print()

    # The determinant audit needs all four endpoint-sign choices.  Reuse its
    # (+,+) entry here instead of evaluating the primary coefficient twice.
    matrix, matrix_tables = connection_matrix(lam, args.base_cutoff, args.levels)
    coefficient = matrix[0][0]
    log_table = matrix_tables[0][0]
    prediction = c_hg * mp.exp(f1_exact * lam)
    print(f"C(+theta0,+theta1; lambda) = {_format(coefficient, digits)}")
    print(f"C_HG exp(f1 lambda)         = {_format(prediction, digits)}")
    print(f"finite-lambda difference    = {_format(coefficient - prediction, digits)}")
    print()
    _print_richardson_table(
        "Backward continued-fraction audit (table stores log(C/C_HG)):",
        log_table,
        table_digits,
        transform=lambda value: c_hg * mp.exp(value),
    )
    print(
        "  transformed last diagonal change = "
        f"{_format(_coefficient_diagonal_change(c_hg, log_table), table_digits)}"
    )
    print()

    determinant = determinant_2x2(matrix)
    determinant_target = -params.theta0 / params.theta1
    print("Wronskian/determinant audit:")
    print("  convention: rows = source signs (+,-), columns = target signs (+,-)")
    if args.show_matrix:
        print(
            "  C = [["
            f"{_format(matrix[0][0], digits)}, {_format(matrix[0][1], digits)}],"
        )
        print(
            "       ["
            f"{_format(matrix[1][0], digits)}, {_format(matrix[1][1], digits)}]]"
        )
    print(f"  det C              = {_format(determinant, digits)}")
    print(f"  -theta0/theta1     = {_format(determinant_target, digits)}")
    print(f"  absolute residual  = {_format(abs(determinant - determinant_target), digits)}")
    print()

    if args.show_scaling:
        print("Weak-coupling scaling (full recurrence versus first-order block prediction):")
        print("  lambda                 full C(lambda)          C_HG exp(f1 lambda)       difference/lambda^2")
        for scale in (2, 1, mp.mpf("0.5"), mp.mpf("0.25")):
            scaled_lambda = scale * lam
            scaled_c, _ = connection_from_backward_recurrence(
                scaled_lambda, params, args.base_cutoff, args.levels
            )
            scaled_prediction = c_hg * mp.exp(f1_exact * scaled_lambda)
            scaled_remainder = (scaled_c - scaled_prediction) / scaled_lambda**2
            print(
                f"  {_format(scaled_lambda, 12):>12}  "
                f"{_format(scaled_c, 22):>25}  "
                f"{_format(scaled_prediction, 22):>25}  "
                f"{_format(scaled_remainder, 16):>22}"
            )
        print()

    if args.check_forward:
        print("Forward/backward recurrence agreement:")
        indices: Iterable[tuple[int, int]]
        indices = ((0, 0), (0, 1), (1, 0), (1, 1)) if args.show_matrix else ((0, 0),)
        signs = (("+", "+"), ("+", "-"), ("-", "+"), ("-", "-"))
        for row, column in indices:
            signed_params = signed_parameter_grid()[row][column]
            forward_table = connection_from_forward_recurrence(
                lam, signed_params, args.base_cutoff, args.levels
            )
            backward_value = matrix[row][column]
            label = signs[2 * row + column]
            print(
                f"  C({label[0]},{label[1]}): backward={_format(backward_value, table_digits)}, "
                f"forward={_format(forward_table.value, table_digits)}, "
                f"|difference|={_format(abs(forward_table.value - backward_value), table_digits)}"
            )
            print(
                f"           forward last diagonal change={_format(forward_table.error_estimate, table_digits)}, "
                "backward C-scale change="
                f"{_format(_coefficient_diagonal_change(hypergeometric_connection(signed_params), matrix_tables[row][column]), table_digits)}"
            )
        print()

    elapsed = perf_counter() - started
    print(
        f"environment: Python {platform.python_version()}, mpmath {mp.__version__}; "
        f"runtime={elapsed:.3f} s"
    )


def main() -> None:
    parser = make_parser()
    args = parser.parse_args()
    try:
        run(args, parser)
    except AuditError as exc:
        parser.error(str(exc))


if __name__ == "__main__":
    main()