Testing Symbolics and Special Functions

Bhuvanesh Bhatt

Wolfram Research

Introduction

While I do some development work (complex analysis, DSolve), my main duties are testing Mathematica's symbolics functions and special functions, with roughly half of my time spent on definite Integrate. Today I will be talking about the development and maintenance of our automated tests for symbolics and special functions.

Writing tests

Verifying results

Analyzing test failures

Improving tests

Reorganizing tests

Background: Structure of a test

Input

Expected output

Expected messages

Unique test ID

Test[
	1/0
	,
	ComplexInfinity
	,
	Power::infy
	,
	TestID->"Unique identifier"
	]

Test sources

MathGroup, other newsgroups

Fixed bugs and implemented suggestions

Functions website (functions.wolfram.com)

MathWorld (mathworld.wolfram.com)

Benchmarks (Kamke, Gradshteyn-Ryzhik, Prudnikov-Brychkov-Marichev, Demidovich, Vavilov, etc.)

Journal articles (e.g. articles from arxiv.org or CA articles)

Wolfram Research staff (including publications)

Developer tests (including conference talks)

QA tests (from test engineer)

Well-known examples

Code coverage (currently over 85% of C code for my areas is covered)

Test sources

MathGroup, other newsgroups

Test[
	"From a November 2003 MathGroup thread. This used to give a wrong result; fixed in 5.1.";
	Integrate[Log[2, 1 + 10*x]*Exp[-x]*x, {x, 0, Infinity}] //Simplify
	,
	(1 - (9*E^(1/10)*ExpIntegralEi[-1/10])/10)/Log[2]
	,
	TestID->"Integrate_MathGroup-20050326-6rub1u"
	]

Test[
	"There was an inconsistency, pointed out by Maxim Rytin on MathGroup, between the formulas
	 used by complex and real Reduce for f/g != h. Both now use f/g != 0 -> f!=0 && g!=0.";
	Reduce[ForAll[x, 1/(x + a) != 0]]
	,
	False
	,
	TestID->"Reduce_MathGroup-20050923-eshz1u"
	]

Test sources

Fixed bugs and implemented suggestions

Test[
	"Stats integral: Variance of beta distribution";
	integrand = (x - alpha/(alpha + beta))^2*((1 - x)^(-1 + beta)*x^(-1 + alpha))/Beta[alpha, beta];
	Integrate[integrand, {x, 0, 1}, Assumptions->alpha>0 && beta>0]
	,
	(alpha*beta)/((alpha + beta)^2*(1 + alpha + beta))
	,
	TestID->"Bugs-59000-20050920-awhr1m-bug-59779"
	]

Test[
	"This limit was one of the causes of bug 61140";
	Limit[CosIntegral[-I - x + I/x], x->Infinity]
	,
	-I Pi
	,
	TestID->"Bugs-61000-20050921-nt721e-bug-61140"
	]

Test sources

MathWorld

Test[
	"Steinmetz solid for two cylinders, CAD to find volume";
	CylindricalDecomposition[x^2 + y^2 < 1 && -2 < z < 2 && y^2 + z^2 < 1 && -2 < x < 2, {x, y, z}]
	,
	-1 < x < 1 && -Sqrt[1 - x^2] < y < Sqrt[1 - x^2] && -Sqrt[1 - y^2] < z < Sqrt[1 - y^2]
	,
	TestID->"CAD_MathWorldExamples-20030506-3u1v6f"
	]

Test sources

Functions website

Test[
	"MathieuCharacteristicB[r, q] == MathieuCharacteristicA[r, q] /; r ∉ Integers && q ∈ Reals";
	"http://functions.wolfram.com/11.06.27.0001.01";
	inst = FindInstance[NotElement[r, Integers] && Element[r|q, Reals], {r, q}, 5];
	Select[inst, N[MathieuCharacteristicB[r, q] == MathieuCharacteristicA[r, q] /. #, 20] =!= True &]
	,
	{}
	,
	TestID->"MathieuCharacteristicB-20050111-6o5f3c"
	]

Test sources

Benchmarks

Test[
	"Kamke 8.3";
	eqns = {Derivative[1][x][t] == a*y[t], Derivative[1][y][t] == b*x[t]};
	sol = DSolve[eqns, {x, y}, t];
	{FreeQ[sol, DSolve|Solve|Integrate], Union[Cases[sol,C[i_Integer],Infinity]], eqns /. sol //Flatten //Simplify}
	,
	{True, {C[1],C[2]}, {True, True}}
	,
	TestID->"Kamke_Chapter8-1-20040519-t4an3f"
	]

Test[
	"Caprasse system, lex order";
	polys = {-2*x + 2*t*x*y - z + y^2*z, 2 + 4*x^2 - 10*t*y + 4*t*x^2*y - 10*y^2 + 2*t*y^3 + 4*x*z - 
		  x^3*z + 4*x*y^2*z, -x + t^2*x - 2*z + 2*t*y*z, 2 - 10*t^2 - 10*t*y + 2*t^3*y + 4*x*z + 
		  4*t^2*x*z + 4*z^2 + 4*t*y*z^2 - x*z^3};
	vars = {x,y,z,t};
	gb = GroebnerBasis[polys, vars];
	expect = {-27 + 1719*t^4 + 8000*t^6 + 11746*t^8 - 11746*t^12 - 8000*t^14 - 1719*t^16 + 27*t^20, 
		  -27*z - 243*t^2*z - 711*t^4*z - 559*t^6*z + 559*t^8*z + 711*t^10*z + 243*t^12*z + 
		  27*t^14*z, -4413 + 29540*t^2 + 143174*t^4 - 12190*t^6 - 186100*t^8 - 34270*t^10 + 
		  46826*t^12 + 17190*t^14 + 513*t^16 - 270*t^18 - 9600*z^2 - 22400*t^2*z^2 + 
		  22400*t^4*z^2 + 9600*t^6*z^2, -34867 + 224387*t^2 + 1063936*t^4 + 1330344*t^6 - 
		  138430*t^8 - 1370330*t^10 - 890936*t^12 - 187344*t^14 + 297*t^16 + 2943*t^18 - 
		  64000*z^2 + 64000*t^4*z^2 - 16000*z^4 + 16000*t^2*z^4, 
		  4233*z + 12098*t^2*z + 14735*t^4*z - 11140*t^6*z - 18537*t^8*z - 6750*t^10*z - 
		  783*t^12*z + 2688*z^3 + 1792*t^2*z^3 + 1664*t^4*z^3 - 384*z^7 + 24*z^9, 
		  1570281*t - 48825575*t^3 - 241936588*t^5 - 311914460*t^7 + 48665150*t^9 + 321282190*t^11 + 
		  198035588*t^13 + 40085460*t^15 - 190431*t^17 - 627615*t^19 + 6144000*y - 5376000*t*z^2 - 
		  3584000*t^3*z^2 - 3328000*t^5*z^2 + 768000*t*z^6 - 48000*t*z^8, 
		  12288*x + 22641*z + 29642*t^2*z + 24591*t^4*z - 28980*t^6*z - 44865*t^8*z - 17334*t^10*z - 
		  2079*t^12*z + 5248*z^3 - 1280*t^2*z^3 - 1920*t^4*z^3 + 1792*z^5 - 128*z^7};
	Union[Abs[Together[gb/expect]]]
	,
	{1}
	,
	TestID->"GroebnerBasis-20050729-b1sn37"
	]

Test sources

Journal articles

Test[
	"Integrate bug mentioned in:
	 Don Page, Average Entropy of a Subsystem, Phys.Rev.Lett. 71 (1993), pp1291-1294";
	res = Integrate[Exp[-x]*x^N*Log[x], {x,0,Infinity}];
	expect = Gamma[1 + N]*PolyGamma[1 + N];
	expcond = Re[N]>-1;
	cond = First[Cases[{res}, If[x_,y__]:>x, Infinity]];
	{!FreeQ[res, If], cond === expcond,
	 (res /. If[x_,y_,z_]:>y) - expect //Simplify,
	 Check[Last[res], "Diverges", Integrate::idiv]}
	,
	{True, True, 0, "Diverges"}
	,
	Integrate::idiv
	,
	TestID->"Articles-20050516-8lrn1l"
	]

Test sources

Wolfram Research staff

Test[
	gb = GroebnerBasis[{x + y + z, x - 2 y + z^3, x^2 - 2y^3 + z}, {x, y, z}];
	expect = {27*z + 12*z^2 + 2*z^3 + 12*z^4 - 6*z^5 + 3*z^6 + 6*z^7 - 2*z^9, 3*y + z - z^3, 3*x + 2*z + z^3};
	Union[Abs[Cancel[gb/expect]]]
	,
	{1}
	,
	TestID->"GroebnerBasis_TrottSymbolicsGuideBook-20050710-3g968r"
	]

Test sources

Developer tests

Test[
	Reduce[x^7 - 3*x*y^6 + 5*y^7 == 3061061 && (x | y) ∈ Integers && 0 < y < 10^100, {x, y}]
	,
	x == 3 && y == 7
	,
	TestID->"ACA2002ConferenceTalk-20030207-992y3d"
	]

Test[
	A = {{0, 0, 0, 0}, {0, 8/9, -4/9, -4/9}, {0, -4/9, 8/9, -4/9}, {0, -4/9, -4/9, 8/9}};
	B = PseudoInverse[A];
	{A.B.A == A, B.A.B == B, A.B == Transpose[A.B], B.A == Transpose[B.A]}
	,
	{True, True, True, True}
	,
	TestID->"PseudoInverse-20030521-2wjn1d-bug-36064"
	]

Test sources

QA tests (from test engineer)

Test[
	MeijerG[{{}, {}}, {{1/2}, {}}, z, 7]
	,
	z^(1/14)/E^z^(1/7)
	,
	TestID->"MeijerGBasic-20050406-0q1n4r"
	]

Test[
	Series[Integrate[Exp[Sin[x]], {x, 0, Log[z]}], {z, 1, 4}]
	,
	SeriesData[z, 1, {1, 0, 0, -1/24}, 1, 5, 1]
	,
	TestID->"SeriesMisc-19980325-r70l3q-bug-58780"
	]

Test sources

Well-known examples

Test[
	Head[Integrate[Sin[x], {x, -Infinity, Infinity}]]
	,
	Integrate
	,
	Integrate::idiv
	,
	TestID->"Definite_Elementary_Divergent-20031022-8n0c1u-bug-56039"
	]

Code coverage

Test[
	GroebnerBasis[{x^2+y^2, x-y}, {x,y}, CoefficientDomain->InexactNumbers[1]] //Head
	,
	GroebnerBasis
	,
	GroebnerBasis::fltgb
	,
	TestID->"GroebnerBasis-20050716-e66u33"
	]

Result verification

Results are verified on a case-by-case basis, using a variety of techniques. We want tests to be as robust as possible, but also easy to understand. When it is not immediately obvious what the expected result should be, sometimes verification is performed before writing the test (e.g. for checking a Reduce result for integer roots of a polynomial), in order to prevent unnecessary computation.

Equal/SameQ (if result is simple enough)

Test[
	Limit[x*Gamma[x], x->Infinity]
	,
	Infinity
	,
	TestID->"Limit_SpecialFunctions-20040831-uuyl3r"
	]

Variant of Simplify (Together/Simplify/FullSimplify) together with a LeafCount bound

Test[
	res = Limit[ArcCot[x], x -> I/3, Direction -> 1];
	expect = ArcCot[I/3];
	{FreeQ[res, Limit], res-expect //FullSimplify, If[LeafCount[res] <= 25, True, LeafCount[res]]}
	,
	{True, 0, True}
	,
	TestID->"Limit_ArcTrigs-19980511-9w1t8r-bug-36007"
	]

Result verification

Random substitution for free variables (Random, FindInstance)

Test[
	"Example from page 331";
	eqn = u[k+1] - k u[k] + Exp[-r] r^k;
	sol = RSolve[eqn == 0, u, k];
	eqn /. sol /. k->Random[Integer,{1,10}] //FullSimplify
	,
	{0}
	,
	TestID->"Milne-Thomson-20010510-pytu4c"
	]

Test[
	integrand = 1/(x Sqrt[(1 - x^2)(1 + x^2/c^2)]);
	res = Integrate[integrand, {x, 1, Infinity}];
	cond = First[Cases[{res}, If[x_,y__]:>x, Infinity]];
	expcond = Inequality[-1, LessEqual, Im[c], LessEqual, 1] || Re[c] != 0;
	inst = FindInstance[cond, c, 10];
	check = Select[inst, Internal`DeactivateMessages[Chop[(res /. #) - 
	NIntegrate[integrand /. #, {x, 1, Infinity}, WorkingPrecision -> 26, MaxRecursion -> 20]]]=!=0&];
	{!FreeQ[res, If], cond === expcond, check}
	,
	{True, True, {}}
	,
	TestID->"Bugs-49000-20040722-yiu819-bug-49983"
	]

Test sources

Substitution with representative set of values

Test[
	exppts = {0, -1, -3, -Sqrt[5], 10/11, Pi/2, Pi, 7, (I/2)*Pi, (-I)*Pi, I, -I, -(-1)^(4/5), 
			 (-11/10)^(2/5), (243*(-1)^(3/5)*3^(1/5))/(100*10^(3/5)), Root[12 - #1 + 2*#1^5 &, 1, 0]};
	Map[
	   Union[Chop[N[
	   		Table[Normal[Series[AiryAi[x], {x, #, 4}]] - AiryAi[x] /. x -> # + (1/100)*Exp[Pi*I*k/6], {k, 0, 11}]
	   		]]] &,
	   exppts
	]
	,
	{{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}
	,
	TestID->"AiryAi-20040924-a2vjd6"
	]

Test sources

Numericizing (to high precision if expression has complicated special functions)

Test[
	"Gives a result involving PFQ. No exceptional points in integration range.";
	res = Integrate[Exp[-Pi*t]*Log[t]^2, {t, 1, 20}];
	{FreeQ[res,Integrate], N[res,20] - NIntegrate[Exp[-Pi*t]*Log[t]^2, {t, 1, 20}] //Chop}
	,
	{True, 0}
	,
	TestID->"Newton-Leibniz-20030919-2oz51w-bug-48497-bug-55825"
	]

Result verification

Using invariants, identities, or properties of the function

Test[
	"Want to give result in terms of elementary functions";
	integrand = Sqrt[x + Sqrt[x^2 + 1]]/x;
	res = Integrate[integrand, x];
	disallowed = {MeijerG, AppellF1, HypergeometricPFQ};
	{FreeQ[res,Integrate], FreeQ[res, Alternatives@@disallowed], D[res, x] - integrand //Simplify, LeafCount[res] < 200}
	,
	{True, True, 0, True}
	,
	TestID->"Indefinite_Algebraic-20050930-tzrn1i"
	]

Test[
	"Gives an implicit solution";
	eqn = {x - y[x] + (3*x + y[x])*y'[x] == 0, y[2] == -1};
	sol = DSolve[eqn, y[x], x];
	solveSol = Solve[D[First[sol], x], y'[x]];
	{FreeQ[sol, DSolve], Simplify[First[eqn] /. solveSol], First[sol] /. x -> 2 /. y[2] -> -1}
	,
	{True, {True}, True}
	,
	{Solve::dinv}
	,
	TestID->"ODE19-19980325-ov2708"
	]

Test[
	HermiteH[100,x] == Expand[(-1)^100*Exp[x^2]*D[Exp[-x^2],{x,100}]]
	,
	True
	,
	TestID->"HermiteH-19980325-48mv84"
	]

Result verification

For conditional results from Integrate, using Refine + pattern matching to extract conditions and corresponding results

Test[
	res = Integrate[Exp[I*p*x]/(1 + p^2), {p, -Infinity, Infinity}];
	cond = First[Cases[{res}, If[x_,y__]:>x, Infinity]];
	{!FreeQ[res, If], cond === Element[x,Reals], Refine[res,cond] - Pi/E^Abs[x] //Simplify,
	 Check[res /. If[x__,y_]:>y, "Diverges", Integrate::idiv]}
	,
	{True, True, 0, "Diverges"}
	,
	{Integrate::idiv}
	,
	TestID->"Bugs-59000-20050926-0dmj1x-bug-59378"
	]

Using a combination of Reduce and a Simplify variant to simplify Integrate conditions
(but also want to keep conditions from becoming too complicated)

Test[
	res = Integrate[1/Sqrt[(1+x)(2+x)(3+x)(4+x)], {x,0,a}];
	cond = First[Cases[{res}, If[x_,y__]:>x, Infinity]];
	simpcond = Reduce[Simplify[cond, Im[a]==0], a];
	expcond = a >= (-5 + Sqrt[5])/2;
	res2 = res /. If[x_,y_,z_]:>y;
	expres = EllipticF[ArcSin[2*Sqrt[2]/3], 3/4] -
			EllipticF[ArcSin[Sqrt[2/3]*Sqrt[(4 + a)/(3 + a)]], 3/4];
	{!FreeQ[res,If], simpcond == expcond, FullSimplify[res2 - expres, simpcond], LeafCount[cond] <= 250}
	,
	{True, True, 0}
	,
	TestID->"Definite_Elliptic-20050930-98dt25-bug-32409-bug-61918"
	]

Test failures

We have nightly testruns on the three major platforms (Windows, OSX, and Linux).

True failure

Report bug -> Add bug number to test -> Investigate -> Fix -> Verify fix -> Add tests

Example

This is from a bug report based on a MathWorld entry (SquareTrianglePicking.html).

This integral should give 11/144, but it used to give a wrong result in our development version.

pw=Boole[
            0 < x ≤ (1/2) ∧ 0 < 1 - y ≤ x] (Piecewise[{{(1/2) (((x - 1))^2) ((((12 (y^4) - 24 (y^3) + 25 (y^2) - 13 y + 1) (x^2) + (-30 (y^4) + 60 (y^3) - 43 (y^2) + 13 y - 1) x + y (3 (y^3) - 6 (y^2) + 4 y - 1))/(x (y - 1) y)) - 
                        3 (x (6 (y^2) - 6 y + 
                                    2) - (((1 - 2 y))^2)) log(
                            (1/x) - 1) + 
                        3 (2 y - 1) (-2 (y^2) + 2 y + 
                              2 x ((y^2) - y + 1) - 1) log(
                            (1/y) - 1)), x ≥ (1/2) ∧ 
                    x ≥ y}, {-(((((x - 1))^2) (-12 (x^2) (y^4) + 
                              30 x (y^4) - 3 (y^4) + 24 (x^2) (y^3) - 
                              60 x (y^3) + 6 (y^3) - 25 (x^2) (y^2) + 
                              43 x (y^2) - 4 (y^2) + 13 (x^2) y - 
                              13 x y + 
                              3 x (y - 
                                    1) (x (6 (y^2) - 6 y + 
                                        2) - (((1 - 2 y))^2)) log(
                                  (1/x) - 1) y + 
                              3 x (2 (y^2) - 3 y + 
                                    1) (-2 (y^2) + 2 y + 
                                    2 x ((y^2) - y + 1) - 1) log(
                                  (y/(1 - y))) y + y - (x^2) + 
                              x))/(2 x (y - 
                              1) y)), x ≥ (1/2) ∧ 
                    y ≤ (1/2)}, {-(((y^2) (12 (y^2) (x^4) + 6 y (x^4) - 15 (x^4) - 
                              24 (y^2) (x^3) - 12 y (x^3) + 30 (x^3) + 
                              25 (y^2) (x^2) - 7 y (x^2) - 14 (x^2) - 
                              13 (y^2) x + 13 y x + 
                              3 (2 (x^2) - 3 x + 1) (y - 
                                    1) (2 ((x^2) - x + 1) y - 
                                    1) log((1/x) - 1) x - 
                              3 (x - 1) (y - 
                                    1) ((6 y - 2) (x^2) + (2 - 
                                        6 y) x + 2 y - 1) log(
                                  (y/(1 - y))) x - x + (y^2) - 
                              y))/(2 (x - 1) x (y - 
                              1))), x ≥ (1/2) ∧ 
                    x + y ≤ 1}, {-(((x^2) (12 (x^2) (y^4) + 6 x (y^4) - 15 (y^4) - 
                              24 (x^2) (y^3) - 12 x (y^3) + 30 (y^3) + 
                              25 (x^2) (y^2) - 7 x (y^2) - 14 (y^2) - 
                              13 (x^2) y + 13 x y - 
                              3 (x - 1) (y - 1) (-2 (y^2) + 
                                    2 y + x (6 (y^2) - 6 y + 2) - 
                                    1) log((x/(1 - x))) y + 
                              3 (2 (y^2) - 3 y + 
                                    1) (2 ((y^2) - y + 
                                        1) (x^2) + (-2 (y^2) + 2 y - 
                                        3) x + 1) log(
                                  (1/y) - 1) y - y + (x^2) - 
                              x))/(2 (x - 1) (y - 
                              1) y)), x + y ≤ 1}, {-(((((y - 1))^2) (-12 (y^2) (x^4) + 
                              30 y (x^4) - 3 (x^4) + 24 (y^2) (x^3) - 
                              60 y (x^3) + 6 (x^3) - 25 (y^2) (x^2) + 
                              43 y (x^2) - 4 (x^2) + 13 (y^2) x - 
                              13 y x + 
                              3 (2 (x^2) - 3 x + 
                                    1) y (2 (y - 1) (x^2) - 
                                    2 (y - 1) x + 2 y - 1) log(
                                  (x/(1 - x))) x + 
                              3 (x - 
                                    1) y ((6 y - 4) (x^2) + (4 - 
                                        6 y) x + 2 y - 1) log(
                                  (1/y) - 1) x + x - (y^2) + 
                              y))/(2 (x - 1) x y)), True}}])

8*
        Integrate[pw,  {x,  0,  1},  {y,  0,  1},  
          GenerateConditions -> True]

Narrowed down to this integral, which is handled by specialized Log integration code:

Integrate[y*Log[1/y  -  1],  {y,  2/3,  1}]

Providing the date of breakage helps narrow down the problem to a commit (or set of commits).

False failure

Fix test using verification techniques outlined above. As is the case when writing tests, the goal is to make the test as robust as possible while keeping it easy to understand.

Conclusion

Test sources

Result verification

Test failures

Questions?