distlink.py

"""Minimum orbit intersection distance and linking coefficients.

This is a source-level Python translation of ``distlink-cpp/distlink.cpp`` by
R. V. Baluev and D. V. Mikryukov (2018-2020).  The original and this port are
distributed under the MIT license reproduced below.

Copyright (c) 2018-2020 R.V. Baluev and D.V. Mikryukov
Copyright (c) 2026 Troy D. Rockwood (Python translation and modifications)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

The public names intentionally follow the C++ API.  Python ``float`` is used
throughout and therefore corresponds to the C++ ``double`` instantiation.
"""

from __future__ import annotations

from dataclasses import dataclass
import cmath
import math
import sys
import time
from typing import Callable, Sequence


DIM = 21
DEG = 16
_PI = math.pi
_CIRC = 2.0 * math.pi
_EPS = sys.float_info.epsilon
_EXPS = tuple(cmath.rect(1.0, 2.0 * j * math.pi / DIM) for j in range(DIM))
_UINT64_MASK = (1 << 64) - 1


class _DeterministicRandom:
    """SplitMix64 stream matching the C++ implementation exactly."""

    __slots__ = ("state",)

    def __init__(self) -> None:
        self.state = 0x243F6A8885A308D3

    def next(self) -> float:
        self.state = (self.state + 0x9E3779B97F4A7C15) & _UINT64_MASK
        value = self.state
        value = ((value ^ (value >> 30)) * 0xBF58476D1CE4E5B9) & _UINT64_MASK
        value = ((value ^ (value >> 27)) * 0x94D049BB133111EB) & _UINT64_MASK
        value ^= value >> 31
        return (value >> 11) * (2.0 ** -53)


def _sqr(x):
    return x * x


def _sign(x: float) -> int:
    return int(x > 0.0) - int(x < 0.0)


def _safe_asin(x: float) -> float:
    return math.asin(max(min(x, 1.0), -1.0))


def _safe_sqrt(x: float) -> float:
    return math.sqrt(max(x, 0.0))


def _safe_acosh(x: float) -> float:
    return math.acosh(max(x, 1.0))


def _sqrt_ieee(x: float) -> float:
    """Return sqrt with the NaN-on-negative behavior of C++ libm."""

    return math.sqrt(x) if x >= 0.0 else math.nan


def _float_divide(numerator: float, denominator: float) -> float:
    """Emulate ordinary IEEE-754 division where Python would raise on zero."""

    if denominator != 0.0:
        return numerator / denominator
    if numerator == 0.0 or math.isnan(numerator):
        return math.nan
    sign = math.copysign(1.0, numerator) * math.copysign(1.0, denominator)
    return math.copysign(math.inf, sign)


def _log_abs(x: float) -> float:
    return -math.inf if x == 0.0 else math.log(abs(x))


def _atan2h(y: float, x: float) -> float:
    return (_log_abs(x + y) - _log_abs(x - y)) / 2.0


def _atan2_smart(y: float, x: float, trig: bool) -> float:
    return math.atan2(y, x) if trig else _atan2h(y, x)


def _dot(a: Sequence[float], b: Sequence[float]) -> float:
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def _cross(a: Sequence[float], b: Sequence[float]) -> tuple[float, float, float]:
    return (
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    )


def _vnorm(a: Sequence[float]) -> float:
    return _dot(a, a)


def _cnorm(z: complex) -> float:
    return z.real * z.real + z.imag * z.imag


def _inverse(z: complex) -> complex:
    return z.conjugate() / _cnorm(z)


def _angle_wrap(x: float) -> float:
    x = math.fmod(x + _PI, _CIRC)
    return x + _PI if x < 0.0 else x - _PI


def detect_suitable_options() -> tuple[float, float, float]:
    """Return ``(max_root_error, min_root_error, max_anom_error)``.

    This is the Python-float equivalent of the C++ template function.
    """

    return math.sqrt(_EPS), 2.0 * _EPS, 1000.0 * _EPS


class COrbitData:
    """Keplerian elements and the associated orthonormal P/Q vectors."""

    __slots__ = ("_a", "_e", "_i", "_w", "_Om", "_P", "_Q")

    def __init__(
        self,
        a: float = 1.0,
        e: float = 0.0,
        i: float = 0.0,
        w: float = 0.0,
        Om: float = 0.0,
    ) -> None:
        self.set_data(a, e, i, w, Om)

    def set_data(self, a: float, e: float, i: float, w: float, Om: float) -> None:
        a, e, i, w, Om = (float(value) for value in (a, e, i, w, Om))
        if not all(math.isfinite(value) for value in (a, e, i, w, Om)):
            raise ValueError("orbital elements must be finite")
        if a == 0.0:
            raise ValueError("semimajor axis must be nonzero")
        if e < 0.0:
            raise ValueError("eccentricity must be nonnegative")
        if e == 1.0:
            raise ValueError("parabolic orbits (e == 1) are unsupported")
        self._a = abs(a) if e < 1.0 else -abs(a)
        self._e, self._i, self._w, self._Om = e, i, w, Om
        self._set_vectors()

    def _set_vectors(self) -> None:
        cw, sw = math.cos(self._w), math.sin(self._w)
        cO, sO = math.cos(self._Om), math.sin(self._Om)
        ci, si = math.cos(self._i), math.sin(self._i)
        self._P = (
            cw * cO - ci * sw * sO,
            cw * sO + ci * sw * cO,
            si * sw,
        )
        self._Q = (
            -sw * cO - ci * cw * sO,
            ci * cw * cO - sw * sO,
            si * cw,
        )

    def get_data(self) -> tuple[float, float, float, float, float]:
        return self._a, self._e, self._i, self._w, self._Om

    def get_a(self) -> float:
        return self._a

    def get_e(self) -> float:
        return self._e

    def get_i(self) -> float:
        return self._i

    def get_w(self) -> float:
        return self._w

    def get_Om(self) -> float:
        return self._Om

    def vectorP(self) -> tuple[float, float, float]:
        return self._P

    def vectorQ(self) -> tuple[float, float, float]:
        return self._Q

    def get_vectors(self) -> tuple[tuple[float, float, float], tuple[float, float, float]]:
        return self._P, self._Q

    def __repr__(self) -> str:
        values = ", ".join(repr(v) for v in self.get_data())
        return f"COrbitData({values})"


def test_peri_apo(O1: COrbitData, O2: COrbitData, limit: float) -> bool:
    a1, e1 = O1.get_a(), O1.get_e()
    a2, e2 = O2.get_a(), O2.get_e()
    rp1, ra1 = a1 * (1.0 - e1), a1 * (1.0 + e1)
    rp2, ra2 = a2 * (1.0 - e2), a2 * (1.0 + e2)
    return rp1 - ra2 < limit and rp2 - ra1 < limit


@dataclass(slots=True)
class SMOIDResult:
    good: bool = True
    distance: float = -1.0
    distance_error: float = 0.0
    u1: float = math.nan
    u1_error: float = 0.0
    u2: float = math.nan
    u2_error: float = 0.0
    root_count: int = 0
    min_delta: float = -1.0
    iter_count: int = 0
    iter_count_2D: int = 0
    time: float = 0.0


@dataclass(slots=True)
class SLCResult:
    I: float = -1.0
    l: float = 0.0
    lmod: float = 0.0
    l2: float = 0.0


class _AuxData:
    __slots__ = (
        "e1", "e2", "a1", "a2", "alpha1", "alpha2", "K",
        "Pp", "Ps", "Sp", "Ss", "p1", "p2", "w1", "w2", "I",
        "abs_w", "P1w", "P2w", "Q1w", "Q2w", "P1", "P2", "Q1", "Q2",
    )

    def __init__(self, O1: COrbitData, O2: COrbitData) -> None:
        self.e1, self.e2 = O1.get_e(), O2.get_e()
        self.a1, self.a2 = O1.get_a(), O2.get_a()
        self.w1, self.w2 = O1.get_w(), O2.get_w()
        self.p1 = self.a1 * (1.0 - self.e1 * self.e1)
        self.p2 = self.a2 * (1.0 - self.e2 * self.e2)
        self.alpha1, self.alpha2 = self.a1 / self.a2, self.a2 / self.a1
        self.K = self.alpha2 * self.e2 * self.e2
        eta1 = math.sqrt(abs(1.0 - self.e1 * self.e1))
        eta2 = math.sqrt(abs(1.0 - self.e2 * self.e2))
        i1, i2 = O1.get_i(), O2.get_i()
        Om1, Om2 = O1.get_Om(), O2.get_Om()
        c1, s1 = math.cos(i1), math.sin(i1)
        c2, s2 = math.cos(i2), math.sin(i2)
        wv = (
            c1 * s2 * math.cos(Om2) - s1 * c2 * math.cos(Om1),
            c1 * s2 * math.sin(Om2) - s1 * c2 * math.sin(Om1),
            s1 * s2 * math.sin(Om2 - Om1),
        )
        self.abs_w = math.sqrt(_vnorm(wv))
        cosI = c1 * c2 + s1 * s2 * math.cos(Om2 - Om1)
        self.I = _safe_asin(self.abs_w) if cosI > 0.0 else _PI - _safe_asin(self.abs_w)

        self.P1, self.P2 = O1.vectorP(), O2.vectorP()
        self.Q1, self.Q2 = O1.vectorQ(), O2.vectorQ()
        self.Pp = _dot(self.P1, self.P2)
        self.Ps = _dot(self.P1, self.Q2) * eta2
        self.Sp = _dot(self.Q1, self.P2) * eta1
        self.Ss = _dot(self.Q1, self.Q2) * eta1 * eta2
        self.P1w, self.P2w = _dot(self.P1, wv), _dot(self.P2, wv)
        self.Q1w, self.Q2w = _dot(self.Q1, wv), _dot(self.Q2, wv)


def _l2(data: _AuxData) -> float:
    f1 = data.p2 * (data.abs_w + data.e1 * data.P1w) - data.p1 * (
        data.abs_w + data.e2 * data.P2w
    )
    f2 = data.p2 * (data.abs_w - data.e1 * data.P1w) - data.p1 * (
        data.abs_w - data.e2 * data.P2w
    )
    return f1 * f2


def _l1(data: _AuxData) -> tuple[float, float]:
    costheta1 = _float_divide(data.P1w, data.abs_w)
    costheta2 = _float_divide(data.P2w, data.abs_w)
    rr = data.p1 / (1.0 + data.e1 * costheta1) - data.p2 / (
        1.0 + data.e2 * costheta2
    )
    RR = data.p1 / (1.0 - data.e1 * costheta1) - data.p2 / (
        1.0 - data.e2 * costheta2
    )
    product = rr * RR
    smaller = RR if abs(RR) < abs(rr) else rr
    return product, smaller * smaller * _sign(product)


def _l3(data: _AuxData) -> float:
    return (
        data.a1 * data.p1
        + data.a2 * data.p2
        - 2.0
        * data.a1
        * data.a2
        * (1.0 - data.e1 * data.e2 * math.cos(data.w1 - data.w2))
    )


def LC(O1: COrbitData, O2: COrbitData, min_mut_incl: float) -> SLCResult:
    if not math.isfinite(min_mut_incl) or min_mut_incl < 0.0:
        raise ValueError("minimum mutual inclination must be finite and nonnegative")
    data = _AuxData(O1, O2)
    result = SLCResult(I=data.I, l2=_l2(data))
    if data.I < min_mut_incl or data.abs_w <= math.sqrt(_EPS):
        result.l = _l3(data)
    else:
        result.l, result.lmod = _l1(data)
    return result


def _radius_vector(
    branch: bool,
    P: Sequence[float],
    Q: Sequence[float],
    a: float,
    e: float,
    u: float,
) -> tuple[float, float, float]:
    if e <= 1.0:
        p, q = math.cos(u), math.sin(u) * math.sqrt(1.0 - e * e)
    else:
        tmp = math.exp(abs(u))
        p = (tmp + 1.0 / tmp) / 2.0
        q = (1.0 / tmp - tmp) / 2.0 * _sign(u) * math.sqrt(e * e - 1.0)
    if not branch:
        p, q = -p, -q
    p -= e
    return tuple(a * (P[k] * p + Q[k] * q) for k in range(3))  # type: ignore[return-value]


def _radius_vector_valder(
    branch: bool,
    P: Sequence[float],
    Q: Sequence[float],
    a: float,
    e: float,
    u: float,
) -> tuple[tuple[float, float, float], tuple[float, float, float]]:
    if e <= 1.0:
        x, y = math.cos(u), math.sin(u)
        eta = math.sqrt(1.0 - e * e)
        h = x * eta
    else:
        tmp = math.exp(abs(u))
        x, y = (tmp + 1.0 / tmp) / 2.0, (1.0 / tmp - tmp) / 2.0 * _sign(u)
        eta = math.sqrt(e * e - 1.0)
        h = -x * eta
    if not branch:
        x, y, h = -x, -y, -h
    p, q = x - e, y * eta
    r = tuple(a * (P[k] * p + Q[k] * q) for k in range(3))
    rd = tuple(-P[k] * y + Q[k] * h for k in range(3))
    return r, rd  # type: ignore[return-value]


def _radius_vector_valder2(
    branch: bool,
    P: Sequence[float],
    Q: Sequence[float],
    a: float,
    e: float,
    u: float,
) -> tuple[
    tuple[float, float, float],
    tuple[float, float, float],
    tuple[float, float, float],
]:
    if e <= 1.0:
        x, y = math.cos(u), math.sin(u)
        eta = math.sqrt(1.0 - e * e)
        h = x * eta
    else:
        tmp = math.exp(abs(u))
        x, y = (tmp + 1.0 / tmp) / 2.0, (1.0 / tmp - tmp) / 2.0 * _sign(u)
        eta = math.sqrt(e * e - 1.0)
        h = -x * eta
    if not branch:
        x, y, h = -x, -y, -h
    p, q = x - e, y * eta
    r = tuple(a * (P[k] * p + Q[k] * q) for k in range(3))
    rd = tuple(-P[k] * y + Q[k] * h for k in range(3))
    rdd = tuple(P[k] * x + Q[k] * q for k in range(3))
    if e <= 1.0:
        rdd = tuple(-v for v in rdd)
    return r, rd, rdd  # type: ignore[return-value]


def _distance_between(data: _AuxData, br1: bool, br2: bool, u1: float, u2: float) -> float:
    r1 = _radius_vector(br1, data.P1, data.Q1, data.a1, data.e1, u1)
    r2 = _radius_vector(br2, data.P2, data.Q2, data.a2, data.e2, u2)
    dr = tuple(r2[k] - r1[k] for k in range(3))
    return math.sqrt(_vnorm(dr))


def _sqdist_vg(
    data: _AuxData, br1: bool, br2: bool, u1: float, u2: float
) -> tuple[float, list[float]]:
    r1, rd1 = _radius_vector_valder(br1, data.P1, data.Q1, data.a1, data.e1, u1)
    r2, rd2 = _radius_vector_valder(br2, data.P2, data.Q2, data.a2, data.e2, u2)
    dr = tuple(r2[k] - r1[k] for k in range(3))
    g = [-_dot(dr, rd1) / data.a2, _dot(dr, rd2) / data.a1]
    return _vnorm(dr) / (2.0 * data.a1 * data.a2), g


def _sqdist_vgH(
    data: _AuxData, br1: bool, br2: bool, u1: float, u2: float
) -> tuple[float, list[float], list[float]]:
    r1, rd1, rdd1 = _radius_vector_valder2(
        br1, data.P1, data.Q1, data.a1, data.e1, u1
    )
    r2, rd2, rdd2 = _radius_vector_valder2(
        br2, data.P2, data.Q2, data.a2, data.e2, u2
    )
    dr = tuple(r2[k] - r1[k] for k in range(3))
    g = [-_dot(dr, rd1) / data.a2, _dot(dr, rd2) / data.a1]
    H = [
        (-_dot(dr, rdd1) + _vnorm(rd1) * data.a1) / data.a2,
        (_dot(dr, rdd2) + _vnorm(rd2) * data.a2) / data.a1,
        -_dot(rd1, rd2),
    ]
    return _vnorm(dr) / (2.0 * data.a1 * data.a2), g, H


def _eliminated_anomaly(
    data: _AuxData, u1: float, coefficient_error: float | None = None
) -> tuple[bool, float, float, float | None]:
    if data.e1 <= 1.0:
        x, y = math.cos(u1), math.sin(u1)
        x_ = x - data.e1
        M = data.Sp * y + data.Pp * x_ + data.alpha2 * data.e2
        N = -data.Ss * y - data.Ps * x_
    else:
        tmp = math.exp(abs(u1))
        x, y = (tmp + 1.0 / tmp) / 2.0, (tmp - 1.0 / tmp) / 2.0 * _sign(u1)
        x_ = x - data.e1
        M = -data.Sp * y + data.Pp * x_ + data.alpha2 * data.e2
        N = data.Ss * y - data.Ps * x_
    A = data.Ps * y - data.Ss * x
    B = data.Pp * y - data.Sp * x
    C = data.e2 * B - data.alpha1 * data.e1 * y * (1.0 - data.e1 * x)
    if data.e2 <= 1.0:
        V, D = A * A + B * B, A * A + B * B - C * C
    else:
        V, D = B * B - A * A, C * C - (B * B - A * A)

    negative_discriminant = D < 0.0
    if negative_discriminant:
        D = 0.0
        if data.e2 <= 1.0:
            u2 = math.atan2(A, B) if C > 0.0 else math.atan2(-A, -B)
        else:
            u2 = _atan2h(A, B)
        u2_other = u2
    else:
        D = math.sqrt(D)
        AC, BD = A * C, B * D
        y1, y2 = (AC - BD) / V, (AC + BD) / V
        if not data.e2 <= 1.0:
            y1, y2 = y2, y1
        BC, AD = B * C, A * D
        x1, x2 = (BC + AD) / V, (BC - AD) / V
        if data.e2 <= 1.0:
            u2, u2_other = math.atan2(y1, x1), math.atan2(y2, x2)
        else:
            if not x1 > 0.0 and x2 > 0.0:
                x1, x2, y1, y2 = x2, x1, y2, y1
            u2, u2_other = _atan2h(y1, x1), _atan2h(y2, x2)
        err = M * y1 + N * x1 - data.K * y1 * x1
        err_other = M * y2 + N * x2 - data.K * y2 * x2
        if not abs(err) < abs(err_other) and (data.e2 <= 1.0 or x2 > 0.0):
            u2, u2_other = u2_other, u2

    if coefficient_error is not None:
        sigma = coefficient_error * math.sqrt(
            1.0 + data.e2 * data.e2 + _sqr(data.alpha1 * data.e1)
        )
        coefficient_error = sigma / math.sqrt(D * D + abs(C) * sigma / 2.0)
    return negative_discriminant, u2, u2_other, coefficient_error


def _min_distance_for(
    data: _AuxData, u1: float, coefficient_error: float | None = None
) -> tuple[float, float, float | None]:
    failed, u2, u2_other, coefficient_error = _eliminated_anomaly(
        data, u1, coefficient_error
    )
    if failed:
        return -_distance_between(data, True, True, u1, u2), u2, coefficient_error
    d = _distance_between(data, True, True, u1, u2)
    d_other = _distance_between(data, True, True, u1, u2_other)
    return (d, u2, coefficient_error) if d < d_other else (
        d_other,
        u2_other,
        coefficient_error,
    )


def _sqdist_2Diter(
    data: _AuxData, br1: bool, br2: bool, u1: float, u2: float
) -> tuple[float, float, float, list[float], list[float], float]:
    rho, g, H = _sqdist_vgH(data, br1, br2, u1, u2)
    detH = H[0] * H[1] - H[2] * H[2]
    if detH == 0.0:
        return u1, u2, 0.0, g, H, detH
    du1 = -(H[1] * g[0] - H[2] * g[1]) / detH
    du2 = -(H[0] * g[1] - H[2] * g[0]) / detH
    return u1 + du1, u2 + du2, du1 * du1 + du2 * du2, g, H, detH


def _newton_sqdist(
    data: _AuxData, u1: float, u2: float, eps: float, maxcount: int
) -> tuple[int, float, float, float, float, float, int, list[float], list[float]]:
    err2 = math.inf
    iter_count = 0
    detH = math.nan
    H = [math.nan, math.nan, math.nan]
    g = [math.nan, math.nan]
    while True:
        u1, u2, next_err2, g, H, detH = _sqdist_2Diter(data, True, True, u1, u2)
        if data.e1 <= 1.0 and abs(u1) > _PI:
            u1 = _angle_wrap(u1)
        if data.e2 <= 1.0 and abs(u2) > _PI:
            u2 = _angle_wrap(u2)
        if detH == 0.0:
            break
        derr2 = err2 - next_err2
        err2 = next_err2
        iter_count += 1
        if not (derr2 > _EPS * err2 and err2 > eps * eps and iter_count <= maxcount):
            break

    rho, g = _sqdist_vg(data, True, True, u1, u2)
    if detH == 0.0:
        u1_error = u2_error = rho_error = 0.0
        if H[0] > 0.0 or H[1] > 0.0:
            sign_H = 1
        elif H[0] < 0.0 or H[1] < 0.0:
            sign_H = -1
        else:
            sign_H = 0
    else:
        u1_error = abs((H[1] * g[0] - H[2] * g[1]) / detH)
        u2_error = abs((H[0] * g[1] - H[2] * g[0]) / detH)
        rho_error = abs(
            (H[1] * g[0] * g[0] + H[0] * g[1] * g[1] - 2.0 * H[2] * g[0] * g[1])
            / (2.0 * detH)
        )
        if detH > 0.0 and H[0] > 0.0:
            sign_H = 2
        elif detH > 0.0 and H[0] < 0.0:
            sign_H = -2
        else:
            sign_H = 0
    return sign_H, u1, u2, u1_error, u2_error, rho, rho_error, iter_count, g, H


def _restrict_search_range_aux(data: _AuxData) -> list[tuple[float, float]]:
    if data.abs_w <= math.sqrt(_EPS):
        if data.e1 < 1.0:
            return [(-_PI, _PI)]
        raise ValueError("cannot bound a coplanar hyperbolic anomaly range")
    theta1, theta2 = math.atan2(data.Q1w, data.P1w), math.atan2(data.Q2w, data.P2w)
    costheta1, costheta2 = math.cos(theta1), math.cos(theta2)
    rp = data.p1 / (1.0 + data.e1 * costheta1)
    rm = data.p1 / (1.0 - data.e1 * costheta1)
    Rp = data.p2 / (1.0 + data.e2 * costheta2)
    Rm = data.p2 / (1.0 - data.e2 * costheta2)
    distances = (abs(rp - Rp), abs(rp + Rm), abs(rm + Rp), abs(rm - Rm))
    flags = (rp > 0.0 and Rp > 0.0, rp > 0.0 and Rm > 0.0,
             rm > 0.0 and Rp > 0.0, rm > 0.0 and Rm > 0.0)
    dO = -1.0
    for valid, distance in zip(flags, distances):
        if valid and (dO > distance or dO < 0.0):
            dO = distance
    A = math.sqrt(abs(1.0 - _sqr(data.e1 * costheta1)))
    k = _float_divide(dO, A * data.abs_w * abs(data.a1))
    sintheta1 = math.sin(theta1)
    eta1 = math.sqrt(abs(1.0 - data.e1 * data.e1))
    phi = _atan2_smart(sintheta1, costheta1 * eta1, data.e1 <= 1.0)
    esth = _float_divide(data.e1 * sintheta1, A)
    if data.e1 <= 1.0:
        if min(abs(esth), 1.0) <= abs(1.0 - k):
            if k < 1.0:
                tmpm, tmpp = _safe_asin(esth - k), _safe_asin(esth + k)
                first = (phi - tmpp, phi - tmpm)
                phi += _PI
                return [first, (phi + tmpm, phi + tmpp)]
            return [(-_PI, _PI)]
        if esth > 0.0:
            tmpm = _safe_asin(esth - k)
            first = [phi + tmpm, phi + _PI - tmpm]
            if first[0] > first[1]:
                first = [phi + _PI / 2.0] * 2
        else:
            tmpp = _safe_asin(esth + k)
            first = [phi - tmpp, phi + _PI + tmpp]
            if first[0] > first[1]:
                first = [phi - _PI / 2.0] * 2
        return [(first[0], first[1])]
    if abs(costheta1) * eta1 <= abs(sintheta1):
        if sintheta1 < 0.0:
            esth = -esth
        tmpm, tmpp = _safe_acosh(esth - k), _safe_acosh(esth + k)
        if tmpm > 0.0:
            return [(-tmpp - phi, -tmpm - phi), (tmpm - phi, tmpp - phi)]
        return [(-tmpp - phi, tmpp - phi)]
    if costheta1 < 0.0:
        esth = -esth
    return [(math.asinh(esth - k) - phi, math.asinh(esth + k) - phi)]


def restrict_search_range(
    O1: COrbitData, O2: COrbitData
) -> tuple[tuple[float, float], ...]:
    """Return the one or two ranges selected by the C++ range restriction."""

    return tuple(_restrict_search_range_aux(_AuxData(O1, O2)))


def _search_in_segment(
    data: _AuxData, a: float, b: float, count: int
) -> tuple[int, float]:
    min_value = -1.0
    min_index = 0
    for index in range(count + 1):
        u1 = (a * (count - index) + b * index) / count
        distance, _, _ = _min_distance_for(data, u1)
        if distance < 0.0:
            continue
        if distance < min_value or min_value < 0.0:
            min_value, min_index = distance, index
    return min_index, min_value


def MOID_direct_search(
    O1: COrbitData,
    O2: COrbitData,
    densities: Sequence[int] | None,
    max_dist_error: float,
    max_anom_error: float,
) -> SMOIDResult:
    """Calculate MOID by the C++ implementation's scan-and-refine method."""

    if not math.isfinite(max_dist_error) or max_dist_error <= 0.0:
        raise ValueError("maximum distance error must be finite and positive")
    if not math.isfinite(max_anom_error) or max_anom_error <= 0.0:
        raise ValueError("maximum anomaly error must be finite and positive")
    started = time.perf_counter()
    results = [SMOIDResult(), SMOIDResult()]
    data = _AuxData(O1, O2)
    other_data = _AuxData(O2, O1)
    if data.abs_w <= math.sqrt(_EPS):
        if data.e1 > 1.0 and data.e2 > 1.0:
            raise ValueError("direct search cannot bound a coplanar hyperbolic pair")
        swapped = data.e1 > 1.0
        if swapped:
            data = other_data
        ranges = _restrict_search_range_aux(data)
    else:
        ranges = _restrict_search_range_aux(data)
        other_ranges = _restrict_search_range_aux(other_data)
        length = sum(max(b - a, 0.0) for a, b in ranges)
        other_length = sum(max(b - a, 0.0) for a, b in other_ranges)
        swapped = other_length < length
        if swapped:
            data, ranges = other_data, other_ranges

    for range_index, (a0, b0) in enumerate(ranges):
        h = (b0 - a0) / 2.0
        a, b = a0, b0
        if h > 0.0:
            min_value = math.nan
            delta = math.inf
            using_densities = densities is not None
            density = 3
            iteration = 0
            while iteration < 1 or h > max_anom_error or abs(delta) > max_dist_error:
                if using_densities:
                    if iteration < len(densities) and densities[iteration] != 0:
                        density = (
                            int(densities[iteration])
                            if iteration > 0
                            else math.ceil(h * densities[iteration] / _PI)
                        )
                    else:
                        using_densities = False
                if density < 3:
                    density = 3
                h = 2.0 * h / density
                min_index, next_min = _search_in_segment(data, a, b, density)
                a += h * min_index
                b, a = a + h, a - h
                if iteration > 0:
                    delta = next_min - min_value
                min_value = next_min
                iteration += 1
            results[range_index].u1_error = h
        else:
            results[range_index].u1_error = 0.0
        results[range_index].u1 = (a + b) / 2.0
        distance, results[range_index].u2, _ = _min_distance_for(
            data, results[range_index].u1
        )
        results[range_index].distance = abs(distance)

    result = results[0] if len(ranges) < 2 or results[0].distance <= results[1].distance else results[1]
    eps = math.sqrt(_EPS)
    distance, u2, _ = _min_distance_for(data, result.u1 + eps)
    distance = abs(distance)
    result.distance_error = abs(distance - result.distance) * _sqr(result.u1_error / eps)
    result.u2_error = abs(u2 - result.u2) * result.u1_error / eps
    if swapped:
        result.u1, result.u2 = result.u2, result.u1
        result.u1_error, result.u2_error = result.u2_error, result.u1_error
    result.good = result.good and result.distance_error < max_dist_error
    result.time = time.perf_counter() - started
    return result


def _polynomial_val(
    degree: int, coefficients: Sequence[complex], z: complex, forward: bool = True
) -> complex:
    if forward:
        value = coefficients[degree]
        for index in range(degree - 1, -1, -1):
            value = value * z + coefficients[index]
    else:
        value = coefficients[0]
        for index in range(1, degree + 1):
            value = value * z + coefficients[index]
    return value


def _polynomial_valder(
    degree: int, coefficients: Sequence[complex], z: complex, forward: bool
) -> tuple[complex, complex]:
    if forward:
        value = coefficients[degree]
        derivative = coefficients[degree] * degree
        for index in range(degree - 1, 1, -1):
            value = value * z + coefficients[index]
            derivative = derivative * z + coefficients[index] * index
        if degree > 1:
            value = value * z + coefficients[1]
            derivative = derivative * z + coefficients[1]
        if degree > 0:
            value = value * z + coefficients[0]
    else:
        value = coefficients[0]
        derivative = coefficients[0] * degree
        for index in range(1, degree - 1):
            value = value * z + coefficients[index]
            derivative = derivative * z + coefficients[index] * (degree - index)
        if degree > 1:
            value = value * z + coefficients[degree - 1]
            derivative = derivative * z + coefficients[degree - 1]
        if degree > 0:
            value = value * z + coefficients[degree]
    return value, derivative


def _polynomial_valder2(
    degree: int, coefficients: Sequence[complex], z: complex, forward: bool
) -> tuple[complex, complex, complex]:
    if forward:
        value = coefficients[degree]
        derivative = coefficients[degree] * degree
        derivative2 = coefficients[degree] * (degree * (degree - 1))
        for index in range(degree - 1, 2, -1):
            value = value * z + coefficients[index]
            derivative = derivative * z + coefficients[index] * index
            derivative2 = derivative2 * z + coefficients[index] * (index * (index - 1))
        if degree > 2:
            value = value * z + coefficients[2]
            derivative = derivative * z + coefficients[2] * 2
            derivative2 = derivative2 * z + coefficients[2] * 2
        if degree > 1:
            value = value * z + coefficients[1]
            derivative = derivative * z + coefficients[1]
        if degree > 0:
            value = value * z + coefficients[0]
    else:
        value = coefficients[0]
        derivative = coefficients[0] * degree
        derivative2 = coefficients[0] * (degree * (degree - 1))
        for index in range(1, degree - 2):
            value = value * z + coefficients[index]
            derivative = derivative * z + coefficients[index] * (degree - index)
            derivative2 = derivative2 * z + coefficients[index] * (
                (degree - index) * (degree - index - 1)
            )
        if degree > 2:
            value = value * z + coefficients[degree - 2]
            derivative = derivative * z + coefficients[degree - 2] * 2
            derivative2 = derivative2 * z + coefficients[degree - 2] * 2
        if degree > 1:
            value = value * z + coefficients[degree - 1]
            derivative = derivative * z + coefficients[degree - 1]
        if degree > 0:
            value = value * z + coefficients[degree]
    return value, derivative, derivative2


def _ratio1(degree: int, coefficients: Sequence[complex], z: complex, forward: bool) -> complex:
    value, derivative = _polynomial_valder(degree, coefficients, z, forward)
    return value / derivative


def _ratio2(coefficients: Sequence[complex], z: complex, forward: bool) -> complex:
    if forward:
        return ((coefficients[2] * z + coefficients[1]) * z + coefficients[0]) / (
            2.0 * coefficients[2] * z + coefficients[1]
        )
    return ((coefficients[0] * z + coefficients[1]) * z + coefficients[2]) / (
        2.0 * coefficients[0] * z + coefficients[1]
    )


def _ratio4(coefficients: Sequence[complex], z: complex, forward: bool) -> complex:
    if forward:
        return (
            (((coefficients[4] * z + coefficients[3]) * z + coefficients[2]) * z + coefficients[1])
            * z
            + coefficients[0]
        ) / (
            ((4.0 * coefficients[4] * z + 3.0 * coefficients[3]) * z + 2.0 * coefficients[2])
            * z
            + coefficients[1]
        )
    return (
        (((coefficients[0] * z + coefficients[1]) * z + coefficients[2]) * z + coefficients[3])
        * z
        + coefficients[4]
    ) / (
        ((4.0 * coefficients[0] * z + 3.0 * coefficients[1]) * z + 2.0 * coefficients[2])
        * z
        + coefficients[3]
    )


def _newton_polynomial(
    degree: int,
    coefficients: Sequence[complex],
    initial: complex,
    max_error: float,
    min_error: float,
    max_count: int,
    restart_count: int,
    random: _DeterministicRandom,
) -> tuple[complex, int]:
    z = initial
    max_error, min_error = abs(max_error), abs(min_error)
    z1 = z2 = z
    if not max_error >= 2.0 * min_error:
        max_error = 2.0 * min_error
    if not max_error >= 2.0 * _EPS:
        max_error = 2.0 * _EPS
    iteration = 0
    error2 = 0.0
    nrm = 0.0
    inner = True
    while True:
        iteration += 1
        cycled = (
            iteration > 2 and _cnorm(z - z2) < _cnorm(z - z1) * 1.0e-6
        ) or iteration % max_count == 0
        multiplier = (
            cmath.rect((2.0 * random.next() + 1.0) / 3.0, random.next() * _CIRC)
            if cycled
            else 1.0 + 0.0j
        )
        nrm = _cnorm(z)
        if inner:
            if nrm > 5.0:
                inner = False
        elif 2.0 * nrm < 1.0:
            inner = True
        z2, z1 = z1, z
        if inner:
            dz = _ratio1(degree, coefficients, z, True)
            z -= dz * multiplier
        else:
            nrm = 1.0 / nrm
            dz = _ratio1(degree, coefficients, z.conjugate() * nrm, False)
            z /= 1.0 - z * dz * multiplier
        error2 = _cnorm(dz)
        if not (error2 > max_error * max_error * nrm and iteration < max_count * restart_count):
            break

    if error2 > min_error * min_error * nrm:
        extra_iterations = 0
        while True:
            nrm = _cnorm(z)
            if inner:
                dz = _ratio1(degree, coefficients, z, True)
                z -= dz
            else:
                nrm = 1.0 / nrm
                dz = _ratio1(degree, coefficients, z.conjugate() * nrm, False)
                z /= 1.0 - z * dz
            next_error2 = _cnorm(dz)
            improvement = error2 - next_error2
            error2 = next_error2
            extra_iterations += 1
            if not (
                improvement > _EPS * error2
                and error2 > min_error * min_error * nrm
                and extra_iterations <= max_count
            ):
                break
        iteration += extra_iterations
    return z, iteration


def _extract_linear_factor(
    degree: int, coefficients: list[complex], root: complex
) -> list[complex]:
    if _cnorm(root) <= 1.0:
        for index in range(degree - 1, 0, -1):
            coefficients[index] += coefficients[index + 1] * root
        return coefficients[1:]
    inverse_root = _inverse(root)
    for index in range(1, degree):
        coefficients[index] += coefficients[index - 1] * inverse_root
    return coefficients


def _solve_quadratic(coefficients: Sequence[complex]) -> list[complex]:
    discriminant = cmath.sqrt(
        coefficients[1] * coefficients[1] - 4.0 * coefficients[2] * coefficients[0]
    )
    same_half_plane = (
        coefficients[1].real * discriminant.real
        + coefficients[1].imag * discriminant.imag
        >= 0.0
    )
    tmp = -coefficients[1] - discriminant if same_half_plane else -coefficients[1] + discriminant
    roots = [tmp / (2.0 * coefficients[2]), (2.0 * coefficients[0]) / tmp]
    for index, root in enumerate(roots):
        if _cnorm(root) <= 1.0:
            roots[index] -= _ratio2(coefficients, root, True)
        else:
            roots[index] /= 1.0 - root * _ratio2(coefficients, _inverse(root), False)
    return roots


def _solve_quartic(c: Sequence[complex]) -> list[complex]:
    r = 4.0 * c[4]
    p = 2.0 * r * c[2] - 3.0 * c[3] * c[3]
    q = 2.0 * c[3] * (c[3] * c[3] - r * c[2]) + r * r * c[1]
    D0 = c[2] * c[2] - 3.0 * c[3] * c[1] + 3.0 * c[0] * r
    D1 = (
        2.0 * c[2] * c[2] * c[2]
        - 9.0 * c[3] * c[2] * c[1]
        + 27.0 * c[3] * c[3] * c[0]
        + 27.0 * c[1] * c[1] * c[4]
        - 18.0 * r * c[2] * c[0]
    )
    tmp = cmath.sqrt(D1 * D1 - 4.0 * D0 * D0 * D0)
    Qc = (D1 + tmp) / 2.0 if D1.real * tmp.real + D1.imag * tmp.imag >= 0.0 else (D1 - tmp) / 2.0
    Qabs, Qarg = abs(Qc) ** (1.0 / 3.0), cmath.phase(Qc) / 3.0
    Q1 = cmath.rect(Qabs, Qarg)
    Q2 = cmath.rect(Qabs, Qarg + 2.0 * _PI / 3.0)
    Q3 = cmath.rect(Qabs, Qarg - 2.0 * _PI / 3.0)
    values = [
        (Q1 + D0 / Q1) * r - p,
        (Q2 + D0 / Q2) * r - p,
        (Q3 + D0 / Q3) * r - p,
    ]
    norms = [_cnorm(value) for value in values]
    largest = 0 if norms[0] > max(norms[1], norms[2]) else (1 if norms[1] > norms[2] else 2)
    tmp = cmath.sqrt(values[largest] / 3.0)
    S = tmp if c[4].real >= 0.0 else -tmp
    tmp = cmath.sqrt(q / S - S * S - p)
    roots = [(-c[3] - S + tmp) / r, (-c[3] - S - tmp) / r]
    tmp = cmath.sqrt(-q / S - S * S - p)
    roots.extend(((-c[3] + S + tmp) / r, (-c[3] + S - tmp) / r))
    for index, root in enumerate(roots):
        if _cnorm(root) <= 1.0:
            roots[index] -= _ratio4(c, root, True)
        else:
            roots[index] /= 1.0 - root * _ratio4(c, _inverse(root), False)
    return roots


def _polynomial_roots(
    degree: int,
    coefficients: Sequence[complex],
    roots: list[complex],
    min_error: float,
    max_error: float,
    trigonometric: bool,
    random: _DeterministicRandom,
) -> int:
    active = list(coefficients)
    iteration_count = 0
    root_index = 0
    active_degree = degree
    while active_degree > 4:
        roots[root_index], used = _newton_polynomial(
            active_degree, active, roots[root_index], max_error, min_error, 300, 10, random
        )
        iteration_count += used
        active = _extract_linear_factor(active_degree, active, roots[root_index])
        active_degree -= 1
        if active_degree > 4 and roots[root_index + 1] == 0.0j:
            if trigonometric:
                roots[root_index + 1] = roots[root_index] / _cnorm(roots[root_index])
            elif abs(roots[root_index].imag) > abs(roots[root_index].real) * 1.0e-3:
                roots[root_index + 1] = roots[root_index].conjugate()
            else:
                roots[root_index + 1] = complex(roots[root_index].imag, roots[root_index].real)
        root_index += 1
    roots[degree - 4 : degree] = _solve_quartic(active)
    return iteration_count + 4


def _func_g_ee_fast(index: int, data: _AuxData) -> float:
    x, y = _EXPS[index].real, _EXPS[index].imag
    A, B = data.Ps * y - data.Ss * x, data.Pp * y - data.Sp * x
    C = data.e2 * B - data.alpha1 * data.e1 * y * (1.0 - data.e1 * x)
    x_ = x - data.e1
    M, N = data.Sp * y + data.Pp * x_ + data.alpha2 * data.e2, -data.Ss * y - data.Ps * x_
    C2, A2, B2 = C * C, A * A, B * B
    Ac, Bc = A2 - C2, B2 - C2
    return data.K * data.K * Ac * Bc + 2.0 * data.K * C * (
        N * A * Ac + M * B * Bc
    ) - (A2 + B2) * (N * N * Ac + M * M * Bc - 2.0 * N * M * A * B)


def _func_g_eh_fast(index: int, data: _AuxData) -> float:
    x, y = _EXPS[index].real, _EXPS[index].imag
    A, B = data.Ps * y - data.Ss * x, data.Pp * y - data.Sp * x
    C = data.e2 * B - data.alpha1 * data.e1 * y * (1.0 - data.e1 * x)
    x_ = x - data.e1
    M, N = data.Sp * y + data.Pp * x_ + data.alpha2 * data.e2, -data.Ss * y - data.Ps * x_
    C2, A2, B2 = C * C, -A * A, B * B
    Ac, Bc = A2 - C2, B2 - C2
    return data.K * data.K * Ac * Bc + 2.0 * data.K * C * (
        -N * A * Ac + M * B * Bc
    ) - (A2 + B2) * (-N * N * Ac + M * M * Bc + 2.0 * N * M * A * B)


def _func_g_he_fast(index: int, data: _AuxData) -> complex:
    x, y = _EXPS[index].real, _EXPS[index].imag
    A, B = complex(data.Ps * y, -data.Ss * x), complex(data.Pp * y, -data.Sp * x)
    C = data.e2 * B - data.alpha1 * data.e1 * y * (1.0 - data.e1 * x)
    x_ = x - data.e1
    M = complex(data.Pp * x_ + data.alpha2 * data.e2, data.Sp * y)
    N = complex(-data.Ps * x_, -data.Ss * y)
    C2, A2, B2 = C * C, A * A, B * B
    Ac, Bc = A2 - C2, B2 - C2
    return data.K * data.K * Ac * Bc + 2.0 * data.K * C * (
        N * A * Ac + M * B * Bc
    ) - (A2 + B2) * (N * N * Ac + M * M * Bc - 2.0 * N * M * A * B)


def _func_g_hh_fast(index: int, data: _AuxData) -> complex:
    x, y = _EXPS[index].real, _EXPS[index].imag
    A, B = complex(data.Ps * y, -data.Ss * x), complex(data.Pp * y, -data.Sp * x)
    C = data.e2 * B - data.alpha1 * data.e1 * y * (1.0 - data.e1 * x)
    x_ = x - data.e1
    M = complex(data.Pp * x_ + data.alpha2 * data.e2, data.Sp * y)
    N = complex(-data.Ps * x_, -data.Ss * y)
    C2, A2, B2 = C * C, -(A * A), B * B
    Ac, Bc = A2 - C2, B2 - C2
    return data.K * data.K * Ac * Bc + 2.0 * data.K * C * (
        -N * A * Ac + M * B * Bc
    ) - (A2 + B2) * (-N * N * Ac + M * M * Bc + 2.0 * N * M * A * B)


def _create_polynomial_elliptic(
    data: _AuxData, function: Callable[[int, _AuxData], float]
) -> tuple[list[complex], float]:
    if not data.e1 <= 1.0:
        return [0.0j] * (DEG + 1), -1.0
    if data.e2 <= 1.0:
        hc = _sqr(data.alpha1 * _sqr(data.e1) / 16.0) * (
            complex(data.Pp - data.Ss - data.e1 * data.e2, data.Sp + data.Ps)
            * complex(data.Pp - data.Ss + data.e1 * data.e2, data.Sp + data.Ps)
            * complex(data.Pp + data.Ss - data.e1 * data.e2, data.Sp - data.Ps)
            * complex(data.Pp + data.Ss + data.e1 * data.e2, data.Sp - data.Ps)
        )
    else:
        hc = _sqr(data.alpha1 * _sqr(data.e1) / 16.0) * (
            complex(data.Pp - data.Ps - data.e1 * data.e2, data.Sp - data.Ss)
            * complex(data.Pp - data.Ps + data.e1 * data.e2, data.Sp - data.Ss)
            * complex(data.Pp + data.Ps - data.e1 * data.e2, data.Sp + data.Ss)
            * complex(data.Pp + data.Ps + data.e1 * data.e2, data.Sp + data.Ss)
        )
    values = [function(j, data) for j in range(DIM)]
    for j in range(1, DIM // 2 + 1):
        total, difference = values[j] + values[DIM - j], values[j] - values[DIM - j]
        values[j], values[DIM - j] = total, difference
    c = [0.0j] * (DEG + 1)
    a = values[0] + sum(values[1 : DIM // 2 + 1])
    c[DEG // 2] = complex(a / DIM)
    for i in range(1, DEG // 2):
        a, b = values[0], 0.0
        for j in range(1, DIM // 2 + 1):
            exponent = _EXPS[(j * i) % DIM]
            a += values[j] * exponent.real
            b += values[DIM - j] * exponent.imag
        a, b = a / DIM, b / DIM
        c[DEG // 2 - i], c[DEG // 2 + i] = complex(a, b), complex(a, -b)
    c[0], c[DEG] = hc, hc.conjugate()
    i = DEG // 2
    a, b = values[0], 0.0
    for j in range(1, DIM // 2 + 1):
        exponent = _EXPS[(j * i) % DIM]
        a += values[j] * exponent.real
        b += values[DIM - j] * exponent.imag
    a, b = a / DIM, b / DIM
    coefficient_error = _cnorm(complex(a, b) - hc)
    for i in range(DEG // 2 + 1, DIM // 2 + 1):
        a, b = values[0], 0.0
        for j in range(1, DIM // 2 + 1):
            exponent = _EXPS[(j * i) % DIM]
            a += values[j] * exponent.real
            b += values[DIM - j] * exponent.imag
        a, b = a / DIM, b / DIM
        coefficient_error += a * a + b * b
    coefficient_error /= (DIM - DEG + 1) // 2
    return c, max(math.sqrt(coefficient_error), _EPS * abs(hc))


def _create_polynomial_hyperbolic(
    data: _AuxData, function: Callable[[int, _AuxData], complex]
) -> tuple[list[complex], float]:
    if not data.e1 > 1.0:
        return [0.0j] * (DEG + 1), -1.0
    d = data.e2 * data.e2
    d1p, d2p = _sqr((data.Pp + data.Sp) / data.e1), _sqr((data.Ps + data.Ss) / data.e1)
    dpp, dpm = d1p + d2p, d1p - d2p
    factor = _sqr(data.alpha1 * _sqr(_sqr(data.e1)) / 16.0)
    h = factor * (dpp * dpp + d * d - 2.0 * d * dpm if data.e2 <= 1.0 else dpm * dpm + d * d - 2.0 * d * dpp)
    d1m, d2m = _sqr((data.Pp - data.Sp) / data.e1), _sqr((data.Ps - data.Ss) / data.e1)
    dmp, dmm = d1m + d2m, d1m - d2m
    hc = factor * (dmp * dmp + d * d - 2.0 * d * dmm if data.e2 <= 1.0 else dmm * dmm + d * d - 2.0 * d * dmp)
    values = [0.0j] + [function(j, data) for j in range(1, DIM // 2 + 1)]
    value0 = function(0, data).real
    c = [0.0j] * (DEG + 1)
    a = sum(value.real for value in values[1:])
    c[DEG // 2] = complex((2.0 * a + value0) / DIM)
    for i in range(1, DEG // 2):
        a = b = 0.0
        for j in range(1, DIM // 2 + 1):
            exponent = _EXPS[(j * i) % DIM]
            a += values[j].real * exponent.real
            b += values[j].imag * exponent.imag
        c[DEG // 2 - i] = complex((2.0 * (a - b) + value0) / DIM)
        c[DEG // 2 + i] = complex((2.0 * (a + b) + value0) / DIM)
    c[0], c[DEG] = complex(hc), complex(h)
    i = DEG // 2
    a = b = 0.0
    for j in range(1, DIM // 2 + 1):
        exponent = _EXPS[(j * i) % DIM]
        a += values[j].real * exponent.real
        b += values[j].imag * exponent.imag
    coefficient_error = _sqr((2.0 * (a - b) + value0) / DIM - hc) + _sqr(
        (2.0 * (a + b) + value0) / DIM - h
    )
    for i in range(DEG // 2 + 1, DIM // 2 + 1):
        a = b = 0.0
        for j in range(1, DIM // 2 + 1):
            exponent = _EXPS[(j * i) % DIM]
            a += values[j].real * exponent.real
            b += values[j].imag * exponent.imag
        coefficient_error += _sqr((2.0 * a + value0) / DIM) + _sqr(2.0 * b / DIM)
    coefficient_error /= DIM - DEG + 1
    return c, max(math.sqrt(coefficient_error), _EPS * max(abs(hc), abs(h)))


def _root_error(
    degree: int, coefficients: Sequence[complex], coefficient_error: float, z: complex
) -> float:
    magnitude2 = _cnorm(z)
    forward = magnitude2 <= 1.0
    if not forward:
        magnitude2 = 1.0 / magnitude2
    function_error2 = 1.0
    for _ in range(1, degree + 1):
        function_error2 = function_error2 * magnitude2 + 1.0
    function_error2 *= coefficient_error * coefficient_error
    value, derivative, derivative2 = _polynomial_valder2(
        degree, coefficients, z if forward else _inverse(z), forward
    )
    discriminant = cmath.sqrt(derivative * derivative - 2.0 * value * derivative2)
    d2 = 4.0 * _cnorm(value) / max(
        _cnorm(derivative + discriminant), _cnorm(derivative - discriminant)
    )
    return math.sqrt((d2 + function_error2 / _cnorm(discriminant)) / magnitude2)


def _max_root_bound(degree: int, coefficients: Sequence[complex], forward: bool) -> float:
    if degree < 1:
        return math.inf
    leading_norm = _cnorm(coefficients[degree if forward else 0])
    result = _cnorm(coefficients[degree - 1 if forward else 1]) / leading_norm
    for index in range(2, degree + 1):
        value = (_cnorm(coefficients[degree - index if forward else index]) / leading_norm) ** (1.0 / index)
        if result < value:
            result = value
    return math.sqrt(result)


def _circular_moid(data: _AuxData, nu: float, started: float) -> SMOIDResult:
    normal1 = _cross(data.P1, data.Q1)
    normal2 = _cross(data.P2, data.Q2)
    line = _cross(normal1, normal2)
    line_norm = _vnorm(line)
    if line_norm <= _EPS:
        direction = data.P1
    else:
        inverse_length = 1.0 / math.sqrt(line_norm)
        direction = tuple(component * inverse_length for component in line)
    result = SMOIDResult(
        distance=abs(data.a1 - data.a2),
        distance_error=(abs(data.a1) + abs(data.a2)) * _EPS * nu,
        u1=math.atan2(_dot(direction, data.Q1), _dot(direction, data.P1)),
        u1_error=_EPS * _PI * nu,
        u2=math.atan2(_dot(direction, data.Q2), _dot(direction, data.P2)),
        u2_error=_EPS * _PI * nu,
        min_delta=math.inf,
    )
    result.time = time.perf_counter() - started
    return result


def MOID_fast(
    O1: COrbitData,
    O2: COrbitData,
    max_root_error: float,
    min_root_error: float,
    nu: float = 1.0,
) -> SMOIDResult:
    """Calculate MOID using the degree-16 polynomial algorithm."""

    started = time.perf_counter()
    result = SMOIDResult()
    if not math.isfinite(max_root_error) or max_root_error <= 0.0:
        raise ValueError("maximum root error must be finite and positive")
    if not math.isfinite(min_root_error) or min_root_error < 0.0:
        raise ValueError("minimum root error must be finite and nonnegative")
    if not math.isfinite(nu) or nu <= 0.0:
        raise ValueError("error scaling must be finite and positive")
    data = _AuxData(O1, O2)
    if data.e1 == 0.0:
        if data.e2 == 0.0:
            return _circular_moid(data, nu, started)
        result = MOID_fast(O2, O1, max_root_error, min_root_error, nu)
        result.u1, result.u2 = result.u2, result.u1
        result.u1_error, result.u2_error = result.u2_error, result.u1_error
        result.time = time.perf_counter() - started
        return result
    if not (data.e1 != 1.0 and data.e2 != 1.0):
        result.good = False
        result.time = time.perf_counter() - started
        return result
    if data.e1 <= 1.0:
        coefficients, coefficient_error = _create_polynomial_elliptic(
            data, _func_g_ee_fast if data.e2 <= 1.0 else _func_g_eh_fast
        )
    else:
        coefficients, coefficient_error = _create_polynomial_hyperbolic(
            data, _func_g_he_fast if data.e2 <= 1.0 else _func_g_hh_fast
        )
    if (
        _cnorm(coefficients[0]) == 0.0
        or _cnorm(coefficients[DEG]) == 0.0
        or not math.isfinite(_cnorm(coefficients[0]))
        or not math.isfinite(_cnorm(coefficients[DEG]))
    ):
        result.good = False
        result.time = time.perf_counter() - started
        return result
    znrm = (
        1.0
        if data.e1 <= 1.0
        else (_cnorm(coefficients[0]) / _cnorm(coefficients[DEG])) ** (1.0 / (2.0 * DEG))
    )
    scale = 1.0
    for index in range(1, DEG // 2 + 1):
        scale *= znrm
        coefficients[DEG // 2 + index] *= scale
        coefficients[DEG // 2 - index] /= scale
    coefficient_error *= max(scale, 1.0 / scale)
    working_coefficients = list(coefficients)
    roots = [0.0j] * DEG
    random = _DeterministicRandom()

    index = 4
    max_abs_root = _max_root_bound(DEG, coefficients, True) + _EPS
    min_abs_root = (
        1.0 / max_abs_root
        if data.e1 <= 1.0
        else 1.0 / _max_root_bound(DEG, coefficients, False) + _EPS
    )
    if data.e1 <= 1.0:
        roots[0] = cmath.rect(min_abs_root, random.next() * _CIRC)
    else:
        roots[0] = cmath.rect(1.0 / znrm, 2.0 * _PI / 5.0)
        roots[2] = cmath.rect(1.0 / znrm, 3.0 * _PI / 5.0 if data.e2 <= 1.0 else _PI / 2.0)

    eta1, eta2 = math.sqrt(abs(1.0 - data.e1 * data.e1)), math.sqrt(abs(1.0 - data.e2 * data.e2))
    wPQ1, wPQ2 = math.hypot(data.P1w, data.Q1w), math.hypot(data.P2w, data.Q2w)
    u1min = 0.0 if data.e1 <= 0.0 else -math.log(max_abs_root * znrm)
    u1max = 0.0 if data.e1 <= 0.0 else -math.log(min_abs_root * znrm)
    u1 = _atan2_smart(data.Q1w * eta1, data.e1 * wPQ1 + data.P1w, data.e1 <= 1.0)
    br1 = wPQ1 + data.e1 * data.P1w >= 0.0
    u2 = _atan2_smart(data.Q2w * eta2, data.e2 * wPQ2 + data.P2w, data.e2 <= 1.0)
    br2 = wPQ2 + data.e2 * data.P2w >= 0.0
    u1_other = _atan2_smart(-data.Q1w * eta1, data.e1 * wPQ1 - data.P1w, data.e1 <= 1.0)
    br1_other = wPQ1 - data.e1 * data.P1w >= 0.0
    u2_other = _atan2_smart(-data.Q2w * eta2, data.e2 * wPQ2 - data.P2w, data.e2 <= 1.0)
    br2_other = wPQ2 - data.e2 * data.P2w >= 0.0

    candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, br1, br2, u1, u2)
    if not data.e1 <= 1.0:
        candidate_u1 = min(max(candidate_u1, u1min), u1max)
    roots[index] = cmath.rect(
        1.0 if data.e1 <= 1.0 else math.exp(-candidate_u1) * (1.0 if br1 else -1.0) / znrm,
        candidate_u1 if data.e1 <= 1.0 else random.next() * 1.0e-2 + 1.0e-3,
    )
    index += 1

    if (data.e1 <= 1.0 or br1_other != br1) and (data.e2 <= 1.0 or br2_other != br2):
        candidate_u1, candidate_u2, *_ = _sqdist_2Diter(
            data, br1_other, br2_other, u1_other, u2_other
        )
        branch_for_root = br1_other
    else:
        candidate_u1 = u1_other if data.e1 <= 1.0 or br1_other != br1 else (u1 + u1_other) / 2.0
        candidate_u2 = u2_other if data.e2 <= 1.0 or br2_other != br2 else (u2 + u2_other) / 2.0
        candidate_u1, candidate_u2, *_ = _sqdist_2Diter(
            data, data.e1 <= 1.0 or not br1, data.e2 <= 1.0 or not br2, candidate_u1, candidate_u2
        )
        branch_for_root = not br1
    if not data.e1 <= 1.0:
        candidate_u1 = min(max(candidate_u1, u1min), u1max)
    roots[index] = cmath.rect(
        1.0 if data.e1 <= 1.0 else math.exp(-candidate_u1) * (1.0 if branch_for_root else -1.0) / znrm,
        candidate_u1 if data.e1 <= 1.0 else random.next() * 1.0e-2 + 1.0e-3,
    )
    index += 1

    if data.e1 <= 1.0 and data.e2 <= 1.0:
        candidate_u1 = math.atan2(-data.P1w * eta1, data.e1 * wPQ1 + data.Q1w)
        candidate_u2 = math.atan2(-data.P2w * eta2, data.e2 * wPQ2 + data.Q2w)
        candidate_u2_other = math.atan2(data.P2w * eta2, data.e2 * wPQ2 - data.Q2w)
        if 2.0 * data.I > _PI:
            candidate_u2, candidate_u2_other = candidate_u2_other, candidate_u2
        candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, True, True, candidate_u1, candidate_u2)
        roots[index] = cmath.rect(1.0, candidate_u1)
        index += 1
        candidate_u1 = math.atan2(data.P1w * eta1, data.e1 * wPQ1 - data.Q1w)
        candidate_u1, candidate_u2_other, *_ = _sqdist_2Diter(
            data, True, True, candidate_u1, candidate_u2_other
        )
        roots[index] = cmath.rect(1.0, candidate_u1)
        index += 1
    else:
        if data.e1 <= 1.0 or br1_other != br1:
            candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, br1_other, br2, u1_other, u2)
            branch_for_root = br1_other
        else:
            candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, not br1, br2, (u1 + u1_other) / 2.0, u2)
            branch_for_root = not br1
        if not data.e1 <= 1.0:
            candidate_u1 = min(max(candidate_u1, u1min), u1max)
        roots[index] = cmath.rect(
            1.0 if data.e1 <= 1.0 else math.exp(-candidate_u1) * (1.0 if branch_for_root else -1.0) / znrm,
            candidate_u1 if data.e1 <= 1.0 else random.next() * 1.0e-2 + 1.0e-3,
        )
        index += 1
        if data.e2 <= 1.0 or br2_other != br2:
            candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, br1, br2_other, u1, u2_other)
        else:
            candidate_u1, candidate_u2, *_ = _sqdist_2Diter(data, br1, not br2, u1, (u2 + u2_other) / 2.0)
        if not data.e1 <= 1.0:
            candidate_u1 = min(max(candidate_u1, u1min), u1max)
        roots[index] = cmath.rect(
            1.0 if data.e1 <= 1.0 else math.exp(-candidate_u1) * (1.0 if br1 else -1.0) / znrm,
            candidate_u1 if data.e1 <= 1.0 else random.next() * 1.0e-2 + 1.0e-3,
        )
        index += 1

    if data.e1 <= 1.0:
        roots[index] = cmath.rect(min_abs_root, random.next() * _CIRC)
    else:
        roots[index] = cmath.rect(1.0 / znrm, 2.0 * _PI / 5.0)
        roots[index + 2] = cmath.rect(
            1.0 / znrm, 3.0 * _PI / 5.0 if data.e2 <= 1.0 else _PI / 2.0
        )

    result.iter_count = _polynomial_roots(
        DEG, working_coefficients, roots, min_root_error, max_root_error,
        data.e1 <= 1.0, random,
    )
    for root in roots:
        root_error = _root_error(DEG, coefficients, coefficient_error, root) * nu
        delta = (
            abs(math.log(_cnorm(root))) / (2.0 * root_error)
            if data.e1 <= 1.0
            else (abs(cmath.phase(root)) if root != 0.0j else 0.0) / root_error
        )
        if delta < 10.0:
            result.good = result.good and root_error < max_root_error
            if delta < 3.0:
                result.root_count += 1
                first_anomaly = (
                    cmath.phase(root)
                    if data.e1 <= 1.0
                    else (-math.log(root.real * znrm) if root.real > 0.0 else math.inf)
                )
                distance, second_anomaly, _ = _min_distance_for(data, first_anomaly)
                distance = abs(distance)
                if distance >= 0.0 and (result.distance < 0.0 or distance < result.distance):
                    result.distance, result.u1, result.u2 = distance, first_anomaly, second_anomaly
        if not delta < 3.0 and (result.min_delta < 0.0 or delta < result.min_delta):
            result.min_delta = delta

    if result.distance < 0.0 or not math.isfinite(result.u1) or not math.isfinite(result.u2):
        result.good = False
        result.time = time.perf_counter() - started
        return result

    (
        hessian_sign,
        result.u1,
        result.u2,
        result.u1_error,
        result.u2_error,
        result.distance,
        result.distance_error,
        result.iter_count_2D,
        gradient,
        hessian,
    ) = _newton_sqdist(data, result.u1, result.u2, min_root_error, 30)
    hessian_sign *= 1 if data.a1 * data.a2 >= 0.0 else -1

    angle_error = _EPS * _PI * nu
    result.u1_error += angle_error
    result.u2_error += angle_error
    scale = 2.0 * data.a1 * data.a2
    x1 = math.cos(result.u1) if data.e1 <= 1.0 else math.cosh(result.u1)
    x2 = math.cos(result.u2) if data.e2 <= 1.0 else math.cosh(result.u2)
    tmp1, tmp2 = data.e1 * x1, data.e2 * x2
    r1, r2 = data.a1 * (1.0 - tmp1), data.a2 * (1.0 - tmp2)
    rd1 = abs(data.a1) * math.sqrt(abs(1.0 - tmp1 * tmp1))
    rd2 = abs(data.a2) * math.sqrt(abs(1.0 - tmp2 * tmp2))
    radial_error = math.hypot(r1, r2) * _EPS * nu
    gradient_error = radial_error * math.hypot(rd1, rd2) / abs(scale)
    max_eigenvalue = abs(hessian[0] + hessian[1]) / 2.0 + math.hypot(
        (hessian[0] - hessian[1]) / 2.0, hessian[2]
    )
    determinant = abs(hessian[0] * hessian[1] - hessian[2] * hessian[2])
    third_anomaly_error = max_eigenvalue / determinant * gradient_error
    result.u1_error += third_anomaly_error
    result.u2_error += third_anomaly_error
    result.distance_error += max_eigenvalue / 2.0 * (
        angle_error * angle_error + gradient_error * gradient_error / determinant
    )
    result.distance *= scale
    result.distance_error *= abs(scale)
    result.distance = _sqrt_ieee(result.distance)
    result.distance_error += 2.0 * result.distance * radial_error + radial_error * radial_error
    result.distance_error /= math.sqrt(
        result.distance * result.distance + result.distance_error / 2.0
    )
    result.good = (
        result.good
        and result.root_count >= (4 if data.e1 <= 1.0 else 2)
        and result.root_count % 2 == 0
        and result.min_delta > 10.0
        and hessian_sign == 2
        and abs(result.u1_error) <= max_root_error
    )
    result.time = time.perf_counter() - started
    return result


__all__ = [
    "COrbitData",
    "SMOIDResult",
    "SLCResult",
    "detect_suitable_options",
    "test_peri_apo",
    "LC",
    "MOID_fast",
    "MOID_direct_search",
    "restrict_search_range",
]