[metapost] Accuracy of "angle"

[Dan Luecking]

[Whoops, forgot the copy to the list]

At 03:48 AM 10/26/2010, Dirk Laurie wrote:
>I am at present revising my package Napoleon which makes geometry
>drawings using Metapost.  I noticed that when the drawing is large,
>the incircle of a triangle fails to touch the sides to an unacceptably
>large extent.  I traced the problem to inaccurate bisectors.  I used
>the following equation for a point z on the bisector of angle z0--z1--z2:
>
>     z = z1 + whatever*((z1-z0) rotated (angle((z2-z1) _zdiv (z0-z1))/2))
>
>where a _zdiv b is complex division,

This not the best way to get the difference in angles,
see below for another.

>I then changed the equation to one based on the theorem that the
>diagonals of a rhombus bisect its angles:
>
>     z = z1 + whatever*(unitvector(z0-z1)+unitvector(z2-z1))
>
>and the incircle is now accurate.
>
>The unitvector computations basically need to do the same as what
>_zdiv does, so I suspect the "angle" function as the source of the
>inaccuracy.


Well, unitvector is a primitive, and presumably optimized
for accuracy. Your _zdiv involves the solution of equations
(to get an inverse). Since solving involves several steps,
inaccuracies can accumulate.

Accuracy of any computation with numbers of finite precision is
a rather deep subject. There are several possible sources of
inaccuracy:
   - coordinates of points are already only known to +- 2^{-16}.
   - coordinates of added/subtracted points might have half that
     precision (+- 2^{-15}).
   - angles of pairs can be no more accurate than the *relative*
     errors in the pairs. That is, roughly, if x is +- e1 and y is +- e2
     then angle z is +- (angle z)*(abs (e1,e2))/(abs z).
   - in solving a set of equations, the determinant of the system
     contributes to inaccuracy. For intersecting lines, if they
     intersect at a small angle, there is greater relative error.
   - Similarly for the inverse of a transformation.

And there are more. I generally try to minimize divisions,
they contribute to inacuracy in two ways: division by a small
denominator magnifies the uncertainty in the numerator, and the
denominator contributes its own uncertainty.

I have experimented with angle to find its accuracy and it
seems to be as accurate as it is possible to be with 2^{16}
accuracy of data.

Here is how I have done this in the past:

% First find the angle from one direction to another.
% Simple, painless, no division:
vardef angle_from_to (expr U,V) =
   angle ( V rotated -angle U)
enddef;

% Get the direction of the bisector between two direction
% vectors. Operator precedence requires the parentheses:
vardef bisector (expr U,V) =
   U rotated (.5*angle_from_to(U,V))
enddef;

% Here is one way to get the incircle:
% First get the center:
vardef incenter (expr A,B,C) =
   save Z; pair Z;
   Z = A + whatever*bisector(B-A,C-A)
     = B + whatever*bisector(C-B,A-B);
   Z
enddef;

% To get the radius I need the distance from a point
% to any of the three lines. For ideosyncratic reasons I
% first define  a command to get the distance from a point
% to a line  through the origin in a given direction
vardef distance_to_line (expr P, U) =
   abs P +-+ abs (P dotprod unitvector U)
enddef;

% Draw the incircle for the triangle with corners z1, z2, z3.
% For best accuracy, one might compute the center 3 times
% (once for each cyclic permutation of z1,z2,z3) and average
% the results. The same for the radius.
beginfig(0);
   z1 = (0,0);
   z2 = (0,72);
   z3 = (18,70);
save C,R; pair C;
   C := incenter(z1,z2,z3);
   R := distance_to_line(C - z1, z2 - z1);
     draw z1--z2--z3--cycle;
     draw fullcircle scaled 2R shifted C;
endfig;
end

Another approach I have used is to use equations to
find the points where the incircle is tangent to the
sides of the triangle (e.g., the points on AB and AC
are equidistant from A). Then mimic the definition of
quartercircle to draw arcs from one such point to
the next having the required angles (180 minus the angle
of the corresponding corners of the triangle). This
guarantees that the resulting curve passes through
the relevant points. Plus, it is visually
indistinguishable from a circle. See grafbase.mp for
the full obfuscated details. Search for incircle (note that
that function expects cyclic path of length three rather
than three points).


Dan

Daniel H. Luecking
Department of Mathematical Sciences
Fayetteville, Arkansas
http://www-cs-faculty.stanford.edu/~knuth/iaq.html