
(* :Name: NumericalMath`Microscope` *)

(* :Title: Microscopic Examination of Roundoff Errors *)

(* :Author: Jerry B. Keiper *)

(* :Summary:
	Functions can be plotted on a microscopic scale to exhibit the
	granularity of machine arithmetic.  Alternatively the actual
	error (measured in upls (units last place)) can be plotted.
*)

(* :Context: NumericalMath`Microscope` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 1992  Wolfram Research, Inc.
	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Revised by Jerry B. Keiper, December, 1991.
	Originally sketched by Jerry B. Keiper, January, 1991
*)

(* :Keywords: machine precision, roundoff error, granularity of machine
	numbers
*)

(* :Source:
*)

(* :Warnings:
*)

(* :Mathematica Version: 2.0 *)

(* :Limitations:  The size of an ulp is defined relative to the size of the
	numbers being considered.  In small neighborhoods of powers of 2
	the size of an ulp changes and it is not clear what to do about
	error plots when this happens.  This is even more of a problem
	when looking at the error near a zero of a function where the
	relative error may even be unbounded.
*)

(* :Discussion:
	This package defines four functions: Ulp[ ], MachineError[ ],
	Microscope[ ] and MicroscopicError[ ].  Ulp[ ] gives the size of an
	ulp for numbers near the argument.  MachineError[ ] gives the error
	involved in evaluating its argument using machine arithmetic.
	Microscope[ ] plots an expression at a magnification such that the
	granularity of machine arithmetic becomes obvious.  MicroscopicError[ ]
	uses high precision (30 digits) to calculate the "true" value and plots
	the difference between the machine precision result and the "true" value
	measured in ulps (units last place).

	Both Microscope[ ] and MicroscopicError[ ] work by evaluating the
	expression at consecutive machine numbers in a neighborhood of
	a point to come up with a set of discrete points.  The option
	PlotJoined can be True (the points are to be simply connected with
	straight line segments), False (the points are not to be connected),
	or Real.  The choice of Real causes the resulting plot to represent
	a mapping from the set of real numbers to the set of machine numbers
	(Microscope[ ]) or the error in such a mapping (MicroscopicError[ ]).
	Note that care must be exercised in interpreting the result when
	PlotJoined -> Real is used.  Functions in computer libraries are
	designed to map the set of machine numbers into the set of machine
	numbers.  It is not fair to blame the function for the error incurred
	in rounding a real number to the nearest machine number prior to the
	function receiving the machine number.
*)

BeginPackage["NumericalMath`Microscope`"]

Ulps::usage =
"Ulps is a unit of error in machine arithmetic.  An ulp is the distance
between two consecutive machine numbers and varies with the size of the
numbers involved."

Ulp::usage =
"Ulp[x] gives the size of an ulp for numbers near x."

MachineError::usage =
"MachineError[f, x -> a] gives the error involved in evaluating f at x == a
using machine arithmetic."

Microscope::usage =
"Microscope[f, {x, a}] plots the expression f (which is supposed to contain
a single variable x) in a small neighborhood of a.  Microscope[f, {x, a, n}]
plots f from a - n ulps to a + n ulps."

MicroscopicError::usage =
"MicroscopicError[f, {x, a}] plots the error incurred by using machine
arithmetic to evaluate the expression f (which is supposed to contain a single
variable x) in a small neighborhood of a.  MicroscopicError[f, {x, a, n}]
plots the error from a - n ulps to a + n ulps.  The scale on the vertical axis
represents ulps in the result."

Options[MicroscopicError] = Options[Microscope] = PlotJoined -> True;
Attributes[MicroscopicError] = Attributes[MachineError] = HoldFirst;

Begin["NumericalMath`Microscope`Private`"]

Ulp[x_] := ulp[x];

MachineError[f_, x_ -> a_] :=
    Module[{aa = N[SetPrecision[a, 30], 30], ff = Hold[f]},
	ff = Release[SetPrecision[f, 30]];
	ff = Re[N[ff /. x -> aa, 30]];
	(SetPrecision[N[f /. x -> a], 30] - ff)/ulp[ff] Ulps
	];

Microscope[e_, {x_Symbol, a_, n_Integer:30}, opt___] :=
    Module[{joined, ans,},
	ans /; (joined = PlotJoined /. {opt} /. (PlotJoined -> True);
		MemberQ[{True, False, Real}, joined] &&
			(ans = mic[e, x, a, n, joined]; True))
	] /; NumberQ[N[a]];

MicroscopicError[ee_, {x_Symbol, a_, n_Integer:30}, opt___] :=
    Module[{joined, ans, e = Hold[ee]},
	ans /; (joined = PlotJoined /. {opt} /. (PlotJoined -> True);
		MemberQ[{True, False, Real}, joined] &&
			(ans = micer[e, x, a, n, joined]; True))
	] /; NumberQ[N[a]];

ulp[_DirectedInfinity] := DirectedInfinity[ ]

ulp[x_] := 
    Module[{t = Abs[N[x]], u},
	If[t < $MinMachineNumber,
		$MinMachineNumber,
		u = N[2^Floor[Log[2, t $MachineEpsilon] - .5]];
		t = t - Release[t - u];
		If[t == 0. || t == 2u, 2u, u]
		]
	];
 
machpts[e_, {x_, a_}, n_] :=
    Module[{h, xtab, ytab, i, na = N[a]},
        h = ulp[na (1-$MachineEpsilon)];
        xtab = Table[i h, {i, -n, n}] + na;
        ytab = e /. x -> xtab;
        {xtab, ytab, na, e /. x -> na}]; 
         
mic[e_, x_, a_, n_, joined_] :=
    Module[{h, xtab, ytab, i, x0, y0, pts, lines, ao}, 
        {xtab, ytab, x0, y0} = machpts[e, {x, a}, n];
	xtab = (xtab - x0)/ulp[x0];
	ytab -= y0;
	ytab *= N[2^-Round[Log[2, Max[Abs[Re[ytab]]]]]];
	ao = {Min[xtab], Min[Re[ytab]]};
	pts = {PointSize[.2/n],
		Table[Point[{xtab[[i]], ytab[[i]]}], {i, Length[xtab]}]};
	If[joined === Real,
            xtab = (Drop[xtab, 1] + Drop[xtab, -1])/2; 
            xtab = Flatten[Table[{xtab[[i]], xtab[[i]]}, {i, Length[xtab]}]];
            ytab = Flatten[Table[{ytab[[i]], ytab[[i]]}, {i, Length[ytab]}]]; 
            ytab = Drop[Drop[ytab, 1], -1]];   
	If[joined === False,
	    lines = {},
	    lines = {Thickness[.001],
		Line[Table[{xtab[[i]], ytab[[i]]}, {i, Length[xtab]}]]}];
	Show[Graphics[{pts, lines}, AxesOrigin -> ao, Axes -> True,
	    PlotRange -> All, Ticks -> {{{0,ToString[x0]}}, {{0,ToString[y0]}}}]]
	];
                 
micer[e_, x_, a_, n_, joined_] :=
    Module[{h, xtab, ytab, yytab, i, x0, y0, pts, lines, ao, sxtab, ee},
        {xtab, ytab, x0, y0} = machpts[Release[e], {x, a}, n];
        xtab = SetPrecision[xtab, 30];
        ytab = SetPrecision[ytab, 30];
	ee = Release[SetPrecision[e, 30]];
        yytab = N[ytab - Re[(ee /. x -> xtab)]]/ulp[y0];
	sxtab = N[xtab - SetPrecision[x0, 30]]/ulp[x0];
	ao = {Min[sxtab], 0};
	pts = {PointSize[.2/n],
		Table[Point[{sxtab[[i]], yytab[[i]]}], {i, Length[sxtab]}]};
	If[joined === Real,
            xtab = (Drop[xtab, 1] + Drop[xtab, -1])/2;
            yytab = ee /. x -> xtab; 
            xtab = Flatten[Table[{xtab[[i]], xtab[[i]]}, {i, Length[xtab]}]];
            yytab = Flatten[Table[{yytab[[i]],yytab[[i]]}, {i,Length[yytab]}]];
            ytab = Flatten[Table[{ytab[[i]], ytab[[i]]}, {i, Length[ytab]}]];
            ytab = Drop[Drop[ytab, 1], -1];
            yytab = N[ytab - yytab]/ulp[y0];
	    sxtab = N[xtab - SetPrecision[x0, 30]]/ulp[x0]
	    ];
	If[joined === False,
	    lines = {},
	    lines = {Thickness[.001],
		Line[Table[{sxtab[[i]], yytab[[i]]}, {i, Length[xtab]}]]}];
	Show[Graphics[{pts, lines}, AxesOrigin -> ao, Axes -> True,
	    PlotRange -> All, Ticks -> {{{0,ToString[x0]}}, Automatic}]]
	];

End[]  (* NumericalMath`Microscope`Private` *)

Protect[Microscope, MicroscopicError];

EndPackage[] (* NumericalMath`Microscope` *)

(* Tests:
Microscope[Log[x], {x, 7}]
Microscope[Log[x], {x, 7, 20}]
Microscope[Log[x], {x, 7}, PlotJoined -> False]
Microscope[Log[x], {x, 7, 20}, PlotJoined -> False]
Microscope[Log[x], {x, 7}, PlotJoined -> True]
Microscope[Log[x], {x, 7, 20}, PlotJoined -> True]
Microscope[Log[x], {x, 7}, PlotJoined -> Real]
Microscope[Log[x], {x, 7, 20}, PlotJoined -> Real]
Microscope[Sqrt[x], {x, 5}]
Microscope[Sqrt[x], {x, 5, 20}]
Microscope[Sqrt[x], {x, 5}, PlotJoined -> False]
Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> False]
Microscope[Sqrt[x], {x, 5}, PlotJoined -> True]
Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> True]
Microscope[Sqrt[x], {x, 5}, PlotJoined -> Real]
Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .5}]
Microscope[ArcTanh[x], {x, .5, 20}]
Microscope[ArcTanh[x], {x, .5}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .5}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .5}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .05}]
Microscope[ArcTanh[x], {x, .05, 20}]
Microscope[ArcTanh[x], {x, .05}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .05}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .05}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .00005}]
Microscope[ArcTanh[x], {x, .00005, 20}]
Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> False]
Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> True]
Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> Real]
Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> Real]

MicroscopicError[Log[x], {x, 7}]
MicroscopicError[Log[x], {x, 7, 20}]
MicroscopicError[Log[x], {x, 7}, PlotJoined -> False]
MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> False]
MicroscopicError[Log[x], {x, 7}, PlotJoined -> True]
MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> True]
MicroscopicError[Log[x], {x, 7}, PlotJoined -> Real]
MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> Real]
MicroscopicError[Sqrt[x], {x, 5}]
MicroscopicError[Sqrt[x], {x, 5, 20}]
MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> False]
MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> False]
MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> True]
MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> True]
MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> Real]
MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .5}]
MicroscopicError[ArcTanh[x], {x, .5, 20}]
MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .05}]
MicroscopicError[ArcTanh[x], {x, .05, 20}]
MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .00005}]
MicroscopicError[ArcTanh[x], {x, .00005, 20}]
MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> False]
MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> True]
MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> Real]
MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> Real]
*)

