### Turning Number and QuickSort

luigi scarso <luigi.scarso <at> gmail.com>

2015-01-21 08:14:27 GMT

It seems that ￼the University College Cork has some problem to deliver emails to the metapost list,

so I forward some notes that

Marc van Dongen has sent to me last days:

"""

I have noticed lots of problems with standard metapost, including

lots of problems related to rounding and an error in turningnumber.

I re-implemented turningnumber (in mpost) and it seems to work for

my purpose. The implementation is quite simple and it avoids rounding

errors. (It's still possible to reduce more rounding but I don't

think it's needed.) Please let me know if you want the implementation.

(I still have to do some proper testing on it.)

"""

The files are in attachment.

--

luigi

% This is a simple implementation to determine whether a path
% is drawn anti clockwise. To do this, it computes the signed
% area of the path consisting of the support points of the path.
% ASSUMPTION: ___pat___ does have a non-zero area.
vardef IsAntiClockwise( expr pat ) =
((cycle pat) and (Area( pat ) > AREA_ACCURACY))
%(turningnumber( pat ) > 0)
enddef;
vardef Area( expr pat ) =
save a, len, sum; numeric a[], len[], sum;
len := length pat;
% First we compute (twice) the signed areas of the triangles
% _(0,0) -- (point i of pat) -- (point (i + 1) of pat) -- cycle_
for i = len - 1 downto 0:
a[ i ] := Det( point i of pat, point (i + 1) of pat );
endfor
% Since adding the sizes may result in large rounding errors,
% we sort them in descending order
SortAbsGeq( a )( len );
% The remaining statements compute the sum of the signed areas
% in a ``robust'' fashion, trying to minimise rounding. A little
% bit more robust would be to first compute the frequencies of
% the entries in ___a___, then multiplying each unique number in
% ___a___ by its frequency, then sorting the products, and finally
% adding the products.
% e.g. [1,2,1,1];
% sort it -> [2,1,1,1];
% compute frequencies -> [2:1,1:3];
% compute products -> [2,3];
% sort -> [2,3];
% add -> 5.
sum := 0;
for i = len - 1 downto 0:
sum := sum + a[ i ];
endfor
% We could have just returned ___sum___ but
% this way we get the actual signed area.
0.5sum
enddef;
endinput

vardef QuickSort( suffix a, cmp )( expr size ) =
QuickSort_( a, cmp )( 0, size - 1 );
enddef;
vardef QuickSort_( suffix a, cmp )( expr lo, hi ) =
if lo < hi:
save sep; numeric sep;
%Swop( a )( lo, floor( (lo + hi) / 2 ) );
Swop( a )( lo, lo + floor( uniformdeviate( hi - lo - 1 ) ) );
sep := Partition( a, cmp, lo, hi );
Swop( a )( lo, sep );
QuickSort_( a, cmp )( lo, sep - 1 );
QuickSort_( a, cmp )( sep + 1, hi );
fi
enddef;
vardef Partition( suffix a, cmp )( expr lo, hi ) =
save dest, pivot; numeric dest, pivot;
dest := lo;
pivot := a[ lo ];
for i = lo + 1 upto hi:
if cmp( a[ i ], pivot ):
dest := dest + 1;
Swop( a )( dest, i );
fi
endfor
dest
enddef;
vardef Swop( suffix a )( expr fst, snd ) =
save tmp; numeric tmp;
tmp := a[ fst ];
a[ fst ] := a[ snd ];
a[ snd ] := tmp;
enddef;
vardef SortAbsGeq( suffix a )( expr size ) =
save CMP__;
vardef CMP__( expr fst, snd ) =
(abs( fst ) >= abs( snd ))
enddef;
QuickSort( a, CMP__ )( size );
enddef;
endinput

--
http://tug.org/metapost/