Multiple-precision arithmetic

Computing moments of phase-type distributions requires Gaussian elimination on the underlying graph. For models with many states or widely varying transition rates, the intermediate multipliers in this elimination can span many orders of magnitude. When the ratio of the largest to smallest multiplier (the condition number) is large, standard 64-bit floating-point arithmetic loses precision and the computed moments become unreliable.

Phasic integrates the MPFR (Multiple Precision Floating-Point Reliable) library to handle these ill-conditioned cases transparently. MPFR provides arbitrary-precision arithmetic so that moment computations remain accurate regardless of condition number.

Key points:

MPFR activates automatically when the condition number exceeds a threshold — no user intervention needed.
Well-conditioned graphs use fast double-precision arithmetic with zero overhead.
Users can force MPFR on, adjust the activation threshold, or set a specific precision.
MPFR applies to expected_waiting_time(), expected_sojourn_time(), expectation(), variance(), and moments().

from phasic import Graph, with_ipv, configure, get_config, get_available_options, set_log_level
import numpy as np
from vscodenb import set_vscode_theme
set_vscode_theme()

Checking MPFR availability

MPFR support is compiled into the C backend when libmpfr and libgmp are present at build time. You can check whether it is available on your system:

opts = get_available_options()
print(f"MPFR available: {opts['mpfr']}")

Automatic activation

By default, phasic monitors the condition number of the elimination graph every time it computes moments. If the condition number exceeds the threshold (default 1e12), MPFR arithmetic is activated automatically for that computation. The precision is chosen adaptively based on the condition number.

This means you do not need to think about MPFR in most cases — it just works. To see it in action, enable info-level logging:

set_log_level('INFO')

nr_samples = 10

@with_ipv([nr_samples] + [0] * (nr_samples - 1))
def coalescent(state):
    transitions = []
    for i in range(state.size):
        for j in range(i, state.size):
            same = int(i == j)
            if same and state[i] < 2:
                continue
            if not same and (state[i] < 1 or state[j] < 1):
                continue
            new = state.copy()
            new[i] -= 1
            new[j] -= 1
            new[i + j + 1] += 1
            transitions.append((new, state[i] * (state[j] - same) / (1 + same)))
    return transitions

graph = Graph(coalescent)
print(f"Graph has {graph.vertices_length()} vertices")

# Moments are computed automatically — MPFR kicks in if needed
ewt = graph.expected_waiting_time()
print(f"Expected waiting time: {ewt[0]:.6f}")

set_log_level('WARNING')

If the condition number is below the threshold, you will see no MPFR-related log messages and the computation uses fast double-precision arithmetic.

Configuring MPFR

The configure() function provides four parameters for controlling MPFR behavior:

Parameter	Default	Description
`force_high_precision`	`False`	Force MPFR for all moment computations regardless of condition number
`mpfr_precision_bits`	`0`	Precision in bits. `0` means auto-determine from condition number
`condition_threshold`	`1e12`	Condition number above which MPFR auto-activates
`enable_condition_warnings`	`True`	Show warnings when ill-conditioning is detected

Forcing MPFR on

If you want maximum accuracy for all moment computations regardless of conditioning:

configure(force_high_precision=True)

ewt = graph.expected_waiting_time()
print(f"Expected waiting time (MPFR forced): {ewt[0]:.6f}")

# Reset to default
configure(force_high_precision=False)

Setting precision

When mpfr_precision_bits=0 (the default), the precision is determined automatically as log₂(condition_number) + 64 bits, clamped to the range [128, 1024]. You can override this with a specific value:

# Use 256-bit precision explicitly
configure(force_high_precision=True, mpfr_precision_bits=256)

ewt = graph.expected_waiting_time()
print(f"Expected waiting time (256-bit): {ewt[0]:.6f}")

# Reset
configure(force_high_precision=False, mpfr_precision_bits=0)

Adjusting the activation threshold

The condition threshold controls when MPFR auto-activates. A lower threshold is more conservative (activates MPFR more often), while a higher threshold tolerates more ill-conditioning before switching.

# More conservative: activate MPFR at lower condition numbers
configure(condition_threshold=1e8)

# Less conservative: only activate for severely ill-conditioned graphs
configure(condition_threshold=1e15)

# Disable auto-activation entirely (only manual force_high_precision works)
configure(condition_threshold=float('inf'))

# Reset to default
configure(condition_threshold=1e12)

Suppressing warnings

If MPFR is not available and a graph is ill-conditioned, phasic logs a warning suggesting you rebuild with MPFR. To silence these:

configure(enable_condition_warnings=False)

# Reset
configure(enable_condition_warnings=True)

Inspecting current configuration

You can read back the active MPFR settings at any time:

config = get_config()
print(f"Force MPFR:          {config.force_high_precision}")
print(f"Precision bits:      {config.mpfr_precision_bits} (0 = auto)")
print(f"Condition threshold: {config.condition_threshold:.0e}")
print(f"Condition warnings:  {config.enable_condition_warnings}")

Environment variables

The same settings can be controlled via environment variables, which is useful in batch scripts or CI pipelines where you cannot modify Python code:

Environment variable	Equivalent `configure()` parameter
`PHASIC_FORCE_MPFR=1`	`force_high_precision=True`
`PHASIC_MPFR_BITS=256`	`mpfr_precision_bits=256`
`PHASIC_CONDITION_THRESHOLD=1e15`	`condition_threshold=1e15`
`PHASIC_DISABLE_CONDITION_WARNINGS=1`	`enable_condition_warnings=False`

# Force 256-bit MPFR for a batch run
PHASIC_FORCE_MPFR=1 PHASIC_MPFR_BITS=256 python my_analysis.py

Calling configure() in Python sets these environment variables internally, so the C backend picks them up.

How it works

When phasic computes moments via Gaussian elimination, the C backend:

Pre-scans the elimination graph to find the largest and smallest multiplier magnitudes.
Computes the condition number as their ratio: max(|multiplier|) / min(|multiplier|).
Checks whether force_high_precision is set or the condition number exceeds condition_threshold.
If MPFR is needed, determines precision either from mpfr_precision_bits or automatically as log₂(condition) + 64 bits (clamped to [128, 1024]).
Re-executes the elimination using MPFR arithmetic at the chosen precision.
Converts the result back to double-precision for return to Python.

The MPFR elimination graph is cached on the graph object, so repeated moment calls on the same graph do not recompute it.

Performance

Scenario	Overhead
Well-conditioned graph (MPFR not activated)	None
MPFR at 128 bits	~5–10× slower than double
MPFR at 256 bits	~8–15× slower
MPFR at 512 bits	~15–25× slower
MPFR at 1024 bits	~25–40× slower

The overhead only applies to the moment computation itself. Graph construction, PDF/CDF evaluation, and sampling are unaffected. In most workflows, moment computation is a small fraction of total runtime, so the practical impact is negligible.

Building with MPFR support

If get_available_options()['mpfr'] returns False, you need to install the MPFR and GMP libraries and rebuild:

pixi add mpfr gmp pkg-config
pixi install
pixi run install-dev

conda install mpfr gmp pkg-config -c conda-forge
pip install --no-build-isolation --force-reinstall --no-deps .

The build system (CMakeLists.txt) detects MPFR automatically and defines the HAVE_MPFR compile flag when it is found.