Animates the convergence of root distributions of expectation polynomials towards their theoretical free probability limits as the dimension/degree $d$ scales up to 300:

• Wigner Semicircle Convergence: GUE expectation polynomials additive convolution $He_d \boxplus_d He_d$ roots converging to the semicircle distribution.
• Marchenko-Pastur Convergence: Laguerre/Wishart expectation polynomials roots converging to the Marchenko-Pastur law.
• Free Jacobi Arcsine Convergence: Jacobi/Legendre polynomials roots converging to the classical Arcsine limit.
• Free Log-Normal Convergence: Multi-precision compound Wishart multiplicative convolutions roots converging to the free Log-normal distribution.

Wigner Semicircle Convergence ($He_{d} \boxplus_{d} He_{d}$ as $d \to 300$)	Marchenko-Pastur Convergence (Laguerre as $d \to 300$)

Interactive companion: Open ensemble convergence to compare sampled spectra and Hermite/Laguerre/Legendre root measures with their limiting laws.

Free Jacobi Arcsine Convergence (Legendre as $d \to 300$)	Free Log-Normal Convergence (Wishart as $m=d \to 300$)

Interactive companion: Open ensemble convergence and select Legendre roots for the arcsine comparison or compound Wishart for a bounded-degree, mean-one free-lognormal comparison (τ = 1).

Showcases GUE spectral universality and eigenvalue spacing statistics:

• CD Bulk Scaling Limit: Christoffel-Darboux kernel bulk scaling limit convergence to the infinite Sine kernel.
• CD Edge Scaling Limit: Christoffel-Darboux kernel edge scaling limit convergence to the infinite Airy kernel.
• Tracy-Widom Edge Convergence: The batch illustration compares HKPV max-eigenvalue samples with a finite-rank continuous Hermite determinant. The interactive companion separately compares finite-rank Hermite CDFs with an Airy-kernel approximation to Tracy–Widom β = 2.
• Nearest-Neighbor Level Repulsion: GUE bulk eigenvalue spacing statistics matching the analytical Wigner Surmise.

CD Bulk Scaling Limit (Sine Kernel)	CD Edge Scaling Limit (Airy Kernel)

Interactive companion: Open Hermite kernels to switch between bulk/sine and soft-edge/Airy sections and vary the projection rank.

Tracy-Widom Edge Convergence ($d=100$)	Nearest-Neighbor Level Repulsion (Bulk Spacings vs Wigner Surmise)

Interactive companion: Open Hermite kernels and select continuous gap / edge CDF. Vary rank and threshold; inspect quadrature and tail-cutoff sensitivity beside the finite-rank and Tracy–Widom β = 2 curves.

Illustrates exact algebraic properties, root interlacing constraints, and geometric distributions of FiniteFree polynomials:

• Root Interlacing Examples: Illustrates interlacing under additive and multiplicative convolution in the displayed examples. Preservation theorems require the corresponding real-rooted and non-negative-root domains; strict interlacing is not an unconditional API guarantee.
• Multivariate Pencil Topography: Visualizes the log-determinant topography and exact eigenvalue boundaries of random symmetric matrix pencils.
• Unitary Eigenvalue Repulsion: Visualizes roots evolving on the complex unit circle $\mathbb{T}$ under unitary Hermite polynomial flows.

Original Roots ($p \prec q$)	Additive Convolution ($p \boxplus_6 r \prec q \boxplus_6 r$)	Multiplicative Convolution ($p \boxtimes_6 r \prec q \boxtimes_6 r$)

Interactive companions: Open convolution interlacing for additive/multiplicative examples, or Open moving line roots for matrix-pencil restrictions and derivative interlacing.

Hyperbolic Matrix Pencil Eigencones & Topography	Unitary Hermite Eigenvalue Repulsion Trajectories

Interactive companions: Open the hyperbolicity cone explorer for linked cone sections, or Open unitary root flow to vary degree and time.

Visualizes the convergence of exact finite free transforms to their continuous free probability limits as the dimension $d$ scales up to 300:

• T-Transform Step Convergence: Fujie-Ueda Finite T-Transform step function convergence to the continuous free Marchenko-Pastur (MP) limit $T(t)=1+t$.
• Asymptotic Decay of Cumulants: Linearizing additive convolution, showing exact high-order finite free cumulants $\kappa_k(d)$ decaying to zero, leaving only $\kappa_2 \to 1$ (the Wigner semicircle limit).
• S-Transform Multiplicativity: S-transforms of Wishart expectation polynomials, demonstrating exact multiplicativity $S_{p \boxtimes_d q}^{(d)} = S_p^{(d)} S_q^{(d)}$ at all discrete nodes and convergence to the free S-transform limit.
• Cauchy Transform Convergence: Domain coloring of the Finite Cauchy Transform $G_d(z) = \frac{1}{d} \frac{p'_d(z)}{p_d(z)}$, showing the roots (poles) and derivative roots (zeros) coalescing into a continuous branch cut $[-2, 2]$ of the semicircle law.

Fujie-Ueda T-Transform Convergence	Asymptotic Decay of Cumulants

Interactive companion: Open finite transforms for the Wishart T-transform and additive-CLT cumulants. The S-transform and Cauchy-domain animations remain separate batch illustrations.

S-Transform Convergence	Cauchy Transform Convergence

• Certified Real-Rooted Polynomial Validation: Verified lazily with exact rational Sturm sequences for small squarefree factors and certified Arb isolation for larger factors. Exact Sturm remains a fallback through degree 30 when Arb cannot obtain a certificate.
• Unitary Circle Geometries ($\mathbb{T}$): UnitaryPolynomial bypasses real-line verification and returns numerical complex roots ordered by angle. The constructor trusts the intended unit-circle domain; it does not certify it. unitary_hermite_polynomial stores its exponential coefficients symbolically.
• Lazy Geometric Domain Properties: After certifying real-rootedness, has_non_negative_roots and has_strictly_positive_roots use an $O(d)$ coefficient sign check to enforce operator domains. Sign alternation alone cannot certify real-rootedness.
• Basic Polynomial Transformations: Supports exact algebraic transformations including variable dilation (dilation), variable shift (shift), positive integer root powers (power), root-reciprocal reversing (reversed_polynomial), derivative (derivative), projection (projection), fractional additive convolution power (additive_power), and the Fujie-Ueda limiting polynomial $\Phi_d$ (phi_d). Noninteger root powers use numerical root isolation and reconstruction.
• Direct Projections: projection(j) computes the monic derivative of order d-j from the leading j+1 coefficients, avoiding intermediate derivative objects. Known Hermite inputs retain their variance, center and numerical root recurrence. Projection to the original degree returns the same object. See the API reference for normalization and benchmark scope.
• High-Degree Scaling & Root Reconstruction:
  • Divide-and-Conquer Polynomial Synthesis: RealRootedPolynomial.from_roots(roots) multiplies rational linear factors using a balanced product tree. Arithmetic costs depend on FLINT multiplication and coefficient bit sizes.
• Working Precision: PrecisionContext(degree, prec=None) sets FLINT's process-wide precision to prec, or max(53, int(2.5*degree)) bits, and restores the previous value on exit. It does not manage memory pools or provide independent thread-local precision.

• Jacobi Polynomials (jacobi_polynomial): Exact $O(n^2)$ recurrence construction of $P_n^{(\alpha, \beta)}(x)$ using flint.fmpq_poly in exact rational arithmetic.
• Hahn Polynomials (hahn_polynomial): Exact $O(n)$ recurrence construction of $Q_n(x; \alpha, \beta, N)$ using sequential running products to prevent redundant Pochhammer factorial operations.
• Jack Polynomials (jack_polynomial): High-performance variables recursion computing symmetric Jack polynomials $J_\lambda^{(\alpha)}(x_1, \dots, x_m)$ strictly over sparse monomial exponent dictionaries, completely bypassing SymPy's slow symbolic expansion engine.
• Chebyshev & Legendre Polynomials (chebyshev_t_polynomial, chebyshev_u_polynomial, legendre_polynomial): Exact recurrence constructions of $T_n(x)$, $U_n(x)$, and $P_n(x)$ over $\mathbb{Q}$ via flint.fmpq_poly.

• Symmetric Additive ($\boxplus_d$): Exact analytical evaluation mapping polynomial convolutions against the symmetric finite combinatorial variance.
• Asymmetric Additive ($\uplus_d$): Exact factorial-weighted coefficient convolution. Optional weights=[a,b] dilates the input roots by a and b; it does not change their degrees or rank.
• Multiplicative ($\boxtimes_d$): Exact discrete multiplicative root transformations via scaled Hadamard projections.

• Finite Cauchy Transform ($G_p^{(d)}$)
• Finite S-Transform ($S_p^{(d)}$): Discrete normalized evaluations bypassing non-linear mapping.
• Finite R-Transform ($R_p^{(d)}$): Computes finite free cumulants ($\kappa_n^{(d)}$) via an exact $O(n^2)$ recursive generating function sequence map over $\mathbb{Q}$, bypassing exponential partition lattice enumeration while verifying exact additivity $\kappa_n^{(d)}(p \boxplus_d q) = \kappa_n^{(d)}(p) + \kappa_n^{(d)}(q)$.
• Finite T-Transform ($T_p^{(d)}$): Step function mapping the right-continuous inverse to the Fujie-Ueda limit $\Phi_d$, evaluated in $O(d)$ algebraically using coefficient sign-alternation validation.
• Symmetric Finite S-Transform: SymmetricFiniteSTransform(p, convention="standard") returns the positive-imaginary square root defined for symmetric real-rooted inputs. The default convention="ratio" preserves the existing ratio of consecutive normalized even coefficients. See the API convention.

• MultivariatePolynomial: Sparse rational polynomials with exact evaluation, rational algebra, explicit variable reordering, simultaneous scalar/polynomial substitution, exact gradient/Hessian polynomials and line restriction. Target symbols and their order are explicit; operations do not certify stability or hyperbolicity. Homogeneous multinomial normalization requires a homogeneous input. See multivariate workflows.
• Sparse Evaluation and Ownership: evaluate uses exact rational coordinates; evaluate_float64, gradient_float64 and hessian_float64 evaluate real batches with explicit float64 rounding/range limits. Numerical derivatives work at zero coordinates and singular determinant points. Numerical work uses bounded point blocks; Hessian output still requires batch size times variable count squared storage. Native inputs, context changes and exports own their data. Costs depend on monomial count, degree, variable count, batch size and coefficient sizes.
• Jacobi SLP & Reverse AD: Generic straight-line programs support reverse-mode gradients. Determinant pencils instead use matrix determinant/inverse and trace identities for gradients and Hessians; exact derivatives require a nonsingular evaluated matrix.
• Prepared Numerical Pencils: Retain an owned float64 snapshot and reuse batched LU factors for slogdet, log-absolute-determinant gradients and Hessian-vector products. NumPy is the default; CuPy is explicitly selected for resident GPU arrays. See the prepared API contract and its memory and singular-domain limits.
• Product-Grid Modular Interpolation (CRT): Evaluates exact determinant polynomials and pencil characteristic polynomials (bypassing symbolic expansion bottlenecks via exact rational interpolation over $\mathbb{Q}[t]$) using C-level nmod_mat solvers and the Chinese Remainder Theorem.
• Zippel Sparse Interpolation: Uses randomized finite-field discovery followed by deterministic coefficient/support verification for sparse determinant evaluations (from_symmetric_matrix_pencil_sparse). An exact balanced-base fallback repairs incomplete discovery. A configurable integer-size limit bounds verification feasibility. See the construction limits.
• Multiplicative Pencils & Diagonal Specialization: Extends exact matrix pencil geometries to generalized asymmetric forms (MultiplicativeMatrixPencil) and computes univariate characteristic projections via optimized 1D Chinese Remainder Theorem loops.
• LMI Cone Verification: Positive definiteness checks ($A(e) \succ 0$) to verify hyperbolic cones.

• to_numpy_poly1d(): Direct coefficient export to NumPy, with a warning above degree 20. It does not rescale coefficients or guarantee finite values; prefer evaluate_roots_float64() for root extraction.
• evaluate_roots_float64(exact=False): Uses scaled symmetric tridiagonal eigenvalues for Hermite, Laguerre, Jacobi, Legendre and Chebyshev polynomials in their orthogonality domains, including GUE/Wishart expectations, affine transforms and proven Hermite additive convolutions. It avoids forming an ill-conditioned monomial companion matrix for these families.
• evaluate_roots_float64(exact=True): Requests the Arb isolation path by default. A high-precision request after a numerical call bypasses the approximate cache. Results are returned as float64; if isolation fails, the general numerical fallbacks can still be used. See the API reference for dispatch, precision and conditioning limits.

• Matrix Samplers: Dense random matrices for GOE (orthogonal invariance, $\beta=1$), GUE (unitary invariance, $\beta=2$), and GSE (symplectic invariance, $\beta=4$). Here $O(d)$ denotes the orthogonal group, not a runtime bound.
• Empirical Validations: Computes theoretical expected characteristic polynomials $\mathbb{E}[\det(xI - M)]$ matching explicit orthogonal sequences.
• Integrated GUE Workflows: Connect sampling, expected polynomials, scaled Hermite kernels and interval-count concentration with reproducible examples and explicit finite-sample uncertainty. See the GUE workflow.
• EmpiricalComparison: Compares all sampled characteristic coefficients, including deterministic ones, with Bonferroni-adjusted two-sided Student-t bands. alpha controls the nominal family significance level; coefficient scaling prevents variance overflow/underflow. Samples must be finite Hermitian matrices of the matching size, with verified eigenvalue pairs for doubled-size representations. See the API reference for sample assumptions and numerical tolerances.

• Kernels: Construct exact discrete DPP kernels and orthogonal polynomial kernels (e.g., Hermite/Laguerre) via the Christoffel-Darboux formula. Exact arguments retain their distinct values. Nearby distinct floating points use the finite basis sum to avoid a cancelled quotient; diagonal calls retain the derivative formula. An empty orthogonal basis gives the zero kernel. See the API reference for conditioning and timing scope.
• Gap Probabilities & Observables: Evaluate exact discrete gap probabilities and approximate continuous Fredholm determinants via Nyström discretization, plus exact $k$-point correlation functions.
• HKPV Sampler: Sample real symmetric DPP correlation kernels using floating-point spectral decomposition and Gram-Schmidt projections. Nonprojection kernels use Bernoulli eigenvector selection. Matrix shape, finiteness, symmetry and spectrum in $[0,1]$ are validated with $10^{-10}$ absolute roundoff tolerance.

You can construct polynomials exactly from sequences of rational coefficients and evaluate their roots using numerical eigensolvers or Arb isolation.

from finitefree.core import RealRootedPolynomial
import sympy as sp

# Initialize exactly via rational/integer coefficients: (x - 1)(x - 2) = x^2 - 3x + 2
p = RealRootedPolynomial([1, -3, 2])

# Dilate, shift, and take positive integer root powers exactly
p_dilated = p.dilation(2)             # x^2 - 6x + 8 (roots scaled by 2)
p_shifted = p.shift(1)                # (x - 2)(x - 3) = x^2 - 5x + 6 (roots + 1)
p_powered = p.power(2)                # (x - 1)(x - 4) = x^2 - 5x + 4 (roots squared)
p_reversed = p.reversed_polynomial()  # x^2 - 1.5x + 0.5 (reciprocal roots)

# Domain properties use coefficient signs after real-rootedness verification
print(p.has_non_negative_roots)       # True
print(p.has_strictly_positive_roots)  # True

# Perform certified complex interval root isolation (Arb)
roots = p.evaluate_roots_float64(exact=True)
print(roots)  # [1.0, 2.0]

# Use stable numerical roots for asymptotic comparisons of known families
from finitefree import gue_expected_poly
gue_roots = gue_expected_poly(100).evaluate_roots_float64(exact=False)
assert len(gue_roots) == 100

Convolutions combine polynomial coefficients exactly. Additive convolution preserves real-rootedness; multiplicative convolution does so when one real-rooted input has non-negative roots. Formal coefficient formulas can also produce polynomials with complex roots, so reconstructed outputs are verified lazily.

from finitefree.core import RealRootedPolynomial
from finitefree.convolutions import symmetric_additive, asymmetric_additive, multiplicative
import sympy as sp

# Instantiate two real-rooted polynomials of degree d=2
p = RealRootedPolynomial([1, -3, 2])  # roots: 1, 2
q = RealRootedPolynomial([1, 0, -4])  # roots: -2, 2

# Symmetric Additive Convolution (p [+]_d q)
res_add = symmetric_additive(p, q, d=2)
print(res_add.coeffs)  # [1, -3, -2]

# Asymmetric Additive Convolution (p [u]_d q), after dilating each input by 1/2
res_asym = asymmetric_additive(p, q, weights=[sp.Rational(1, 2), sp.Rational(1, 2)], d=2)
print(res_asym.coeffs)

# Multiplicative Convolution (p [*]_d q)
res_mult = multiplicative(p, q, d=2)
print(res_mult.coeffs)  # Hadamard-like projection

Finite free probability transforms compute expected spectral properties and algebraic limits without combinatorial partition search.

from finitefree import FiniteRTransform, FiniteTTransform, SymmetricFiniteSTransform
from finitefree.orthogonal import laguerre_polynomial

# Initialize a polynomial (must be non-negative rooted for T-transform)
poly = laguerre_polynomial(n=3, alpha=1)

# Compute Finite Free Cumulants exactly via generating function recurrences (O(d^2))
# Additivity holds: kappa(p [+] q) = kappa(p) + kappa(q)
r_transform = FiniteRTransform(poly)
cumulant_3 = r_transform[2]
print(f"3rd Finite Free Cumulant: {cumulant_3}")

# Map inverse limit points using the Fujie-Ueda Finite T-Transform
t_transform = FiniteTTransform(poly)
print(t_transform(0.5))

Finite free cumulants use the normalization in Definition 2.14 of Arizmendi et al., $\kappa_n^{(d)}=(-d)^{n-1}c_n/(n-1)!$, where $c_n$ is the classical cumulant of the normalized coefficient sequence. additive_power uses the matching inverse. Orders above the ambient dimension return zero as an API convention; they are outside the finite cumulant definition.

normalized_coeffs(d) divides out the leading coefficient before applying the binomial normalization, so e_0=1 even when a polynomial was stored with monic=False. Finite cumulants and coefficient convolutions therefore depend on the roots rather than a scalar multiple of the polynomial. The stored coefficients and polynomial evaluations retain that scalar; reconstructed/convolved results are monic. Larger ambient dimensions pad roots with zeros.

FiniteRTransform(..., numerical=True, prec=256) centers the requested coefficient prefix in rational arithmetic, computes the recurrence with Arb, and restores the first cumulant (the mean). This avoids cancellation from large translations without shifting the entire polynomial or computing all cumulants exactly. Working arithmetic uses the requested precision, but outputs are Python floats rather than arbitrary-precision values or error bounds. High orders can require more precision, and float range/rounding still apply.

FiniteTTransform requires positive degree and non-negative roots, and evaluates a right-continuous step function for finite real inputs in $(0,1)$. Rational inputs retain their exact position, including values arbitrarily close to an endpoint. Floats select intervals using their stored binary values; use sympy.Rational(k, d) when an exact grid boundary is intended.

Reconstruction from normalized coefficients produces a formal polynomial with lazy real-rootedness validation. Root extraction and positive-root domain checks reject complex-rooted inputs; formal coefficient convolutions remain available.

Jacobi and Laguerre polynomials constructed outside their orthogonality domains are also validated lazily. Formal construction remains available, but a root request rejects a complex-rooted result such as laguerre_polynomial(2, -3).

FiniteFree builds classical and symmetric orthogonal systems exactly via optimized recursive relations.

from finitefree.orthogonal import (
    jacobi_polynomial,
    hahn_polynomial,
    jack_polynomial,
    hermite_polynomial,
)
from finitefree.ensembles import sample_gue, gue_expected_poly

# Construct Jacobi, Hahn, and Hermite recurrence relations exactly over Q
h_prob = hermite_polynomial(n=4, physicist=False)  # Probabilist Hermite He_4
jacobi = jacobi_polynomial(n=3, alpha=1, beta=1)     # Jacobi P_3^(1, 1)
hahn   = hahn_polynomial(n=2, alpha=1, beta=1, N=5)  # Hahn Q_2

# High-performance multivariate Jack polynomials using sparse exponent dicts
jack = jack_polynomial(m=3, partition=[2, 1], alpha=2)
print(jack.expr)

# Compare random matrix characteristic polynomials with theoretical sequences
M = sample_gue(d=4)
expected_poly = gue_expected_poly(d=4)
print(expected_poly.coeffs)

Evaluate homogeneous determinants $\det(x_1 A_1 + \dots + x_m A_m)$ exactly via modular matrix interpolation. General rational polynomials also support explicit context changes and batched numerical derivatives:

import numpy as np
import sympy as sp
from finitefree.multivariate import MultivariatePolynomial

x, y, u, v = sp.symbols("x y u v")
p = MultivariatePolynomial(x**2 + 3*x*y + y**2, [x, y])
assert p.reorder_variables([y, x]).evaluate([2, 1]) == p.evaluate([1, 2])
q = p.substitute({x: u+v, y: u-v}, variables=[u, v])
assert q.expr == 5*u**2 - v**2
np.testing.assert_array_equal(q.gradient_float64([[2, 1], [0, 0]]), [[20, -2], [0, 0]])
np.testing.assert_array_equal(q.hessian_float64([0, 0]), [[10, 0], [0, -2]])

Matrix-pencil constructors copy integer/rational entries before preparing a separate float64 numerical view. Characteristic polynomials, rational SLP derivatives and multivariate determinants use the copied values, including integers above $2^{53}$. Float inputs preserve their supplied binary values; precision lost before construction cannot be recovered. Matrices must be square, share one shape, contain finite real rational values that fit the finite numerical view, and be symmetric on the supplied values for SymmetricMatrixPencil. Construct a new pencil to change its entries.

from finitefree.hyperbolic import SymmetricMatrixPencil
from finitefree.multivariate import MultivariatePolynomial
import numpy as np

# Define a matrix pencil
A1 = np.eye(3)
A2 = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])
pencil = SymmetricMatrixPencil([A1, A2])

# Construct the multivariate polynomial exactly using modular determinants and CRT
poly_crt = MultivariatePolynomial.from_symmetric_matrix_pencil_interpolated(pencil)
print(poly_crt.expr)

# Or use sparse discovery with deterministic verification and exact fallback
poly_sparse = MultivariatePolynomial.from_symmetric_matrix_pencil_sparse(pencil)
print(poly_sparse.expr)

Construct correlation kernels, evaluate k-point joint intensities exactly, compute Fredholm determinant gap probabilities, and sample point configurations numerically.

from finitefree import (
    DiscreteFiniteKernel,
    OrthogonalPolynomialKernel,
    hermite_polynomial,
    gap_probability_discrete,
    gap_probability_continuous,
    sample_discrete,
)
import numpy as np
import math
import sympy as sp

# --- 1. Discrete DPP Kernel & Exact Gap Probability ---
# 3x3 nonprojection correlation kernel with eigenvalues 1, 1, 1/3
K_mat = [
    [sp.Rational(2, 3), sp.Rational(1, 3), 0],
    [sp.Rational(1, 3), sp.Rational(2, 3), 0],
    [0, 0, 1],
]
discrete_kernel = DiscreteFiniteKernel(K_mat)

# Exact gap probability on states {0, 1}: det(I - K_{[0,1]})
prob_discrete = gap_probability_discrete(discrete_kernel, [0, 1])
print(f"Exact Gap Probability on {{0,1}}: {prob_discrete}")  # 0

# --- 2. Orthogonal Polynomial Kernel & CD Formula ---
# Probabilist's Hermite polynomials for GUE / Wigner ensemble DPP
norms = [1, 1, 2] # h_0, h_1, h_2
polys = [
    hermite_polynomial(0, physicist=False),
    hermite_polynomial(1, physicist=False),
    hermite_polynomial(2, physicist=False),
    hermite_polynomial(3, physicist=False),
]
hermite_kernel = OrthogonalPolynomialKernel(polys, norms)

# Evaluate kernel off-diagonal and diagonal (CD derivatives)
print(f"CD Kernel K(0.5, 1.2): {hermite_kernel(0.5, 1.2)}")
print(f"CD Kernel Diagonal K(0.5, 0.5): {hermite_kernel(0.5, 0.5)}")

# --- 3. Continuous Fredholm Determinant ---
# Hermite orthogonal polynomial weight function
def weight_func(x):
    return math.exp(-x**2 / 2.0) / math.sqrt(2.0 * math.pi)

gap_prob = gap_probability_continuous(hermite_kernel, -1.0, 1.0, n_points=20, weight_func=weight_func)
print(f"Continuous Gap Probability over [-1.0, 1.0]: {gap_prob}")

# --- 4. HKPV Discrete Sampling ---
sampled_states = sample_discrete(discrete_kernel, state_space=[0, 1, 2])
print(f"HKPV Sampled point configuration: {sampled_states}")