Shape *s1,*s2;

void epa()
{
	static std::vector<Vec2D> s;
		
	s.clear();

        //a,b,c are points from gjk simplex
	s.push_back(a);
	s.push_back(b);
	s.push_back(c);

	for (int i = 0; i < 32; i++)
	{
		Vec2D norm = Vec2D::Zero();
		float dist = 0;
		int index = -1;

		epaClosestEdge(s, index, norm, dist);

		Vec2D p = support(norm);

		float d = p * norm;

		if (d - dist < 1e-6)
		{
			collisionNormal = norm;
			penetration = d;
			return;
		}
		else
		{
			s.insert(s.begin() + index, p);
		}
	}
}

void epaClosestEdge(std::vector<Vec2D>& s, int& edgeIndex, Vec2D& normal, float& dist)
{
	dist = INFINITY;

	for (int i = 0; i < s.size(); i++)
	{
		int j = i + 1 == s.size()? 0 : i + 1;

		Vec2D a = s[i];
		Vec2D b = s[j];

		Vec2D e = b - a;

		Vec2D oa = a;

		Vec2D n = tripleProduct(e, oa, e);

		n.normalize();

		float d = n * a;

		if (d < dist)
		{
			dist = d;
			normal = n;
			edgeIndex = j;
		}
	}
}

Vec2D support(Vec2D& n)
{
	Vec2D a = getfar(s1, n);
	Vec2D b = getfar(s2, -n);

	return a-b;
}

Vec2D getfar(Shape* s, Vec2D& n)
{
	float maxDot = -INFINITY;
	Vec2D best;

	for (int i = 0; i < s->vsCount(); i++)
	{
		Vec2D v = s->vs(i);

		float d = v * n;

		if (d > maxDot)
		{
			maxDot = d;
			best = v;
		}
	}

	return best;
}

Plain copy-paste-adaptation to cpp from dyn4j.

1. What about simplex' winding? Should it be CW, CCW or arbitrary given?

2. Any mistakes in my implementation? Any situations when this implementation fails?

3. How can I get contatct points? Why is it said that EPA gives only one contact point? All i get are penetration and normal.