When I first started learning about rule-based programming with Mathematica, I tried to translate this algorithm for computing the convex hull of a set of 2-D points in $O(n \log(n))$ time, to use rule-based replacement. For convenience, I'll paste the pseudocode of the algorithm below:

Input: a list P of points in the plane.

  • Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate).

  • Initialize U and L as empty lists. The lists will hold the vertices of upper and lower hulls respectively.

  • for i = 1, 2, ..., n:

    while L contains at least two points and the sequence of last two points of L and the point P[i] does not make a counter-clockwise turn:
    remove the last point from L
    append P[i] to L

  • for i = n, n-1, ..., 1:

    while U contains at least two points and the sequence of last two points of U and the point P[i] does not make a counter-clockwise turn:
    remove the last point from U
    append P[i] to U

  • Remove the last point of each list (it's the same as the first point of the other list).

  • Concatenate L and U to obtain the convex hull of P. Points in the result will be listed in counter-clockwise order.

(Actually, I used a version of this algorithm from the book "Computational Geometry" by de Berg & others, but it's pretty much the same thing except for the "counter-clockwise turn condition" that uses determinant instead of cross-product, which isn't really relevant to this question.)

This was my implementation, which has the same end result as the stated algorithm, but as you can see doesn't operate quite the same way:

convexHullPM[pts_] := Module[{li, isRightTurn}, 
  isRightTurn[p1_, p2_, p3_] := 
     Sign[Det@MapThread[Prepend, {{p1, p2, p3}, {1, 1, 1}}]] == -1;
  li = Sort[pts];
  li = Join[li, Reverse@Most@li];
  li //. {pre___, a_List, b_List, c_List, post___} :>
         {pre, a, c, post} /; Not[isRightTurn[a, b, c]]]

It works, giving the same sequence of points as Mathematica's own ConvexHull function from the Computational Geometry Package.

However it's slow (referring to its time complexity, as the input size increases); timings with increasingly larger sets of input points seem to indicate $O(n^2)$ complexity, which seems about right - because the way I constructed the rule in Mathematica, it starts matching from the start every time a replacement is made. Whereas the listed algorithm has the notion of a current point, and points are removed from the set (if need be) from the most recent points.

So what would be a better rule-based implementation?

  • $\begingroup$ Probably you need ReplacePart and FixedPoint. $\endgroup$
    – swish
    Commented May 25, 2013 at 10:33
  • $\begingroup$ @swish ReplacePart won't help, since it copies the entire list even when only one or a few points are modified. Besides, it does not delete points, so Delete will likely be more appropriate (but will have the same problem). $\endgroup$ Commented May 25, 2013 at 13:45
  • 1
    $\begingroup$ A few pointers: 1. you can compactly implement your test, like so (though I have already incorporated the negation): notRightTurn = NonNegative[Det[PadLeft[{##}, {3, 3}, 1]]] &;; 2. Here's a shorter way to generate li: n = Length[pts]; li = ArrayPad[Sort[pts], {{0, n - 1}}, "Reflected"] $\endgroup$ Commented May 25, 2013 at 15:14
  • $\begingroup$ @J.M. Thanks. Mathematica seems to have more functions than I have brain cells! $\endgroup$
    – Aky
    Commented May 25, 2013 at 15:22

1 Answer 1



I suggest to use linked lists and recursion. Note that this won't be the fastest possible solution, but it will be idiomatic, rule-based, will have the right asymptotic complexity, and will have decent performance. I want to stress that I am here after these features, rather than raw performance, as those seems to be at the core of your question.

Linked lists

Here is the linked list API we will need:

ClearAll[ll, toLinkedList, fromLinkedList];
SetAttributes[ll, HoldAllComplete];
toLinkedList[l_List] := Fold[ll[#2, #1] &, ll[], Reverse@l];
fromLinkedList[lst_ll] := List @@ Flatten[lst, Infinity, ll];

This has been explained in more details in the linked post. We will need an additional function to reverse a linked list. One possible (idiomatic) implementation can be:

llrev[lst_ll] := llrev[lst, ll[]];
llrev[ll[], accum_] := accum;
llrev[ll[head_, tail_], accum_] := llrev[tail, ll[head, accum]];

(of course, one could also use Composition[toLinkedList, Reverse, fromLinkedList], which would also be somewhat faster, but I like the recursive version).

Some simple examples how this is used:

tst = toLinkedList[Range[5]]

(* ll[1, ll[2, ll[3, ll[4, ll[5, ll[]]]]]] *)


(* {1, 2, 3, 4, 5} *)


(* ll[5, ll[4, ll[3, ll[2, ll[1, ll[]]]]]] *)


Now we are ready to rewrite your algorithm. One subtlety here is that we'd either need to consider 4 adjacent points in this approach, or use FixedPoint (many runs). I would pick the FixedPoint route. Here is the code:

isRightTurn[p1_, p2_, p3_] := 
   Sign[Det@MapThread[Prepend, {{p1, p2, p3}, {1, 1, 1}}]] == -1;

and the main function:

convexHullLL[pts_List] :=
   convexHullLL[toLinkedList[Join[#, Reverse@Most@#] &[Sort[pts]]]];

convexHullLL[pts_ll] :=
    llrev[convexHullLL[#, ll[]]] &,

convexHullLL[ll[], accum_] := accum; 
convexHullLL[ll[a_, ll[b_, tail : ll[c_, _]]], accum_]/;Not[isRightTurn[a, b, c]] :=
   convexHullLL[ll[a, tail], accum];
convexHullLL[ll[head_, tail_], accum_] := convexHullLL[tail, ll[head, accum]];

The main recursive body here is typical for the solutions based on linked lists, and the ideas behind this technique were explained in more details in the linked post.

Examples and benchmarks

We can test that it gives the same results as your code:

rpts = RandomInteger[{-10,10},{100,2}];

(* True *)

But it has a different complexity:

rptsLarge = RandomReal[{-5,5},{2*10^3,2}];
(cnv1 = convexHullPM[rptsLarge]);//AbsoluteTiming
(cnv2 = convexHullLL[rptsLarge]);//AbsoluteTiming

  • $\begingroup$ Thanks for your reply, Leonid. I'll go through it carefully, along with your linked list discussion, when my mind is fresh. You're right that my question is more about idiomaticity (if that's even a word) rather than raw performance. $\endgroup$
    – Aky
    Commented May 25, 2013 at 14:12
  • $\begingroup$ @Aky Sure, was glad to help. Good if it works for you. $\endgroup$ Commented May 25, 2013 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.