home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / NUMERICA.PAK / MICROSCO.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  11.9 KB  |  294 lines
  1.  
  2. (* :Name: NumericalMath`Microscope` *)
  3.  
  4. (* :Title: Microscopic Examination of Roundoff Errors *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Functions can be plotted on a microscopic scale to exhibit the
  10.     granularity of machine arithmetic.  Alternatively the actual
  11.     error (measured in upls (units last place)) can be plotted.
  12. *)
  13.  
  14. (* :Context: NumericalMath`Microscope` *)
  15.  
  16. (* :Package Version: 1.0 *)
  17.  
  18. (* :Copyright: Copyright 1992  Wolfram Research, Inc.
  19.     Permission is hereby granted to modify and/or make copies of
  20.     this file for any purpose other than direct profit, or as part
  21.     of a commercial product, provided this copyright notice is left
  22.     intact.  Sale, other than for the cost of media, is prohibited.
  23.  
  24.     Permission is hereby granted to reproduce part or all of
  25.     this file, provided that the source is acknowledged.
  26. *)
  27.  
  28. (* :History:
  29.     Revised by Jerry B. Keiper, December, 1991.
  30.     Originally sketched by Jerry B. Keiper, January, 1991
  31. *)
  32.  
  33. (* :Keywords: machine precision, roundoff error, granularity of machine
  34.     numbers
  35. *)
  36.  
  37. (* :Source:
  38. *)
  39.  
  40. (* :Warnings:
  41. *)
  42.  
  43. (* :Mathematica Version: 2.0 *)
  44.  
  45. (* :Limitations:  The size of an ulp is defined relative to the size of the
  46.     numbers being considered.  In small neighborhoods of powers of 2
  47.     the size of an ulp changes and it is not clear what to do about
  48.     error plots when this happens.  This is even more of a problem
  49.     when looking at the error near a zero of a function where the
  50.     relative error may even be unbounded.
  51. *)
  52.  
  53. (* :Discussion:
  54.     This package defines four functions: Ulp[ ], MachineError[ ],
  55.     Microscope[ ] and MicroscopicError[ ].  Ulp[ ] gives the size of an
  56.     ulp for numbers near the argument.  MachineError[ ] gives the error
  57.     involved in evaluating its argument using machine arithmetic.
  58.     Microscope[ ] plots an expression at a magnification such that the
  59.     granularity of machine arithmetic becomes obvious.  MicroscopicError[ ]
  60.     uses high precision (30 digits) to calculate the "true" value and plots
  61.     the difference between the machine precision result and the "true" value
  62.     measured in ulps (units last place).
  63.  
  64.     Both Microscope[ ] and MicroscopicError[ ] work by evaluating the
  65.     expression at consecutive machine numbers in a neighborhood of
  66.     a point to come up with a set of discrete points.  The option
  67.     PlotJoined can be True (the points are to be simply connected with
  68.     straight line segments), False (the points are not to be connected),
  69.     or Real.  The choice of Real causes the resulting plot to represent
  70.     a mapping from the set of real numbers to the set of machine numbers
  71.     (Microscope[ ]) or the error in such a mapping (MicroscopicError[ ]).
  72.     Note that care must be exercised in interpreting the result when
  73.     PlotJoined -> Real is used.  Functions in computer libraries are
  74.     designed to map the set of machine numbers into the set of machine
  75.     numbers.  It is not fair to blame the function for the error incurred
  76.     in rounding a real number to the nearest machine number prior to the
  77.     function receiving the machine number.
  78. *)
  79.  
  80. BeginPackage["NumericalMath`Microscope`"]
  81.  
  82. Ulps::usage =
  83. "Ulps is a unit of error in machine arithmetic.  An ulp is the distance
  84. between two consecutive machine numbers and varies with the size of the
  85. numbers involved."
  86.  
  87. Ulp::usage =
  88. "Ulp[x] gives the size of an ulp for numbers near x."
  89.  
  90. MachineError::usage =
  91. "MachineError[f, x -> a] gives the error involved in evaluating f at x == a
  92. using machine arithmetic."
  93.  
  94. Microscope::usage =
  95. "Microscope[f, {x, a}] plots the expression f (which is supposed to contain
  96. a single variable x) in a small neighborhood of a.  Microscope[f, {x, a, n}]
  97. plots f from a - n ulps to a + n ulps."
  98.  
  99. MicroscopicError::usage =
  100. "MicroscopicError[f, {x, a}] plots the error incurred by using machine
  101. arithmetic to evaluate the expression f (which is supposed to contain a single
  102. variable x) in a small neighborhood of a.  MicroscopicError[f, {x, a, n}]
  103. plots the error from a - n ulps to a + n ulps.  The scale on the vertical axis
  104. represents ulps in the result."
  105.  
  106. Options[MicroscopicError] = Options[Microscope] = PlotJoined -> True;
  107. Attributes[MicroscopicError] = Attributes[MachineError] = HoldFirst;
  108.  
  109. Begin["NumericalMath`Microscope`Private`"]
  110.  
  111. Ulp[x_] := ulp[x];
  112.  
  113. MachineError[f_, x_ -> a_] :=
  114.     Module[{aa = N[SetPrecision[a, 30], 30], ff = Hold[f]},
  115.     ff = Release[SetPrecision[f, 30]];
  116.     ff = Re[N[ff /. x -> aa, 30]];
  117.     (SetPrecision[N[f /. x -> a], 30] - ff)/ulp[ff] Ulps
  118.     ];
  119.  
  120. Microscope[e_, {x_Symbol, a_, n_Integer:30}, opt___] :=
  121.     Module[{joined, ans,},
  122.     ans /; (joined = PlotJoined /. {opt} /. (PlotJoined -> True);
  123.         MemberQ[{True, False, Real}, joined] &&
  124.             (ans = mic[e, x, a, n, joined]; True))
  125.     ] /; NumberQ[N[a]];
  126.  
  127. MicroscopicError[ee_, {x_Symbol, a_, n_Integer:30}, opt___] :=
  128.     Module[{joined, ans, e = Hold[ee]},
  129.     ans /; (joined = PlotJoined /. {opt} /. (PlotJoined -> True);
  130.         MemberQ[{True, False, Real}, joined] &&
  131.             (ans = micer[e, x, a, n, joined]; True))
  132.     ] /; NumberQ[N[a]];
  133.  
  134. ulp[_DirectedInfinity] := DirectedInfinity[ ]
  135.  
  136. ulp[x_] := 
  137.     Module[{t = Abs[N[x]], u},
  138.     If[t < $MinMachineNumber,
  139.         $MinMachineNumber,
  140.         u = N[2^Floor[Log[2, t $MachineEpsilon] - .5]];
  141.         t = t - Release[t - u];
  142.         If[t == 0. || t == 2u, 2u, u]
  143.         ]
  144.     ];
  145.  
  146. machpts[e_, {x_, a_}, n_] :=
  147.     Module[{h, xtab, ytab, i, na = N[a]},
  148.         h = ulp[na (1-$MachineEpsilon)];
  149.         xtab = Table[i h, {i, -n, n}] + na;
  150.         ytab = e /. x -> xtab;
  151.         {xtab, ytab, na, e /. x -> na}]; 
  152.          
  153. mic[e_, x_, a_, n_, joined_] :=
  154.     Module[{h, xtab, ytab, i, x0, y0, pts, lines, ao}, 
  155.         {xtab, ytab, x0, y0} = machpts[e, {x, a}, n];
  156.     xtab = (xtab - x0)/ulp[x0];
  157.     ytab -= y0;
  158.     ytab *= N[2^-Round[Log[2, Max[Abs[Re[ytab]]]]]];
  159.     ao = {Min[xtab], Min[Re[ytab]]};
  160.     pts = {PointSize[.2/n],
  161.         Table[Point[{xtab[[i]], ytab[[i]]}], {i, Length[xtab]}]};
  162.     If[joined === Real,
  163.             xtab = (Drop[xtab, 1] + Drop[xtab, -1])/2; 
  164.             xtab = Flatten[Table[{xtab[[i]], xtab[[i]]}, {i, Length[xtab]}]];
  165.             ytab = Flatten[Table[{ytab[[i]], ytab[[i]]}, {i, Length[ytab]}]]; 
  166.             ytab = Drop[Drop[ytab, 1], -1]];   
  167.     If[joined === False,
  168.         lines = {},
  169.         lines = {Thickness[.001],
  170.         Line[Table[{xtab[[i]], ytab[[i]]}, {i, Length[xtab]}]]}];
  171.     Show[Graphics[{pts, lines}, AxesOrigin -> ao, Axes -> True,
  172.         PlotRange -> All, Ticks -> {{{0,ToString[x0]}}, {{0,ToString[y0]}}}]]
  173.     ];
  174.                  
  175. micer[e_, x_, a_, n_, joined_] :=
  176.     Module[{h, xtab, ytab, yytab, i, x0, y0, pts, lines, ao, sxtab, ee},
  177.         {xtab, ytab, x0, y0} = machpts[Release[e], {x, a}, n];
  178.         xtab = SetPrecision[xtab, 30];
  179.         ytab = SetPrecision[ytab, 30];
  180.     ee = Release[SetPrecision[e, 30]];
  181.         yytab = N[ytab - Re[(ee /. x -> xtab)]]/ulp[y0];
  182.     sxtab = N[xtab - SetPrecision[x0, 30]]/ulp[x0];
  183.     ao = {Min[sxtab], 0};
  184.     pts = {PointSize[.2/n],
  185.         Table[Point[{sxtab[[i]], yytab[[i]]}], {i, Length[sxtab]}]};
  186.     If[joined === Real,
  187.             xtab = (Drop[xtab, 1] + Drop[xtab, -1])/2;
  188.             yytab = ee /. x -> xtab; 
  189.             xtab = Flatten[Table[{xtab[[i]], xtab[[i]]}, {i, Length[xtab]}]];
  190.             yytab = Flatten[Table[{yytab[[i]],yytab[[i]]}, {i,Length[yytab]}]];
  191.             ytab = Flatten[Table[{ytab[[i]], ytab[[i]]}, {i, Length[ytab]}]];
  192.             ytab = Drop[Drop[ytab, 1], -1];
  193.             yytab = N[ytab - yytab]/ulp[y0];
  194.         sxtab = N[xtab - SetPrecision[x0, 30]]/ulp[x0]
  195.         ];
  196.     If[joined === False,
  197.         lines = {},
  198.         lines = {Thickness[.001],
  199.         Line[Table[{sxtab[[i]], yytab[[i]]}, {i, Length[xtab]}]]}];
  200.     Show[Graphics[{pts, lines}, AxesOrigin -> ao, Axes -> True,
  201.         PlotRange -> All, Ticks -> {{{0,ToString[x0]}}, Automatic}]]
  202.     ];
  203.  
  204. End[]  (* NumericalMath`Microscope`Private` *)
  205.  
  206. Protect[Microscope, MicroscopicError];
  207.  
  208. EndPackage[] (* NumericalMath`Microscope` *)
  209.  
  210. (* Tests:
  211. Microscope[Log[x], {x, 7}]
  212. Microscope[Log[x], {x, 7, 20}]
  213. Microscope[Log[x], {x, 7}, PlotJoined -> False]
  214. Microscope[Log[x], {x, 7, 20}, PlotJoined -> False]
  215. Microscope[Log[x], {x, 7}, PlotJoined -> True]
  216. Microscope[Log[x], {x, 7, 20}, PlotJoined -> True]
  217. Microscope[Log[x], {x, 7}, PlotJoined -> Real]
  218. Microscope[Log[x], {x, 7, 20}, PlotJoined -> Real]
  219. Microscope[Sqrt[x], {x, 5}]
  220. Microscope[Sqrt[x], {x, 5, 20}]
  221. Microscope[Sqrt[x], {x, 5}, PlotJoined -> False]
  222. Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> False]
  223. Microscope[Sqrt[x], {x, 5}, PlotJoined -> True]
  224. Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> True]
  225. Microscope[Sqrt[x], {x, 5}, PlotJoined -> Real]
  226. Microscope[Sqrt[x], {x, 5, 20}, PlotJoined -> Real]
  227. Microscope[ArcTanh[x], {x, .5}]
  228. Microscope[ArcTanh[x], {x, .5, 20}]
  229. Microscope[ArcTanh[x], {x, .5}, PlotJoined -> False]
  230. Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> False]
  231. Microscope[ArcTanh[x], {x, .5}, PlotJoined -> True]
  232. Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> True]
  233. Microscope[ArcTanh[x], {x, .5}, PlotJoined -> Real]
  234. Microscope[ArcTanh[x], {x, .5, 20}, PlotJoined -> Real]
  235. Microscope[ArcTanh[x], {x, .05}]
  236. Microscope[ArcTanh[x], {x, .05, 20}]
  237. Microscope[ArcTanh[x], {x, .05}, PlotJoined -> False]
  238. Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> False]
  239. Microscope[ArcTanh[x], {x, .05}, PlotJoined -> True]
  240. Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> True]
  241. Microscope[ArcTanh[x], {x, .05}, PlotJoined -> Real]
  242. Microscope[ArcTanh[x], {x, .05, 20}, PlotJoined -> Real]
  243. Microscope[ArcTanh[x], {x, .00005}]
  244. Microscope[ArcTanh[x], {x, .00005, 20}]
  245. Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> False]
  246. Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> False]
  247. Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> True]
  248. Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> True]
  249. Microscope[ArcTanh[x], {x, .00005}, PlotJoined -> Real]
  250. Microscope[ArcTanh[x], {x, .00005, 20}, PlotJoined -> Real]
  251.  
  252. MicroscopicError[Log[x], {x, 7}]
  253. MicroscopicError[Log[x], {x, 7, 20}]
  254. MicroscopicError[Log[x], {x, 7}, PlotJoined -> False]
  255. MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> False]
  256. MicroscopicError[Log[x], {x, 7}, PlotJoined -> True]
  257. MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> True]
  258. MicroscopicError[Log[x], {x, 7}, PlotJoined -> Real]
  259. MicroscopicError[Log[x], {x, 7, 20}, PlotJoined -> Real]
  260. MicroscopicError[Sqrt[x], {x, 5}]
  261. MicroscopicError[Sqrt[x], {x, 5, 20}]
  262. MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> False]
  263. MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> False]
  264. MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> True]
  265. MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> True]
  266. MicroscopicError[Sqrt[x], {x, 5}, PlotJoined -> Real]
  267. MicroscopicError[Sqrt[x], {x, 5, 20}, PlotJoined -> Real]
  268. MicroscopicError[ArcTanh[x], {x, .5}]
  269. MicroscopicError[ArcTanh[x], {x, .5, 20}]
  270. MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> False]
  271. MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> False]
  272. MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> True]
  273. MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> True]
  274. MicroscopicError[ArcTanh[x], {x, .5}, PlotJoined -> Real]
  275. MicroscopicError[ArcTanh[x], {x, .5, 20}, PlotJoined -> Real]
  276. MicroscopicError[ArcTanh[x], {x, .05}]
  277. MicroscopicError[ArcTanh[x], {x, .05, 20}]
  278. MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> False]
  279. MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> False]
  280. MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> True]
  281. MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> True]
  282. MicroscopicError[ArcTanh[x], {x, .05}, PlotJoined -> Real]
  283. MicroscopicError[ArcTanh[x], {x, .05, 20}, PlotJoined -> Real]
  284. MicroscopicError[ArcTanh[x], {x, .00005}]
  285. MicroscopicError[ArcTanh[x], {x, .00005, 20}]
  286. MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> False]
  287. MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> False]
  288. MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> True]
  289. MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> True]
  290. MicroscopicError[ArcTanh[x], {x, .00005}, PlotJoined -> Real]
  291. MicroscopicError[ArcTanh[x], {x, .00005, 20}, PlotJoined -> Real]
  292. *)
  293.  
  294.