21 Jan 09:14 2015
Turning Number and QuickSort
luigi scarso <luigi.scarso <at> gmail.com>
2015-01-21 08:14:27 GMT
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.
% 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