Source code is available at: GitHub and binaries (executable) here: Download OuelletOnlineConvexHull.zip
Introduction
This article has 2 major intends: Explain the newly developed “Online” (dynamic add point one by one) part of Ouellet Convex Hull algorithm and show some usage and how flexible it is.
This article is the 3^{rd} of a series of 3:
What a Convex Hull is?
A quick definition: When X is a bounded subset of the plane, the convex hull may be visualized as the shape enclosed by a rubber band stretched around X.
A more detailed definition could be found in my first article: A Convex Hull Algorithm and its implementation in O(n log h) or directly in Wikipedia at: Convex Hull and/or Convex Hull Algorithms.
Greeen nails is the Convex Hull or all nails
Image from: Brilliant.org
Finding the Convex hull is an old problem that started in the 60’. Many algorithms to find it appeared since then. Ouellet Convex Hull is a pretty recent Convex Hull algorithm. Its first publication occurs in 2007 as a new algorithm for computing the Convex Hull of a planar point set by Liu and Chen in Journal of Zhejiang UniversitySCIENCE A. Liu and Chen article is pretty hidden and very few articles refer to it, at least in 2014. Liu and Chen article does not contain any code or any details about its implementation. No comparison test was available to prove its real efficiency and nobody seems to have discovered its real power.
Although not the first, I created it totally by myself without having any idea of Liu and Chen paper existence. I discovered all by myself the same algorithm and created an implementation of it. Then I tested it against Chan algorithm because it appeared as being the best one at that moment. In 2014, both had very close performance, either if Ouellet was made in “C#” and Chan in “C” language. “C” has a reputation to be faster than “C#”. In order to know its exact performance advantage, I decided to write my algorithm in “C”. The final Ouellet algorithm in “C” was around 4x faster than Chan in general cases (Circle and Throw away random generator). But it had a major drawback discovered recently. The issue was a performance problem in a specific case where “h” is close to “n”. Although this case should happen rarely, it makes the algorithm not stable in all cases. The problem came from the collection container that was used to keep Convex Hull points. The container was a list which is using an array as its inner container. In Ouellet algorithm, Hull points should always stay sorted. To insert into an array and stay sorted, points need to be moved on one side of the insertion point, usually to the right. The algorithm works fine when Hull points (named “h”) stay very low, either when the source points (named “n”) count is extremely high. But when “h” become too big, around 100 000 points, the performance decrease quickly. In order to fix that issue, the "List" has been replaced by an AVL tree. At the same time, many optimizations have also been added. An article had been written to explain those new implementations and optimizations (see Fast and improved 2D Convex Hull algorithm and its implementation in O(n log h) ).
Today’s article is to explain the latest main addition in the algorithm: the possibility to add a point one by one dynamically. Being able to process one point at a time refer to “online” as per Wikipedia / Preparata and Shamos. The algorithm can, not only process one point at a time, it can do it in O(log h) per point. It could also be initialized from a batch call (like most standard Convex Hull usage) prior to add points one by one.
“h”
 When speaking of performance/complexity, “h” refer to the count of result count which is the Convex Hull points.

“n”
 When speaking of performance/complexity, “n” refer to the count of input points.

Batch
 It is the conventional way to call a Convex Hull algorithm giving a bunch of source points to an algorithm and do one pass to calculate the Convex Hull. It is used in this article often in opposition to “Online”.

Convex Hull
 See: Wikipedia – Convex Hull or Brilliant – Convex Hull or University of Glasgow – pdf document for excellent definitions

General use
(General cases)
 “General use” is often used in this article. It means that in most tested cases, “h” stay extremely low. In code, 2 random generators are mainly used as representing real life cases : Circle and Throw Away. As reference, in those 2 cases, when “n” = 10 000 000 points, “h” stay usually under 1000 points.

Online
 “Online” stand for dynamically add a point one by one, having a valid Convex Hull at any iteration. It is opposed to “Batch” which is applied on a collection of Point and calculate the Convex Hull in one single shot. Wikipedia  "Convex Hull algorithms" has a good definition of “Online” Convex hull.

What the "Online" Ouellet Convex Hull is
Verfy fast Convex Hull Algorithm that enable insertion of points by batch or one by one.
 OR 
Ouellet Convex Hull Online = Ouellet Convex Hull Avl v2 + Online feature = Ouellet Convex Hull Avl v3
Ouellet Convex Hull “Online” algorithm is “Ouellet Convex Hull” algorithm based on the “AVL v2” version with added online functionalities and a little bit more. It is called “Ouellet Convex Hull AVL v3”.
Ouellet Convex Hull algorithm (basic, first version, not AVL, not online) is pretty well explained in:
Ouellet Convex Hull AVL and AVL v2 is pretty well explained in:
Ouellet Convex Hull AVL v3 is based on Ouellet Convex Hull AVL v2. Differences, come from the addition of online functionality and added helper functions in order to simplify querying results in some circumstances. Those added functionalities could help getting results faster in some scenarios.
Added functions are listed in the section Online Convex Hull source code  Added functions.
Some little adjustments have been made in order to add those functions but the impact on performance is insignificant.
In order to understand, you need to know the inner working of Ouellet Convex Hull. See what the Ouellet Convex Hull “Online” is.
The ‘Online’ function has 3 steps executed sequentially:
 It checks if the “Init” has been done and if not, it only does it, otherwise it continues. “Init” initialize each quadrant with the first point and return. “Init” is done only once per Convex Hull calculation on the first call (either Online or Batch).
 It verify if the new point becomes a new boundary (Min/Max X or Min/Max Y) of existing Convex hull points. If so, it update the boundary or each affected quadrant, it then invalidates appropriate neighbors points, finally it return.
 If the new point does not define a new boundary, it does a regular add point as if it was calling in a standard way. The way we add the point is as the quadrant were not disjoint because it is the most general case – which mean we are losing a potential little performance in case where quadrants are disjoint. But evaluating if quadrant is disjoint would be more time consuming than simply call the function for not disjoint quadrants.
 O(1) per point when the point is a Hull Point, otherwise O(0). It can be considered as ~O(h).
There are few advantages about using “Standard/Batch” mode over “Online” mode where few optimizations are only possible in that mode due to the disjoint quadrant check which cost too much to be done in “online” mode. Disjoint quadrant check has those 2 advantages (when apply):
 An optimization enable a quicker verification of quadrant belonging.
 An optimization enables to skip more quadrant belonging verification.
Pros
 The algorithm has a complexity of O(n log h), which is currently the bestknown performance possible.
 According to our test, it is actually the fastest Convex Hull algorithm.
 Based on Ouellet Convex Hull AVL v2 which has been pretty well optimized and tested.
 The same algorithm could be used for both usages: “Standard/Batch” and “Online”. The same algorithm enable the use or standard/Batch call followed by any “Online” call. It could also be used exclusively as an “Online” algorithm from the beginning.
 Has similar performance as standard/Batch Ouellet Convex Hull, although a little bit slower.
 The input of the standard/batch functionalities can be any IEnumerable of Point. That means that it is a lot less restrictive than any algorithm which requires an array as input. Most of them would be able to be fed by an IList in C#. But IList is still more restrictive than IEnumerable.
 Another positive side effect on the low input requirement of IEnumerable is that you can use any object, or a combination of many objects. If they support IEnumerable, it would be enough. The iteratoradapter helper class would provide a single IEnumerable interface that iterate over many different IEnumerable objects.
 Resulting Convex Hull is not in place. That means that the original point order stay the same. In some circumstance it could be a requirement to keep original sources point ordered. For inplace algorithm, when original source order has to stay the same, it would be necessary to make a full copy of source points before calling the algorithm.
 The size required to keep Convex Hull points is in order of “h”. It could be a little higher in general usage, due to possible neighbor invalidation of newly added Convex Hull point, but it should stay very close to “h”. In worse case, it could be close to “n” either if h is extremely low. The worse case could occur when “h” ~= “n” but the 4 last points to add to the Convex Hull are Hull point, one for each quadrant that would invalidate all previous ones. That should statically almost never happen in general case.
 The algorithm is divided in 3. The second part is the one who takes the most time. But it is also the part which is the easiest to multithreaded.
 The algorithm 3^{rd} part is the merge of individual Quadrant result into an array. That part could be avoided depending on scenario requirements.
 Space required stays in O(h). Although not inplace, the footprint should stay very low in general use.
 It appears to be the most reliable algorithm/implementation to use in any circumstance. Best performance in general usage and excellent performance either when h ~= n. Graphics results below show that stability.
 It is not “In Place”. The space used to store Convex Hull points is outside of the source points which mean that it requires additional space. Space required is always in O(h).
Actually, there is 3 known Convex Hull Algorithm which is in O(n log h): Chan, the ultimate planar convex hull algorithm and Ouellet Convex Hull.
The ultimate Convex Hull is not included and has not been tested. No implementation has been found. But there is a very high probability to be quite slow because it requires the calculation of the median. Chan algorithm implementation has been taken as a reference in performance tests. Performance results come from the code supplied at GitHub. Both, Chan and Ouellet algorithms should be the 2 algorithms with best performance. That’s why only both are shown in most performance graphics below.
Graphics shown in this section are obtained from tests with the provided benchmark. For each call to calculate the Convex Hull result, its time taken is represented by a point. For each timing results per algorithm, a line is drawn to represent the linear regression of every result per algorithm. It helps see the trend, if it stays linear or not.
Ouellet algorithm suffix
(Ouellet C# Avl v3…)
 Valid for
 Description

(standard/batch)
 All speed and
validity tests
 Regular/batch algorithm usage. Input: Array, Output: Array

(one by one in a loop)
 All speed and
validity tests
 Simulate a regular/batch algorithm usage.
It receives an array and iterate over all points and passes them to the algorithm like it was used as a regular/batch algorithm. The final Convex Hull result is queried only once when all points have been inserted. It demonstrates how fast an online algorithm could be. Actually, Ouellet Online is almost as fast as the regular/batch part of Ouellet algorithm.

(** Only for Online performance test)
 Only for: Speed Test ‘Online’
 This choice should only be selected only in the “Online” speed test. Its purpose is to see the advantage of being able to add point dynamically vs algorithm than are not “Online”.

Type
 Description

Speed test
 Standard/batch mode. Call the algorithm passing an array of points and get an array of Convex Hull points as result.

Speed test ‘online’
 Call the algorithm one point at a time for “online” algorithm (currently only supported by Ouellet Avl v3).
For all other algorithm, the previous Convex Hull result is resized (array enlarged by one) and add the new point to be evaluated at the end of this enlarged copy. Because for each iteration the array grow by one, that mean a new full copy of the array has to be created. Call the algorithm again with the newly enlarged array. Loop until every point is processed.
NOTE: This test has a special parameter which applies only to Ouellet Avl v3. It could save some time by preventing to ask for a full copy of the result each time (which is required on regular/batch algorithm). Sometimes, context could only require to know if the point had been added or not as a Convex Hull point. Sometimes, only knowing the neighbors of the added point can suffice.

Figure 1
Overview of general performance of many Convex Hull implementations found on the web and included in the benchmark. All O(n log h) algorithms are faster than MI Convex Hull (Delaunay/Voronoi). They are all grouped and sits below MI Convex Hull (Delaunay/Voronoi).
It shows performance advantages of O(n log h) algorithms vs others.
Figure 2
Figure 3
Only algorithms that are in O(n log h) are chosen in order to compare fastest algorithms.
Those figures show performances of Convex Hull algorithm implementations when called in the standard batch way. Both tests show the speed of algorithms for general usage. The difference between the 2 pictures are the type of random generator used (shown in the graphic title).
Figure 2 (Circle generator) shows results for points randomly generated in a circle.
Figure 3 (Throw away) show results for a generator that start from a point and generate the next one in random direction at random distance. That usually creates many distinct clusters of points.
Adding a point one by one has similar performance than doing the same thing in batch mode.
The main difference that exists between the 2 pictures is caused by the random number generator. The distribution of points for “Throw away” random generator has many neighbors closer together while "Circle" has points totally random. Having throw away generator distribution does make one performance optimization pretty efficient. This optimization only exists in batch mode. That is why “Ouellet C# Avl v3 standard/batch” (add all points in batch) has a bigger difference than “Ouellet C# Avl v3 (one by one in a loop)” with "Throw away" generator.
“Ouellet C# Avl v2”, is slightly slower than “Ouellet C# Avl v3” (either when both are in batch mode), although the second one is based on the first one for which they should be extremely close. That is due to 2 things:
 The fact that I put back an optimization in “v3” version that I commented out a long time ago in “v2”. I preferred to let the optimization commented out in “v2” because it makes comparison easier. The optimization is “Quick Discard” (see here for more details).
 The second has an opposite effect, it slow “v3”. It is a little change done on the Online version (“v3”) which slow it down, but so little that it is pretty hard to see. The change is to obtain the result of adding one point: “added” or “not added” as part of the Convex Hull. That modification implied code that adds a very light delay for the “Online” version. To prevent this new delay, I could have duplicate those 2 functions (one for standard/batch, which would have stayed the same, and another for online) that does the exact same thing but because the difference in time is so minor, I kept only one function in order to keep the code smaller and easier to maintain.
Figure 4
Figure 5
Figures 4 and 5 show the performance test for “Speed test Online” usage. The test is the same as regular “Speed test” but it adds points one at a time. For online algorithm, it adds points one by one normally as the interface is done for that. For all other algorithms, which are not “online” implementation, the resulting Convex Hull previously calculated is enlarged by one, the new point to evaluate is stored as the last point in the new enlarged array. The algorithm is then invoked with the new enlarged array containing the additional point. Then we loop until all points are added.
In picture 5, the 2 slower implementation (Multithreaded) has been removed to show better resolution of the 4 others.
As shown on figure 4, it could be observed that multithreaded (mainly "MT" but also for "4T") version have a pretty long steady offset. That offset is due to the fact that they start threads and wait for them to complete before continuing but that step takes quite a long time when compared to the very short time required to insert one point. When a thread synchronization is implied (waits thread before continuing), it implies a thread switch of the main thread waiting for the other to complete. That thread has to wait for the OS to reschedule CPU to continue its work.
We can see that Chan implementation speed difference with any Ouellet used in Batch mode (all but the blue one) is quite smaller. The difference should probably come from the fact that it is done in “C” where no instantiation of objects is required.
The Ouellet C# Avl v3 used in "Online" mode (blue one) has a clear advantage because its instance and its states always stay alive. There is no need to recreate its object and their dependent objects. Also, due to Ouellet algorithm design, which is based on quadrants, the online version benefits of a lot of information which enable very quick evaluation and insertion of Convex Hull point.
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Those figures show the difference when using different approaches to get results of the online Ouellet Convex Hull Algorithm. Usage is described in the title of each figure. Be careful with Figure 9 where the count of point is 10x less due to very low response time. Figures 9 and 10 use the “Arc” random generator where all points are Convex Hull points.
When “online” context is required, there are various ways to get Convex Hull results from Ouellet Convex Hull Avl v3. Choosing the best strategy for your context could make the difference between superfast and slow, between milliseconds and hours. Some rule of thumbs:
 Do not ask any Convex Hull information if you only want to know if the newly evaluated point is part of the Convex Hull or not.
 Only ask for neighbors for newly added Convex Hull point if it is sufficient in your context.
 Ask for the full Convex Hull point array only when a new evaluated point has been added.
When the context permits it, it is far quicker to query neighbors of the newly added point instead of the full Convex Hull result. Looking at the Y axis of Figures 610 shows clearly that not asking a full copy of the points forming the Convex Hull save a lot of time. Choosing the proper function, with less impact on performance, could lead to a good amount of time saved, mainly when the Convex Hull is formed by lots of points.
Figures 9 and 10 purpose is to show the impact on the function used to query Convex Hull information which does not appear between figures 7 and 8. In fact, it is surprising to see how asking for a full copy of points of the Convex Hull require less time than asking for neighbors (Figure 7 vs Figure 8). That weird result should be due to a very low count of Convex Hull result. But when the number of points forming the Convex Hull is pretty large, the advantage of only querying neighbors become very clear (Figures 9 and 10). Although querying for neighbors of a newly added Convex Hull point is not always sufficient, it seems highly recommended to do so when the context permits it.
Figure 11
Figure 12
Figure 13
All figures are based on the “Arc” point generator where the Convex Hull result has “h” = “n”. That means that all Convex Hull result is made of all source points. The figure 11 took days to do.
In figure 11, Ouellet v3 is used in 2 difference ways:
 as a standard/batch algorithm
 as an online algorithm
Figures 12 and 13 only included Ouellet Convex Hull Online to better show its own result in more details.
Figure 13 is the same test as in figure 12 but instead of querying the full results, only the neighbors of the added point are queried. Figure 12 is based on 100 000 points (too slow) while Figure 13 is based on 1 000 000 points.
Having “h” = “n” largely amplify the advantage of online algorithm. But the main idea was to show that difference and according to results, the goal is reached.
Although there is a very low risk to happen in real life, all those tests easily show the main advantage of having an ‘Online’ Convex Hull algorithm when to count of Convex Hull points become very large.
Figure 12 is likely to become in O(n^{2}) because for “n” we query the full results which are “n” in this test. In the case of 1000000 points or more, days would have been required to get all results.
Figure 13 stays pretty much linear and very fast. In fact it stays in O(n log n) because we only get the information about the neighbors of the added point.
“Ouellet Convex Hull AVL v3” has a pretty constant behavior in most circumstances. When “h” = “n”, it is a lot more efficient than any other algorithm. It offers the best performance* in all general cases for conventional/batch use and also when “Online” usage is required. According to current tests, Ouellet Convex Hull AVL v3 offer excellent performance and is very stable whichever the circumstance. It is also the most flexible algorithm enabling “online” use, with or without being previously initialized once by a conventional/batch use.
*Here, Ouellet CPP and all Ouellet algorithm with array based containers are not taken into account because of their problem when “h” is close to “n”.
Code summary
The code contains:
 3 performance tests to compare algorithms speed.
 3 validity tests to ensure Convex Hull algorithms return proper results.
 Many different algorithm implementation (Convex Hull and Smallest Enclosing Circle) found on the web.
 Many Ouellet Convex Hull algorithm implementations in order to see evolution and compare implementation details like optimizations, language and more.
The executable (binary) is provided here. A link is also available at the top of this article.
The executable provides many tests either for performance estimation or for testing validity of algorithms. Each test is selfexplained in the user interface.
The source code is provided on GitHub. It includes everything in order to build the executable. It also includes each and every source code of each algorithm implementation offers as possible tests in the GUI. It also includes a little sample showing different way to call the algorithm, also described a little further in the next section.
Using the code
Added features (functions) for Online use
Added features  Description 
EnumConvexHullPoint  Enum type that is the result of calling either "Evaluate" or "TryAddOnePoint" functions. Possible values are:
(NotConvexHullPoint , AlreadyExists, ConvexHullPoint)

bool IsExists(Point pt)
 Verify if a point is already part of the Convex Hull.

EnumConvexHullPoint Evaluate(Point pt)
 Evaluate if a point would be a Convex Hull point if try to add, without adding it.

EnumConvexHullPoint TryAddOnePoint(Point pt)
 Try to add a point as a Convex Hull.

Point GetNextPoint(Point pt)
 Return the next point of a point

Point GetPreviousPoint(Point pt)
 Return the previous point of a point

Tuple<Point, Point> GetNeighbors(Point pt)
 Return previous and next point of a point

public void CallOnlineConvexHullDifferentWays()
{
GeneratePoints(); // Points are generated in "Point[] _points";
ConvexHull convexHull = new ConvexHull();
// First way to call: Standard way to call (standard/batch)
convexHull.CalcConvexHull(_points);
// Usage of IEnumerable wrapper to loop over each convex hull points
// that are node in avl tree of each quadrant.
foreach (Point pt in convexHull) { Debug.Print(pt.ToString()); }
//Or can get an array copy of the convex hull points
Point[] pointsStandardCall = convexHull.GetResultsAsArrayOfPoint();
convexHull = new ConvexHull();
// Second way to call: Adding one point at a time (online).
foreach (Point pt in _points)
{
convexHull.TryAddOnePoint(pt);
}
Point[] pointsOnlineCall = convexHull.GetResultsAsArrayOfPoint();
DifferencesInPath diffs = ConvexHullUtil.GetPathDifferences(nameof(ConvexHull), _points, pointsStandardCall, pointsOnlineCall);
Debug.Assert(diffs.HasErrors == false);
convexHull = new ConvexHull();
Point[] allPoints = new Point[_points.Length * 2];
Array.Copy(_points, allPoints, _points.Length);
// Third way to call: Standard/Batch then Online
convexHull.CalcConvexHull(_points);
GeneratePoints(_points.Length); // Points are generated in "Point[] _points";
Array.Copy(_points, 0, allPoints, _points.Length, _points.Length);
foreach (Point pt in _points)
{
if (convexHull.TryAddOnePoint(pt) == EnumConvexHullPoint.ConvexHullPoint)
{
// Do nothing if only knowing if Hull point or not is enough
// Get Only neighbors if enough
var neighbors = convexHull.GetNeighbors(pt);
// Query full result as an array
var convexHullPoints = convexHull.GetResultsAsArrayOfPoint();
}
}
// HERE: Verifying previous result
OuelletConvexHull.ConvexHull convexHullOnline2 = new OuelletConvexHull.ConvexHull(allPoints); // Original First Convex Hull (highly tested)
convexHullOnline2.CalcConvexHull();
diffs = ConvexHullUtil.GetPathDifferences(nameof(ConvexHull), _points,
convexHull.GetResultsAsArrayOfPoint(), convexHullOnline2.GetResultsAsArrayOfPoint());
Debug.Assert(diffs.HasErrors == false);
}
Requirement to add full dynamic feature (managing add/remove) in Ouellet Convex Hull :
Requirement
 Reason

Fast search (derived requirement: Keep a sorted copy of source points)
 Removing an existing Convex Hull point could make any existing source points a possible new Convex Hull point. Any points, as many as all of source points could become a new Convex Hull point when only one actual Convex Hull point is removed. In order to not have to look for each and every source point, we need to keep them sorted.

Fast Insert/Remove
 Adding or removing a point requires to update the sorted copy of points. In order to save time, a collection with fast access could improve performance.

* Almost all points could be invalidated if points are in only one quadrant where a point invalidates all points, but the first and last one.
To ensure 3 functionalities:
 To keep source point sorted
 Have fast search
 Support fast insert/remove
A tree appears to be the best compromise as a collection to manage points. Because of a well balance between insert/remove and read, it appears to give an advantage to RedBlack tree over AVLTree.
Also, in order to keep fast access to point tree node between source points and Convex Hull points, each node could include a reference to its counterpart or being the same node if ever possible. That should help keeping a better performance.
Considering every aspect discuss here, it could be inferred that a dynamic Convex Hull algorithm could be in O(log n + log h) = O(log n) per point (add or remove). By using cross reference between both collections (all points and hull points), we can probably remove the "log h" from the equation. Although not really important in the final big O, it could help achieve better performance. I suspect that it is the best possible reachable performance to find Dynamic Convex Hull. Sounds like it could be better than what it is written in Wikipedia:
According to “Preparata, Shamos, Computational Geometry, Chapter "Convex Hulls: Basic Algorithms”: “The online version may be handled with O(log n) per point, which is asymptotically optimal”. The dynamic version may be handled with O(log^{2} n) per operation.”
According to our tests, Ouellet Convex Hull algorithm does it in O(Log h) per point in online which is better than O(log n). Also, it is highly probable that a future dynamic version of Ouellet Convex Hull could achieve O(Log n + Log h) per point using the same principle as the one used and adding another tree to keep source points sorted. That would also be better than O(log^{2} n).
 Why not using RedBlack tree for keeping Convex Hull point. I chose AVL tree over RedBlack tree because there is a lot less point added as Convex Hull point than read operation in regular operation. For 1000000 source points, the count of points in the final Convex Hull result is less than 1000 on every test. That means that most points are rejected while a very small part or ever added as Convex Hull point. That is based on specific random generator but I’m pretty confident that it should reflect most general usage.
 Why not having done the fully dynamic Convex Hull? The question is partly answered in the Future section. But the mains reason is due to complexity of dynamic vs online. While “Online” is compatible with initial batch Convex Hull, the full dynamic Convex Hull is a lot harder and takes lots more time to code/test. It should be really fun to do but should take a lot of time.
 Why result from getting the Convex Hull array (GetResultsAsArrayOfPoint) is not exactly the same as the iterator. The iterator only return Convex Hull points without closing the path. GetResultsAsArrayOfPoint can return the path with or without a closing point. GetResultsAsArrayOfPoint start usually from the second point of Q1 and the only reason is for performance. On the iterator side, it starts at the first element of Q1 because it was also faster this way. In short it was for performance and/or keep code easier to read.
 Why there is no direct access to either the trees or tree nodes. It is because it could lead to wrong utilization and it is highly errorprone because a good undertsanding of the algorithm is required to not make any mistake.
This Article History
 20180228, Creation of this article
 20180304, Upated added features for Online use to reflect recent code changes. Update code to fix warning, fix C++ compilation/reference. Add MSVCRT dlls in code to simplify usage and makes code independant. Fixed, added and corrected some text to improve comprehension.
 20180514, Fixed information about the first appearance of the algorithm article. Added some clarifications,