<figure><img src="image_1.png" /></figure>

HAL
open science

# Certified Algebraic Path Tracking with Algpath

## Alexandre Guillemot

► To cite this version:

Alexandre Guillemot. Certified Algebraic Path Tracking with Algpath. 2025. hal-05308166

HAL Id: hal-05308166

https://inria.hal.science/hal-05308166v1

Preprint submitted on 10 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Certified Algebraic Path Tracking with Algpath

Alexandre Guillemot

Inria

Université Paris-Saclay

Palaiseau, France, 91120

alexandre.guillemot@inria.fr

## Abstract

Algpath is a certified homotopy continuation software. We upgrade the previous fixed-precision Rust implementation by incorporating mixed, adaptive precision with minimal overhead. This allows us to tackle problems on which the initial implementation fails due to the inability to increase precision, and where uncertified methods may fail or path jump.

## 1 Introduction

Let $F: \mathbb{C} \times \mathbb{C}^n \rightarrow \mathbb{C}^n$ be a polynomial map (we denote by $F_t: \mathbb{C}^n \rightarrow \mathbb{C}^n$ the specialization of $F$ in its first variable), let $\zeta \in \mathbb{C}^n$ be a zero of $F_0$. We assume that there exists a unique continuation $\zeta_t$ of $\zeta$, meaning that $\zeta_t$ depends continuously on $t$, $\zeta_0 = \zeta$ and $F_t(\zeta_t) = 0$ for all $t \in [0, 1]$. Given $F$ and $\zeta$, homotopy continuation algorithms “track” $\zeta_t$, by initializing $z = \zeta$ at $t = 0$, then incrementing $t$ in small steps while updating $z$ to stay close to $\zeta_t$, ultimately stopping when $t$ reaches 1. When applied to polynomial system solving, heuristic tracking methods combined with a posteriori certification may suffice. Although various tools exist for this purpose [1, 3, 5, 18], they are inadequate if we care about the continuations themselves. The algorithm we implement in this work is certified, meaning that it outputs an interval path $\tilde{\zeta}_t$ such that for all $t \in [0, 1]$, $\zeta_t \in \tilde{\zeta}_t$ and is the only root of $F_t$ contained within. Previous theoretical work and implementations on this topic include [2, 6, 7, 9, 11, 13, 19]. Certified homotopy continuation finds applications in computations of Galois groups [8], braids [15] and certified algebraic curve projections [4].

In [7] we describe a corrector-predictor loop based on interval arithmetic, where termination is ensured through precision management. We previously implemented a fixed-precision version of the algorithm in Rust, which aborted computations whenever more precision was required. In the present work we extend this to mixed adaptive precision: not only can it increase precision instead of aborting computations, using Arb [10] as a backend (an arbitrary precision interval arithmetic library based on balls $[x \pm r]$), but when possible it falls back to the fixed-precision interval arithmetic used in the initial implementation (based on intervals of the form $[a, b]$ where $a$ and $b$ are IEEE 754 64-bit floating point numbers, implemented following [12]), thereby benefiting from the performance advantages of fixed precision. Managing two interval arithmetics simultaneously combined with the use of Rust's generics to avoid code duplication presents some implementation challenges.
## 2 Adaptive corrector-predictor loop

The algorithm in [7] introduces precision increases to ensure termination. While this is sufficient for termination, incorporating precision decreases improves practical performance, though it requires careful handling of certified data. The algorithm relies on *Moore boxes*, which are zero-isolating regions for polynomial systems. Decreasing the precision of some data representing a Moore box could compromise its isolating property, and in such a case an additional validation step is necessary. *Precision management*. If precision is part of (interval-based) data and managed globally, precision decreases modify the data itself, which we want to avoid. To ensure correctness, we adopt a less restrictive model where operations and functions depend on precision, while precision is not encoded directly in interval data: changing precision changes operations, while (certified) data stays unchanged. In this model, the precision of the data is indirectly changed by performing operations on it. This is the model we chose in [7], so we can lower precision by one bit before a round of corrector-predictor; this will not hinder correction nor termination. Arb [10] also follows this model, in contrast with MPFI [14], an adaptive precision interval arithmetic library that encodes precision within intervals, and operations are performed at the precision of the result operand. Moreover, we observed that Arb outperforms MPFI by a factor of two when evaluating polynomials. For these reasons, our implementation uses Arb.

## 3 Incorporating mixed precision

While we are able to increase and decrease precision using Arb, we lose performance compared to the initial implementation when computations can be performed using double precision only. To remedy this, we introduce mixed precision. In this section, “double precision” refers to the double precision interval arithmetic present in the initial fixed-precision implementation.

*Exact conversions.* We cannot assume exact conversions between Arb and double precision, even from double precision. For example, consider $[-2^{-50}, 2] \subseteq \mathbb{R}$. This corresponds to the ball $[(1 - 2^{-51}) \pm (1 + 2^{-51})]$, but $1 + 2^{-51}$ cannot be represented by a mag_t (radius type for Arb balls) since this type has a fixed mantissa of 30 bits. We could extend the idea of the previous section: we build a union arithmetic for which intervals are either Arb balls or double precision intervals. Operations would determine the appropriate arithmetic for the output depending on the precision; for example adding two double precision intervals would return an Arb ball if the precision is over 53. This significantly complexifies the implementation when trying to minimize overhead over running the algorithm entirely in double precision. Not only would we need to implement elementary operations of the union arithmetic, but also large portions of the code where precision does not change, to avoid unnecessary precision checks. Falling back to data conversions, we notice that the only certified data in our algorithm are Moore boxes $m$ for a certain time $t$, which are elements of $\mathbb{C}^n \times \mathbb{R}_> \times \mathbb{C}^{n \times n}$, and are, in practice, represented by singleton intervals (so is $t$). Converting singletons behaves as expected: any double precision interval $[a, a]$ can be converted to a ball $[a \pm 0]$ and conversely if the precision of $a$ does not exceed 53 bits, and its exponent fits in the exponent of a double. So before a round of corrector-predictor, we try to convert $m$ and $t$ to double precision (this might fail) and we try to perform a double precision corrector-predictor round (which may also fail). If either of these steps raise an error, we convert $m$ and $t$ to Arb and perform an adaptive precision corrector-predictor round. In practice, this approach allows us to preserve the performances of the implementation presented in [7] on problems that can be solved using fixed precision only.
<table><caption>Table 1: Comparison of *HomotopyContinuation.jl* and *Algpath*. *dim.*: number of variables; *deg.*: total degree of the system; *time*: total time to track all paths; *fail.*: number of reported failures; *jump.*: path jumping occured; *med.*: median number of iterations used over all paths; *max.*: maximum number of iterations in a single path; *prec.*: maximum number of precision used in a single path. * A hundred randomly picked starting zeros of a total degree homotopy. Except the last three rows, the problems we test are linear homotopies $F_t(x) = tf(x) + \gamma(1-t)g(x)$ where $\gamma \in_R [0, 1] + i[0, 1]$ ($i \in_R$ stands for picked randomly uniformly in) and $g_i = x_i^{\deg f_i} - 1$. Names correspond to the following: *dense*: $f$ is a random dense polynomial, meaning that $f = \sum_{i=0}^d a_i x^i$ where $a_i \in_R [0, 1] + i[0, 1]$; *katsura*: classical benchmark system, generated in SageMath using `sage.rings.ideal.Katsura`; *resultants*: we pick two random dense polynomials $h_1, h_2 \in \mathbb{C}[x_1, \dots, x_n, y]$, set $f_1$ to $\text{Res}_y(h_1, h_2)$, and $f_2, \dots, f_n$ are random dense polynomials; *structured*: $f_i = \pm 1 + \sum_{k=1}^{10} \left( \sum_{l=1}^3 a_{ikl} x_l \right)^d$, kept factorized, where $a_{ikl} \in_R [-1, 1]$; $W_d$: $f$ is the degree $d$ Wilkinson polynomial $W_d = \prod_{i=1}^d (x - i)$ in expanded form. For $W_{20} + 1$, we keep $W_{20}$ factorized. The last three lines are one-dimensional homotopies where $F_{1/2}$ has clustered roots so that path jumping is likely to happen. For a precise construction, we refer to [16, section 5.4] (the parameters we used are $(n_c, CS, \alpha) = (5, 5, 10)$ and $(10, 5, 10)$).</caption><thead><tr><th rowspan="2">name</th><th rowspan="2">dim.</th><th rowspan="2">deg.</th><th colspan="4">HomotopyContinuation.jl</th><th colspan="4">Algpath</th></tr><tr><th>time (s)</th><th>fail.</th><th>jump.</th><th>med.</th><th>max.</th><th>time (s)</th><th>prec.</th><th>med.</th><th>max.</th></tr></thead><tbody><tr><td>dense</td><td>1</td><td>1000</td><td>11</td><td></td><td></td><td>15</td><td>100</td><td>27 min</td><td>59</td><td>76</td><td>24 k</td></tr><tr><td>dense</td><td>1</td><td>2000</td><td>48</td><td>3</td><td></td><td>24</td><td>79</td><td>2 h</td><td>62</td><td>95</td><td>71 k</td></tr><tr><td>katsura</td><td>21</td><td>2<sup>20</sup></td><td>5 h</td><td></td><td></td><td>209</td><td>468</td><td>89 h</td><td>65</td><td>427</td><td>13 k</td></tr><tr><td>resultants</td><td>3</td><td>16 · 4<sup>2</sup></td><td>7.6</td><td></td><td>×</td><td>54</td><td>128</td><td>160</td><td>58</td><td>440</td><td>3510</td></tr><tr><td>resultants</td><td>2</td><td>40 · 5</td><td>4.2</td><td>200</td><td></td><td></td><td></td><td>13 min</td><td>69</td><td>458</td><td>2205</td></tr><tr><td>structured *</td><td>3</td><td>20<sup>3</sup></td><td>3.7</td><td>12</td><td>×</td><td>107</td><td>164</td><td>6.3</td><td>56</td><td>271</td><td>634</td></tr><tr><td>structured *</td><td>3</td><td>30<sup>3</sup></td><td>3.5</td><td>92</td><td></td><td>111</td><td>133</td><td>97</td><td>71</td><td>383</td><td>820</td></tr><tr><td>W<sub>19</sub></td><td>1</td><td>19</td><td>2.0</td><td>17</td><td></td><td>97</td><td>101</td><td>1 h</td><td>111</td><td>249 k</td><td>1006 k</td></tr><tr><td>W<sub>20</sub></td><td>1</td><td>20</td><td>2.4</td><td>20</td><td></td><td></td><td></td><td>3 h</td><td>118</td><td>559 k</td><td>2316 k</td></tr><tr><td>W<sub>20</sub> + 1</td><td>1</td><td>20</td><td>2.7</td><td>20</td><td></td><td></td><td></td><td>0.6</td><td>63</td><td>309</td><td>402</td></tr><tr><td>clustered (5, 5, 10)</td><td>1</td><td>25</td><td>3.6</td><td></td><td>×</td><td>80</td><td>151</td><td>0.4</td><td>53</td><td>261</td><td>834</td></tr><tr><td>clustered (10, 5, 10)</td><td>1</td><td>50</td><td>3.8</td><td></td><td>×</td><td>80</td><td>153</td><td>153</td><td>70</td><td>992</td><td>17 k</td></tr></tbody></table>

# 4 Implementation and benchmarks

*Implementation details.* Many parts are written generically not only with respect to real intervals, but also to complex intervals, interval matrices and Taylor models, as for example when using Arb we prefer using `acb mat t` (a type provided by Arb for matrices with complex interval coefficients) instead of generic matrices over Arb complex intervals. Rust offers interface-based generics, so handling many type parameters requires a lot of work on interfaces. Ultimately, the setup required when generically implementing high level function tends to be very heavy just for a few lines of code. As of now only two arithmetics are usable, so duplicating the code would have been reasonable, though the generic approach allows for an easier addition of new interval arithmetics.

*Benchmarks.* We compare *Algpath* (code available at https://gitlab.inria.fr/numag/algpath) against *HomotopyContinuation.jl* [3], a software for uncertified homotopy continuation. We used an Intel Xeon E3-1220v3 CPU and 16GB of RAM (on one thread, even though the implementation is multi-threaded). The raw results are available at https://gitlab.inria.fr/numag/algpath-bench. On some examples (e.g. $W_{19}$, $W_{20}$) the number of iterations required by *Algpath* is particularly high, and while HC.jl cannot solve those problems, more robust methods can suffice with much better performances (for example, PHCpack solves $W_{19}$ performing at most 83 steps [16]).
References

[1] D. J. Bates et al. *Numerically Solving Polynomial Systems with Bertini*. Software, Environments, and Tools. Society for Industrial and Applied Mathematics, 2013. 372 pp.

[2] C. Beltrán and A. Leykin. "Certified Numerical Homotopy Tracking". In: *Experimental Mathematics* 21.1 (2012), pp. 69–83.

[3] P. Breiding and S. Timme. "HomotopyContinuation.Jl: A Package for Homotopy Continuation in Julia". In: *Mathematical Software – ICMS 2018*. Ed. by J. H. Davenport et al. Lecture Notes in Computer Science. Springer International Publishing, 2018, pp. 458–465.

[4] M. Burr, M. Byrd, and K. Lee. "Certified Algebraic Curve Projections by Path Tracking". In: (2025). Pre-published.

[5] T. Chen, T.-L. Lee, and T.-Y. Li. "Hom4PS-3: A Parallel Numerical Solver for Systems of Polynomial Equations Based on Polyhedral Homotopy Continuation Methods". In: *Mathematical Software – ICMS 2014*. Ed. by H. Hong and C. Yap. Springer, 2014, pp. 183–190.

[6] T. Duff and K. Lee. "Certified Homotopy Tracking Using the Krawczyk Method". In: *Proc. ISSAC 2024*. Association for Computing Machinery, 2024, pp. 274–282.

[7] A. Guillemot and P. Lairez. "Validated Numerics for Algebraic Path Tracking". In: *Proc. ISSAC 2024*. Association for Computing Machinery, 2024, pp. 36–45.

[8] J. D. Hauenstein, J. I. Rodriguez, and F. Sottile. "Numerical Computation of Galois Groups". In: *Found Comput Math* 18.4 (2018), pp. 867–890.

[9] J. van der Hoeven. "Reliable Homotopy Continuation". In: (2015).

[10] F. Johansson. "Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic". In: *IEEE Transactions on Computers* 66.8 (2017), pp. 1281–1292.

[11] R. B. Kearfott and Z. Xing. "An Interval Step Control for Continuation Methods". In: *SIAM J. Numer. Anal.* 31.3 (1994), pp. 892–914.

[12] B. Lambov. "Interval Arithmetic Using SSE-2". In: *Reliab. Implement. Real Number Algorithms*. Ed. by P. Hertling et al. Lecture Notes in Computer Science. Springer, 2008, pp. 102–113.

[13] M. Á. Marco-Buzunariz and M. Rodríguez. "SIROCCO: A Library for Certified Polynomial Root Continuation". In: *Mathematical Software – ICMS 2016*. Ed. by G.-M. Greuel et al. Vol. 9725. Springer International Publishing, 2016, pp. 191–197.

[14] N. Revol and F. Rouillier. "MPFI". In: (1999).

[15] J. I. Rodriguez and B. Wang. "Numerical Computation of Braid Groups". In: (2017). Pre-published.

[16] S. Telen, M. V. Barel, and J. Verschelde. "A Robust Numerical Path Tracking Algorithm for Polynomial Homotopy Continuation". In: *SIAM J. on Sci. Comput.* 42.6 (2020), A3610–A3637.

[17] S. Timme. "Mixed Precision Path Tracking for Polynomial Homotopy Continuation". In: *Adv Comput Math* 47.5 (2021), p. 75.

[18] J. Verschelde. "Algorithm 795: PHCpack: A General-Purpose Solver for Polynomial Systems by Homotopy Continuation". In: *ACM Trans. Math. Softw.* 25.2 (1999), pp. 251–276.

[19] J. Xu, M. Burr, and C. Yap. "An Approach for Certifying Homotopy Continuation Paths: Univariate Case". In: *Proc. ISSAC 2018*. Association for Computing Machinery, 2018, pp. 399–406.