New algebraic functionality in Mathematica 6

Adam Strzebonski

Wolfram Research, Inc.

Abstract

In the upcoming version of Mathematica many functions formerly available through add-on packages have been redesigned, extended, often reimplemented, and are included in the kernel. In my talk I will discuss new functionality in polynomial algebra, algebraic number and number field computations, equation an inequality solving and exact global optimization.

Equations and Inequalities

Cylindrical algebraic decomposition improvements

(Affected functions: CylindricalDecomposition, Reduce, Resolve, FindInstance, Minimize, Maximize, all functions using assumptions.)

Several improvements have been added in the Mathematica implementation of the CAD algorithm, resulting in substantial speed-up for some examples.

This example (from H.Hong, R.Liska, S.Steinberg, Testing Stability by Quantifier Elimination, Journal of Symbolic Computation, 24, 1997) takes ~9000 seconds in Mathematica 5.2.

Reduce[(Subscript[∀, {a,b,c},(a|b|c)∈Reals&&0≤a≤1&&0≤b≤1])((-1+a-b) (((a-b))^2) (1+a-b) (c^4)≤0&&(a^2) (-1+(a^2))+4 (-1+a) (a^2) b c+2 a b (1-2 a-2 b+3 a b) (c^2)+4 a (-1+b) (b^2) (c^3)+(b^2) (-1+(b^2)) (c^4)≤0&&(4 a-3 (a^2)-4 (a^3)+4 (a^4)-4 b+2 a b+8 (a^2) b-8 (a^3) b+5 (b^2)-12 a (b^2)+8 (a^2) (b^2)+2 ((a^2)-2 (a^3)+2 a b-2 (a^2) b+4 (a^3) b+(b^2)-2 a (b^2)-4 (a^2) (b^2)-2 (b^3)+4 a (b^3)) c+(-4 a+5 (a^2)+4 b+2 a b-12 (a^2) b-3 (b^2)+8 a (b^2)+8 (a^2) (b^2)-4 (b^3)-8 a (b^3)+4 (b^4)) (c^2)≤0||2 a (-a+(a^3)+b-2 (a^2) b-2 (b^2)+3 a (b^2)) (c^2)+4 a b (-a+(a^2)-b+(b^2)) (c^3)+2 b (a-2 (a^2)-b+3 (a^2) b-2 a (b^2)+(b^3)) (c^4)≤0))]//Timing

This example (from P. Dorato, W. Yang, C. Abdallah, Robust Multi-Objective Feedback Design by Quantifier Elimination, Journal of Symbolic Computation, 24, 1997) couldn't be done with Mathematica 5.2 in 24 hours. The paper authors report being able to check only special cases with fixed nbyd variable.

CylindricalDecomposition[(Subscript[∃, {q1,q2},q1>1&&q2>0])((Subscript[∀, w])(nbyd>0&&(nbyd-(q1^2)) (w^4)+(nbyd ((((1+q1))^2)-2 q2)-((q2^2)+(q1^2))) (w^2)+(nbyd-1) (q2^2)≥0&&(nbyd-(q1^2)) (w^4)+(nbyd ((((-1+q1))^2)-2 q2)-((q2^2)+(q1^2))) (w^2)+(nbyd-1) (q2^2)≥0)),nbyd]//Timing

Exact global optimization improvements

Parametric optimization

Minimize and Maximize can handle polynomial optimization problems involving parameters. The minimum value and the minimum point are usually piecewise functions of the parameters. The method used is based on the cylindrical algebraic decomposition algorithm.

Minimize[a (x^2)+b x+c,x]

Maximize[{x y,a (x^2)+b (y^2)==1},{x,y}]

Linear fractional programming

Optimization problems with linear constraints and a rational objective function with linear numerator and denominator can be reduced to linear programming problems. Without the reduction this example was solved by CAD and took 71 seconds.

Minimize[{(11+(Subscript[x, 1])+2 (Subscript[x, 2])+3 (Subscript[x, 3])+4 (Subscript[x, 4])+5 (Subscript[x, 5])+6 (Subscript[x, 6])+7 (Subscript[x, 7])+8 (Subscript[x, 8])+9 (Subscript[x, 9])+10 (Subscript[x, 10]))/(111+101 (Subscript[x, 1])+102 (Subscript[x, 2])+103 (Subscript[x, 3])+104 (Subscript[x, 4])+105 (Subscript[x, 5])+106 (Subscript[x, 6])+107 (Subscript[x, 7])+108 (Subscript[x, 8])+109 (Subscript[x, 9])+110 (Subscript[x, 10])),(Subscript[x, 1])≥1&&(Subscript[x, 2])≥2&&(Subscript[x, 3])≥3&&(Subscript[x, 4])≥4&&(Subscript[x, 5])≥5&&(Subscript[x, 6])≥6&&(Subscript[x, 7])≥7&&(Subscript[x, 8])≥8&&(Subscript[x, 9])≥9&&(Subscript[x, 10])≥10&&(Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])==100},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8]),(Subscript[x, 9]),(Subscript[x, 10])}]//Timing

Integer linear programming

A heuristic has been added to solve integer linear programming problems with unbounded constraints for which the real minimum is finite.

Minimize[{((46 (Subscript[x, 1]))/5)+((44 (Subscript[x, 2]))/5)+((44 (Subscript[x, 3]))/5)+((42 (Subscript[x, 4]))/5)+8 (Subscript[x, 5])+8 (Subscript[x, 6])+((42 (Subscript[x, 7]))/5)+((44 (Subscript[x, 8]))/5)+((44 (Subscript[x, 9]))/5)+((46 (Subscript[x, 10]))/5)+((48 (Subscript[x, 11]))/5)+((48 (Subscript[x, 12]))/5),(Subscript[x, 1])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 12])≥3&&(Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 10])+(Subscript[x, 11])≥3&&(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 11])+(Subscript[x, 12])≥8&&(Subscript[x, 1])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 12])≥8&&(Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 4])+(Subscript[x, 5])≥10&&(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 5])+(Subscript[x, 6])≥10&&(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 6])+(Subscript[x, 7])≥8&&(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 7])+(Subscript[x, 8])≥8&&(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 8])+(Subscript[x, 9])≥14&&(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 9])+(Subscript[x, 10])≥14&&(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 10])+(Subscript[x, 11])≥5&&(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 11])+(Subscript[x, 12])≥5&&((Subscript[x, 1])|(Subscript[x, 2])|(Subscript[x, 3])|(Subscript[x, 4])|(Subscript[x, 5])|(Subscript[x, 6])|(Subscript[x, 7])|(Subscript[x, 8])|(Subscript[x, 9])|(Subscript[x, 10])|(Subscript[x, 11])|(Subscript[x, 12]))∈Integers},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8]),(Subscript[x, 9]),(Subscript[x, 10]),(Subscript[x, 11]),(Subscript[x, 12])}]//Timing

GenericCylindricalDecomposition

(Formerly Experimental`GenericCylindricalAlgebraicDecomposition.)

GenericCylindricalDecomposition runs a simplified version of the CAD algorithms, finding the solution set of inequalities correct up to lower dimensional components.

ineqs=(x^2)+(y^2)≤1&&(((x-(1/2)))^2)+(((y-(1/2)))^2)≤1;

CylindricalDecomposition[ineqs,{x,y}]

GenericCylindricalDecomposition[ineqs,{x,y}]

Show[RegionPlot[ineqs,{x,-1,1},{y,-1,1}],ContourPlot[#,{x,-1,1},{y,-1,1}]&/@List@@%[[2]]]

Now let us try a more complicated example.

ineqs=(x^2)+(y^2)+(z^2)-3 x y z<1&&(x^2)+(y^2)≥2 (z^3)&&(x^4)+(y^4)+(z^4)≤3;

RegionPlot3D[ineqs,{x,-2,2},{y,-2,2},{z,-2,2}]

CylindricalDecomposition[ineqs,{x,y,z}];//Timing

GenericCylindricalDecomposition[ineqs,{x,y,z}];//Timing

SemialgebraicComponents

(Replaces and extends SemialgebraicComponents from the Algebra`AlgebraicInequalities` package.)

SemialgebraicComponents gives at least one point in each connected component of a semialgebraic set. The package function worked only for solutions of strict inequalities, the new kernel function works for arbitrary semialgebraic sets.

eq=(x^3)+(((2 z y))^2)-x-(1/4)==0;
pts=SemialgebraicComponents[eq&&-2<x<2&&-2<y<2&&-2<z<2,{x,y,z}];

Show[{ContourPlot3D[Evaluate[eq],{x,-2,2},{y,-2,2},{z,-2,2}],Graphics3D[{Red,PointSize[0.025],Point[{x,y,z}/.pts]}]}]

ineqs=(x^2)+(y^2)+(z^2)-7 x y z<1&&(x^2)+(y^2)>2 (z^3)+(5/4)&&(x^2)+(y^2)+(z^2)<3;
pts=SemialgebraicComponents[ineqs,{x,y,z}];

And@@(ineqs/.pts)

Show[{RegionPlot3D[ineqs,{x,-2,2},{y,-2,2},{z,-2,2}],Graphics3D[{Red,PointSize[0.015],Point[{x,y,z}/.pts]}]}]

CountRoots

(Replaces CountRoots from the Algebra`RootIsolation` package.)

CountRoots now counts roots in closed intervals (rectangles).

poly=(x^17)-1;

Graphics[{{GrayLevel[0.7],Rectangle[{0,0},{1,1}]},{Red,PointSize[0.02],Point[{Re[x],Im[x]}/.NSolve[poly==0,x]]}}]

CountRoots[poly,{x,0,1+I}]

RootIntervals

(Replaces RealRootIntervals and ComplexRootIntervals from the Algebra`RootIsolation` package.)

This finds isolating intervals for real roots of polynomials. The second list indicates which root belongs to which polynomials. For instance the second element {1, 2, 2} means the root in the second interval ({0, 2}) is a simple root of the first polynomial and a double root of the second polynomial.

RootIntervals[{(x^3)-1,(((x-1))^2),(x^3)-2 x+2}]

This finds isolating intervals for all complex roots.

RootIntervals[{(x^3)-1,(((x-1))^2),(x^3)-2 x+2},Complexes]

IsolatingInterval

(Replaces and extends ContractInterval from the Algebra`RootIsolation` package.)

Gives an isolating interval for an algebraic number, i.e. an interval which does not contain other roots of the minimal polynomial.

IsolatingInterval[(Sqrt[2])]

This gives an isolating interval of width at most (10^(-10)).

IsolatingInterval[(Sqrt[2]),(1/(10^10))]

FrobeniusNumber

(Replaces FrobeniusF from the NumberTheory`Frobenius` package.)

the Frobenius number of (Subscript[a, 1]),...,(Subscript[a, n]) is the largest integer k for which the Frobenius equation

(Subscript[a, 1])(Subscript[x, 1])+...+(Subscript[a, n])(Subscript[x, n]) ==k

has no nonnegative integer solutions.

This gives the highest amount of postage (in cents) that cannot be paid with postcard, letter and international stamps.

FrobeniusNumber[{24,39,80}]

FrobeniusNumber can handle quite large inputs: million-digit coefficients for n<=3, thousand-digit coefficients for n<=10, large n with coefficients <= 10^6

N[FrobeniusNumber[Table[RandomInteger[{1,(10^(10^6))}],{3}]]]//Timing

N[FrobeniusNumber[Table[RandomInteger[{1,(10^(10^3))}],{5}]]]//Timing

FrobeniusNumber[Table[RandomInteger[{1,(10^6)}],{1000}]]//Timing

FrobeniusSolve

(Replaces FrobeniusSolve and FrobeniusInstance from the NumberTheory`Frobenius` package.)

This gives all ways to pay $3.99 postage with postcard, letter and international stamps.

FrobeniusSolve[{24,39,80},399]

Algebraic Numbers

AlgebraicNumber

(Replaces and extends Algebraic from the NumberTheory`AlgebraicNumberFields` package.)

AlgebraicNumber[θ,{(Subscript[c, 0]),(Subscript[c, 1]),…,(Subscript[c, n])}]represents the algebraic number if the field ℚ[θ] given by (Subscript[c, 0])+(Subscript[c, 1])θ +…+(Subscript[c, n]) (θ^n).

AlgebraicNumber[(Sqrt[2]),{1,2,3,4}]

RootReduce[%]

The field generator returned is always an algebraic integer.

a=AlgebraicNumber[Root[2 ((#1)^5)+3 #1-7&,1],{1,2,3,4,5}]

Arithmetic with AlgebraicNumber objects with the same field generators is performed automatically.

b=AlgebraicNumber[a[[1]],{9,8,7,6,5}]

(((a^2)+3 a (b^3)-11)/((a^5)-(b^5)))

AlgebraicNumber objects are numeric.

N[%,20]

RootReduce converts AlgebraicNumber objects to Root objects.

RootReduce[%%]

MinimalPolynomial

(Replaces and extends MinimalPolynomial and NumberFieldCharacteristicPolynomial from the NumberTheory`AlgebraicNumberFields` package.)

Gives the minimal polynomial of an algebraic number.

MinimalPolynomial[(Sqrt[2]),x]

MinimalPolynomial[(Sqrt[2]),x,Extension->(Sqrt[2])+(Sqrt[3])]

AlgebraicNumberNorm

(Replaces and extends AlgebraicNumberNorm and NumberFieldNorm from the NumberTheory`AlgebraicNumberFields` package.)

Gives the norm of an algebraic number.

AlgebraicNumberNorm[(Sqrt[2])]

Norm of an algebraic number relative to a specified number field.

AlgebraicNumberNorm[(Sqrt[2]),Extension->(Sqrt[2])+(Sqrt[3])]

AlgebraicNumberTrace

(Replaces and extends AlgebraicNumberTrace and NumberFieldTrace from the NumberTheory`AlgebraicNumberFields` package.)

Gives the trace of an algebraic number.

AlgebraicNumberTrace[Root[((#1)^5)-3 ((#1)^4)+11&,1]]

AlgebraicIntegerQ

(Replaces and extends AlgebraicIntegerQ from the NumberTheory`AlgebraicNumberFields` package.)

AlgebraicIntegerQ[GoldenRatio]

AlgebraicNumberDenominator

(Replaces and extends AlgebraicNumberDenominator from the NumberTheory`AlgebraicNumberFields` package.)

AlgebraicNumberDenominator[(1/4) (1+(Sqrt[5]))]

AlgebraicUnitQ

(Replaces and extends AlgebraicUnitQ from the NumberTheory`AlgebraicNumberFields` package.)

AlgebraicUnitQ[GoldenRatio]

RootOfUnityQ

(Replaces and extends RootOfUnityQ from the NumberTheory`AlgebraicNumberFields` package.)

RootOfUnityQ[GoldenRatio]

RootOfUnityQ[-(1/2)+((I (Sqrt[3]))/2)]

RootApproximant

(New, extends the functionality of Recognize from the NumberTheory`Recognize` package.)

Finds an algebraic number close to a given inexact number.

N[(Sqrt[2])+(Sqrt[3])]

RootApproximant[%,4]

RootApproximant[%%]

RootReduce[(Sqrt[2])+(Sqrt[3])]

RootApproximant[N[π],5]

N[%-π]

Algebraic Number Fields

ToNumberField

(Replaces and extends ToCommonField and ToNumberFieldElement of NumberTheory`AlgebraicNumberFields`.)

ToNumberField[{(Sqrt[2]),(Sqrt[3])}]

ToNumberField[(Sqrt[2]),(2^(1/8))]

NumberFieldIntegralBasis

(Replaces IntegralBasis from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldIntegralBasis[(Sqrt[2])+(Sqrt[3])]

NumberFieldRootsOfUnity

(Replaces RootsOfUnity from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldRootsOfUnity[(Sqrt[3])+I]

RootReduce[%]

ToRadicals[%]

ToRadicals[%%%]

NumberFieldFundamentalUnits

(Replaces FundamentalUnits from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldFundamentalUnits[(Sqrt[3])+I]

NumberFieldNormRepresentatives

(Replaces NormRepresentatives from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldNormRepresentatives[Root[((#1)^3)-2 #1+7&,1],3]

AlgebraicNumberNorm[%,Extension->Automatic]

NumberFieldSignature

(Replaces NumberFieldSignature from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldSignature[(Sqrt[2])+(Sqrt[3])]

NumberFieldSignature[((#1)^7)-11 #1+9&]

NumberFieldDiscriminant

(Replaces NumberFieldSignature from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldDiscriminant[(Sqrt[2])+(Sqrt[3])]

NumberFieldDiscriminant[((#1)^7)-11 #1+9&]

b=NumberFieldIntegralBasis[Root[((#1)^7)-11 #1+9&,1]]

Det[Table[AlgebraicNumberTrace[b[[i]] b[[j]],Extension->Root[((#1)^7)-11 #1+9&,1]],{i,7},{j,7}]]

Since here ℤ[θ]==(Subscript[ℤ, ℚ[θ]]), we can compute the discriminant of ℚ[θ] as

Discriminant[(((#1)^7)-11 #1+9&)[x],x]

NumberFieldRegulator

(Replaces NumberFieldRegulator from NumberTheory`AlgebraicNumberFields` package.)

NumberFieldRegulator[(Sqrt[2])+(Sqrt[3])]

Polynomial Algebra

Discriminant

(New.)

The discriminant is zero iff the polynomial has multiple roots.

a=.;b=.;c=.;

Discriminant[a (x^2)+b x+c,x]

R(f,f')==(((-1))^(n(n-1)/2))(Subscript[a, n])D(f).

f=3 (x^4)-4 (x^3)-2 (x^2)-x+77;

d=Discriminant[f,x]

r=Resultant[f,D[f,x],x]

r/d

D(f)==((Subscript[a, n])^(2n-2))(∏_(1≤i<j≤n))((((Subscript[r, i])-(Subscript[r, j])))^2).

rts=x/.{ToRules[Reduce[f==0,x]]};

RootReduce[Product[(rts[[i]]-rts[[j]])^2,{i,4},{j,i-1}]]

(3^(2 4-2)) %==d

SymmetricPolynomial

(Formerly in the Algebra`SymmetricPolynomials` package.)

SymmetricPolynomial[3,{x,y,z,t}]

SymmetricReduction

(Formerly in the Algebra`SymmetricPolynomials` package.)

SymmetricReduction[(x^3)+(y^3)+3 x y z,{x,y,z}]

Plus@@%//Expand

SymmetricReduction[(x^3)+(y^3)+3 x y z,{x,y,z},{(Subscript[s, 1]),(Subscript[s, 2]),(Subscript[s, 3])}]

PolynomialExtendedGCD

(Replaces PolynomialExtendedGCD from the Algebra`PolynomialExtendedGCD` package.)

PolynomialExtendedGCD[2 (x^5)-2 x,((((x^2)-1))^2),x]

%[[2]].{2 (x^5)-2 x,((((x^2)-1))^2)}//Expand

PolynomialPowerMod

(Replaces PolynomialPowerMod from the Algebra`PolynomialPowerMod` package.)

PolynomialPowerMod[(x^2)-2 x+5,7,{2 (x^3)-3 x+9,11}]

PolynomialMod[PolynomialRemainder[((((x^2)-2 x+5))^7),2 (x^3)-3 x+9,x],11]

PolynomialPowerMod can handle very high powers.

PolynomialPowerMod[(x^2)-2 x+5,1000000,{2 (x^3)-3 x+9,11}]

TimeConstrained[PolynomialMod[PolynomialRemainder[((((x^2)-2 x+5))^1000000),2 (x^3)-3 x+9,x],11],60]

PolynomialQuotientRemainder

(Formerly Developer`PolynomialDivision.)

f=(x^5)-3 x+9;
g=(x^3)-7 x+11;
{q,r}=PolynomialQuotientRemainder[f,g,x]

Expand[f==q g+r]

PolynomialQuotientRemainder (as well as PolynomialQuotient and PolynomialRemainder) has been extended to handle rational functions, according to the following definition.

PolynomialQuotientRemainder[(f/(d h)), ((d g)/k), x] := {((k q)/((d^2) h)),(r/d)}, such that  (f/(d h))==(((d g) (k q))/(((d^2) h) k))+(r/d), and deg r<deg g.If gcd (f, d h) == gcd (d g, k) == gcd (g h) == 1, q and r are unique.

This somewhat technical definition has a useful property that if gcd (f, g) == 1, PolynomialRemainder[(1/f), g, x] gives the inverse of f modulo g.

f=(x^3)+3 x+4;
g=(x^3)-5 x+7;
invf=PolynomialRemainder[(1/f),g,x]

PolynomialRemainder[f invf,g,x]