GJK+EPA infinite loop

I am implementing custom physics in my game and having a go at gjk+epa for collision detection and response. The GJK part seems to work correctly (in 2D for now), but I am running into an infinite loop on my epa code and I do not understand why.

My algorithm looks like this:

     public static Vector2 GJK(List<Vector2> a, List<Vector2> b)
     {
         //GJK
         List<Vector2> s = new List<Vector2>();
         Vector2 dir = Vector2.left;
         s.Add(support(a, dir) - support(b, -dir));
         dir = -s[0];
 
         while (s.Count < 3) //3 because we are in 2D, so an enclosing simplex is a triangle
         {
             //add the vertex furthest along in the search direction onto the simplex
             Vector2 p = support(a, dir) - support(b, -dir);
             //if that vertex is not on the other side of the origin, report no collision
             if (Vector2.Dot(p, dir) < 0) return Vector2.zero;
             s.Add(p);
             //keep only the piece of the simplex that the origin is closest to
             do_simplex(ref s, ref dir);
         }
 
         //EPA
         int     min_i = 0;
         float   min_d = float.PositiveInfinity;
         Vector2 min_n = Vector2.zero;
 
         while (min_d == float.PositiveInfinity)
         {
             //get closest edge to origin
             for (int i = 0; i < s.Count; ++i)
             {
                 Vector2 e = s[(i+1) % s.Count] - s[i];
                 Vector2 n = new Vector2(e.y, -e.x).normalized;
                 float e_d = Vector2.Dot(n, s[i]);
                 if (e_d < 0) { n *= -1; e_d *= -1;}
 
                 if (e_d < min_d)
                 {
                     min_d = e_d;
                     min_n = n;
                     min_i = (i + 1) % s.Count;
                 }
             }
 
             //search for verticies in direction normal to the min edge
             Vector2 p = support(a, min_n) - support(b, -min_n);
             float p_d = Vector2.Dot(p, min_n);
 
             //if the polytope has been expanded, we have more faces to search
             if (Mathf.Abs(p_d - min_d) > 0.001f)
             {
                 min_d = float.PositiveInfinity;
                 //inserting the vertex in the correct order should take care of
                 //the polygon automatically in 2D. (why is this still convex?)
                 s.Insert(min_i, p);
             }
         }
 
         return min_n * (min_d + 0.001f);
     }

The problem is that if I run it on almost any vertex list, e.g. a = b = {(-0.2, 0.6), (0.8, 0.6), (0.8, 1.6), (-0.2, 1.6)}, I get stuck in the second while loop forever, as points are cyclically added into the polytope.

I feel like there is just something I fundamentally don't understand about the algorithm and all of the readings about it online are extremely terse, so I'm hoping someone here who has implemented it before can see the error.