Lecture 10 : Exact vs. Inexact Arithmetic, Packages and Namespaces

Mathematica tracks the precision (number of significant digits) of the numerical quantities being manipulated. Typically, with each operation, the precision of the result is less than that of the input, so when the precision of the output is important, you should use as high precision input as possible. One exception is that of machine precision-- any calculations done with machine precision return results in machine precision, regardless of whether or not Mathematica could prove that the result indeed is of machine precision. That is, for computational efficiency, Mathematica doesn't track the precision of calculations 'contaminated' with machine precision quantities. 

Enter a number with a given precision either by enumerating all the significant digits, or using the numbermark notation. The precision of numbers can be checked with the Precision command, and set with SetPrecision (zeros are inserted as needed). One can also use N to evaluate an expression to a desired precision,  and Rationalize, Round, Ceiling, Floor and other such functions to convert expressions to infinite precision. $MachinePrecision (and MachinePrecision) store the number of digits of precision on your computer. MachineNumberQ tests an expression to see if it has MachinePrecision. see tutorial/MachinePrecisionNumbers for more information on machine precision calculations.

Numerical operators such as NSolve, NIntegrate, and NDSolve typically accept the directives WorkingPrecision and PrecisionGoal which specify, respectively, the maximum precision with which to carry out internal calculations and the precision which you desire in your final results. To be effective, the inputs (if numerical) have to have precision higher than WorkingPrecision and PrecisionGoal. If left unspecified, assume that Mathematica will work in machine precision (for speed). Note that the higher precision the input is, the longer the calculations take. 

Mathematica uses $MaxExtraPrecision extra digits of precision in its internal numerical calculations. For example, when you call N to evaluate a quantity, $MaxExtraPrecision affects the intermediate calculations. There are also some nonobvious instances where $MaxExtraPrecison comes into play; for example, when comparing equality of symbolic quantities, Mathematica will sometimes evaluate them numerically before comparison. This can lead to problems sometimes, e.g. if you compare two complex expressions that are in fact equal, then Mathematica will compare them numerically, and no amount of precision will allow it to determine which is the largest. If you need to, you can use a Block to locally change the value of $MaxExtraPrecision.

For more discussion of precision in Mathematica, see tutorial/ArbitraryPrecisionNumbers

Note there's an associated notion of accuracy: while precision is the number of significant digits in a quantity, accuracy is the number of significant digits after the decimal point. Pretty much the same comments above apply to accuracy: corresponding to Precision there's Accuracy, to SetPrecision there's SetAccuracy, there's AccuracyGoal, etc. Accuracy is important, but you can usually control it by controlling precision, so we'll concentrate on that.

Basic Examples

1/3 //Precision
$MachinePrecision
N[1/3, MachinePrecision]
3.33 // Precision (* defaults to MachinePrecision, because it fits *)
MachineNumberQ[3.33333333333333333333333333333333333333]  (* doesn't fit MachinePrecision *)

3.33333333333333333333333333333333333333 // Precision 
3.33333333333333333333333333333333333333 // First@RealDigits[#]& // Length
3.33333333333333333333333333333333333333 // Accuracy
3.33333333333333333333333333333333333333`20 //Precision
3.33333333333333333333333333333333333333`50 //Precision
3.3`∞ //Precision
3.333`1
3.333`1 //FullForm
3.333`1 == 3333/1000 (* even though Mathematica only displays the sig digits, it uses the rest in the comparison *)

NSum[(1/(i^2)),{i,1,2000}] (* gives machine precision by default *)
NSum[(1/(i^2)),{i,1,2000}, PrecisionGoal->100] (* fails-- look at the poor error estimate-- maybe because we aren't using enough precision on intermediate steps *)
NSum[(1/(i^2)),{i,1,2000}, PrecisionGoal->100, WorkingPrecision->1000]
% // FullForm (* shows the precision with the numbermark formatting *)

Inheritance of precision

(* example where N can't give the precision asked for, because you of inherited machine precision; example taken from http://dirac.org/apps/mathematica/numerics.php *)
x={2.0,3.0,5.0};
A={{6.0,2.0,1.0},{2.0,3.0,1.0},{1.0,1.0,1.0}};
For[k=0,k<15,++k,lambda=x.A.x/(x.x);
y=LinearSolve[A,x];
x=y/Norm[y,Infinity];]
N[lambda,30]

x={2,3,5};
A={{6,2,1},{2,3,1},{1,1,1}};
For[k=0,k<15,++k,lambda=x.A.x/(x.x);
y=LinearSolve[A,x];
x=y/Norm[y,Infinity];]
N[lambda,30]

(* or equivalently (could also use Round, Ceiling, Floor to do the same thing)) *)
x=Rationalize@{2.0,3.0,5.0};
A=Rationalize@{{6.0,2.0,1.0},{2.0,3.0,1.0},{1.0,1.0,1.0}};
For[k=0,k<15,++k,lambda=x.A.x/(x.x);
y=LinearSolve[A,x];
x=y/Norm[y,Infinity];]
N[lambda,30]

(* you could start off with higher, but finite, precision numbers *)
x=SetPrecision[{2.0,3.0,5.0}, 5000];
A=SetPrecision[{{6.0,2.0,1.0},{2.0,3.0,1.0},{1.0,1.0,1.0}}, 5000];
For[k=0,k<15,++k,lambda=x.A.x/(x.x);
y=LinearSolve[A,x];
x=y/Norm[y,Infinity];]
N[lambda,30]

Cost of Extra Precision

(* on average, exact operations are the most expensive, then arbitrary precision operations, then machine precision operations *)
exacttimes = {};
inexacttimes = {};
precisiontimes = {};
For[ n=1, n< 100, n++,
Remove["A"];
A = Table[(1/(j+k-1)),{j,1,10}, {k,1,10}];
(* exact numbers *)
exacttimes = Append[exacttimes, First@AbsoluteTiming[ A //Eigenvalues]];
(* inexact numbers: machine precision by default *)
inexacttimes = Append[inexacttimes, First@AbsoluteTiming[A //N//Eigenvalues]];
(* inexact numbers: specify a precision *)
precisiontimes = Append[precisiontimes, First@AbsoluteTiming[A // N[#, 40] & // Eigenvalues]];
];
ListPlot[{exacttimes, inexacttimes, precisiontimes}, Filling->Axis, PlotStyle->{Red, Green, Yellow}]

(* inexact numbers: compare to Matlab *)
b = RandomReal[{0,1}, 5000];
A = RandomReal[{0,1}, {5000,5000}];
First@Timing[x = LinearSolve[A,b]]
Norm[A . x - b,∞]

(* some computations should only be done on numerical data! *)
A = HilbertMatrix[50];
(* Timing@Short@Eigenvalues@A (* a list of 50 Root objects, takes 32 seconds to compute -- Mathematica calculated the characteristic poly of A!*) *)

A = N[HilbertMatrix[50]]; (* usually you only want approximations to the eigenvalues *)
A // Eigenvalues (* almost immediate *)

(* the advantage over Matlab is I can ask for more than MachinePrecision in the eigenvalues *)
A = N[HilbertMatrix[50], 100];
A // Eigenvalues //Short
% // Precision

Packages and Contexts

Contexts are used in Mathematica to help organize collections of symbols. They are analogous to namespaces in C++, etc. Packages are coherent collections of code, which usually define useful symbols and functions. For instance, Mathematica ships with a Histograms package which has useful functions for graphing histograms, and a VectorAnalysis package which defines the Grad, Curl, Div, Laplacian, etc. operations in several different coordinates. There are tons more packages; Combinatorica is a pretty neat one, which provides many convenience codes for dealing with discrete mathematics and graph theory.

Each symbol has a full name, which consists of a context and a short name. Contexts are hierarchial, so a symbol could for instance have the full name context`subcontext`shortname (the backquote separates contexts and short names). Generally, you refer to symbols by their short name and Mathematica uses the current $Context and then $ContextPath to determine the full name of the variable you're referring to. If you consider contexts to be like the folders in a file system, then $Context would be the current directory and $ContextPath would be like your path. If Mathematica can't find the symbol you referenced, it creates it in the current context. By default the current context is Global`.

You import a package before calling the functions it provides. Do so with the Needs command, Needs["context`"]. This will load the package if it has not yet been loaded. Alternatively, you can use the Get command, abbreviated <<, e.g. <<Combinatorica` or <<VectorAnalysis`. Note that in both cases, you tell Mathematica that you're loading a context by using the context mark. What Mathematica then does is determine the filename corresponding to that context, and executes the commands in the file. (Needs calls Get)

Be careful that you don't call the commands in a package before loading it, because that will generate symbols in the current context that will shadow the package definitions after you load it; this can sometimes cause a problem. If you do this by accident, use Remove to remove the shadowing symbol.

$Packages gives a list of the currently loaded packages (really, their contexts). Contexts gives a list of currently defined contexts, and Names gives a list of symbols matching a string pattern. Mathematica built in objects are located in the System` context, and user defined variables are by default in Global`, so you can delete all your variables by e.g. Remove["Global`*"]. You can use Context to tell the context a symbol is in. There doesn't seem to be a formal command for unloading a package, but you could try e.g. Remove["PackageContext`*"]

see tutorial/Contexts and tutorial/ContextsAndPackages for more details.

$Packages
Contexts["Alg*"]
Names["Global`*"]

Grad[f] (* grad not defined yet, so this created Global`Grad *)
Names["Global`*"]
Needs["VectorAnalysis`"] (* I get a warning *)
Remove["Global`Grad"]; (* erase the problem def *)
Grad[Xx+Yy+(Zz^3)] (* uses the correct definition of Grad *)

Needs["Combinatorica`"];
ShowGraph[RankedEmbedding[GridGraph[5,5],{13}]] (* example taken from Combinatorica/tutorial/Combinatorica *)

Needs["GraphUtilities`"];
ExpressionTreePlot[Hold@Table[(1/(Sqrt[6]))(E^(I π j k/3)),{j,1,6}, {k,1,6}]] (* check out the Tooltips that show the subexpressions under each node! *)
ExpressionTreePlot[Table[(1/(Sqrt[6]))(E^(I π j k/3)),{j,1,6}, {k,1,6}]]