Symbolic integration is a core component of computer algebra systems such as Maple, where the selection of an appropriate integration sub-algorithm has traditionally relied on rule-based heuristics. Because different sub-algorithms are suited to different classes of integrands, accurate algorithm selection significantly affects both the quality and efficiency of the integration result. Prior work has applied LSTM and TreeLSTM models to this task, with TreeLSTM achieving stronger performance. We first reproduce the TreeLSTM experimental framework and then introduce a Transformer model for comparison. Multiple Transformer configurations are evaluated on the original dataset. Results show that one Transformer variant matches the reproduced TreeLSTM in exact optimal matching accuracy but lags behind under 5
GCDHEU is a heuristic algorithm for computing greatest common divisors of multivariate polynomials over the integers and remains available in several computer algebra systems. In this paper, we revisit the GCDHEU algorithm from both theoretical and practical viewpoints. On the theoretical side, we show that the standard correctness argument leaves a gap in the genuinely multivariate setting, and we close this gap by establishing a multivariate Cauchy-type bound. This yields a rigorous proof of correctness for the algorithm. On the practical side, we observe that, for multivariate inputs, the order of variable specialization can substantially affect coefficient growth in the induced integer GCD computations. Motivated by this observation, we propose a simple degree-aware variable-ordering heuristic that helps control such growth. A small modification to SymPy’s implementation of GCDHEU yields an average speedup of about 13.7% on our benchmark set.
Motivated by algorithms manipulating parametric polyhedral sets, we discuss the computation of convex hulls of polyhedra. We extend the algorithm proposed by Fukuda, Liebling and Lütolf, by relaxing its assumptions. In particular, we show how to compute the H-representation (resp. V-representation) of the convex hull (resp. the closure of the convex hull) of polyhedral sets P_1, … , P_k , given by H-representation and such that some P_i ’s may not be of maximum dimension (resp. may not be bounded).
New finite element method (FEM) schemes have been developed for solving self-adjoint boundary value problems (BVPs) for a system of N coupled elliptic equations with second-order mixed partial derivatives (PDEs) on multidimensional domains using symbolic-numerical algorithms implemented in Wolfram Mathematica. Basis functions are constructed using multidimensional Hermite interpolation polynomials (HIPs) on hyperparallelepiped meshes that preserve the continuity of derivatives up to a given order of FEM solutions on the boundaries of adjacent finite elements. The Helmholtz problem for hyperellipsoids mapped in generalized non-orthogonal hyperspherical coordinates to hyperparallelepipeds was chosen as a test example. Benchmark calculations of the lower part of the spectra of BVPs confirm the declared order of accuracy of the finite element method schemes.
The Jacobian of a hyperelliptic curve is the subject of many conjectures and open problems in algebraic geometry and number theory. Efficient Jacobian arithmetic plays a crucial role in the experimental investigation of these problems. For curves of low genus, explicit formulas effectively realize this arithmetic, but manually developing such formulas is not feasible in higher genera. We introduce the first software package for automating the process of generating and verifying explicit formulas for curves of arbitrary genus. Our formulas outperform generic algorithms for Jacobian arithmetic and exhibit comparable operation counts against manually fine-tuned formulas for small genera.
We investigate the interest of integral equations with fractions for the treatment of systems of parametric nonlinear differential equations. We construct a new parametric system that admits a nontrivial integral input-output equation involving nested division and integration operators. We show the practical interest of such an equation by performing the identifiability analysis and parameter estimation.
This survey is the review of the author’s papers dedicated to the use of the computer algebra package NORT for constructing approximations of periodic solutions families of nonlinear dynamical systems, as well as studying the integrability of systems of this type. A brief description of the package’s structure and key concepts are provided with some applications examples. Here we discuss several topics of practical application of the resonance normal form method, based on A.D. Bruno’s fundamental theorem on the convergence of a normalizing transformation.
Polynomial evaluation is ubiquitous in scientific computing, yet simultaneously achieving high performance and numerical accuracy remains challenging. Traditional methods such as Horner’s scheme are efficient but inherently sequential, while the Estrin scheme offers a parallel-friendly evaluation tree structure but does not inherently improve accuracy. This paper introduces a compensated Estrin scheme that combines the low-depth evaluation tree of the Estrin method with error compensation techniques based on error-free transformations. The proposed algorithm dynamically tracks and compensates for rounding errors at each stage of the evaluation, delivering accuracy close to double-double precision arithmetic while using only binary64 error-free transformations. A rigorous rounding error analysis proves that the relative forward error is bounded by u + cond(p,x)· O(u^2) , where u is the unit roundoff and cond(p,x) is the condition number of the polynomial evaluation. Numerical experiments on ill-conditioned polynomials demonstrate that the compensated Estrin scheme maintains high accuracy over a much larger range than the standard Estrin scheme, up to the regime where the condition number is too large for any binary64-output twice-working-precision-like method to remain accurate, providing a practical and efficient solution for accurate polynomial evaluation on modern architectures.
Let K(t)[x;σ ,δ ] be a skew polynomial ring over the field of rational function K(t). We first explore the possible forms of the involutions on K(t)[x], K(t)[x;δ ] and K(t)[x;σ ] , respectively. Then, we give some necessary and sufficient conditions for an involution ϕ of K(t) to be extended on K(t)[x], K(t)[x;δ ] , and K(t)[x;σ ] as involutions, respectively. Finally, we prove that the involutions ϕ on K(t) cannot be extended to K(t)[x;σ ,δ ] if σ 1 and δ 0 .
In this paper, we present a method for solving the inverse kinematics (IK) problem of a 7-Degree-of-Freedom (7-DOF) robot manipulator using Comprehensive Gröbner Systems (CGS). Due to kinematic redundancy in 7-DOF manipulators, additional constraints are required to obtain finite solutions. In our approach, the angle of the final joint is predetermined to satisfy task-specific constraints such as obstacle avoidance, thereby reducing the problem to a 6-DOF IK formulation. The remaining joint angles are derived systematically in stages: first, determining the elbow joint angle, then solving for wrist joint angles using a CGS-based method for polynomial systems, and finally computing the shoulder joint angles via polynomial equations. By precomputing the CGS for parametric polynomial systems, the proposed method avoids repeated Gröbner basis computations for different end-effector configurations, enabling more efficient handling of trajectory planning problems. The algorithm is implemented in a computer algebra system and evaluated through numerical experiments, demonstrating that the proposed method achieves high accuracy and practical computational performance. Furthermore, the results confirm that the number of valid IK solutions is bounded and that the method can effectively identify feasible configurations among multiple candidates.
A multilinear Boolean polynomial f is a polynomial over GF(2) in which each variable has degree at most 1. Such polynomials arise in Boolean circuit optimization, yet their efficient factorization remains challenging. We present two Monte-Carlo algorithms that advance this problem. Our first factorization algorithm assumes f is given in the sparse representation. It has algebraic complexity O(n^2t) over a suitable extension field GF(2^k) , where n is the number of variables and t is the number of terms of f. Our C implementation achieves substantial speedups over both the FD and GCD algorithms of Emelyanov–Ponomaryov [6]. Our second algorithm assumes f is given by a black box B for its evaluation. Here we apply our recently developed black box factorization algorithm CMBBSHL which is implemented in Maple and C. The black box representation allows us to reduce t (the number of terms of the input polynomial f) to T = s_max( ∑#f_i + ℭ( probe B) ) (which can be ≪ t ), where the f_i are the irreducible factors of f, s_max is the maximum number of terms in any coefficient (in the highest ranking variable) in any factor, ∑#f_i is the total number of terms in all irreducible factors, and ℭ( probe B) is the cost of a single black box probe. This yields an overall algebraic complexity O( n^2 T ) . Our Maple and C implementation of our second algorithm is the fastest algorithm when T ≪ t .
Tameness and bifurcation sets of complex polynomial mappings are investigated from the viewpoint of symbolic computation. Based on results of S. Broughton, effective algorithms are developed for testing the tameness of complex polynomial mappings and for computing the bifurcation sets of tame complex polynomial mappings. These algorithms rely on comprehensive Gröbner systems as their main computational tool. Several examples illustrating the computation of bifurcation sets of tame polynomial mappings are presented.
Montgomery curves are an important family of elliptic curves in computational arithmetic, especially because of their efficiency in cryptographic applications. For an ordinary elliptic curve over a finite field, the unit root is a fundamental arithmetic invariant. In this paper, we give a formula, called the unit-root formula, expressing the unit root as a value of a p-adic analytic function. More precisely, for an ordinary Montgomery curve, we express the unit root as a product of values of a p-adic hypergeometric function. Since the Frobenius trace, and hence the number of rational points, can be recovered from the unit root, our formula also gives a connection from the arithmetic of Montgomery curves to the computation of special values of p-adic hypergeometric functions. We further discuss how these values can be computed using Asakura’s algorithm, which computes such values to precision p^N with bit complexity O(N^4(log N)^3) for fixed p, and briefly indicate possible applications to point-counting methods.
Given a finite number of sparse polynomials, we present a new algorithm, along with its C++ implementation, to compute an efficient straight-line program which is able to simultaneously evaluate these polynomials at a given point.
This survey article provides an overview of the fundamental principles used to simplify multi-sums into indefinite nested sums over hypergeometric products in the setting of difference rings. We place special emphasis on the algorithmic translation between hypergeometric sums and the formal difference ring setting. Furthermore, we detail the core summation paradigms of telescoping, creative telescoping, and recurrence solving within difference rings, illustrating these techniques and their underlying algorithms with concrete examples.
Let p be a prime and =ℤ/p^2ℤ the ring of integers modulo p^2 . Any A∈^n× n is unimodularly equivalent to its Smith form S=diag (1,… ,1_r_0, p,… ,p_r_1, 0,… ,0_r_2 ) ∈^n× n, i.e., there exist U,V∈^n× n such that UAV=S , with U, V∈^* (where ^* is the set of units in , elements not equivalent to 0 p ). Our goal in this paper is to determine r_0,r_1,r_2 quickly when A is sparse or structured. By “sparse” we mean A is given by a black box such that for any v∈^n× 1 we can compute v↦ Av with operations in , which captures having few nonzero elements or a multiplicative structure (e.g., Hankel or Toeplitz matrices). We present a randomized algorithm which requires an expected number of to compute the Smith form, where ω is the exponent of dense matrix multiplication. Using standard cubic matrix multiplication ( ω =3 ) our algorithm thus requires operations in , while using the current asymptotically fastest matrix multiplication, with ω <2.371339 , our algorithm requires operations in . Our algorithm is probabilistic of the Monte Carlo type, meaning it fails on any invocation with controllably small probability. We employ iterative block-Wiedemann-style matrix techniques and structured preconditioners. To our knowledge, this is the first algorithm to compute the modular Smith Normal Form modulo p^2 requiring fewer than operations in , i.e., faster than any dense algorithm.
We study algorithms for factorization in the quantum plane of (dilation) skew polynomials over a function field of a cyclotomic field: =(t)[x;σ ], =ℚ(ω ), σ (t)=ω t, where ω∈ℂ is a primitive m -th root of unity. We start with the established approach through central elements and factor the central left multiples, staying in characteristic zero, to obtain a partial decomposition. A two-level modular approach is proposed: specialize a central parameter to good algebraic values, study the resulting cyclic algebras over number fields, and then reduce further at good inert primes so that fast finite-field skew-factorization algorithms apply. A prototype SageMath implementation is provided to experiment with the algorithms. We then look at the effect of extending the field of constants from ℚ(ω ) to ℚ , an algebraic closure of ℚ , and factoring over ℚ(t)[x;σ ] . In this case we show factorization is decidable in the exact algebraic model based on finite extensions.
We present a novel S-pair selection strategy called Homogeneity Entropy, for deciding the sequence of S-polynomials to construct in Buchberger’s algorithm to compute a Gröbner basis. The strategy uses an information-theoretic measure derived from the distribution of degrees among the monomials of the S-polynomial: a very different approach to the classical heuristics such as Degree, Normal and Sugar, or indeed the more recent machine learning approaches to the problem. We implement this strategy and evaluate it on two different datasets: (1) variations of randomly generated polynomial systems with controlled numbers of variables, degrees, densities and number of polynomials per system; and (2) the PHCpack benchmark dataset sourced from real world problems. The Homogeneity Entropy strategy significantly outperforms classical strategies on random polynomial datasets, but on the PHCpack dataset the classical strategies perform better. This suggests the right strategy varies with the shape of the data and we explore this in several experiments. The new strategy offers practically meaningful gains on certain distributions, and represents the first use of such information-theoretic guidance in the optimisation of symbolic computation algorithms.
This paper presents the issues arising in implementing a fast integer division algorithm on general purpose GPUs. The algorithm uses a Newton iteration based on the shifted inverse operation, keeping all arithmetic in the integer domain and relying on data-parallel operators. The principal contribution is an efficient GPU/cuda implementation for integer precisions from 2^15 to 2^18 – sizes not supported by cgbn division. We propose algorithmic refinements, define a cost model in terms of multiplications, build on prefix sums and previous work on multi-precision multiplication, and present an evaluation showing near-optimal performance relative to the model for the target precision.
Tagging of generic dynamic values is important in symbolic-computation and dynamic-language systems, but the trade-offs change as machine architectures and workloads evolve. In particular, old folklore about boxed values, immediate values, and type tags must be recalibrated from time to time. We revisit the performance of badged object headers, low-bit tagging, and two NaN-boxing layouts on a range of platforms in use today, including AArch64 and x86-64 architectures from different manufacturers. The experiments isolate two distinct effects: the cost avoided by not heap-allocating common scalar values, and the cost avoided by obtaining tag information from the value word rather than by performing a heap read. The results show that several local bit operations are often cheaper than opening a heap object to obtain a tag or small value. Low-bit tagging remains the simplest and usually fastest choice for mostly symbolic workloads, while NaN-boxing is close in access cost and avoids the time and space of heap allocation for ordinary floating-point values.
