[metapost] precision of division

[Dan Luecking]

At 10:15 AM 10/10/2006, you wrote:
>Dear users of metapost,
>
>   I'm not an expert but I have used metapost for some drawings and, 
> recently,
>I've faced a problem when I try to build cycles involving points on a circle.
>You can see it on the following pdf file:
>http://www-sop.inria.fr/everest/Clement.Hurlin/figure.pdf
>
>I've divided the quarter circle from p0 to p2 into six equal paths: q0 -- q1,
>q1 -- q2, ..., by doing:
>
>   pair q[]; % calcul des points d'approximation
>   for $:=0 upto i:
>     q[$] = point $ * length pt[1]/i of pt[1];
>
>          (1)

These won't necessarily be "equal" parts. To get that you would use
"p0 rotated (90/6)". What you have is "equal time" subdivisions.
Anyway, MP's quartercircle (if that's how you got pt[1]) has two segments,
so length pt[1] is 2, and pt[1]/i = 1/3 which MP stores internally
as 21845/65536 (actually as the integer 21845). That happens to be
accurate to 1 part in 196608, i.e., plus or minus .000005.

In general, simple calculations are rounded to 16 binary places, and 
so are accurate to plus or minus 2^{-17}

>  I try to fill the different parts of the camembert by doing:
>   % calcul des parts
>   for $:=0 upto i - 1:
>     show $;
>     qt[$] := q[$] -- p[3] -- q[$+1];
>     pt[3] := pt[1] cutbefore q[$] cutafter q[$+1];

However, path intersection is a different issue, and may be off by
4 * 2^{-16) (in the time variable). MP may decide that paths do
intersect which actually do not, and vice versa.

>     show point infinity of qt[$];
>     show point infinity of (pt[3]);
>
>     show point 0 of  qt[$];
>     show point 0 of (pt[3]);
>     qt[$] := buildcycle(qt[$], pt[3]);

Like cutbefore/after, buildcycle relies on detecting path intersection.

>The problem is that, sometimes the buildcycle method fails to do the job (as
>you can see on the sample for the second part, which isn't filled) because, I
>think, of insufficient precision in the division (1).
>I think so because the points where paths should intersect exhibit a small
>difference when it fails: (you can see values outputted by the multiple show
>commands of the last loop mentionned, when compiling the sample .mp file):
>
> >> 0
> >> (-136.87585,-36.78265)
> >> (-136.87585,-36.78265)
> >> (-141.73225,0)
> >> (-141.73225,0)
> >> 1
> >> (-122.79346,-70.77658)
> >> (-122.79332,-70.77682)
> >> (-136.87585,-36.78265)
> >> (-136.87585,-36.78265)
> >> 2
>....
>
>You can see that the first buildcycle succeed because extremities of the two
>paths are the same.
>Yet for the second one it fails, because the first coordinate of the first
>extremity are not the same: -122.79346 != -122.79332
>Plus by changing i and k and the ratio of my figure, it leads to hazardous
>results (sometimes works fine, sometimes many parts of the camembert are not
>filled despite they should be).

One point is calculated from the division, the other from the path
intersection operation. Both are slightly off, but the latter usually
more than the former.

One should never actually rely on MP to detect that paths intersect
except in one circumstance: if the paths intersect transversally
(i.e., not tangentially) at a points on each that are strictly between
nodes (i.e., at non-integer times). Even then the points of
intersection that one obtains need not be strictly equal.


>  does it make sense and how can I fix it ?

Does what make sense?

Here is what I would do. Assuming that you really want the points
   q[$] := point $*length pt[1]/i of pt[1];
I would save those times instead of the points:
   t[$] := $*length pt[1]/i;
then if I need to fill a wedge, I would write:
   fill p3-- subpath(t[$],t[$+1]) of p --cycle;


Dan


Daniel H. Luecking
Department of Mathematical Sciences
University of Arkansas
"I reject your reality, and substitute my own." -- Adam Savage