Transform recursion into loop

Question

The C# code below computes the smallest enclosing circle of a set of points but causes stack overflow exceptions due to it's recursive implementation.

It was extracted from a paper which specified the algorithm as recursive.

The algorithm is not tail recursive and it is not trivially convertible to a loop structure.

The recursive function is

public static Circle MiniDiskImpl(IReadOnlyList<Point2D> points, List<Point2D> boundary)
{
    if (!points.Any() || boundary.Count == 3)
    {
        if (boundary.Count == 0)
            return null;

        if (boundary.Count == 1)
            return null;

        if (boundary.Count == 2)
        {
            var radius = boundary[0].DistanceTo(boundary[1])/2;
            var center = (boundary[0] + boundary[1])*0.5;
            return new Circle(Plane.XY, center, radius);
        }

        return new Circle(Plane.XY,boundary[0],boundary[1],boundary[2]);
    }

    var p = points[0];
    var Q = points.GetRange(1,points.Count - 1);
    var D = MiniDiskImpl(Q, boundary);

    if (D==null || D.Center.DistanceTo(p) >= D.Radius)
    {
        D = MiniDiskImpl(Q, boundary.Concat(p).ToList());
    }

    return D;
}

The user calls this function.

public static Circle MiniDisk(List<Point2D> points)
{
    points = points.Slice(0); // Clone the list so we can manipulate in place
    points.Shuffle(); // sorted points cause poor algorithm performance
    return MiniDiskImpl(points, new List<Point2D>());
}

The question is, what is the transformation required to make the algorithm non recursive?

Answer 1

When you use recursion, you make use of the call stack to perform something, you can always transform a recursive algorithm into it's iterative version by dropping your own objects to an actual stack, see DFS example on wiki: http://en.wikipedia.org/wiki/Depth-first_search

Answer 2

I just amazed myself with the following.

/// <summary>
/// Find the smallest enclosing circle. See reference
/// http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf
/// "Smallest enclosing disks (balls and ellipsoids) EMO WELZL
/// </summary>
/// <returns></returns>
public static class Circles 
{

    private static Return<Circle> Left
        (List<Point2D> points, List<Point2D> boundary)
    {

        if ( !points.Any() || boundary.Count == 3 )
        {
            if ( boundary.Count <= 1 )
            {
                return Return.Value<Circle>(null);
            }

            if ( boundary.Count == 2 )
            {
                var radius = boundary[0].DistanceTo(boundary[1]) / 2;
                var center = ( boundary[0] + boundary[1] ) * 0.5;
                return Return.Value(new Circle(Plane.XY, center, radius));
            }

            return Return.Value(new Circle(Plane.XY, boundary[0], boundary[1], boundary[2]));
        }

        var p = points[0];
        var q = points.GetRange(1, points.Count - 1);

        return from circle in Return.Call(() => Left(q, boundary))
               from result in 
                   (circle == null || circle.Center.DistanceTo(p) >= circle.Radius)
                           ? Return.Call(() => Left(q, boundary.Concat(p).ToList()))
                           : Return.Value(circle)
               select result;
    }

    public static Circle MiniDisk( List<Point2D> points )
    {
        points = points.Slice(0);
        points.Shuffle();
        return TrampolineRecursion.Do(()=>LeftLeft(points, new List<Point2D>()));
    }

}

using the following implementation of the Trampoline monad I just magicked up.

namespace Weingartner.Numerics.Recursion
{
    public class Return<T>
    {
        #region return value mode
        public bool HasReturnValue;
        public T ReturnValue;
        #endregion

        #region call and continue mode
        public Func<Return<T>> Call;
        public Func<T, Return<T>> Continue;
        #endregion

        /// <summary>
        /// To avoid GC pressure we use a thread static variable
        /// to handle the return value from functions.
        /// </summary>
        [ThreadStatic] private static Return<T> _ReturnRegister;

        public static Return<T> R()
        {
            if (_ReturnRegister == null)
            {
                _ReturnRegister = new Return<T>();
            }
            return _ReturnRegister;
        }

        public Return<T> Bind(Func<T, Return<T>> fn)
        {
            if (HasReturnValue)
                return Return.Call(() => Return.Value(ReturnValue), fn);

            if(Continue!=null)
                return Return.Call(Call, t => Continue(t).SelectMany(fn));

            return Return.Call(Call, fn);
        }

        public Return<T> SelectMany(Func<T, Return<T>> fn)
        {
            return Bind(fn);
        }
        public Return<T> SelectMany( Func<T, Return<T>> f, Func<T, T, T> g)
        {
            return Bind(x => f(x).Bind(y => Return.Value(g(x, y))));
        }

        public Return<T> Select(Func<T, T> fn)
        {
            if (HasReturnValue)
                return Return.Value(fn(ReturnValue));

            if(Continue!=null)
                return Return.Call(Call, t => Continue(t).Select(fn));

            return Return.Call(Call, t => Continue(t).Select(fn));
        }

        public T Run()
        {

            var stack = new Stack<Func<T,Return<T>>>();
            if(Continue!=null)
                stack.Push(Continue);
            stack.Push(r => Call());

            var @return = Return.Value(default(T));

            while ( stack.Count > 0 )
            {
                @return = stack.Peek()(@return.ReturnValue);
                stack.Pop();
                if (!@return.HasReturnValue)
                {
                    if(@return.Continue!=null)
                        stack.Push(@return.Continue);
                    var t = @return;
                    stack.Push(r => t.Call());
                }
            }

            return ReturnValue = @return.ReturnValue;
        }
    }

    public class Return
    {
        public static Return<T> Value<T>( T t )
        {
            var r = Return<T>.R();
            r.HasReturnValue = true;
            r.ReturnValue = t;
            return r;
        }

        public static Return<T> Call<T>
            (Func<Return<T>> call, Func<T, Return<T>> @continue = null)
        {
            var r = Return<T>.R();
            r.HasReturnValue = false;
            r.Call = call;
            r.Continue = @continue;
            return r;
        }

    }

    public static class TrampolineRecursion
    {
        public static T Do<T>(Func<Return<T>> call)
        {
            return Return.Call(call).Run();
        }
    }
}