Algorithm

Sparse interpolation

Cost: With a known support of $t$ monomials, $t$ evaluations and $O(t^2)$ field operations to solve for the coefficients; with an unknown support, $2t$ evaluations, $O(t^2)$ operations for Berlekamp–Massey, root finding of a degree- $t$ polynomial and $t$ discrete logarithms.

The sparse GCD and factoring algorithms of FLINT reduce a problem in many variables to many problems in one or two variables, whose solutions are the values of the unknown polynomial (a GCD, a cofactor, a lifted factor) at chosen points. Sparse interpolation reconstructs a polynomial with $t$ terms from $O(t)$ such values, independently of the number of variables and of the degrees, which would determine the cost of dense interpolation. The two classical approaches are Zippel's, which needs the support (set of monomials) of the result and finds it variable by variable [Zip1979; Zip1981], and Ben-Or and Tiwari's, which finds support and coefficients together from twice as many values [BenOrTiw1988; KalLak1988]. Both work over a finite field $\mathbb{F}_p$ with $p$ large enough for the points chosen to be good with high probability; over $\mathbb{Z}$ the coefficients are then lifted or combined over several primes.

All methods need strict degree bounds $\deg_{x_j} f < r_j$ . The number of values of $f$ (probes) they use, for $f$ with $t$ terms in $m$ variables, compares as follows; the presentation below, in particular of the Berlekamp–Massey algorithm as a half-gcd computation, follows Schultz's notes [Sch-mpoly].

Method	Probes	Main difficulties
dense, variable by variable (Newton)	$\prod_j r_j$	enough evaluation points
Zippel (Zippel's variable-by-variable method)	about $t \sum_j r_j$	enough points; inconsistent or underdetermined systems when the shape assumption fails; the linear algebra
Ben-Or–Tiwari (Unknown support: Ben-Or–Tiwari)	$2t$	an element of order at least $\prod_j r_j$ ; the linear algebra, root finding and discrete logarithms

For Zippel's method the count comes from adding one variable $x_j$ at a time: each of the $r_j$ images needed for the new variable is found from as many probes as the current support has terms, at most $t$ .

Throughout, $f = \sum_{i=1}^{t} c_i x^{e_i} \in \mathbb{F}_p[x_1, \ldots, x_m]$ with distinct $e_i$ and nonzero $c_i$ . For a point $\alpha \in (\mathbb{F}_p^*)^m$ write $\mu_i = \alpha^{e_i}$ for the value of the $i$ -th monomial.

Known support: transposed Vandermonde systems#

If the $e_i$ are known, evaluating $f$ at the powers $\alpha, \alpha^2, \ldots, \alpha^t$ (componentwise) gives

v_k = f(\alpha^k) = \sum_{i=1}^{t} c_i \mu_i^k, \qquad k = 1, \ldots, t,

a linear system for the $c_i$ whose matrix $(\mu_i^k)$ is a transposed Vandermonde matrix scaled by the $\mu_i$ . It is nonsingular if and only if the $\mu_i$ are distinct, and it can be solved in quadratic time.

Lemma 1 (Transposed Vandermonde solution). Let $M(z) = \prod_{i=1}^{t} (z - \mu_i)$ with distinct nonzero $\mu_i$ , and let $q_i(z) = M(z)/(z - \mu_i) = \sum_{k=0}^{t-1} q_{i,k} z^k$ . Then

c_i = \frac{\sum_{k=0}^{t-1} q_{i,k} \, v_{k+1}}{\mu_i \, q_i(\mu_i)} .

All $c_i$ can be computed with $O(t^2)$ operations from the $v_k$ and the coefficients of $M$ .

Proof. $\sum_k q_{i,k} v_{k+1} = \sum_j c_j \mu_j \sum_k q_{i,k} \mu_j^k = \sum_j c_j \mu_j q_i(\mu_j)$ , and $q_i(\mu_j) = 0$ for $j \ne i$ , while $q_i(\mu_i) = M'(\mu_i) \ne 0$ . Given $M$ , the coefficients of $q_i$ follow by synthetic division in $O(t)$ operations, and the numerator and denominator are accumulated during that division.

_nmod_zip_vand_solve implements this. With more than $t$ values it also checks the remaining equations, which detects a wrong support: it returns $-1$ if two $\mu_i$ coincide (choose another point) and $0$ if the system is inconsistent (the assumed support was wrong).

Zippel's variable-by-variable method#

Zippel's algorithm applies this inside a recursion on the variables. Suppose $f$ has been interpolated in $x_1, \ldots, x_{j}$ with the remaining variables at random values; write it as $\sum_\ell f_\ell(x_1, \ldots, x_j) \, x_{j+1}^{\ell}$ , the shape. To add $x_{j+1}$ , Zippel assumes that each coefficient of $f$ as a polynomial in $x_{j+1}, \ldots$ has the same monomials in $x_1, \ldots, x_j$ as in the image; new images at other values of $x_{j+1}$ are then found by Lemma 1 instead of a full recursive computation, and the dependence on $x_{j+1}$ is recovered by dense univariate interpolation.

Lemma 2 (Probability of a wrong shape). Let $f \in \mathbb{F}_p[x_1, \ldots, x_m]$ and let $f = \sum_{u} f_u(x_{j+1}, \ldots, x_m) \, x^{u}$ over monomials $x^u$ in $x_1, \ldots, x_j$ . If the values $\beta$ of $x_{j+1}, \ldots, x_m$ are chosen uniformly at random from $S \subseteq \mathbb{F}_p$ , the image $f(x_1, \ldots, x_j, \beta)$ has fewer monomials than $f$ has distinct $x^u$ with probability at most $t \cdot \deg f / |S|$ .

Proof. The image lacks $x^u$ exactly when $f_u(\beta) = 0$ ; each $f_u$ is a nonzero polynomial of degree at most $\deg f$ , so this has probability at most $\deg f / |S|$ by the Schwartz–Zippel lemma (Lemma 1 of GCDs over small finite fields), and there are at most $t$ of them.

The method needs $O(t d)$ values per variable for degree $d$ , against $d^m$ for dense interpolation. A wrong shape is detected later, because the linear systems become inconsistent or the final verification (a division) fails, and the computation is restarted with new points.

Scaled images#

In GCD computations the values are not values of $f$ itself: an image $\gcd(A(x_0, \alpha^k), B(x_0, \alpha^k))$ is only determined up to a nonzero scalar $s_k$ unless the leading coefficient of the GCD is known. If the GCD is $\sum_\ell x_0^\ell \sum_i c_{\ell i} x^{e_{\ell i}}$ and the images are $\sum_\ell x_0^\ell g_{\ell k}$ , then

\sum_i c_{\ell i} \mu_{\ell i}^{k} = s_k \, g_{\ell k} \quad \text{for all } \ell, k,

a linear system in the $c_{\ell i}$ and the $s_k$ , homogeneous, and normalised by fixing one coefficient. De Kleine, Monagan and Wittkopf [dKMW2005] solve it by eliminating the $c_{\ell i}$ with the structure of Lemma 1 for each $\ell$ , which leaves a small system for the scalars. nmod_mpolyl_gcds_zippel does this (with a shortcut when the leading coefficient is a monomial, which makes the images monic).

Unknown support: Ben-Or–Tiwari#

Theorem 3 (Ben-Or–Tiwari). With $v_k = f(\alpha^k) = \sum_{i=1}^{t} c_i \mu_i^k$ for $k \ge 0$ and distinct $\mu_i$ , the sequence $(v_k)$ satisfies the linear recurrence with characteristic polynomial $\Lambda(z) = \prod_{i=1}^t (z - \mu_i) = z^t + \lambda_{t-1} z^{t-1} + \cdots + \lambda_0$ , that is, $v_{k+t} + \sum_{j<t} \lambda_j v_{k+j} = 0$ for all $k \ge 0$ , and no recurrence of smaller order. $\Lambda$ is determined by $v_0, \ldots, v_{2t-1}$ and can be computed from them by the Berlekamp–Massey algorithm in $O(t^2)$ operations.

Proof. $v_{k+t} + \sum_j \lambda_j v_{k+j} = \sum_i c_i \mu_i^k \Lambda(\mu_i) = 0$ . A recurrence of order $s < t$ would give $\sum_i c_i \mu_i^k P(\mu_i) = 0$ for a nonzero $P$ of degree $s$ and all $k$ , that is, a vanishing combination of the geometric sequences $(\mu_i^k)$ with coefficients $c_i P(\mu_i)$ , not all zero since $P$ has fewer than $t$ roots; this contradicts the independence of geometric sequences with distinct ratios (a Vandermonde determinant). The minimal polynomial of a linearly recurrent sequence of order $t$ is determined by its first $2t$ terms, and the Berlekamp–Massey algorithm [Mas1969; ADL2006] finds it with $O(t^2)$ operations.

The roots of $\Lambda$ are the monomial values $\mu_i$ , from which the exponents must be recovered, and then the coefficients follow by Lemma 1. Two devices make the recovery cheap:

Kronecker substitution: with degree bounds $d_j > \deg_{x_j} f$ , the substitution $x_j \mapsto z^{D_j}$ , $D_j = \prod_{l > j} d_l$ , maps the monomials of $f$ injectively to powers of one variable $z$ , with exponents $E_i < D = \prod_j d_j$ (Lemma 4 of Sparse and dense multiplication).
Discrete logarithms: choose a prime $p > D$ for which $p - 1$ has only small prime factors, and a generator $\omega$ of $\mathbb{F}_p^*$ . Evaluating at $z = \omega^k$ gives $\mu_i = \omega^{E_i}$ , the $\mu_i$ are distinct because $E_i < D < p - 1$ , and $E_i$ is the discrete logarithm of $\mu_i$ , computed by the Pohlig–Hellman algorithm [PohHel1978] in time polynomial in the largest prime factor of $p - 1$ .

In practice the number of terms $t$ is unknown: the Berlekamp–Massey algorithm is run incrementally as values arrive, and the candidate $\Lambda$ is accepted when it stops changing, its roots are all in $\mathbb{F}_p$ and distinct, the discrete logarithms are below $D$ , and the resulting polynomial passes a check at a random point. Compared with Zippel's method this needs about $2t$ values in total instead of $O(t)$ per variable, but the prime must be larger than $D$ , which grows with the degrees and the number of variables.

fmpz_mpolyl_gcd_zippel2 uses this for the coefficients of a GCD with respect to its two main variables (Sparse interpolation with Berlekamp–Massey): fmpz_next_smooth_prime chooses the prime, nmod_berlekamp_massey_t (or fmpz_mod_berlekamp_massey_t for primes beyond a word) runs the recurrence, the roots are found with nmod_poly root finding, and nmod_discrete_log_pohlig_hellman_run gives the exponents. Images modulo further primes reuse the support and only solve transposed Vandermonde systems.

Berlekamp–Massey as a partial Euclidean algorithm#

FLINT computes the recurrence of Theorem 3 by rational reconstruction with the extended Euclidean algorithm, in a form that accepts the values in batches and reuses the previous work, with the half-gcd for large steps [ThullYap1990]. For values $a_1, a_2, \ldots \in K$ ( $K$ a field) write

S_n = a_1 z^{n-1} + a_2 z^{n-2} + \cdots + a_n \in K[z], \qquad \sum_{i \ge 1} a_i z^{-i} \in K[[z^{-1}]],

so that $S_n / z^n$ is the truncation of the series. If $a_i = f(\alpha^{s + i - 1})$ with $f$ , $\alpha$ , $\mu_k$ as above and a shift $s \in \mathbb{Z}$ , then

\sum_{i \ge 1} a_i z^{-i} = \sum_{k=1}^{t} \frac{c_k \mu_k^{s}}{z - \mu_k},

a rational function whose reduced denominator is $\Lambda = \prod_k (z - \mu_k)$ when the $\mu_k$ are distinct, since every numerator $c_k \mu_k^s$ is nonzero. Finding $\Lambda$ is therefore a rational reconstruction problem.

Theorem 4 (Rational reconstruction). Suppose that $\sum_{i \ge 1} a_i z^{-i} = -\bar u / \bar v$ with $\deg \bar u < \deg \bar v \le n/2$ , and that $u, v, r \in K[z]$ satisfy

u z^n + v S_n = r, \qquad \deg u < \deg v \le n/2, \qquad \deg r < n/2, \qquad \deg r < \deg v .

Then $u / v = \bar u / \bar v$ .

Proof. Dividing by $v z^n$ gives $u/v + S_n/z^n = r/(v z^n)$ . Here $S_n / z^n = -\bar u/\bar v + O(z^{-n-1})$ , and $r/(v z^n) = O(z^{-n-1})$ because $\deg r < \deg v$ . Hence $u/v - \bar u/\bar v = (u \bar v - \bar u v)/(v \bar v)$ is $O(z^{-n-1})$ . A nonzero numerator would make it at least of order $z^{-\deg v - \deg \bar v}$ , and $\deg v + \deg \bar v \le n$ , so the numerator is zero.

The Euclidean remainder sequence $r_0 = z^n$ , $r_1 = S_n$ , $r_2, \ldots$ of Half-GCD comes with Bézout coefficients $u_j z^n + v_j S_n = r_j$ , $\deg v_j = n - \deg r_{j-1}$ . At the index $j$ of the half-gcd, $2 \deg r_{j-1} \ge n > 2 \deg r_j$ , so $\deg v_j \le n/2$ and $\deg r_j < n/2$ : if also $\deg r_j < \deg v_j$ , the theorem applies to $(u_j, v_j, r_j)$ . Conversely, if the series has a denominator of degree $t \le n/2$ , then $\bar u z^n + \bar v S_n$ has degree less than $t$ , and by the uniqueness of rational reconstruction [GG2003, §5.7] $(\bar u, \bar v)$ is a multiple of the row $(u_j, v_j)$ , a scalar multiple because $\bar u$ and $\bar v$ are coprime. So $2t$ values determine $\Lambda$ , as the half-gcd row $v_j$ made monic.

Since $t$ is not known, the computation is repeated as values are added, and $v_j$ is accepted once further values no longer change it. The next theorem shows how the half-gcd for $n$ values is extended to $n + m$ values. We use the notation of Half-GCD: $(a, b)^T = M (A, B)^T$ with $M$ a product of quotient matrices, and $\deg a$ the degree of the first entry.

Theorem 5 (Extending the half-gcd). Let $M$ be the matrix of the half-gcd of $(z^n, S_n)$ , with remainders $(R_0, R_1)$ , $2 \deg R_0 \ge n > 2 \deg R_1$ , and let $S_{n+m} = z^m S_n + s$ with $\deg s < m$ . Define

\begin{pmatrix} A \\ B \end{pmatrix} = M^{-1} \begin{pmatrix} z^{n+m} \\ S_{n+m} \end{pmatrix} = z^m \begin{pmatrix} R_0 \\ R_1 \end{pmatrix} + M^{-1} \begin{pmatrix} 0 \\ s \end{pmatrix}.

Then $A$ , $B$ are consecutive remainders of $(z^{n+m}, S_{n+m})$ with $2 \deg A \ge n + m$ and $\deg B < m + n/2$ . If $2 \deg B < n + m$ , then $M$ is the matrix of the half-gcd of $(z^{n+m}, S_{n+m})$ . Otherwise let $C = B$ , $D = A \bmod B$ , $q = A \operatorname{quo} B$ , $l = \deg C$ and $k = n + m - l$ . Then $0 < k \le l$ , and if $H$ is the matrix of the half-gcd of $(C \operatorname{quo} z^k, D \operatorname{quo} z^k)$ , the matrix of the half-gcd of $(z^{n+m}, S_{n+m})$ is $M Q H$ with $Q = \bigl( \begin{smallmatrix} q & 1 \\ 1 & 0 \end{smallmatrix} \bigr)$ .

Proof. The entries of $M$ , and hence of $M^{-1} = \pm \operatorname{adj} M$ , have degree at most $n - \deg R_0 \le n/2$ . So the second term has degree less than $m + n/2$ , while $z^m R_0$ has degree $m + \deg R_0 \ge m + n/2$ : thus $\deg A = m + \deg R_0 \ge (n + m)/2$ , and $\deg B < m + n/2$ because $z^m R_1$ also has degree less than $m + n/2$ . As $\deg A > \deg B$ , Lemma 1 of Half-GCD shows that $A, B$ are consecutive remainders of $(z^{n+m}, S_{n+m})$ (this is the lifting lemma Lemma 3 of Half-GCD with $a_0 = 0$ ). If $2 \deg B < n + m$ they are the half-gcd remainders. Otherwise $2l \ge n + m$ gives $k \le l$ , and $l < m + n/2 \le n + m$ gives $k > 0$ . The pair $(C, D)$ consists of the next consecutive remainders. By the lifting lemma applied to $(C, D)$ split at $z^k$ , the remainders $(C', D')^T = H^{-1}(C, D)^T$ are consecutive with $\deg C' = k + \deg C^* \ge k + (l - k)/2 = (n + m)/2$ , where $C^*$ is the first half-gcd remainder of $(C \operatorname{quo} z^k, D \operatorname{quo} z^k)$ , whose first entry has degree $l - k$ . The bound on the degrees of the entries of $H^{-1}$ used above gives, in the same way, $\deg D' < k + (l - k)/2 = (n + m)/2$ . So $(C', D')$ are the half-gcd remainders of $(z^{n+m}, S_{n+m})$ .

The algorithm keeps the remainders $R_0, R_1$ and the second Bézout coefficients $V_0, V_1$ (the first ones are not needed), with $U_i z^n + V_i S_n = R_i$ ; multiplying $(R_0, R_1)$ and $(V_0, V_1)$ by the same matrix on the left preserves these relations.

Algorithm 1: Incremental Berlekamp–Massey

State: $n$ ; $R_0, R_1, V_0, V_1$ as above, initially $n = 0$ , $R_0 = 1$ , $R_1 = 0$ , $V_0 = 0$ , $V_1 = 1$ . Input: new values $a_{n+1}, \ldots, a_{n+m}$ . Output: the updated state; $V_1$ is the candidate denominator.

Let $s = \sum_{i=1}^{m} a_{n+i} z^{m-i}$ . Set $R_i \gets z^m R_i + V_i s$ for $i = 0, 1$ and $n \gets n + m$ .
If $2 \deg R_1 < n$ , return: $V_1$ has not changed.
Do one Euclidean step: $q = R_0 \operatorname{quo} R_1$ , $(R_0, R_1) \gets (R_1, R_0 - q R_1)$ , $(V_0, V_1) \gets (V_1, V_0 - q V_1)$ .
Let $l = \deg R_0$ and $k = n - l$ . If $l - k$ is small, do Euclidean steps until $2 \deg R_1 < n$ . Otherwise compute the half-gcd matrix $H$ of $(R_0 \operatorname{quo} z^k, R_1 \operatorname{quo} z^k)$ and replace $(R_0, R_1)^T$ and $(V_0, V_1)^T$ by $H^{-1} (R_0, R_1)^T$ and $H^{-1}(V_0, V_1)^T$ .

By Theorem 5 the state after step 4 is the half-gcd of $(z^n, S_n)$ with its Bézout coefficients, and by Theorem 4 $V_1$ is the denominator of the series as soon as $n \ge 2t$ . When the new values do not change the reconstruction, only step 1 is executed, at the cost of two products of polynomials of degrees at most $n/2$ and less than $m$ . When they do, the half-gcd works on polynomials of degree $l - k = 2l - n$ , about twice the number of Euclidean steps still needed, rather than on polynomials of degree $n$ .

nmod_berlekamp_massey_reduce implements Algorithm 1, using Euclidean steps when $l - k < 10$ and nmod_poly_hgcd otherwise; it returns whether $V_1$ changed. fmpz_mod_berlekamp_massey_reduce has the same structure but always uses Euclidean steps. In fmpz_mpolyl_gcd_zippel2 each coefficient of the bivariate GCD has its own instance; the values $a_i$ are the images at the successive powers $\omega^{s}, \omega^{s+1}, \ldots$ of the evaluation point, the reduction is run after every eighth value, and a reduction that leaves a set of candidates unchanged triggers the root finding of Unknown support: Ben-Or–Tiwari and the verification. The coefficients $c_k$ then follow from the roots and the first $t$ values by Lemma 1, with the shift $\mu_k^{s}$ divided out.

Sparse interpolation in Hensel lifting#

In sparse Hensel lifting [Zip1981; Kal1985; MonTun2020] the unknowns are the coefficients of the lifted factors as polynomials in the new variable $x_j$ : each coefficient of $x_j^\ell$ is assumed to have the same support in the previous variables as the factor itself, and its values are obtained from bivariate factorisations at the points $(\beta_2^k, \ldots, \beta_{j-1}^k)$ for $k = 1, 2, \ldots$ . The values of all factors at the same points are only determined up to scalars whose product is known, and the leading coefficients are imposed as in the GCD case. nmod_mpoly_hlift_zippel (over $\mathbb{F}_p$ ) and the $p$ -adic lifting inside fmpz_mpoly_factor_irred_zippel work this way; see Sparse Hensel lifting.

Use in FLINT#

Transposed Vandermonde systems (Lemma 1): _nmod_zip_vand_solve, _n_fq_zip_vand_solve, _n_fqp_zip_vand_solve, with the evaluation helpers _nmod_zip_eval_step and n_polyun_zip_start, n_polyun_zip_solve.
Zippel's GCD with scaled images (Scaled images): nmod_mpolyl_gcds_zippel, nmod_mpolyl_gcdp_zippel_smprime, used by fmpz_mpolyl_gcd_zippel and nmod_mpoly_gcd_zippel.
Ben-Or–Tiwari (Theorem 3): fmpz_mpolyl_gcd_zippel2, with nmod_berlekamp_massey_t, fmpz_mod_berlekamp_massey_t, nmod_discrete_log_pohlig_hellman_t and fmpz_mod_discrete_log_pohlig_hellman_t.
Incremental Berlekamp–Massey (Algorithm 1): nmod_berlekamp_massey_reduce (with nmod_poly_hgcd), fmpz_mod_berlekamp_massey_reduce.
Sparse Hensel lifting: nmod_mpoly_hlift_zippel, fmpz_mpoly_factor_irred_zippel, nmod_mpoly_factor_irred_smprime_zippel.

int fmpz_mpolyl_gcd_zippel(fmpz_mpoly_t G, fmpz_mpoly_t Abar, fmpz_mpoly_t Bbar, const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx, flint_rand_t randstate)

int fmpz_mpolyl_gcd_zippel2(fmpz_mpoly_t G, fmpz_mpoly_t Abar, fmpz_mpoly_t Bbar, const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_t Gamma, const fmpz_mpoly_ctx_t ctx)

int nmod_mpolyl_gcds_zippel(nmod_mpoly_t G, const ulong * Gmarks, slong Gmarkslen, nmod_mpoly_t A, nmod_mpoly_t B, slong *perm, slong l, slong var, const nmod_mpoly_ctx_t ctx, flint_rand_t state, slong * Gdegbound, n_poly_t Amarks, n_poly_t Bmarks)

int nmod_mpoly_hlift_zippel(slong m, nmod_mpoly_struct * B, slong r, const ulong * alpha, const nmod_mpoly_t A, const slong * degs, const nmod_mpoly_ctx_t ctx, flint_rand_t state)

int _nmod_zip_vand_solve(ulong * coeffs, const ulong * monomials, slong mlength, const ulong * evals, slong elength, const ulong * master, ulong * scratch, nmod_t ctx)

Referenced from#

Greatest common divisors (algorithm)
Multivariate Hensel lifting (algorithm)
GCDs over small finite fields (algorithm)
Multivariate polynomials: representations and algorithms (overview)
nmod_berlekamp_massey_reduce

References#

[Zip1979] R. Zippel. Probabilistic algorithms for sparse polynomials. Symbolic and Algebraic Computation (EUROSAM 1979), Lecture Notes in Computer Science 72, Springer, 1979, 216–226. https://doi.org/10.1007/3-540-09519-5_73
[Zip1981] R. Zippel. Newton's iteration and the sparse Hensel algorithm. Proceedings of SYMSAC 1981, ACM, 68–72. https://doi.org/10.1145/800206.806373
[BenOrTiw1988] M. Ben-Or and P. Tiwari. A deterministic algorithm for sparse multivariate polynomial interpolation. Proceedings of STOC 1988, ACM, 301–309. https://doi.org/10.1145/62212.62241
[KalLak1988] E. Kaltofen and Y. N. Lakshman. Improved sparse multivariate polynomial interpolation algorithms. Symbolic and Algebraic Computation (ISSAC '88), Lecture Notes in Computer Science 358, Springer, 1989, 467–474. https://doi.org/10.1007/3-540-51084-2_44
[Sch-mpoly] D. Schultz. Algorithms for multivariate polynomials. Notes in the FLINT source tree, src/mpoly/doc/MPolyAlgorithms.tex.
[dKMW2005] J. de Kleine, M. Monagan and A. Wittkopf. Algorithms for the non-monic case of the sparse modular GCD algorithm. Proceedings of ISSAC 2005, ACM, 124–131. https://doi.org/10.1145/1073884.1073903
[Mas1969] J. L. Massey. Shift-register synthesis and BCH decoding. IEEE Transactions on Information Theory 15 (1969), 122–127.
[ADL2006] N. B. Atti, G. M. Diaz-Toca and H. Lombardi. The Berlekamp–Massey algorithm revisited. Applicable Algebra in Engineering, Communication and Computing 17 (2006), 75–82. https://doi.org/10.1007/s00200-005-0190-z
[PohHel1978] S. C. Pohlig and M. E. Hellman, "An improved algorithm for computing logarithms over GF(p) and its cryptographic significance", IEEE Transactions on Information Theory 24 (1978), 106–110, https://doi.org/10.1109/TIT.1978.1055817
[ThullYap1990] Thull, K. and Yap, C. : A Unified Approach to HGCD Algorithms for Polynomials and Integers, (1990)
[GG2003] J. von zur Gathen and J. Gerhard, Modern Computer Algebra, second edition, Cambridge University Press (2003)
[Kal1985] E. Kaltofen. Sparse Hensel lifting. EUROCAL '85, Lecture Notes in Computer Science 204, Springer, 1985, 4–17. https://doi.org/10.1007/3-540-15984-3_230
[MonTun2020] M. Monagan and B. Tuncer. The complexity of sparse Hensel lifting and sparse polynomial factorization. Journal of Symbolic Computation 99 (2020), 189–230. https://doi.org/10.1016/j.jsc.2019.05.001

Cite this page

The FLINT team, “Sparse interpolation,” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/algorithms/multivariate/sparse-interpolation.html (generated 2026-11-02).

@misc{flint-doc-algorithms-multivariate-sparse-interpolation,
  author       = {{The FLINT team}},
  title        = {{Sparse interpolation}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/multivariate/sparse-interpolation.html},
  note         = {Generated 2026-11-02}
}

Edit this page · PDF version · Built with fdoc from 0a4634bf95