LACE for Maple

The Maple package LACE computes, with certified results:

• the real roots of univariate polynomials, isolated in intervals, with constraints or through a parameterization;
• the Rational Univariate Representation (RUR) and the real solutions of zero-dimensional polynomial systems;
• the resultant, the decomposition and the RUR of bivariate systems, and points of plane curves;
• interval Newton iterations and Kantorovich tests on square systems;
• for discriminant varieties of parametric systems, over Q: elimination ideals, ideals at infinity and saturations by a polynomial, Groebner bases included.

Its functions are exported by the module LACE, which is also available under the names RealSolving and RS4, so that code written for those packages runs unchanged. Maple’s own fgbrs package, used by RootFinding, is left untouched. Every function of the Julia package LACE.jl is there too, with the same name, arguments, options and results (see The interface of LACE.jl).

Installation

Download the archive for Maple 2026 and your system: LACE-Maple2026-Darwin_arm64.tar.gz (macOS 12 or later, Apple silicon), LACE-Maple2026-Linux_x86_64.tar.gz (Linux x86-64) or LACE-Maple2026-Windows_x86_64.zip (Windows x86-64); other versions of Maple are not available yet. It holds a directory LACE with README.md and lib/ (LACE.mla, the Maple code and its help pages LACE.help, and liblace_mpl.so, the compiled library). It needs nothing but Maple: on macOS and Linux it uses the GMP, MPFR and FLINT of Maple and includes MPFI (LGPL); on Windows GMP, MPFR, MPFI and FLINT (LGPL) are included.

Quit Maple, then install the package as a Maple toolbox, in a terminal, from the directory that holds the archive:

mkdir -p ~/maple/toolbox/2026
tar -C ~/maple/toolbox/2026 -xzf LACE-Maple2026-Darwin_arm64.tar.gz
xattr -cr ~/maple/toolbox/2026/LACE

On Linux, the same without xattr. On Windows, extract the zip archive into %USERPROFILE%\maple\toolbox\2026 (the folder LACE directly under it).

xattr removes the quarantine that macOS puts on files received from the Internet or by AirDrop; without it, Maple may be refused the library. Maple adds ~/maple/toolbox/2026/LACE/lib to libname, and reads the help pages, when it starts: start Maple and check:

LACE:-library();     # the file .../maple/toolbox/2026/LACE/lib/liblace_mpl.so
LACE:-backend();     # "native"
?LACE

To update, extract the new archive over the old one, in a closed Maple session. To uninstall, remove the directory ~/maple/toolbox/2026/LACE.

The archive must match the Maple that loads it: an archive for Maple 2026 does not load in another version of Maple.

Troubleshooting

• LACE:-library() returns only a file name, or the first call fails with an error about an external library: the toolbox is not on libname. Check libname in a new session and the location of the directory LACE (the year in the path must be the version of Maple).
• The library is found but refused: remove the quarantine with xattr -cr as above.
• ?LACE finds nothing: check that LACE.help is in ~/maple/toolbox/2026/LACE/lib, then quit Maple entirely and start it again; Maple reads the help databases when it starts.
• The library can be taken elsewhere: LACEBinDir := "/some/dir": before the first call, or the environment variable MAPLELDPATH, take precedence over the toolbox.
• Call LACE:-finish() before quit in scripts: Maple does not always release the objects held by the library at exit.

First steps

f := x^3 - 7*x + 6:
LACE:-rs_isolate_uni_classical(f, precision = 64);         # the roots -3, 1, 2 in intervals
LACE:-rur([x^2 + y^2 - 5, x - y - 1], [x, y]);             # a RUR of the system
LACE:-real_solutions([x^2 + y^2 - 5, x - y - 1], [x, y]);  # its two real solutions
LACE:-elim([a*x^2 + a + 1, 2*a*x], [x], [a]);              # [a + 1]: the critical values (a double root in x)
LACE:-w_infty([a*x^2 + a + 1], [x], [a]);                  # [true, [[a]]]: where a root in x goes to infinity
with(LACE):                                                # the short names

Long computations can be interrupted as usual in Maple.

Conventions

• Polynomials are ordinary Maple polynomials with rational coefficients; each one is scaled to integer coefficients by a positive factor before it is passed to the library (except by the Newton functions).
• An interval is a list [lo, hi] of rationals. Two equal bounds denote an exact value; otherwise the open interval (lo, hi) contains exactly one root, the bounds excluded (a bound may be another root, given as its own exact value). The functions that take intervals (real_solution_at, rur_coords, signs_at) read them the same way.
• Intervals are returned in increasing order of the root (of the eliminant, for systems).
• precision is a number of significant bits (a relative precision), except for isolate_roots and refine_roots, where it is an absolute width, and for the Newton functions, where it is the precision of the interval arithmetic.
• threads is the number of workers; 0, the default, takes the library’s default: ACE_NB_THREADS (read when the library is loaded), else OMP_NUM_THREADS, else the performance cores of the machine. The exact results (RUR, ideals, polynomials) do not depend on it, nor do the primes the lifts over Q use (last_nprimes, last_biv_nprimes; for elim, w_infty, saturate up to a rare change of learn prime). The intervals do: the isolation of a univariate polynomial steals work between the threads, so that its isolating intervals, and the intervals and boxes refined from them, depend on the number of threads and, from two threads on, may differ from one run to the next. What they enclose does not: the number of roots or solutions, their order, the exact roots; refined to the same width, the intervals of the same rank overlap. Use threads = 1 for reproducible intervals.
• Errors of the library are raised as Maple errors.
• The functions of LACE.jl take and return Maple objects: polynomials rather than coefficient vectors, an interval [lo, hi] of rationals rather than a Root or a BigInterval, a list rather than a tuple or a named tuple (its fields in the same order), FAIL for nothing, a name (sparse) for a symbol (:sparse). A Julia name ending in ! loses it, with set_ in front for a setter: refine! is refine, gb_backend! is set_gb_backend. A Julia keyword whose name is also an export of LACE (pmv, refine) is to be quoted after with(LACE): 'pmv' = true.

The functions of LACE.jl and their Maple names:

LACE.jl	Maple
zdim_parameterization(sys; ...)	zdim_parameterization(sys, vars, ...)
rur_modp(sys; p, ...), rur_core(sys)	rur_modp(sys, vars, p, ...), rur_core(sys, vars, p)
real_solutions(...) (four forms)	real_solutions
real_solution_at, rur_coords, signs_at, eval_bfi	the same, see below
guess_infos, guess_lowest_input_precision	the same, with vars
biv_parameterization, biv_decomposition, biv_psc, biv_resultant, pmv	the same, with vars = [x_1, x_2]
biv_decomposition_modp(sys; p, ...)	biv_decomposition_modp(sys, vars, p, ...)
isolate, Root, lo, hi, width, isexact	isolate and the same on [lo, hi]
isolation, refine!, isolated_root, refine_root!, set_root!	isolation, refine, isolated_root, refine_root, set_root
squarefree_part, dyadic_outward, exact_rational	the same
elim(F, vars), w_infty, saturate	the same (params optional)
tci_points, tci_columns, tci_track	the same, with vars
interval_newton, interval_classic_newton, kantorovich_test	the same, with vars
gb_backend, gb_backend!, rur_engine!, gb_time, gb_time_reset!, real_times_reset!	gb_backend, set_gb_backend, set_rur_engine, gb_time, gb_time_reset, real_times_reset
last_search, last_nprimes, last_biv_nprimes, nb_threads, set_verbose	the same
RS.jl: rs_isolate, rs_rur, rs_tci_points, rs_newton	rs_isolate and the same
RS4.jl: bz_biv_*, set_nb_threads, get_nb_threads	the bz_biv_* above (the forms of Maple’s RS4), bz_biv_rur_all, set_nb_threads, get_nb_threads

The random strategies of the separating search draw their forms from a generator of their own: the forms they find may differ from those of Julia, not their validity.

Module functions

function	effect
LACE:-backend()	the arithmetic the library was built with ("native", "flint")
LACE:-library()	the file the external functions are bound to
LACE:-set_verbose(v)	0 or false: silent; 1 or true: a summary of each computation; 2: details
LACE:-reinit()	initialises the library if needed
LACE:-finish()	releases the objects held by the library; call it before quit
LACE:-nb_threads()	the default number of workers
LACE:-set_nb_threads(n)	sets it; n <= 0: the performance cores
LACE:-perf_cores()	the physical performance cores the process may use
LACE:-memory_budget()	the memory, in MB, the multi-modular loops may occupy (0: the available memory)
LACE:-set_memory_budget(mb)	sets it
LACE:-rur_engine()	the engine of the RUR modulo a prime, table (default) or fglm
LACE:-set_rur_engine(e)	sets it
LACE:-real_times()	the seconds of the real solving since the last reset: [square-free part, isolation, coordinates]
LACE:-real_times_reset()	resets them
LACE:-last_nprimes()	the primes of the last rur, real_solutions, elim, w_infty or saturate (-1 if none)
LACE:-last_biv_nprimes()	the primes of the last bivariate decomposition (-1 if none)
LACE:-last_search()	what the last separating search found: [quotient, minpoly, sep_lin, nvars] (all -1 before any)
LACE:-gb_backend()	the engine of the Groebner bases: always lace, the F4 of the library
LACE:-set_gb_backend(b)	lace (or f4c); groebner, Groebner.jl, exists in Julia only
LACE:-gb_time()	[seconds, calls] of the Groebner bases computed modulo a prime since the last reset
LACE:-gb_time_reset()	resets them
LACE:-get_nb_threads()	the same as nb_threads() (the name of RS4.jl)
LACE:-KANTO_HESSIAN, KANTO_PRECOND, KANTO_SCHULZ	1, 2, 4: the flags of the Kantorovich test

Univariate polynomials

rs_isolate_uni_classical

LACE:-rs_isolate_uni_classical(f, precision = 2, intervalLeft = -infinity,
                               intervalRight = infinity, constraints = [],
                               verbose = 0, threads = 0, squarefree = false)

The real roots of the univariate polynomial f, each in an interval of precision significant bits. The square-free part of f is isolated unless squarefree = true says that f is square-free; a multiple root is returned once.

• intervalLeft, intervalRight keep the roots of that range only (a root lying on a bound is returned as the point [q, q]).
• constraints is a list of polynomials in the variable of f. The intervals are then such that each constraint has a constant sign on them, and the result is the pair [roots, values], values[i][j] an interval containing the value of constraint j at root i (the point [0, 0] when the constraint vanishes at that root).

f := x^3 - 7*x + 6:                                   # (x - 1)(x - 2)(x + 3)
LACE:-rs_isolate_uni_classical(f, precision = 64);    # three intervals, around -3, 1, 2
LACE:-rs_isolate_uni_classical(f, intervalLeft = 0, precision = 64);   # the roots 1 and 2
roots, vals := op(LACE:-rs_isolate_uni_classical(f, precision = 64,
                                                 constraints = [x - 1]));
# vals[i][1] encloses the value of x - 1 at root i; it is [0, 0] at x = 1

rs_isolate_rur

LACE:-rs_isolate_rur(rext, rden, [num_1, ..., num_m], precision = 64,
                     exactinterval = FAIL, verbose = 0, threads = 0,
                     squarefree = false)

The real points of a parameterization in one variable T: rext(T) = 0, x_k = num_k(T) / rden(T), where rden is coprime with rext (a constant is allowed). The result has one element per real root of rext, the list of the intervals of x_1, ..., x_m, each with precision significant bits. exactinterval = [lo, hi] keeps only the root of rext that this interval isolates.

LACE:-rs_isolate_rur(T^2 - 2, 1, [T, T^2], precision = 64);
# two points: x_1 around -sqrt(2) and sqrt(2), x_2 an interval containing 2
LACE:-rs_isolate_rur(T^2 - 2, 1, [T, T^2], exactinterval = [1, 2]);
# the point with x_1 = sqrt(2) only

isolate_roots

LACE:-isolate_roots(f, precision = 0, threads = 0, max_roots = 0,
                    newton = true, squarefree = false)

The isolating intervals of the real roots of f, with dyadic bounds num / 2^k, refined to a width below 2^(-precision) when precision > 0. max_roots > 0 stops after that many roots; newton = false turns off the Newton steps of the isolation, which help on clustered roots and may cost on small polynomials.

refine_roots

LACE:-refine_roots(f, roots, precision, threads = 0, squarefree = false)

The intervals roots of f, as isolate_roots returns them, refined to a width below 2^(-precision).

f := expand((x^2 - 2)^2 * (x - 1)):
r := LACE:-isolate_roots(f):                 # three intervals
LACE:-refine_roots(f, r, 100)[3];            # sqrt(2) to 100 bits

squarefree_part

LACE:-squarefree_part(f)

The primitive square-free part of f, with integer coefficients.

LACE:-squarefree_part(expand((x^2 - 2)^2 * (x - 1)));   # x^3 - x^2 - 2*x + 2

Zero-dimensional systems

rur

LACE:-rur(sys, vars, tvar = _T, sep = 0, strategy = 0, nbits = 28,
          maxprimes = 0, threads = 0, multiplicities = false)

The Rational Univariate Representation of a zero-dimensional system sys with rational coefficients in the variables vars. The Groebner bases are computed by the library. The result is a list of components

[f, g, [n_1, ..., n_n], tvar = form, mult]

with vars[k] = n_k(tvar) / g(tvar) at the roots of f(tvar), form the separating linear form in vars and mult the multiplicity. The result is [] when the system has no solution; a system that is not zero-dimensional raises an error.

option	meaning
tvar	the name of the parameter
sep	0: the library chooses the separating element; a name of vars, or its rank, imposes that variable
strategy	the search of the separating element, a name or its number (sparse = 0, the default), see zdim_parameterization
nbits	the primes are taken below 2^nbits
maxprimes	the lifting fails beyond that many primes; 0: 200000
multiplicities	true: the multiplicities of the solutions as well

LACE:-rur([x^2 + y^2 - 5, x - y - 1], [x, y]);
# [[_T^2 + _T - 2, 2*_T + 1, [_T + 5, 2*_T^2 + _T], _T = y, 1]]

that is _T = y, _T^2 + _T - 2 = 0, x = (_T + 5) / (2 _T + 1): the points (2, 1) and (-1, -2).

When no variable separates the solutions, the form is a genuine linear combination:

R := LACE:-rur([x^2 - 1, y^2 - 1], [x, y]):
R[1][4];                                      # _T = -2*x + y

sep imposes a variable:

LACE:-rur([x^2 - 2, y - x - 1/3, z - x*y], [x, y, z], sep = z)[1][4];   # _T = z

With multiplicities = true the RUR does not change and each component gets a sixth element [[j, h_j], ...]: h_j(tvar) is the monic polynomial whose roots are the values of the separating form at the solutions of multiplicity exactly j; the h_j are pairwise coprime and their product is f made monic.

R := LACE:-rur([(x^2 - 2)^2 * (x - 1), y^2 - x], [x, y], multiplicities = true):
R[1][4], R[1][6];
# _T = y, [[1, _T^2 - 1], [2, _T^4 - 2]]

real_solutions

LACE:-real_solutions(sys, vars, prec = 64, threads = 0, nn = 28, verbose = false,
                     sep = 0, strategy = 0, maxprimes = 0)
LACE:-real_solutions(rur, prec = 64, sep = FAIL, threads = 0, squarefree = false)
LACE:-real_solutions(rext, rden, nums, prec = 64, sep = FAIL, threads = 0, squarefree = false)
LACE:-real_solutions(f, prec = 64, constraints = [], constraint_values = false,
                     ensure_derivative_nonzero = false, search_interval = FAIL,
                     threads = 0, squarefree = false)

The real solutions of the zero-dimensional system sys: one list of intervals [lo, hi] per solution, in the order of vars, each with prec significant bits; [] when there is none. nn is the bit size of the primes; the other options are those of rur, which is computed first. precision and nbits are accepted for prec and nn.

The other forms are those of LACE.jl:

• rur = [f, n_1, ..., n_n], a parameterization in the format of zdim_parameterization: x_k = n_k(T) / f'(T) at the roots of f, rational coefficients. sep, the separating vector, may be given: when it is the unit vector of a variable, that variable is T itself and its numerator is not evaluated. squarefree = true says that f is square-free.
• rext, rden, nums: rext(T) = 0, x_k = nums[k](T) / rden(T), rden any polynomial prime to rext.
• f univariate (RS.jl’s rs_isolate): every root of f as a one-element list [box], on which every polynomial of constraints (in the variable of f) has a constant sign; with constraint_values = true, the pair [roots, values], values[i][j] an interval of constraint j at root i (the point [0, 0] where it vanishes), the constraints being scaled to integer coefficients by a positive factor. ensure_derivative_nonzero appends f' to the constraints. search_interval = [l, r] keeps the roots of [l, r] (rational or infinite bounds; a root on a bound is the point); with finite dyadic bounds and no constraint they are isolated in [l, r] only.

B := LACE:-real_solutions([x^2 + y^2 - 5, x - y - 1], [x, y], precision = 64):
nops(B);                                      # 2
map(s -> map(b -> evalf((b[1] + b[2]) / 2), s), B);
# [[-1., -2.], [2., 1.]]

rur_modp

LACE:-rur_modp(sys, vars, p, sep = 0, search = true, strategy = sparse, tvar = _T,
               multiplicities = false)

The RUR of sys modulo the prime p (below 2^30), the Groebner basis being computed by the library: [flag, [f, n_1, ...], sep_lin, sep_var], the polynomials in tvar with coefficients in 0 .. p-1. f is the eliminant, then come the numerators of every variable after a search of the separating form (search = true), of the variables other than x_sep otherwise. sep_lin is the separating form, sep_var the separating variable (from 1, 0 for a genuine form) and flag is 1 on success. strategy is a strategy of zdim_parameterization, by its name or its number. With multiplicities = true the result has a fifth element, the multiplicity factors [[j, h_j], ...] of rur, modulo p.

guess_infos

LACE:-guess_infos(sys, vars, nn = 28, verbose = true, search_strategy = sparse)

The structure of sys modulo one prime of nn bits (nbits, strategy accepted for nn, search_strategy; the strategies of zdim_parameterization): [found, type, quotient, minpoly, sep_lin]. type is shape when the last variable separates and the ideal is radical, cyclic when the ideal is radical but another form is needed, general otherwise; quotient is the dimension of the quotient and minpoly the degree of the minimal polynomial of the separating form sep_lin. The answer is right for all but finitely many primes; it is not certified.

guess_lowest_input_precision

LACE:-guess_lowest_input_precision(sys, vars, nn = 28, verbose = true, lowest_prec = 2)

The fewest decimal digits the coefficients of sys can be rounded to without changing quotient and minpoly of guess_infos: the pair [rounded system, digits].

Discriminant varieties

sys is a list of polynomials with rational coefficients in the unknowns vars (eliminated) and the parameters params. Each result is a reduced Groebner basis for the degree reverse lexicographic order on params, in that order, made of primitive integer polynomials: [1] for the whole ring, [] for the zero ideal. params may be left out: the other names of the polynomials, in lexicographic order. last_nprimes() gives the primes used.

elim

LACE:-elim(sys, vars, params = <the other names>, threads = 0)

The elimination ideal <sys> ∩ Q[params]. For the critical values, add to sys the maximal minors of its Jacobian matrix with respect to vars.

F := [a*x^2 + a + 1]:
LACE:-elim([op(F), diff(F[1], x)], [x], [a]);
# [a + 1]

w_infty

LACE:-w_infty(sys, vars, params = <the other names>, threads = 0)

[zerodim, [I_1, ..., I_k]]: whether sys is generically zero-dimensional in vars and, when it is, for each vars[j], the ideal of Q[params] generated by the leading coefficients of the elements of a Groebner basis for DRL(vars) > DRL(params) whose leading monomial in vars is a power of vars[j]. The list is empty when zerodim is false.

LACE:-w_infty([a*x^6 + b*y^2 - 1, x^2 - a*y - b], [x, y], [a, b]);
# [true, [[1], [a^4]]]

saturate

LACE:-saturate(sys, g, vars, params = <the other names>, threads = 0)

(<sys> : g^∞) ∩ Q[params]; vars may be [] (the saturation itself, in params).

G := [(a - 1)*(s - a), (a - 1)*(a - 2)*s]:
LACE:-elim(G, [s], [a]);            # [a^3 - 3*a^2 + 2*a]
LACE:-saturate(G, a - 1, [s], [a]); # [a^2 - 2*a]

Plane curves

tci_points

LACE:-tci_points(curve, [x, y], a1, b1, a2, b2, nb, prec = 63, precision = 0,
                 threads = 0, verbose = 0)

Points of the curve curve(x, y) = 0 in the box [a1, b1] x [a2, b2] (x in [a1, b1]), computed on a grid of nb lines per direction: the pair [X, Y] of the lists of their coordinates, as floats. Each point is exact to the width of the grid, enough to draw the curve, or to 2^(-precision) when that is finer. The eighth argument prec, as in LACE.jl, is accepted and ignored.

P := LACE:-tci_points(x^2 + y^2 - 1, [x, y], -2, 2, -2, 2, 50):
plots:-pointplot(zip((a, b) -> [a, b], P[1], P[2]));

Interval Newton and the Kantorovich test

These functions take a square system sys in the variables vars, with rational or float coefficients (floats are read exactly and the coefficients are not scaled), and a starting point x0: a list of numbers, or a list of intervals [lo, hi] for a box; one equation and one number for one variable. Their arguments are those of LACE.jl, after vars: input_precision (accepted: the coefficients are read exactly and made intervals of work_precision bits), work_precision, the number of bits of the interval arithmetic, output_precision, the significant bits the results are rounded outward to, and for the iterations target_precision and max_nb_loops. verbose sets the messages of the library.

The Kantorovich test certifies that Newton’s method converges from x0 to a zero of sys, unique in the ball of radius 2 eta around x0. Its classical bound depends on how the equations are scaled; flags, a set among {hessian, precond, schulz} (or the integer sum of KANTO_HESSIAN = 1, KANTO_PRECOND = 2, KANTO_SCHULZ = 4), selects bounds that depend less, or not at all, on that scaling:

flag	bound
hessian	a bound independent of a diagonal scaling of the equations
precond	a bound on a preconditioned system, independent of the scaling of the equations
schulz	with precond, a better preconditioner when the first one is poor

The constants are [2 eta, mu] without flags, else [2 eta, mu, mu1, mu2, delta], each an interval: mu is the best bound computed (the point is certified when mu < 1), mu1 the one of hessian, mu2 the one of precond (infinity when not computed), and delta the defect of the preconditioner (-1 without precond).

kantorovich_test

LACE:-kantorovich_test(sys, vars, x0, input_precision = 63, work_precision = 63,
                       output_precision = 2, verbose = 1, flags = 0)

The constants, rounded outward to output_precision bits; the point is certified when mu < 1, that is when constants[2][2] < 1. precision is accepted for work_precision.

S := [x^2 + y, x^2 + y + x/1024]:
LACE:-kantorovich_test(S, [x, y], [1/256, 1/256], 63, 128, 128);
# not certified: mu = 64
LACE:-kantorovich_test(S, [x, y], [1/256, 1/256], 63, 128, 128, flags = {hessian, precond});
# certified: mu1 = 32, mu2 = 1/64

interval_newton

LACE:-interval_newton(sys, vars, x0, input_precision = 63, work_precision = 63,
                      target_precision = 2, output_precision = 2, max_nb_loops = 10,
                      refine = 1, verbose = 1, flags = 0)

The interval Newton method from x0, after the Kantorovich test: [status, box, constants], box a list of intervals, one per variable, rounded outward to output_precision bits. Each step is intersected with the current box (refine = 1), or only when the step contains 0 (refine = 0). status is 0 for a singular Jacobian, 1 when a zero is proved to lie in the box and every interval is narrower than 2^(-target_precision), 2 when the box no longer shrinks, 3 after max_nb_loops steps, 4 when the starting box is proved to hold no zero, 5 when the width is reached without a proof that the box holds a zero. From a point (x0 a list of numbers), the point is first inflated into a box proved to hold a zero (epsilon-inflation): the status 1 then means that the box returned holds a zero, and any other status gives no enclosure. From a box, the status 1 also proves a zero in the box returned. The names of the previous versions are accepted: precision, target, maxsteps for work_precision, target_precision, max_nb_loops, and always_intersect = true for refine = 1.

LACE:-interval_newton([x^2 - 2, y - x], [x, y], [1.4, 1.4], 63, 128, 100, 128);

interval_classic_newton

LACE:-interval_classic_newton(sys, vars, x0, input_precision = 63, work_precision = 63,
                              target_precision = 2, output_precision = 2,
                              max_nb_loops = 10, verbose = 1, flags = 0)

The same without intersections, for a start that may not enclose the solution.

Bivariate systems

The functions below take two polynomials p1, p2 in exactly two variables and the variable v1 to eliminate; the other variable, v2, is found from the polynomials. The polynomials in the results are polynomials in v2 (or in the parameter).

Each of them, and biv_parameterization, biv_decomposition and biv_psc below, takes certified = false: the fast version, whose lift stops when it reduces to the image of a witness prime. certified = true lifts to the bound and checks the decomposition over the integers by exact divisions, the primes all drawn again when a check fails (rs_biv_set_certify).

bz_biv_resultant

LACE:-bz_biv_resultant(p1, p2, v1, threads = 0)

The resultant of p1 and p2 with respect to v1, a polynomial in v2 (defined up to a constant factor).

bz_biv_decomposition

res, asymp, res_gen, res_non_gen := LACE:-bz_biv_decomposition(p1, p2, v1,
                                        nbits = 62, threads = 0)

The solutions split by level, the level k of a value u of v2 being the degree of gcd(p1(v1, u), p2(v1, u)). With S_k = s_0 + s_1 v1 + ... + s_k v1^k the subresultant of index k:

• res: the resultant with respect to v1;
• asymp: the polynomial of the values of v2 above which there are solutions at infinity, 1 when there are none;
• res_gen[k] = [s_(k-1) + k s_k v1, univ]: at the roots of univ there is one solution, v1 being the root of the first polynomial, which is linear in v1;
• res_non_gen[k] = [S_k, rest]: at the roots of rest, the values of v1 are the roots of S_k.

An empty component is [0, 0].

c := g^2 + o^2 - 1:
res, asymp, gen, nongen := LACE:-bz_biv_decomposition(c, diff(c, g), g):
# gen[1] = [2*g, o^2 - 1] up to constant factors: g = 0 above o = -1 and o = 1

bz_biv_is_generic(res_non_gen) is true when every non-generic component is empty; bz_biv_is_solved(res, res_gen, res_non_gen, v2) is true when, in addition, the generic components describe the solutions.

bz_biv_rur

LACE:-bz_biv_rur(p1, p2, v1, nbits = 62, threads = 0)

A RUR with parameter v2: a list of components [k, rext, rden, num_v1, num_v2], k the multiplicity, rext(v2) = 0, v1 = num_v1 / rden and num_v2 = v2 * rden. An error is raised when v2 does not separate the solutions; use bz_biv_rur_shear then.

bz_biv_rur_shear

LACE:-bz_biv_rur_shear(p1, p2, v1, tvar = _T, nbits = 62, threads = 0,
                       sep = 0, a = FAIL)

A RUR with parameter tvar = v2 + a v1, the integer a being chosen by the library. The result is [a, tvar, comps], each component [mult, rext, rden, num_v1, num_v2] with rext(tvar) = 0, v1 = num_v1 / rden, v2 = num_v2 / rden. a = 0 imposes tvar = v2; another imposed value of a is refused.

c2 := g^2 - o^3:
a, T, comps := op(LACE:-bz_biv_rur_shear(c2, diff(c2, g), g, tvar = _T)):

bz_biv_decomposition_shear

a, t, res, asymp, res_gen, res_non_gen :=
    LACE:-bz_biv_decomposition_shear(p1, p2, v1, a = FAIL, tvar = T,
                                     nbits = 62, threads = 0)

The decomposition in the coordinates (v1, t), v2 = t - a v1. Without an imposed a: when bz_biv_decomposition already solves the system, a = 0 and t = v2; otherwise a is the shear of bz_biv_rur_shear and t = tvar.

The interface of LACE.jl

The functions of the Julia package LACE.jl, with its arguments, options and results, the conventions above turning its objects into Maple ones. A system is given with its variables, which Julia takes from the ring of the polynomials; the polynomials of a result are in tvar (_T by default) when they are univariate in a new variable.

zdim_parameterization

LACE:-zdim_parameterization(sys, vars, nn = 28, verbose = true, threads = 0,
                            get_separating_element = false, max_primes = 200000,
                            search_strategy = sparse, eliminate_linear = true,
                            parallelism = default, composite = 4,
                            multiplicities = false, characteristic = 0, tvar = _T)

The RUR of a zero-dimensional system in the format of RationalUnivariateRepresentation.jl: the list [f, n_1, ..., n_n], f the eliminant and n_k the numerator of vars[k] (T f'(T) for the separating variable), vars[k] = n_k(T) / f'(T) at the roots of f, with rational coefficients; with characteristic = p, the same modulo the prime p, coefficients in 0 .. p-1. nn is the bit size of the primes, max_primes bounds their number. get_separating_element = true returns [rur, lin], lin the separating linear form.

search_strategy, a name or its number: sparse (0, the free variables, then the sparse escalation on them), random_free (1, the free variables, then random forms on them), variables (2, the free variables only), incremental (3, T_1 = y_1, then T_k = T_(k-1) + a y_k with a = 0, 1, -1, 2, -2, ... until y_k is a function of T_k), adaptive (4, the free variables, then random forms with small coefficients), random_last (5, the last variable, then the forms of random_free), random_direct (6, the forms of random_free only), l0_norm / mron_0l (7 / 8, the free variables, then supports of increasing size, x_1 / x_n first), deterministic (9, the last variable, then the forms (1, i, i^2, ...) on the free variables). All the forms are on the free variables, every other linear form taking the same values at the solutions as one of them; the whole search runs in the library. random, the former name of random_free, is still accepted, with a warning. eliminate_linear = false keeps the linear variables in the engine; parallelism = serial forces one thread; composite is accepted and ignored.

multiplicities = true asks for the multiplicities of the solutions as well: the RUR does not change and the result becomes [rur, mult] ([rur, lin, mult] with get_separating_element), mult the list of the pairs [j, h_j] of rur. A system without solution raises an error.

R := LACE:-zdim_parameterization([x^2 + y^2 - 5, x*y - 2], [x, y], verbose = false);
# [_T^4 - 5*_T^2 + 4, ..., 4*_T^4 - 10*_T^2]: T = y, x = R[2] / diff(R[1], _T)

rur_core

LACE:-rur_core(sys, vars, p, verbose = true, tvar = _T)

The RUR modulo p with the last variable as separating element: [flag, [f, n_1, ...]], as rur_modp(..., search = false).

biv_parameterization

LACE:-biv_parameterization([p1, p2], [x_1, x_2], nbits = 62, threads = 0, sep = 0,
                           verbose = true, tvar = _T)

The RUR of two polynomials in two variables by the subresultants over Z (no Groebner basis): one component per multiplicity, each the list [rext, rden, num_X, num_Y, mult, lin], rext(T) = 0, x_1 = num_Y / rden, x_2 = num_X / rden at the roots of rext, T = x_2 + lin[1] x_1 (sep = 2 imposes T = x_2).

biv_decomposition

LACE:-biv_decomposition([p1, p2], [x_1, x_2], nbits = 62, threads = 0,
                        verbose = true, psc = false)

The decomposition of bz_biv_decomposition, as a list [res, asymp, res_gen, res_non_gen], x_1 eliminated; a constant resultant raises an error. With psc = true a fifth element, the signed principal subresultant coefficients [sRes_0, ..., sRes_p] of the pair with respect to x_1, polynomials in x_2; for (P, diff(P, x_1)) it is the Sturm-Habicht sequence of P: above a value u of x_2 where sRes_p(u) <> 0, the number of distinct real roots of P(x_1, u) is pmv([sRes_p(u), ..., sRes_0(u)]). Without psc = true these coefficients are not computed. biv_psc computes that sequence alone.

biv_decomposition_modp

LACE:-biv_decomposition_modp([p1, p2], [x_1, x_2], p, verbose = true, psc = false)

The decomposition of biv_decomposition modulo the prime p < 2^63, in the same format, the coefficients in [0, p): the resultant and the S_k are those of the pair reduced modulo p, the other factors monic. p must exceed the degree in x_1 and that of the resultant, and leave enough evaluation points for the subresultants. Where two roots of the resultant meet modulo p, or a leading coefficient vanishes at one of them, the result is the decomposition of the reduced pair, not the reduction of the decomposition over Q.

LACE:-biv_decomposition_modp([x^2 + y^2 - 5, x*y - 2], [x, y], 101);
       [y^4 + 96*y^2 + 4, 1, [[x*y + 99, y^4 + 96*y^2 + 4]], [[0, 0]]]

biv_psc

LACE:-biv_psc([p1, p2], [x_1, x_2], nbits = 62, threads = 0, verbose = true)

The signed principal subresultant coefficients [sRes_0, ..., sRes_p] of p1, p2 with respect to x_1, polynomials in x_2, lifted multi-modularly without the decomposition. For [P, diff(P, x_1)] it is the Sturm-Habicht sequence of P.

SH := LACE:-biv_psc([x^2 + y^2 - 1, 2*x], [x, y]);      # [4 - 4*y^2, 2, 1]
LACE:-pmv(ListTools:-Reverse(eval(SH, y = 0)));        # 2 real roots in x above y = 0

biv_resultant

LACE:-biv_resultant([p1, p2], [x_1, x_2], threads = 0)

The resultant with respect to x_1, primitive, a polynomial in x_2.

pmv

LACE:-pmv(s)

Permanences minus variations of s = [s_p, ..., s_0] (decreasing indices, s_p <> 0): the sum, over the consecutive nonzero entries s_i, s_j (i > j, zeros between) with i - j odd, of eps_(i-j) signum(s_i s_j), eps_k = (-1)^(k(k-1)/2).

isolate

LACE:-isolate(f, prec = 0, threads = 0, max_roots = 0, newton = true, trunc = true,
              squarefree = false)
LACE:-isolate(f, [l, r], prec = 0, threads = 0, newton = true, trunc = true,
              squarefree = false)
LACE:-isolate([[lo_0, hi_0], ..., [lo_n, hi_n]], [l, r], precision = 64)
LACE:-isolate(iso)

The real roots of the univariate f, sorted increasing, each an interval [lo, hi] with dyadic bounds (an exact root when lo = hi), refined to a width below 2^(-prec) when prec > 0. The square-free part is isolated unless squarefree = true; max_roots > 0 stops after that many roots (a budget, not the smallest ones); newton, trunc are the Newton jump and the truncated frame of the isolation.

• With [l, r] (dyadic or infinite bounds), the roots of the closed interval [l, r] only, isolated there; a root on a bound is the point.
• A list of intervals is a polynomial with interval coefficients (constant term first, the bounds rounded outward to precision bits): [roots, nb_undecided], the roots in the open interval (l, r), simple, of every polynomial of the box, and what the width of the coefficients left undecided. An error when the leading coefficient holds 0.
• On an isolation, its current intervals.

Root(num_lo, den_lo, num_hi, den_hi) is [num_lo / 2^den_lo, num_hi / 2^den_hi]; lo(r), hi(r), width(r) and isexact(r) are its bounds, its width and lo = hi. dyadic_outward(q, up, k) is q when it is dyadic, else q rounded up (up = true) or down to a multiple of 2^(-k); exact_rational(x) is the rational of the number x, every digit of a float kept.

LACE:-isolate(6 - 11*x + 6*x^2 - x^3);          # [[0, 2], [2, 2], [2, 4]]
LACE:-isolate(6 - 11*x + 6*x^2 - x^3, [3/2, 3]); # [[2, 2], [3, 3]]

isolation

LACE:-isolation(f, threads = 0, max_roots = 0, newton = true, trunc = true,
                squarefree = false)

An isolation kept to be tightened later, as a record: sf the square-free part isolated, nb the number of roots, prec the width reached (0: isolation only), reached whether the last refinement reached its precision on every interval, threads, t_squarefree, t_isolation the seconds of the square-free part and of the isolation, t_refinements the pairs [bits, seconds] of the refinements. Its intervals live in the library until finish().

• refine(iso, prec, threads = iso:-threads) tightens every interval below 2^(-prec), in place, from where it stands, and returns the seconds; a coarser request costs nothing.
• isolate(iso) gives the intervals, isolated_root(iso, k) interval k.
• refine_root(iso, k, prec) refines interval k alone and returns it.
• set_root(iso, k, lo, hi) replaces interval k (dyadic bounds, not checked); the next refinement starts from there.
• reached_width(iso): the bits reached at least by the widest interval, infinity when they are all exact.

iso := LACE:-isolation(x^2 - 2):
LACE:-refine(iso, 40):
LACE:-isolate(iso);                 # the two roots to 2^-40
LACE:-refine_root(iso, 2, 200);     # sqrt(2) to 2^-200

real_solution_at

LACE:-real_solution_at(f, g, nums, [lo, hi], prec = 64)

The values nums[k](T) / g(T) at the one root of the square-free f in the interval [lo, hi] (dyadic bounds; open when lo < hi, a point when lo = hi), each with prec significant bits, a vanishing one the point 0: [root, values], the root refined, or FAIL when [lo, hi] does not hold exactly one root of f. g = 1 is allowed.

rur_coords

LACE:-rur_coords(f, g, nums, intervals, prec = 64, threads = 0)

The same at the roots of f isolated by intervals, a list of [lo, hi] with dyadic bounds, exactly one root in each (open when lo < hi, the bounds excluded; a point when lo = hi; an interval where f does not change sign is an error): one [root, values] per interval, in their order. The roots are refined, not isolated again; for a parameterization with a separating form t, the interval t(B) of a box B will do.

signs_at

LACE:-signs_at(f, polys, intervals, nonzero = [], pmv = false, threads = 0)

The exact signs of polys at the roots of the square-free f isolated by intervals: [roots, signs], signs[i][k] the sign of polys[k] at root i, each root refined only as far as the signs need (roots[i] its interval then). A zero is decided by a gcd; nonzero lists the k known not to vanish. With pmv = true an entry between two nonzero ones of opposite signs is left undecided, as 0: it changes nothing to pmv of the sequence.

LACE:-signs_at(x^2 - 2, [x - 1, x^2 - 2, 3], [[-2, -1], [1, 2]])[2];
# [[-1, 0, 1], [1, 0, 1]]

eval_bfi

LACE:-eval_bfi(polys, [lo, hi], prec = 64)
LACE:-eval_bfi([[lo_0, hi_0], ..., [lo_n, hi_n]], [lo, hi], prec = 64)

The univariate polys at the interval [lo, hi] (its bounds rounded outward to prec bits), each enclosed at precision prec by the evaluation of the library: hardware intervals when they suffice, MPFI otherwise; or the polynomial with interval coefficients at [lo, hi].

LACE:-eval_bfi([1 + 2*x + 3*x^2, 5*x^3 - 1], [1/4, 1/2]);   # [[27/16, 11/4], [-59/64, -3/8]]

tci_columns

LACE:-tci_columns(P, [x_1, x_2], a, b, nb, prec, ymin = -infinity, ymax = infinity,
                  squarefree = false, threads = 0)

The roots of P(x_1, X_i) in [ymin, ymax] for X_i = a + i (b - a) / nb, i = 0 .. nb - 1 (x_2 the abscissa, a, b dyadic), refined to a width at most 2^(-prec), and the number of roots under ymin above each X_i: [roots, below]. ymin, ymax must be doubles exactly; squarefree = true says every P(x_1, X_i) is square-free. prec is required.

tci_track

LACE:-tci_track(P, [x_1, x_2], arcs, crit = 0, window, eps, tol, ycrit = [],
                maxsteps = 0, threads = 0)

The arcs of the curve P = 0 followed step by step, x_2 the abscissa, x_1 the ordinate. Each arc is [x0, ylo, yhi, dir, kind, e1, e2, e3, e4, inf] (dyadic numbers and integers): it starts at the root of P(x_1, x0) isolated by [ylo, yhi] (refined as the other points: the interval need only isolate it) and goes in the direction dir (1 or -1) up to e1 (kind = 0, a side of the window) or up to the root of crit (a polynomial in x_2) isolated by [e1, e2] (kind = 1), where it ends at the point of ordinate in [e3, e4] (inf = 0), goes to -infinity, infinity (inf = -1, 1), or crosses a simple root of P (inf = 2: the arc ends at the middle of [e1, e2], at the root of P there). crit and ycrit come from the caller: tci_track neither computes nor checks them. window = [xmin, xmax, ymin, ymax]; the boxes have a diagonal below eps, the points a width below tol; ycrit the intervals of the abscissae of the horizontal tangents. Returns, for each arc, [x, y, certified, steps, last]: its points as floats, whether every piece is certified, the steps made, and the abscissa where it was given up (FAIL otherwise).

rs_isolate

LACE:-rs_isolate(f, output_precision = 2, constraints = [], constraint_values = false,
                 search_interval = FAIL, ensure_derivative_nonzero = false, threads = 0,
                 verbose = 0, hybrid_limit = 2^30)
LACE:-rs_isolate(sys, vars, output_precision = 2, threads = 0, nn = 28, verbose = 0)
LACE:-rs_isolate(rext, rden, nums, sepvector, output_precision = 2, threads = 0)
LACE:-rs_isolate(rur, sepvector, output_precision = 2, threads = 0)

RS.jl’s isolation: the forms of real_solutions with output_precision for prec; sepvector is checked against the numerators before it is trusted, hybrid_limit is ignored. rs_rur(sys, vars, verbose = 0, threads = 0, nn = 28) is RS.jl’s RUR, [ext, den, coords, sep_vect] with integer coefficients, den a multiple of ext'; rs_tci_points is tci_points; rs_get_last_sep_vector() the separating form of the last search; rs_newton raises an error (use kantorovich_test and interval_newton).

bz_biv_rur_all

LACE:-bz_biv_rur_all(p1, p2, v1, nbits = 62, threads = 0)

RS4.jl’s form of bz_biv_rur: every component as [rext, rden, num_X, num_Y, mult], num_X the numerator of v1, num_Y = T rden.