25154

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.

<img alt="enter image description here" class="b-lazy" data-src="https://i.stack.imgur.com/UI2A1.png" data-original="https://i.stack.imgur.com/UI2A1.png" src="https://etrip.eimg.top/images/2019/05/07/timg.gif" />

It was extracted from <a href="http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf" rel="nofollow">a paper</a> which specified the algorithm as recursive.

The algorithm is <strong>not</strong> 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?

Answer1:

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(); } } }

Answer2:

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: <a href="http://en.wikipedia.org/wiki/Depth-first_search" rel="nofollow">http://en.wikipedia.org/wiki/Depth-first_search</a>

Recommend

  • How to fix this floating point square root algorithm
  • Graph matplotlib to show total count in the histogram bins
  • Suspending event listeners
  • ObjectMaterialize in EF not firing on first level query
  • Visual Studio 2010 - 2015 does not use ymm* registers for AVX optimization
  • Update data in d3.js group
  • AngularJS : transclude ng-repeat inside directive
  • Can I customize a Jackson ObjectMapper by adding a module?
  • Why doesn't CGContextSetStrokeColorWithColor set text colour to black?
  • How to bind comma separated list of values to List
  • Stitching 2 images (OpenCV)
  • Simplify where clause with repeated associated type restrictions
  • How to open html table in xls on click of a button
  • Setting WPF Window Background to Resource Dictionary Brush User Setting
  • NHibernate manually control fetching
  • How to estimate the Kalman Filter with 'KFAS' R package, with an AR(1) transition equation
  • How can Delete be both a DDL and a DML statement
  • How can I include If-None-Match header in HttpRequestMessage
  • Calculating ratio of reciprocated ties for each node in igraph
  • How to create CGPath from a SKSpriteNode in SWIFT
  • Generate random number from custom distribution
  • Inline R code in YAML for rmarkdown doesn't run
  • Mysterious problem with floating point in LISP - time axis generation
  • dc-js disable selecting slices on click for pie chart
  • Unity3D & Android: Difference between “UnityMain” and “main” threads?
  • Does CUDA 5 support STL or THRUST inside the device code?
  • PHP: When would you need the self:: keyword?
  • JTable with a ScrollPane misbehaving
  • How to get Windows thread pool to call class member function?
  • Authorize attributes not working in MVC 4
  • Add sale price programmatically to product variations
  • unknown Exception android
  • Django query for large number of relationships
  • Busy indicator not showing up in wpf window [duplicate]
  • failed to connect to specific WiFi in android programmatically
  • Unable to use reactive element in my shiny app
  • Python/Django TangoWithDjango Models and Databases
  • Net Present Value in Excel for Grouped Recurring CF
  • How can I use threading to 'tick' a timer to be accessed by other threads?
  • How do I use LINQ to get all the Items that have a particular SubItem?