MetaPost, Plotting, and numerical precision

So, to write up diagrams in LaTeX, you need to use Metapost. But Metapost doesn't use floating point arithmetic.

As Claudio Beccari's "Floating point numbers and METAFONT, METAPOST, TEX, and PostScript Type 1 fonts" (TUGboat [pdf]) notes, 32 bit integers are used to represent real numbers. The first 16 bits form the fractional part, 14 bits the integer part, 1 bit for the sign, and 1 bit for special purposes.

So, that means we have 16 log(2)/log(10) digits of precision, or about 4 digits. Now lets remember:

1 PS point = 1.00375 points
1 pica = 12 PS points
1 inch = 72 PS points = 72.27 points = 6 pica
1 cm = 28.3464567 PS points = 2.36220472 pica

We have precision of 2^-16 points, or about 0.000868055556 inches, or 0.00220486111 centimeters.

That's decent for output but not for intermediate computations. For example, if we were to plot x^x, we may lose a lot of precision.

Plots in Metapost

Lets consider a simple plot of $f(x)=x^{2}$.

numeric u;
u := 1pc; % units

vardef f(expr x) = x*x enddef;

beginfig(0)
  % draw the axes
  drawdblarrow (-3u-ahlength,0)--(3u+ahlength,0);
  drawdblarrow (0,0-2ahlength)--(0,9u+ahlength);

  % plot the function
  draw (-3u,f(-3)*u)
    for i=-3+0.05 step 0.05 until 3:
    ..(i*u,f(i)*u)
  endfor;
endfig;
end;

Remember that ahlength is the length of the arrow head.

This basic scheme can be generalized if we add numerics x0 and x1 which control where the plot begins and ends (respectively), as well as the step size dx which is taken to be "small enough".

Revising our code:

numeric u;
numeric dx;
u := 1pc; % units
dx := 0.05; 

vardef f(expr x) = x*x enddef;

beginfig(0)
  numeric x[];

  x0 := -3; % start plotting at x=-3
  x1 := 3; % stop plotting at x=+3

  % draw the axes
  drawdblarrow (x0*u-ahlength,0)--(x1*u+ahlength,0); % x-axis
  drawdblarrow (0,0-2ahlength)--(0,f(x1)*u+ahlength); % y-axis

  % plot the function
  draw (x0*u,f(x0)*u)
    for i=x0+dx step dx until 3:
    ..(i*u,f(i)*u)
  endfor;
endfig;
end;

This makes things a little complicated. What we are doing is computing the pairs (x,y) and then scaling them, then plotting.

The dx is the change in x before scaling. The points plotted have a change in x that amounts to dx*u=0.6pt approximately.

But we can do more! If we specify how big we want this plot to be, i.e. it has to fit within X inches, then we can determine the scale u by this.

Specifically, u := X/(x[1]-x[0]) is the scale factor definition.

If we demand that dx*u=0.6pt hold, which is "sufficiently good" for practical purposes, then we also define dx := (3pt)/(5*u).

The interested reader may want to read Learning Metapost by Doing.