[metapost] Turning Number and QuickSort

[luigi scarso]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://tug.org/pipermail/metapost/attachments/20150121/562008e0/attachment.html>
-------------- next part --------------
% 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
-------------- next part --------------
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