Skip to main content

omnibias.core.verified.riesz

0.5.0rc1Release candidatePython >=3.10Apache-2.0

Coefficient recurrences, activation contracts, jet algebra, and verified arithmetic shared by the numerical backends.

import omnibias.core.verified.riesz
Module contract & mathematical conventions

Verified 2-D Riesz transform & Leray projection on a closed-form blob basis.

The 1-D self-similar machinery in omnibias.core.verified.line works in the basis where the line Hilbert transform acts *exactly* (the Poisson / conjugate-Poisson pair). Two-dimensional model equations (2-D Euler / SQG) need the planar analogue: the **Riesz transforms** R_j (Fourier multiplier -i\xi_j/|\xi|) and the **Leray projection** P = I - \nabla\Delta^{-1}\nabla\cdot (multiplier \delta_{ik} - \xi_i\xi_k/|\xi|^2) that removes the divergence of a vector field.

The trick that makes this *closed form* (hence amenable to interval arithmetic) is to work on a basis of **radial blobs that are the Laplacian of an explicit Newtonian potential**. For a scale a > 0 set D = x^2 + y^2 + a^2 and

N_a(x) = \tfrac{1}{4\pi}\,\ln D \quad(\text{the smoothed Newtonian potential}),\qquad f_a = \Delta N_a = \frac{a^2}{\pi\,D^2}\ (\ge 0),

so N_a = \Delta^{-1} f_a is known in closed form. The *second-order* Riesz composite that the Leray projection is built from is then **elementary**:

R_i R_k = -\partial_i\partial_k\Delta^{-1} \;\Longrightarrow\; (R_i R_k f_a)(x) = -\partial_i\partial_k N_a(x) = \frac{1}{2\pi}\,\frac{2 x_i x_k - \delta_{ik} D}{D^2}.

(The *single* Riesz transform of a radial blob is not elementary -- it carries a half-Laplacian |\xi|^{-1} -- but the Leray projection only ever needs the composite R_iR_k, which is.) Two consequences are used as built-in checks:

  • R_{11} f + R_{22} f = -f exactly (the multiplier identity (\xi_1^2+\xi_2^2)/|\xi|^2 = 1);
  • for a vector blob v = (c_1, c_2) f_a the Leray projection (Pv)_i = c_i f_a + \sum_k (R_iR_k f_a)\,c_k is **divergence-free by an algebraic cancellation**, so \nabla\cdot(Pv) \equiv 0 -- here certified by a closed-form interval residual that encloses 0.

Every quantity is an outward-rounded Interval; the only transcendental is the rigorous ln (the potential) and exp/ln in the far-field tail bound. The closed forms are *proved* (Fourier multipliers) and *checked* in the tests against an independent mpmath derivative of N_a.

Shared formulas do not imply universal floating-point identity. Verify dtype, shape, derivative-order support, and the stated numerical assumptions for the path you use.

FUNCTION

blob​

omnibias.core.verified.blobSource ↗
blob(x: float, y: float, a: float) -> Interval

Verified radial blob f_a = \Delta N_a = a^2 / (\pi D^2) (>= 0).

Parameters

NameType & defaultContract
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ReturnsInterval
Explicit validation & failure cases
  • When a == 0.0, raises ValueError('blob scale a must be non-zero').

These are explicit source guards. Backend operations may also raise for invalid shapes, dtypes, or unsupported devices.

Executed regression · input validation

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_input_validation() -> None:
    with pytest.raises(ValueError):
        blob(0.1, 0.2, 0.0)  # zero scale
    with pytest.raises(ValueError):
        riesz_double_blob(2, 0, 0.1, 0.2, 0.7)  # bad axis
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_input_validation" -q
Regression example · blob is laplacian of potential

This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_blob_is_laplacian_of_potential() -> None:
    for x, y in _PTS:
        hxx = mp.diff(lambda t: _N(t, y, _A), x, 2)  # noqa: B023
        hyy = mp.diff(lambda t: _N(x, t, _A), y, 2)  # noqa: B023
        lap = float(hxx + hyy)
        fb = blob(x, y, _A)
        assert fb.lo <= lap <= fb.hi
python -m pytest packages/omnibias-core/tests/verified/test_riesz.py -k test_blob_is_laplacian_of_potential -q
Read the implementation
def blob(x: float, y: float, a: float) -> Interval:
    r"""Verified radial blob ``f_a = \Delta N_a = a^2 / (\pi D^2)`` (``>= 0``)."""
    if a == 0.0:
        raise ValueError("blob scale a must be non-zero")
    d = _denom(x, y, a)
    return Interval.point(a).pow_int(2) * (_PI * d.pow_int(2)).reciprocal()
Other public import paths
  • omnibias.core.verified.blob
  • omnibias.core.verified.riesz.blob
FUNCTION

blob_gradient​

omnibias.core.verified.blob_gradientSource ↗
blob_gradient(x: float, y: float, a: float) -> tuple[Interval, Interval]

Verified gradient of the radial blob \nabla f_a.

From f_a = a^2/(\pi D^2) with D = x^2+y^2+a^2 one differentiates to the closed form

\nabla f_a = -\frac{4 a^2}{\pi\,D^3}\,(x, y),

a purely *radial* field (parallel to (x, y)). This is the \nabla\omega of a 2-D Euler vorticity blob; paired with the tangential Biot--Savart velocity u = \nabla^\perp N_a \parallel (-y, x) it makes the steady-state residual u\cdot\nabla\omega vanish by exact perpendicularity.

Parameters

NameType & defaultContract
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
Returnstuple[Interval, Interval]
Explicit validation & failure cases
  • When a == 0.0, raises ValueError('blob scale a must be non-zero').

These are explicit source guards. Backend operations may also raise for invalid shapes, dtypes, or unsupported devices.

Executed regression · blob gradient zero scale rejected

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_blob_gradient_zero_scale_rejected() -> None:
    with pytest.raises(ValueError):
        blob_gradient(0.1, 0.2, 0.0)
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_blob_gradient_zero_scale_rejected" -q
Regression example · blob gradient matches mpmath

This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_blob_gradient_matches_mpmath() -> None:
    for x, y in _PTS:
        fx = float(mp.diff(lambda t: _f(t, y, _A), x))  # noqa: B023
        fy = float(mp.diff(lambda t: _f(x, t, _A), y))  # noqa: B023
        gx, gy = blob_gradient(x, y, _A)
        assert gx.lo <= fx <= gx.hi
        assert gy.lo <= fy <= gy.hi
python -m pytest packages/omnibias-core/tests/verified/test_riesz.py -k test_blob_gradient_matches_mpmath -q
Read the implementation
def blob_gradient(x: float, y: float, a: float) -> tuple[Interval, Interval]:
    r"""Verified gradient of the radial blob ``\nabla f_a``.

    From ``f_a = a^2/(\pi D^2)`` with ``D = x^2+y^2+a^2`` one differentiates to the
    closed form

    .. math::

        \nabla f_a = -\frac{4 a^2}{\pi\,D^3}\,(x, y),

    a purely *radial* field (parallel to ``(x, y)``).  This is the
    ``\nabla\omega`` of a 2-D Euler vorticity blob; paired with the tangential
    Biot--Savart velocity ``u = \nabla^\perp N_a \parallel (-y, x)`` it makes the
    steady-state residual ``u\cdot\nabla\omega`` vanish by exact perpendicularity.
    """
    if a == 0.0:
        raise ValueError("blob scale a must be non-zero")
    d = _denom(x, y, a)
    fac = Interval.point(-4.0) * Interval.point(a).pow_int(2) * (_PI * d.pow_int(3)).reciprocal()
    return Interval.point(x) * fac, Interval.point(y) * fac
Other public import paths
  • omnibias.core.verified.blob_gradient
  • omnibias.core.verified.riesz.blob_gradient
FUNCTION

leray_blob_field​

omnibias.core.verified.leray_blob_fieldSource ↗
leray_blob_field(c1: float, c2: float, x: float, y: float, a: float) -> tuple[Interval, Interval]

Verified Leray projection of the vector blob v = (c_1, c_2) f_a.

(Pv)_i = c_i f_a + \sum_k (R_iR_k f_a) c_k collapses to the closed forms

(Pv)_1 = \frac{c_1(a^2+x^2-y^2) + 2 c_2 x y}{2\pi D^2},\qquad (Pv)_2 = \frac{c_2(a^2+y^2-x^2) + 2 c_1 x y}{2\pi D^2},

which are divergence-free (see leray_divergence_residual).

Parameters

NameType & defaultContract
c1positional-or-keywordfloatRequiredSee the operation contract and implementation below.
c2positional-or-keywordfloatRequiredSee the operation contract and implementation below.
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
Returnstuple[Interval, Interval]
Executed regression · leray field matches mpmath closed form

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_leray_field_matches_mpmath_closed_form() -> None:
    c1, c2 = 1.3, -0.8
    for x, y in _PTS:
        d = x * x + y * y + _A * _A
        denom = 2 * mp.pi * d * d
        pv1_true = (c1 * (_A**2 + x * x - y * y) + 2 * c2 * x * y) / denom
        pv2_true = (c2 * (_A**2 + y * y - x * x) + 2 * c1 * x * y) / denom
        pv1, pv2 = leray_blob_field(c1, c2, x, y, _A)
        assert pv1.lo <= float(pv1_true) <= pv1.hi
        assert pv2.lo <= float(pv2_true) <= pv2.hi
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_leray_field_matches_mpmath_closed_form" -q
Read the implementation
def leray_blob_field(
    c1: float, c2: float, x: float, y: float, a: float
) -> tuple[Interval, Interval]:
    r"""Verified Leray projection of the vector blob ``v = (c_1, c_2) f_a``.

    ``(Pv)_i = c_i f_a + \sum_k (R_iR_k f_a) c_k`` collapses to the closed forms

    .. math::

        (Pv)_1 = \frac{c_1(a^2+x^2-y^2) + 2 c_2 x y}{2\pi D^2},\qquad
        (Pv)_2 = \frac{c_2(a^2+y^2-x^2) + 2 c_1 x y}{2\pi D^2},

    which are divergence-free (see :func:`leray_divergence_residual`).
    """
    d = _denom(x, y, a)
    inv = (_TWO_PI * d.pow_int(2)).reciprocal()
    a2 = Interval.point(a).pow_int(2)
    x2 = Interval.point(x).pow_int(2)
    y2 = Interval.point(y).pow_int(2)
    xy = Interval.point(x) * Interval.point(y)
    c1i, c2i = Interval.point(c1), Interval.point(c2)
    pv1 = (c1i * (a2 + x2 - y2) + Interval.point(2.0) * c2i * xy) * inv
    pv2 = (c2i * (a2 + y2 - x2) + Interval.point(2.0) * c1i * xy) * inv
    return pv1, pv2
Other public import paths
  • omnibias.core.verified.leray_blob_field
  • omnibias.core.verified.riesz.leray_blob_field
FUNCTION

leray_divergence_residual​

omnibias.core.verified.leray_divergence_residualSource ↗
leray_divergence_residual(c1: float, c2: float, x: float, y: float, a: float) -> Interval

Closed-form enclosure of \nabla\cdot(Pv) for the blob field above.

Computes \partial_x (Pv)_1 + \partial_y (Pv)_2 from the *separate* component derivatives (each a non-trivial rational function); the sum is analytically 0, so the returned interval is a certificate that the two independently-evaluated rationals cancel -- it must enclose 0 with a width reflecting only outward rounding.

Parameters

NameType & defaultContract
c1positional-or-keywordfloatRequiredSee the operation contract and implementation below.
c2positional-or-keywordfloatRequiredSee the operation contract and implementation below.
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ReturnsInterval
Executed regression · leray is divergence free residual encloses zero

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_leray_is_divergence_free_residual_encloses_zero() -> None:
    for c1, c2 in [(1.3, -0.8), (0.0, 1.0), (-2.1, 0.5)]:
        for x, y in _PTS:
            res = leray_divergence_residual(c1, c2, x, y, _A)
            assert res.lo <= 0.0 <= res.hi
            assert res.width < 1e-10
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_leray_is_divergence_free_residual_encloses_zero" -q
Read the implementation
def leray_divergence_residual(
    c1: float, c2: float, x: float, y: float, a: float
) -> Interval:
    r"""Closed-form enclosure of ``\nabla\cdot(Pv)`` for the blob field above.

    Computes ``\partial_x (Pv)_1 + \partial_y (Pv)_2`` from the *separate*
    component derivatives (each a non-trivial rational function); the sum is
    analytically ``0``, so the returned interval is a certificate that the two
    independently-evaluated rationals cancel -- it must enclose ``0`` with a width
    reflecting only outward rounding.
    """
    d = _denom(x, y, a)
    inv3 = (_TWO_PI * d.pow_int(3)).reciprocal()
    a2 = Interval.point(a).pow_int(2)
    x2 = Interval.point(x).pow_int(2)
    y2 = Interval.point(y).pow_int(2)
    xi, yi = Interval.point(x), Interval.point(y)
    xy = xi * yi
    c1i, c2i = Interval.point(c1), Interval.point(c2)
    two = Interval.point(2.0)
    four = Interval.point(4.0)
    a1 = c1i * (a2 + x2 - y2) + two * c2i * xy  # numerator of 2pi D^2 (Pv)_1
    a2_num = c2i * (a2 + y2 - x2) + two * c1i * xy  # numerator of 2pi D^2 (Pv)_2
    grad_s = two * c1i * xi + two * c2i * yi  # d/dx and d/dy of the linear parts
    dx_pv1 = (grad_s * d - four * xi * a1) * inv3
    dy_pv2 = (grad_s * d - four * yi * a2_num) * inv3
    return dx_pv1 + dy_pv2
Other public import paths
  • omnibias.core.verified.leray_divergence_residual
  • omnibias.core.verified.riesz.leray_divergence_residual
FUNCTION

newtonian_potential​

omnibias.core.verified.newtonian_potentialSource ↗
newtonian_potential(x: float, y: float, a: float) -> Interval

Verified N_a = \Delta^{-1} f_a = (1/4\pi)\ln(x^2+y^2+a^2).

Parameters

NameType & defaultContract
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ReturnsInterval
Executed regression · potential matches mpmath

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_potential_matches_mpmath() -> None:
    for x, y in _PTS:
        enc = newtonian_potential(x, y, _A)
        assert enc.lo <= float(_N(x, y, _A)) <= enc.hi
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_potential_matches_mpmath" -q
Read the implementation
def newtonian_potential(x: float, y: float, a: float) -> Interval:
    r"""Verified ``N_a = \Delta^{-1} f_a = (1/4\pi)\ln(x^2+y^2+a^2)``."""
    return ln_iv(_denom(x, y, a)) * _FOUR_PI.reciprocal()
Other public import paths
  • omnibias.core.verified.newtonian_potential
  • omnibias.core.verified.riesz.newtonian_potential
FUNCTION

potential_gradient​

omnibias.core.verified.potential_gradientSource ↗
potential_gradient(x: float, y: float, a: float) -> tuple[Interval, Interval]

Verified \nabla N_a = (1/2\pi)(x, y)/D.

Parameters

NameType & defaultContract
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
Returnstuple[Interval, Interval]
Read the implementation
def potential_gradient(x: float, y: float, a: float) -> tuple[Interval, Interval]:
    r"""Verified ``\nabla N_a = (1/2\pi)(x, y)/D``."""
    d = _denom(x, y, a)
    inv = (_TWO_PI * d).reciprocal()
    return Interval.point(x) * inv, Interval.point(y) * inv
Other public import paths
  • omnibias.core.verified.potential_gradient
  • omnibias.core.verified.riesz.potential_gradient
FUNCTION

riesz_double_blob​

omnibias.core.verified.riesz_double_blobSource ↗
riesz_double_blob(i: int, k: int, x: float, y: float, a: float) -> Interval

Verified second-order Riesz composite (R_i R_k f_a)(x).

Equals -\partial_i\partial_k N_a = (1/2\pi)(2 x_i x_k - \delta_{ik} D)/D^2; i, k index the axes (0 = x, 1 = y). This is the closed-form Calderon--Zygmund building block of the Leray projection.

Parameters

NameType & defaultContract
ipositional-or-keywordintRequiredSee the operation contract and implementation below.
kpositional-or-keywordintRequiredSee the operation contract and implementation below.
xpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ypositional-or-keywordfloatRequiredSee the operation contract and implementation below.
apositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ReturnsInterval
Executed regression · input validation

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_input_validation() -> None:
    with pytest.raises(ValueError):
        blob(0.1, 0.2, 0.0)  # zero scale
    with pytest.raises(ValueError):
        riesz_double_blob(2, 0, 0.1, 0.2, 0.7)  # bad axis
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_input_validation" -q
Regression example · riesz trace identity equals minus blob

This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_riesz_trace_identity_equals_minus_blob() -> None:
    # R_11 f + R_22 f = -f  (the (xi_1^2 + xi_2^2)/|xi|^2 = 1 multiplier identity).
    for x, y in _PTS:
        r11 = riesz_double_blob(0, 0, x, y, _A)
        r22 = riesz_double_blob(1, 1, x, y, _A)
        fb = blob(x, y, _A)
        s = r11 + r22
        assert s.lo <= -fb.mid <= s.hi
        # tight: the trace sum and -blob agree to rounding
        assert (s + fb).abs().hi < 1e-12
python -m pytest packages/omnibias-core/tests/verified/test_riesz.py -k test_riesz_trace_identity_equals_minus_blob -q
Read the implementation
def riesz_double_blob(i: int, k: int, x: float, y: float, a: float) -> Interval:
    r"""Verified second-order Riesz composite ``(R_i R_k f_a)(x)``.

    Equals ``-\partial_i\partial_k N_a = (1/2\pi)(2 x_i x_k - \delta_{ik} D)/D^2``;
    ``i, k`` index the axes (``0`` = x, ``1`` = y).  This is the closed-form
    Calderon--Zygmund building block of the Leray projection.
    """
    _require_axis("i", i)
    _require_axis("k", k)
    comp = (x, y)
    d = _denom(x, y, a)
    num = Interval.point(2.0) * Interval.point(comp[i]) * Interval.point(comp[k])
    if i == k:
        num = num - d
    return num * (_TWO_PI * d.pow_int(2)).reciprocal()
Other public import paths
  • omnibias.core.verified.riesz_double_blob
  • omnibias.core.verified.riesz.riesz_double_blob
FUNCTION

riesz_tail_bound​

omnibias.core.verified.riesz_tail_boundSource ↗
riesz_tail_bound(kernel_const: float, decay_const: float, decay_power: float, x_trunc: float, core_radius: float) -> Interval

Rigorous far-field bound on a 2-D Calderon--Zygmund (Riesz) convolution.

The second-order Riesz kernel obeys |K(x)| <= C_K / |x|^2. If a profile g satisfies |g(t)| <= C |t|^{-p} for |t| >= X (p > 0) and the evaluation point has |x_0| <= \rho < X, then the contribution of the region |t| >= X to (K * g)(x_0) is bounded by

|\text{tail}| \le \frac{2\pi\,C_K\,C\,X^{-p}}{p\,(1 - \rho/X)^2},

from |x_0 - t| \ge |t|(1 - \rho/X) and the planar polar integral \int_{|t|\ge X} |t|^{-p-2}\,dt = 2\pi X^{-p}/p. The returned symmetric interval [-B, B] is an outward-rounded enclosure of that bound, so a finite-basis Leray/Riesz evaluation can be upgraded to a rigorous full-plane statement by adding it.

Parameters mirror the inequality: kernel_const = C_K (a bound on |K(x)|\,|x|^2), decay_const = C, decay_power = p, x_trunc = X, core_radius = \rho.

Parameters

NameType & defaultContract
kernel_constpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
decay_constpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
decay_powerpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
x_truncpositional-or-keywordfloatRequiredSee the operation contract and implementation below.
core_radiuspositional-or-keywordfloatRequiredSee the operation contract and implementation below.
ReturnsInterval
Explicit validation & failure cases
  • When decay_power <= 0.0, raises ValueError('decay_power p must be > 0 for the tail to converge').
  • When not (x_trunc > 0.0 and 0.0 <= core_radius < x_trunc), raises ValueError('require 0 <= core_radius < x_trunc and x_trunc > 0').
  • When decay_const < 0.0 or kernel_const < 0.0, raises ValueError('decay_const C and kernel_const C_K must be non-negative').

These are explicit source guards. Backend operations may also raise for invalid shapes, dtypes, or unsupported devices.

Executed regression · riesz tail bound grows with core radius

Recorded coverage confirms this test executed the API body. This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_riesz_tail_bound_grows_with_core_radius() -> None:
    base = riesz_tail_bound(1.0, 1.0, 3.0, 2.0, 0.0).hi
    closer = riesz_tail_bound(1.0, 1.0, 3.0, 2.0, 1.5).hi
    assert closer > base  # shrinking the gap X - rho inflates the bound
python -m pytest "packages/omnibias-core/tests/verified/test_riesz.py::test_riesz_tail_bound_grows_with_core_radius" -q
Regression example · riesz tail bound validation

This example belongs to the pinned upstream test suite. It may use fixtures and imports from its full file. Run it from that checkout with its test dependencies installed.

def test_riesz_tail_bound_validation() -> None:
    with pytest.raises(ValueError):
        riesz_tail_bound(1.0, 1.0, 0.0, 2.0, 0.0)  # p <= 0
    with pytest.raises(ValueError):
        riesz_tail_bound(1.0, 1.0, 3.0, 2.0, 2.0)  # rho >= X
    with pytest.raises(ValueError):
        riesz_tail_bound(-1.0, 1.0, 3.0, 2.0, 0.0)  # negative kernel const
python -m pytest packages/omnibias-core/tests/verified/test_riesz.py -k test_riesz_tail_bound_validation -q
Read the implementation
def riesz_tail_bound(
    kernel_const: float,
    decay_const: float,
    decay_power: float,
    x_trunc: float,
    core_radius: float,
) -> Interval:
    r"""Rigorous far-field bound on a 2-D Calderon--Zygmund (Riesz) convolution.

    The second-order Riesz kernel obeys ``|K(x)| <= C_K / |x|^2``.  If a profile
    ``g`` satisfies ``|g(t)| <= C |t|^{-p}`` for ``|t| >= X`` (``p > 0``) and the
    evaluation point has ``|x_0| <= \rho < X``, then the contribution of the
    region ``|t| >= X`` to ``(K * g)(x_0)`` is bounded by

    .. math::

        |\text{tail}| \le
        \frac{2\pi\,C_K\,C\,X^{-p}}{p\,(1 - \rho/X)^2},

    from ``|x_0 - t| \ge |t|(1 - \rho/X)`` and the planar polar integral
    ``\int_{|t|\ge X} |t|^{-p-2}\,dt = 2\pi X^{-p}/p``.  The returned symmetric
    interval ``[-B, B]`` is an outward-rounded enclosure of that bound, so a
    finite-basis Leray/Riesz evaluation can be upgraded to a rigorous full-plane
    statement by adding it.

    Parameters mirror the inequality: ``kernel_const = C_K`` (a bound on
    ``|K(x)|\,|x|^2``), ``decay_const = C``, ``decay_power = p``,
    ``x_trunc = X``, ``core_radius = \rho``.
    """
    if decay_power <= 0.0:
        raise ValueError("decay_power p must be > 0 for the tail to converge")
    if not (x_trunc > 0.0 and 0.0 <= core_radius < x_trunc):
        raise ValueError("require 0 <= core_radius < x_trunc and x_trunc > 0")
    if decay_const < 0.0 or kernel_const < 0.0:
        raise ValueError("decay_const C and kernel_const C_K must be non-negative")
    x_iv = Interval.point(x_trunc)
    x_neg_p = exp_iv(Interval.point(-decay_power) * ln_iv(x_iv))  # X^{-p}
    one_minus = Interval.point(1.0) - Interval.point(core_radius) * x_iv.reciprocal()
    denom = Interval.point(decay_power) * one_minus.pow_int(2)
    bound = (
        _TWO_PI
        * Interval.point(kernel_const)
        * Interval.point(decay_const)
        * x_neg_p
        * denom.reciprocal()
    )
    b = bound.hi
    return Interval(-b, b)
Other public import paths
  • omnibias.core.verified.riesz_tail_bound
  • omnibias.core.verified.riesz.riesz_tail_bound