(*	benchmarks.m (obtained from Wolfram Research 10/1/90)
 		     (modified by Lyle Wright 10/2/90 and 10/3/90)

	This file contains a set of Mathematica expressions for use as
	timing benchmarks.

	To use it, just type <<benchmarks.m at a Mathematica input prompt.

	Output is written to a file called math.timings in the current
	working directory.

	The output file contains only the timings, in seconds.

	This output file has been updated by Lyle Wright (10/2/90) to display
     	the operation followed by the time that it takes to perform the
	operation.

	This version is set up for Unix workstations only.  To use on 
	another operating system, change $Display to write to a temporary
	file.

	The benchmarks in the middle section are derived from a review of
	Mathematica by Eugene A. Herman, which appeared in the November
	1988 issue of Notices of the American Mathematical Society.
	(Notices of the AMS, vol. 35, p. 1333)

	The final set of benchmarks are the standard benchmarks distributed
	with MatLab.

	Benchmarks in the first section are used locally by Wolfram Research,
	Inc. for timing tests.

	This benchmark suite was developed by Kate Ebneter, Software Quality
	Assurance Manager, Wolfram Research, Inc.

*)

(* Set things up. *)

OpenWrite["math.timings", FormatType->OutputForm]
printIt[expr_] := Write["math.timings", expr]
$Display = "/dev/null"

printIt["This benchmark were performed on ptnext.claremont.edu"]
printIt["which is a 68040 NeXTstation running NeXTStep 2.0."]
printIt[""]

(*  WRI benchmarks  *)

printIt["The following set of benchmarks are those used by"]
printIt["Wolfram Research, Inc. for timing tests."]
printIt[""]

printIt["Factor[Expand[(1+x)^100]]"]
printIt[Timing[Factor[Expand[(1+x)^100]];]]
printIt[""]

printIt["Factor[Expand[(1+x)^20]]"]
printIt[Timing[Factor[Expand[(1+x)^20]];]]
printIt[""]

printIt["1989^1989"]
printIt[Timing[1989^1989;]]
printIt[""]

printIt["Eigenvalues[Table[Random[], {40}, {40} ]]"]
printIt[Timing[Eigenvalues[Table[Random[], {40}, {40} ]];]]
printIt[""]

printIt["<<Graphics/Polyhedra.m"]
<<Graphics/Polyhedra.m
printIt["Show[Graphics3D[Stellate[Stellate[Icosahedron[] ]]]]"]
printIt[Timing[Show[Graphics3D[Stellate[Stellate[Icosahedron[] ]]]]]]
printIt[""]

printIt["N[Pi, 1000]"]
printIt[Timing[N[Pi, 1000];]]
printIt[""]

printIt["N[E, 1000]"]
printIt[Timing[N[E, 1000];]]
printIt[""]

printIt["N[EulerGamma,400]"]
printIt["This takes too long so I am not going to evaluate it."]
(* printIt[Timing[N[EulerGamma,400];]] *)
printIt[""]

printIt["BernoulliB[60]"]
printIt[Timing[BernoulliB[60];]]
printIt[""]

printIt["1000!"]
printIt[Timing[1000!;]]
printIt[""]

printIt["FactorInteger[2^2^6+1]"]
printIt[Timing[FactorInteger[2^2^6+1];]]
printIt[""]

(*  Benchmarks from Notices of the AMS review  *)

printIt["The following set of benchmarks were obtained from"]
printIt["Notices of the AMS review."]
printIt[""]

printIt["Solve[x^3 - x^2 - x - 1 == 0, x]"]
printIt[Timing[Solve[x^3 - x^2 - x - 1 == 0, x];]]
printIt[""]

printIt["NRoots[x^3 - x^2 - x - 1 == 0, x, 10]"]
printIt[Timing[NRoots[x^3 - x^2 - x - 1 == 0, x, 10];]]
printIt[""]

printIt["NRoots[x^3 - x^2 - x - 1 == 0, x, 20]"]
printIt[Timing[NRoots[x^3 - x^2 - x - 1 == 0, x, 20];]]
printIt[""]

printIt["The function definition was a little complicated"]
printIt["and I did not feel like displaying it."]
printIt[""]

H[n_] := Table[1/(i + j - 1), {i, n}, {j, n}]
ERR[n_] :=
        Block[{error, h, hinv, eye},
                h = N[H[n], 15];
                hinv = Inverse[h];
                eye = IdentityMatrix[n];
                error = Max[Abs[eye - h.hinv], Abs[eye - hinv.h]]
        ]

printIt[""]
printIt["ERR[11]"]
printIt[Timing[ERR[11];]]
printIt[""]

printIt["ERR[12]"]
printIt[Timing[ERR[12];]]
printIt[""]

printIt["ERR[13]"]
printIt[Timing[ERR[13];]]
printIt[""]

printIt["ERR[14]"]
printIt[Timing[ERR[14];]]
printIt[""]

printIt["Clear[H];"]
Clear[H];
printIt[""]

printIt["Clear[ERR];"]
Clear[ERR];
printIt[""]

printIt["f[x_, y_, z_] := Sin[x] + y^2 + Log[z]"]
f[x_, y_, z_] := Sin[x] + y^2 + Log[z]
printIt[""]

printIt["g[x_, y_, z_] := 3 x + 2^y - z^3"]
g[x_, y_, z_] := 3 x + 2^y - z^3
printIt[""]

printIt["h[x_, y_, z_] := x + y + z"]
h[x_, y_, z_] := x + y + z
printIt[""]

printIt["FindRoot[{f[x, y, z] == 7, g[x, y, z] == -1, h[x, y, z] == 5},"]
printIt["{x, 0}, {y, 1}, {z, 1}, AccuracyGoal->13, WorkingPrecision->26]"]
printIt[Timing[FindRoot[{f[x, y, z] == 7, g[x, y, z] == -1, h[x, y, z] == 5},
	{x, 0}, {y, 1}, {z, 1}, AccuracyGoal->13, WorkingPrecision->26];]]
printIt[""]


printIt["FindRoot[{f[x, y, z] == 7, g[x, y, z] == -1, h[x, y, z] == 5},"]
printIt["{x, 0}, {y, 0}, {z, 1}, AccuracyGoal->13, WorkingPrecision->26]"]
printIt[Timing[FindRoot[{f[x, y, z] == 7, g[x, y, z] == -1, h[x, y, z] == 5},
	{x, 0}, {y, 0}, {z, 1}, AccuracyGoal->13, WorkingPrecision->26];]]
printIt[""]

printIt["Clear[h];"]
Clear[h];
printIt[""]

printIt["Clear[g];"]
Clear[g];
printIt[""]

printIt["Clear[f];"]
Clear[f];
printIt[""]


printIt["f20[x_] := N[Sum[Sin[(2 k - 1) x]/(2 k - 1), {k, 20}]]"]
f20[x_] := N[Sum[Sin[(2 k - 1) x]/(2 k - 1), {k, 20}]]
printIt[""]

printIt["Plot[f20[x], {x, 0, Pi/2}]"]
printIt[Timing[Plot[f20[x], {x, 0, Pi/2}]]]
printIt[""]


printIt["NIntegrate[f20[x], {x, 1, 2}, AccuracyGoal->10, 
WorkingPrecision->20];"]
 printIt[Timing[NIntegrate[f20[x], {x, 1, 2}, AccuracyGoal->10,
WorkingPrecision->20];]] 
printIt[""]

printIt["N[Integrate[f20[x], {x, 1, 2}], 10];"]
printIt[Timing[N[Integrate[f20[x], {x, 1, 2}], 10];]]
printIt[""]

printIt["-FindMinimum[-f20[x], {x, 0, 0, Pi/2}];"]
printIt[Timing[-FindMinimum[-f20[x], {x, 0, 0, Pi/2}];]]
printIt[""]

printIt["Clear[f20];"]
Clear[f20];
printIt[""]

printIt["f[x_, n_Integer] := N[Sum[Sin[(2 k - 1) x]/(2 k - 1), {k, n}]]"]
f[x_, n_Integer] := N[Sum[Sin[(2 k - 1) x]/(2 k - 1), {k, n}]]
printIt[""]

printIt["Plot[{f[x, 8], f[x, 9], f[x, 10]}, {x, 0, Pi}]"]
printIt[Timing[Plot[{f[x, 8], f[x, 9], f[x, 10]}, {x, 0, Pi}]]]
printIt[""]

printIt["f[1, 500];"]
printIt[Timing[f[1, 500];]]
printIt[""]

printIt["f[1, 10000];"]
printIt[Timing[f[1, 10000];]]
printIt[""]

printIt["FactorInteger[266382004787];"]
printIt[Timing[FactorInteger[266382004787];]]
printIt[""]

printIt["Clear[f];"]
Clear[f];
printIt[""]

printIt["<<Examples/RungeKutta.m"]
<<Examples/RungeKutta.m
printIt["tmp = RungeKutta[1 - 2 x y, {x, y}, {-1, 0}, 0.1, 30];"]
printIt[Timing[tmp = RungeKutta[1 - 2 x y, {x, y}, {-1, 0}, 0.1, 30];]]
printIt[""]

printIt["ListPlot[tmp, PlotJoined->True]"]
printIt[Timing[ListPlot[tmp, PlotJoined->True];]]
printIt[""]

printIt["f[t_] := Integrate[Exp[x^2], {x, -1, t}]"]
f[t_] := Integrate[Exp[x^2], {x, -1, t}]
printIt[""]


printIt["Plot[f[y], {y, -1, 2}]"]
printIt[Timing[Plot[f[y], {y, -1, 2}]]]
printIt[""]

printIt["Clear[f];"]
Clear[f];
printIt[""]

printIt["tmp = Integrate[(1 + x)/(1 + x^2 + x^4 + x^6), x];"]
printIt[Timing[tmp = Integrate[(1 + x)/(1 + x^2 + x^4 + x^6), x];]]
printIt[""]

printIt["tmp = D[tmp, x];"]
printIt[Timing[tmp = D[tmp, x];]]
printIt[""]

printIt["Simplify[ tmp ];"]
printIt[Timing[ Simplify[ tmp ]; ]]
printIt[""]


printIt["y = Series[f[x], {x, 0, 6}];
y = y /. {f[0] -> 0, f'[0] -> 1, Derivative[n_][f][0] -> a[n]};
eqn = D[y, {x, 2}] + Series[Sin[x], {x, 0, 6}] y == 0;
var = Table[a[k], {k, 2, 6}];
Solve[eqn, var];"]

y = Series[f[x], {x, 0, 6}];
y = y /. {f[0] -> 0, f'[0] -> 1, Derivative[n_][f][0] -> a[n]};
eqn = D[y, {x, 2}] + Series[Sin[x], {x, 0, 6}] y == 0;
var = Table[a[k], {k, 2, 6}];
printIt[Timing[Solve[eqn, var];]]
printIt[""]

printIt["y =.
eqn =.
var =."]
y =.
eqn =.
var =.
printIt[""]

(*  Tests from MatLab benchmark suite  *)

printIt["The following set of benchmarks are those from the"]
printIt["MatLab benchmark suite."]
printIt[""]

printIt["a = Table[Random[], {50}, {50}];"]
a = Table[Random[], {50}, {50}];
printIt[""]

printIt["printIt[a.a;]"]
printIt[Timing[a.a;]]
printIt[""]

printIt["a = Table[Random[], {50}, {50}];"]
a = Table[Random[], {50}, {50}];
printIt[""]

printIt["Inverse[a];"]
printIt[Timing[Inverse[a];]]
printIt[""]

printIt["a = Table[Random[], {25}, {25}];"]
a = Table[Random[], {25}, {25}];
printIt[""]

printIt["Eigenvalues[a];"]
printIt[Timing[Eigenvalues[a];]]
printIt[""]

printIt["a = Table[Random[], {i, 4096}] + I;"]
a = Table[Random[], {i, 4096}] + I;
printIt[""]

printIt["Fourier[a];"]
printIt[Timing[Fourier[a];]]
printIt[""]
         
printIt["a = Table[Random[], {100}, {100}];"]
a = Table[Random[], {100}, {100}];
printIt[""]

printIt["b = Table[Random[], {100}];"]
b = Table[Random[], {100}];
printIt[""]

printIt["x = LinearSolve[a, b];"]
printIt[Timing[x = LinearSolve[a, b];]]
printIt[""]

printIt["a = Table[0, {i, 1000}];"]
a = Table[0, {i, 1000}];
printIt[""]

printIt["Do[a[[i]] = 1, {i, 1000}];"]
printIt[Timing[Do[a[[i]] = 1, {i, 1000}];]]
printIt[""]

printIt["a = IdentityMatrix[25];"]
a = IdentityMatrix[25];
printIt[""]

printIt["ListPlot3D[a];"]
printIt[Timing[ListPlot3D[a];]]
printIt[""]

printIt["This concludes the benchmark tests that are"]
printIt["included in this package."]
