A.9.5 Algebra and Calculus

Polynomial manipulation

For univariate polynomials, Factor uses a variant of the Cantor-Zassenhaus algorithm to factor modulo a prime, then uses Hensel lifting and recombination to build up factors over the integers.

Factoring over algebraic number fields is done by finding a primitive element over the rationals and then using Trager's algorithm.

For multivariate polynomials Factor works by substituting appropriate choices of integers for all but one variable, then factoring the resulting univariate polynomials, and reconstructing multivariate factors using Wang's algorithm.

The internal code for Factor exclusive of general polynomial manipulation is about 250 pages long.

FactorSquareFree works by finding a derivative and then iteratively computing GCDs.

Resultant uses either explicit subresultant polynomial remainder sequences or modular sequences accompanied by the Chinese Remainder Theorem.

Apart uses either a version of the Padé technique or the method of undetermined coefficients.

PolynomialGCD and Together usually use modular algorithms, including Zippel's sparse modular algorithm, but in some cases use subresultant polynomial remainder sequences.

For multivariate polynomials the Chinese Remainder Theorem together with sparse interpolation are also used.

Symbolic linear algebra

RowReduce, LinearSolve, NullSpace and MatrixRank are based on Gaussian elimination.

Inverse uses cofactor expansion and row reduction. Pivots are chosen heuristically by looking for simple expressions.

Det uses direct cofactor expansion for small matrices, and Gaussian elimination for larger ones.

MatrixExp finds eigenvalues and then uses Putzer's method.

Zero testing for various functions is done using symbolic transformations and interval-based numerical approximations after random numerical values have been substituted for variables.

Exact equation solving and reduction

For linear equations Gaussian elimination and other methods of linear algebra are used.

Root objects representing algebraic numbers are usually isolated and manipulated using validated numerical methods. With ExactRootIsolation->True, Root uses for real roots a continued fraction version of an algorithm based on Descartes' rule of signs, and for complex roots the Collins-Krandick algorithm.

For single polynomial equations, Solve uses explicit formulas up to degree four, attempts to reduce polynomials using Factor and Decompose, and recognizes cyclotomic and other special polynomials.

For systems of polynomial equations, Solve constructs a Gröbner basis.

Solve and GroebnerBasis use an efficient version of the Buchberger algorithm.

For non-polynomial equations, Solve attempts to change variables and add polynomial side conditions.

The code for Solve and related functions is about 500 pages long.

For polynomial systems Reduce uses cylindrical algebraic decomposition for real domains and Gröbner basis methods for complex domains.

With algebraic functions, Reduce constructs equivalent purely polynomial systems. With transcendental functions, Reduce generates polynomial systems composed with transcendental conditions, then reduces these using functional relations and a database of inverse image information.

CylindricalDecomposition uses the Collins-Hong algorithm with Brown-McCallum projection for well-oriented sets and Hong projection for other sets. CAD construction is done by Strzebonski's genealogy-based method using validated numerics backed up by exact algebraic number computation. For zero-dimensional systems Gröbner basis methods are used.

For Diophantine systems, Reduce solves linear equations using Hermite normal form, and linear inequalities using Contejean-Devie methods. For univariate polynomial equations it uses an improved Cucker-Koiran-Smale method, while for bivariate quadratic equations, it uses Hardy-Muskat-Williams methods for ellipses, and classical techniques for Pell and other cases. Reduce includes specialized methods for about 25 classes of Diophantine equations, including the Tzanakis-de Weger algorithm for Thue equations.

With prime moduli, Reduce uses linear algebra for linear equations and Gröbner bases over prime fields for polynomial equations. For composite moduli, it uses Hermite normal form and Gröbner bases over integers.

Resolve mainly uses an optimized subset of the methods from Reduce.

Reduce and related functions use about 350 pages of Mathematica code and 1400 pages of C code.

Exact optimization

For linear cases, Minimize and Maximize use exact linear programming methods. For polynomial cases they use cylindrical algebraic decomposition.

Simplification

FullSimplify automatically applies about 40 types of general algebraic transformations, as well as about 400 types of rules for specific mathematical functions.

Generalized hypergeometric functions are simplified using about 70 pages of Mathematica transformation rules. These functions are fundamental to many calculus operations in Mathematica.

FunctionExpand uses an extension of Gauss's algorithm to expand trigonometric functions with arguments that are rational multiples of .

Simplify and FullSimplify cache results when appropriate.

When assumptions specify that variables are real, polynomial constraints are handled by cylindrical algebraic decomposition, while linear constraints are handled by the simplex algorithm or Loos-Weispfenning linear quantifier elimination. For strict polynomial inequalities, Strzebonski's generic CAD algorithm is used.

When assumptions involve equations among polynomials, Gröbner basis methods are used.

For non-algebraic functions, a database of relations is used to determine the domains of function values from the domains of their arguments. Polynomial-oriented algorithms are used whenever the resulting domains correspond to semi-algebraic sets.

For integer functions, several hundred theorems of number theory are used in the form of Mathematica rules.

Differentiation and integration

Differentiation uses caching to avoid recomputing partial results.

For indefinite integrals, an extended version of the Risch algorithm is used whenever both the integrand and integral can be expressed in terms of elementary functions, exponential integral functions, polylogarithms and other related functions.

For other indefinite integrals, heuristic simplification followed by pattern matching is used.

The algorithms in Mathematica cover all of the indefinite integrals in standard reference books such as Gradshteyn-Ryzhik.

Definite integrals that involve no singularities are mostly done by taking limits of the indefinite integrals.

Many other definite integrals are done using Marichev-Adamchik Mellin transform methods. The results are often initially expressed in terms of Meijer G functions, which are converted into hypergeometric functions using Slater's Theorem and then simplified.

Integrate uses about 500 pages of Mathematica code and 600 pages of C code.

Differential equations

Systems of linear equations with constant coefficients are solved using matrix exponentiation.

Second-order linear equations with variable coefficients whose solutions can be expressed in terms of elementary functions and their integrals are solved using the Kovacic algorithm.

Higher-order linear equations are solved using Abramov and Bronstein algorithms.

Systems of linear equations with rational function coefficients whose solutions can be given as rational functions are solved using Abramov-Bronstein elimination algorithms.

Linear equations with polynomial coefficients are solved in terms of special functions by using Mellin transforms.

When possible, nonlinear equations are solved by symmetry reduction techniques. For first-order equations classical techniques are used; for second-order equations and systems integrating factor and Bocharov techniques are used.

The algorithms in Mathematica cover most of the ordinary differential equations in standard reference books such as Kamke.

For partial differential equations, separation of variables and symmetry reduction are used.

For differential-algebraic equations, a method based on isolating singular parts by core nilpotent decomposition is used.

DSolve uses about 300 pages of Mathematica code and 200 pages of C code.

Sums and products

Polynomial series are summed using Bernoulli and Euler polynomials.

Series involving rational and factorial functions are summed using Adamchik techniques in terms of generalized hypergeometric functions, which are then simplified.

Series involving polygamma functions are summed using integral representations.

Dirichlet and related series are summed using pattern matching.

For infinite series, d'Alembert and Raabe convergence tests are used.

The algorithms in Mathematica cover at least 90% of the sums in standard reference books such as Gradshteyn-Ryzhik.

Products are done primarily using pattern matching.

Sum and Product use about 100 pages of Mathematica code.

Series and limits

Series works by recursively composing series expansions of functions with series expansions of their arguments.

Limits are found from series and using other methods.

Recurrence equations

RSolve solves systems of linear equations with constant coefficients using matrix powers.

Linear equations with polynomial coefficients whose solutions can be given as hypergeometric terms are solved using van Hoeij algorithms.

Systems of linear equations with rational function coefficients whose solutions can be given as rational functions are solved using Abramov-Bronstein elimination algorithms.

Nonlinear equations are solved by transformation of variables, Göktas symmetry reduction methods or Germundsson trigonometric power methods.

The algorithms in Mathematica cover most of the ordinary and -difference equations ever discussed in the mathematical literature.

For difference-algebraic equations, a method based on isolating singular parts by core nilpotent decomposition is used.