Recent Symbolic Computation Developments in Mathematica

Adam Strzeboński, Wolfram Research, Inc .

Symbolic computation Developments in 2020’s

Algebra

• Finite fields (FiniteField, FiniteFieldElement, FiniteFieldEmbedding, FrobeniusAutomorphism, FiniteFieldElementTrace, FiniteFieldElementNorm, FiniteFieldElementPrimitiveQ, MinimalPolynomial, MultiplicativeOrder)

• Performance improvements in polynomial algebra (PolynomialGCD, PolynomialLCM, Together, Cancel, Factor, FactorSquareFree, Resultant, Discriminant, Subresultants)

• Multivariate polynomial factorization over (Subscript[ℤ, p])

• Linear algebra for univariate rational function matrices over ℚ and (Subscript[ℤ, p]) (MatrixPower, Det, Inverse, LinearSolve, RowReduce, NullSpace, MatrixRank)

• Structured matrices (DiagonalMatrix, PermutationMatrix, LowerTriangularMatrix, UpperTriangularMatrix, BlockDiagonalMatrix, BlockLowerTriangularMatrix, BlockUpperTriangularMatrix, VandermondeMatrix, CauchyMatrix, ToeplitzMatrix, HankelMatrix, FourierMatrix, HilbertMatrix, TargetStructure)

• New linear algebra functions (Adjugate, CoreNilpotentDecomposition, DrazinInverse)

Solving Equations and inequalities

• System Root objects

• Solving systems of analytic function equations in bounded regions

• Finding roots of meromorphic functions in bounded regions

• Efficient representation of cylindrical algebraic formulas (CylindricalDecompositionFunction)

• Representing polynomials as sums of squares (PolynomialSumOfSquaresList)

• Performance improvements in polynomial real root isolation

• Assumptions in Solve

• SolveValues

Exact Global optimization

• Exact convex optimization

• Linear (MPLP) and quadratic (MPQP) parametric optimization

Calculus

• Fractional calculus (FractionalD, CaputoD, DSolve, AsymptoticDSolveValue)

• Implicit differentiation (ImplicitD)

• Curve and surface integration (LineIntegrate, SurfaceIntegrate)

• Contour integration (ContourIntegrate)

• Integral transforms (BilateralLaplaceTransform, InverseBilateralLaplaceTransform, BilateralZTransform, InverseBilateralZTransform)

• Sum of residues (ResidueSum)

• Change of variables (IntegrateChangeVariables, DSolveChangeVariables)

• Improvements in ODE, PDE, and OΔE solving (DSolve, DSolveValue, CompleteIntegral, RSolve)

• Improvements in integration (Integrate)

• Improvements in DifferentialRoot and DifferenceRoot

• Improvements in integral transforms (LaplaceTransform, InverseLaplaceTransform)

Asymptotics

• Asymptotics of functions and sequences (Asymptotic, DiscreteAsymptotic)

• Asymptotic statistics (AsymptoticProbability, AsymptoticExpectation)

Function Properties

• Set-theoretical properties (FunctionInjective, FunctionSurjective, FunctionBijective)

• Topological and analytic properties (FunctionContinuous, FunctionAnalytic, FunctionMeromorphic, FunctionDiscontinuities, FunctionSingularities, FunctionPoles)

• Inequality-related properties of real functions (FunctionSign, FunctionMonotonicity, FunctionConvexity)

Special functions

• HeunB, HeunBPrime, HeunC, HeunCPrime, HeunD, HeunDPrime, HeunG, HeunGPrime, HeunT, HeunTPrime

• JacobiEpsilon, JacobiZN

• LameC, LameCPrime, LameS,  LameSPrime, LameEigenvalueA, LameEigenvalueB

• FoxH and FoxHReduce

• CarlsonRC, CarlsonRD, CarlsonRE, CarlsonRF, CarlsonRG, CarlsonRJ, CarlsonRK, CarlsonRM

• CoulombF,CoulombG, CoulombH1, CoulombH2

• MittagLefflerE

• AppellF2, AppellF3, AppellF4

• CenteredInterval evaluation for 161 mathematical functions

roots of systems of equations

This finds an approximate root of a system of equations.

FindRoot[Gamma[x+y-2]==Sin[x]+Beta[x,y]+3&&SinIntegral[x y]==x-2,{{x,3.4},{y,1.9}},WorkingPrecision->20]

Can we prove that the result is an approximation of an actual root? The numeric method does not guarantee that.

FindRoot[(2x-y)^2+E^(-1000y)==0&&y^2==2,{{x,1},{y,2}},WorkingPrecision->20]

Center-radius interval arithmetic

({{, }, {real interval, complex interval}})

f[(Subscript[Δ, 1]),…,(Subscript[Δ, n])] yields a center-radius interval Δ which contains f[(Subscript[a, 1]),…,(Subscript[a, n])] for any (Subscript[a, i])∈(Subscript[Δ, i]).

{a,b}={CenteredInterval[1,10^-10],CenteredInterval[2+3I,10^-10+10^-10I]}

BesselJ[Log[a]+1,a/b^2]

CenteredInterval evaluation works for 161 mathematical functions. The implementation uses the Arb library.

Krawczyk operator

Let 𝕂 be ℝ or ℂ, let F:(𝕂^n)⟶(𝕂^n) be a (𝒞^1) mapping, let A be an invertible matrix in (𝕂^(n×n)), let x be an interval in 𝕂, and let (Subscript[x, 0])∈x. Define

K((Subscript[x, 0]),x,F)=(Subscript[x, 0])-A F((Subscript[x, 0]))+(I-A F'(x))(x-(Subscript[x, 0]))

where I is the identity matrix in (𝕂^(n×n)).

Theorem: 
(1) If y∈x and F(y)=0 then y∈K((Subscript[x, 0]),x,F).
(2) If K((Subscript[x, 0]),x,F)⊆x and (Subscript[‖I-A F'(x)‖, ∞])<α then F has a unique root in x, where α=1 if 𝕂 is ℝ and α=(1/(Sqrt[2])) if 𝕂 is ℂ .

In practice we choose (Subscript[x, 0]) to be the center of x, A to be an approximation of ((F'((Subscript[x, 0])))^(-1)), and use F((Subscript[x, 0])) instead of F((Subscript[x, 0])) where (Subscript[x, 0]) is an interval centered at (Subscript[x, 0]) with a radius much smaller than the radius of x.

Representing roots of systems of equations

Root[{{(Subscript[f, 1]),…,(Subscript[f, n])},{(Subscript[c, 1]),…,(Subscript[c, n])}},j] represents j-th coordinate of the unique root of the system of equations {(Subscript[f, 1])[(Subscript[x, 1]),…,(Subscript[x, n])]==0,…,(Subscript[f, n])[(Subscript[x, 1]),…,(Subscript[x, n])]==0} in the product of intervals/rectangles represented by {(Subscript[c, 1]),…,(Subscript[c, n])}.

f=Gamma[#1+#2-2]-Sin[#1]-Beta[#1,#2]-3&;
g=SinIntegral[#1 #2]-#1+2&;
{c,d}={x,y}/.FindRoot[f[x,y]==0&&g[x,y]==0,{{x,3.4},{y,1.9}},WorkingPrecision->20]

This uses interval arithmetic to prove existence and uniqueness of the roots.

r={Root[{{f,g},{c,d}},1],Root[{{f,g},{c,d}},2]}

The Root objects represent exact numbers.

Precision[r]

N[r,50]

Now let’s try an approximation that does not correspond to an actual root.

f=(2#1-#2)^2+E^(-1000#2)&;
g=#2^2-2&;
{c,d}={x,y}/.FindRoot[f[x,y]==0&&g[x,y]==0,{{x,1},{y,2}},WorkingPrecision->20]

Root[{{f,g},{c,d}},1]

Finding solution instances

FindInstance uses a method which finds an approximate solution using FindRoot and then proves that a box around the approximate solution contains a unique root using interval arithmetic.

A real solution:

FindInstance[Gamma[x+y+z]==Sin[x]+Beta[x,y]+3&&SinIntegral[x y]==x-Cos[z]+2&&BesselJ[2,x-y]≥3 z,{x,y,z},Reals]

A complex solution:

FindInstance[Sinh[x y]-Exp[I y]==0&&2 x Log[y]-y Log[x]==1,{x,y}]

Solving systems of equations in bounded regions

Interval arithmetic provides criteria for proving that a box contains a unique root and for proving that a box contains no roots. Hence if a system contains a finite number of isolated simple roots in a box, then the roots can be found using a subdivision method.

A real system:

Solve[Gamma[x+y+1]-Sin[x y]==1&&Erf[x^2-y]-E^y-x+4==0&&0<x<3&&0<y<3,{x,y}, Reals]

A complex system:

Reduce[SinhIntegral[x+y]-y AiryAi[x]^2==1+I&&x^2-BesselJ[3,x-y] -Cos[3x y]==-3+I&&0≤Re[x]≤1&&0≤Im[x]≤1&&0≤Re[y]≤1&&0≤Im[y]≤1,{x,y}]//AbsoluteTiming

roots of univariate meromorphic functions

The algorithm

The following method finds the roots of a univariate meromorphic function in a bounded box:
1) Find approximate roots using a numeric method.
2) Prove existence and uniqueness of each root.
3) Prove that all roots have been found by counting the number of roots in the box.

Counting the number of roots can be done by computing the winding number of the function along the boundary of the box and by computing the number of poles.

The winding number can be computed using interval arithmetic.

Example

rts=Solve[Gamma[x^2]==Sin[Cos[x]]&&-3/2<=Re[x]<=3/2&&-3/2<=Im[x]<=3/2,x]

ComplexPlot[Gamma[x^2]-Sin[Cos[x]],{x,-1.5-1.5I,1.5+1.5I},Epilog->{Blue,Point[ReIm[x]/.rts]}]

FunctionPoles[{Gamma[x^2]-Sin[Cos[x]],-3/2<=Re[x]<=3/2&&-3/2<=Im[x]<=3/2},x]

Exact convex optimization

Convex optimization problems

Find the global minimum of a convex function f over the set

({{({{C={x∈}})(ℝ^n): (Subscript[g, 1])(x)≤0∧…∧(Subscript[g, m])(x)≤0∧(Subscript[h, 1])(x)==0∧…∧(Subscript[h, p])(x)==0}}})

where  (Subscript[g, 1]),…,(Subscript[g, m]) are convex functions and (Subscript[h, 1]),…,(Subscript[h, p]) are affine functions.

KKT conditions for convex optimization problems

Theorem (Sufficiency of KKT conditions)

Suppose that f,(Subscript[g, 1]),…,(Subscript[g, m]) are continuously differentiable over (ℝ^n), 
and there exist a∈C, (Subscript[λ, 1]),…,(Subscript[λ, m])≥0 , and (Subscript[μ, 1]),…,(Subscript[μ, p])∈ℝ such that

({{({{({{ (Subscript[∇, x])f(a)+(∑_(i=1)^m)(Subscript[λ, i])(Subscript[∇, x])(Subscript[g, i])(a)+(∑_(j=1)^p)(Subscript[μ, j])(Subscript[∇, x])(Subscript[h, j])(a)==0 }})
(Subscript[∀, 1≤i≤m])(Subscript[λ, i])(Subscript[g, i])(a)==0  }})}})

Then f(a)  is the global minimum of f on C.

Solving the KKT conditions

Functions f,(Subscript[g, 1]),…,(Subscript[g, m]) need not to be polynomial, so we cannot solve the KKT conditions directly using symbolic methods. However, we can solve convex optimization problems numerically. The method finds an exact minimum point by using an approximate minimum point as a hint and then proving that the obtained exact point satisfies the KKT conditions. The algorithm proceeds as follows.

Find an approximate minimum point (OverTilde[a]).

From the approximate solution guess the set I⊆{1,…,m} such that (Subscript[∀, i∈I])(Subscript[g, i])(a)==0.

Find an exact real solution (a,λ,μ) with a near (OverTilde[a]) of the square system of equations

({{({{({{ (Subscript[∇, x])f(x)+(∑_(i∈I))(Subscript[λ, i])(Subscript[∇, x])(Subscript[g, i])(x)+(∑_(j=1)^p)(Subscript[μ, j])(Subscript[∇, x])(Subscript[h, j])(x)==0 }})
(Subscript[∀, i∈I])(Subscript[g, i])(x)==0 
(Subscript[∀, 1≤j≤p])(Subscript[h, j])(x)==0 }})}})

Verify that (Subscript[g, i])(a)≤0 for i∉I and (Subscript[λ, i])≥0 for i∈I.

Example

Find an exact solution to a convex optimization problem:

{min,pt}=Minimize[{Log[(E^(x+y))+(E^(2 x-y))],{x,y,x+y+1}(Underscript[\[VectorGreaterEqual], (ExponentialConeString)])0},{x,y}]

Show[Plot3D[Log[(E^(x+y))+(E^(2 x-y))],{x,-2,(5/4)},{y,-(1/4),1},MeshFunctions->{#3&},RegionFunction->Function[{x,y},{x,y,x+y+1}(Underscript[\[VectorGreaterEqual], (ExponentialConeString)])0]],Graphics3D[{Red,PointSize[0.03],Point[N[{x,y,min}/. pt]]}]]

Finite fieldS

Representation

• FiniteField

• FiniteFieldElement

Element Operations

• Arithmetic (Plus, Times, integer Power)

• FiniteFieldElementTrace

• FiniteFieldElementNorm

• FiniteFieldElementPrimitiveQ

• MinimalPolynomial

• MultiplicativeOrder

Finite field Mappings

• FrobeniusAutomorphism

• FiniteFieldEmbedding

Polynomials over Finite fields

PolynomialGCD, PolynomialLCM, PolynomialExtendedGCD, Expand, Factor, FactorList, FactorSquareFree, FactorSquareFreeList, IrreduciblePolynomialQ, PolynomialQuotientRemainder, PolynomialQuotient, PolynomialRemainder, Together, Cancel, Resultant, Discriminant

COMING Soon

• Linear algebra over finite fields

• Solving systems of equations and inequations

Finite field Representation

A finite field ff is represented by FiniteField[p, f, rep], where p is the characteristic of ff, f is a polynomial function irreducible over (Subscript[ℤ, p]), and rep is “Polynomial” or “Exponential”.

ff has (p^d) elements, where d is the degree of f, indexed by integers between 0 and (p^d)-1.

This constructs a finite field with (17^3) elements:

ff=FiniteField[17,3]

This gives the element with index 123:

a=ff[123]

Field elements are represented by atomic FiniteFieldElement objects:

{Head[a],AtomQ[a]}

Elements of ff can be expressed as polynomials of degree at most d-1 in the field generator θ, where f(θ)=0. This specifies a field element through a coefficient list of a polynomial in θ:

b=ff[{1,2,3}]

Extract information about field and field elements:

Information[ff,{"Characteristic","FieldIrreducible","FieldSize"}]

Information[a,{"Index","Field"}]

a["Index"]

“Polynomial” Representation

In the default "Polynomial" representation the index of the field element c=(Subscript[c, 0])+(Subscript[c, 1])θ+⋯+(Subscript[c, d-1])(θ^(d-1)) is ind = (Subscript[c, 0])+(Subscript[c, 1])p+⋯+(Subscript[c, d-1])(p^(d-1)):

c=ff[{3,5,7}]

ind = {3,5,7}.(17^Range[0,2])

Find the coefficients from the index:

Reverse[IntegerDigits[ind,17]]

“Exponential” Representation

This constructs a finite field with (23^4) elements in the “Exponential” representation:

ff=FiniteField[23,4,"Exponential"]

In the "Exponential" representation the field irreducible is chosen to be a primitive polynomial and hence each nonzero element of the field is a power of the field generator. The index of the field element (θ^k), 0≤k≤(p^d)-2, is k+1. The index of the field zero is 0:

θ=ff[{0,1}]

θ^123

θ^5+θ^7

Each element still can be represented as a polynomial in θ of degree at most d. Find the coefficients:

ff[1234]["Coefficients"]

%.θ^Range[0,3]

The "Exponential" representation provides somewhat faster finite field element arithmetic than the "Polynomial" representation, at the cost of a pre-computation necessary when the field is created. The precomputed table has size proportional to the number of field elements, hence the “Exponential” representation can be used only for relatively small fields ((p^d)≤(10^9)).

FiniteField[11,10,"Exponential"]

FiniteField[11,10]

Finite field Operations

ff=FiniteField[7,5]

{a,b} ={ff[123],ff[456]}

Do arithmetic in the finite field:

(a^2+2b)/(5a b-7)

Integers and rational numbers valid mod p are automatically treated as field elements:

a+3

Compute finite field element properties:

FiniteFieldElementTrace[a]

MinimalPolynomial[a,x]

MultiplicativeOrder[b]

Polynomials oveR Finite fieldS

Factor[x^5+2x+11,Extension->ff]

PolynomialRemainder[x^7-1,a x^3+b x+ a b,x]

Discriminant[a x^3+ a^2 b x y+a-b,x]

Cancel[(x-a)/MinimalPolynomial[a, x]]

Finite field Mappings

The Frobenius automorphism is given by ϕ:ff∋a⟶(a^p)∈ff. All automorphisms of ff are functional powers of the Frobenius automorphism (the first d powers are different).

ff=FiniteField[23,4]

{a,b} ={ff[123],ff[456]}

Automorphisms preserve field operations:

FrobeniusAutomorphism[a]+FrobeniusAutomorphism[b]==FrobeniusAutomorphism[a+b]

FrobeniusAutomorphism[a] FrobeniusAutomorphism[b]==FrobeniusAutomorphism[a b]

Find all conjugates of a field element:

Table[FrobeniusAutomorphism[a,k],{k,0,4}]

SameQ@@(MinimalPolynomial/@%)

Construct a larger field:

gg=FiniteField[23,12]

ff is isomorphic to a subfield of gg, but not automatically considered a subfield of gg, because there is no “standard” embedding of ff in gg:

a+gg[123]

Construct an embedding of ff in gg:

emb=FiniteFieldEmbedding[ff,gg]

emb[a]

%+gg[123]

Embeddings preserve field operations:

{emb[a]+emb[b]==emb[a+b], emb[a] emb[b]==emb[a b]}

Polynomial Algebra Performance improvements

The FLINT library provides efficient implementations of polynomial algebra algorithms within a fixed polynomial ring. The available rings are univariate polynomials (dense representation) and multivariate polynomials (sparse distributed representation) over the integers, the rationals, and the integers modulo p (separate implementations for machine-sized integers and big integers). The library also provides linear algebra functions for matrices of univariate polynomials over the integers and the integers modulo p for machine-sized p.

Direct improvements

These functions recognize cases that can be handled by FLINT library functions (more efficiently than by the Mathematica internal code) and call FLINT to do the computation.

• PolynomialGCD

• Resultant

• FactorSquareFree

• Factor

• MatrixPower

• Det

• Inverse

• LinearSolve

• RowReduce

• NullSpace

• MatrixRank

Indirect improvements

These functions use the directly improved functions for a large part of their computation or call FLINT functions to do some parts of the computation.

• PolynomialLCM

• Together

• Discriminant

• Subresultants (uses FLINT for the coefficient arithmetic)

• Polynomial real root isolation (uses FLINT Taylor shift implementation)

• All symbolic computation functions use polynomial algebra

performance improvements

Example

Generate a 10x10 matrix of randomly generated univariate polynomials of degree 100:

SeedRandom[1234];
rpoly[n_]:=RandomInteger[{-2^10,2^10},{n+1}].x^Range[0,n]
m=Table[rpoly[100],{10},{10}];

Compute the determinant of this matrix in Mathematica 13.3:

((Subscript[d, new])=Det[m]);//AbsoluteTiming

LeafCount[(Subscript[d, new])]

Compare with the results in Mathematica 13.1:

((Subscript[d, old])=Det[m]);//AbsoluteTiming

LeafCount[(Subscript[d, old])]

Benchmark plots

PolynomialGCD[p*r, q*r], with p, q, r randomly generated dense univariate polynomials of degree n with integer coefficients between -1000 and 1000.

Factor[p*q], with p, q randomly generated dense univariate polynomials of degree n with integer coefficients between -1000 and 1000.

LinearSolve[m, b], where m is a n x n matrix and b is a n vector of randomly generated univariate polynomials degree 10 with integer coefficients between -1000 and 1000.

Polynomial Algebra Functionality extensions

Factor

Factor now works for multivariate polynomials over (Subscript[ℤ, p]).

Expand[(x+y+x y)(x^2+x y+1),Modulus->2]

Factor[%,Modulus->2]

Linear algebra

Det, Inverse, LinearSolve, RowReduce, NullSpace, and MatrixRank now work for matrices of univariate rational functions over (Subscript[ℤ, p]) (for machine-sized prime p).

LinearSolve[{{2x+3,x^3-7},{x/(x-1),(x+5)/(x^2-3x+2)}},{x^2+1,x^3+2},Modulus->17]

Together[{{2x+3,x^3-7},{x/(x-1),(x+5)/(x^2-3x+2)}}.%-{x^2+1,x^3+2},Modulus->17]

Det now works for matrices of multivariate rational functions over (Subscript[ℤ, p]) .

Det[{{x+y,x-y},{x^2-5y,y^2-7x}},Modulus->1009]

Thank you for listening!