Triangular Gregory patch

I found an implementation of the Triangular Gregory Patch here on CGSE. Their question is about the normal vectors, which I don't need at the moment because I only need the geometry I was thinking of trying the code on the CPU to test the algorithm, but it seems that the control points (p400, p310, p220...,p004) in the tessellation control shader are calculated incorrectly.

Vector3D p0 = tescOut[gl_InvocationID].p0 = gl_in[gl_InvocationID].gl_Position.xyz();
Vector3D n = tescOut[gl_InvocationID].n = vertOut[gl_InvocationID].normal;

const int nextInvID = gl_InvocationID < 2 ? gl_InvocationID + 1 : 0;
Vector3D edge = gl_in[nextInvID].gl_Position.xyz() - p0;
Vector3D nNext = vertOut[nextInvID].normal;
float d = edge.length();
float a = n.dot(nNext);
Vector3D gama = edge / d;
float a0 = n.dot(gama);
float a1 = nNext.dot(gama);
float ro = 6.0 * (2.0 * a0 + a * a1) / (4.0 - a * a);
float sigma = 6.0 * (2.0 * a1 + a * a0) / (4.0 - a * a);

Vector3D v[4];
v[0] = p0;
v[1] = p0 + (gama * 6.0 - n * ro * 2.0 + nNext * sigma) * (d / 18.0);
v[2] = gl_in[nextInvID].gl_Position.xyz() - (gama * 6.0 + n * ro - nNext * 2.0 * sigma) * (d / 18.0);
v[3] = gl_in[nextInvID].gl_Position.xyz();

edge = gl_in[nextInvID].gl_Position.xyz();

Vector3D w[3];
w[0] = v[1] - v[0];
w[1] = v[2] - v[1];
w[2] = v[3] - v[2];

Vector3D A[3];
A[0] = n.cross(w[0].normalized());
A[2] = nNext.cross(w[2].normalized());

A[1] = (A[0] + A[2]).normalized();

Vector3D l[5];
l[0] = v[0];
l[1] = (v[0] + v[1] * 3.0) * 0.25;
l[2] = (v[1] * 2.0 + v[2] * 2.0) * 0.25;
l[3] = (v[2] * 3.0 + v[3]) * 0.25;
l[4] = v[3];

Vector3D p1 = l[1];
Vector3D p2 = l[2];
Vector3D p3 = l[3];

tescOut[gl_InvocationID].p1 = l[1];
tescOut[gl_InvocationID].p2 = l[2];
tescOut[gl_InvocationID].p3 = l[3];


The input parameters are: gl_in[gl_InvocationID].gl_Position in there is the vertex position in model space. vertOut[nextInvID].normal in there is the vertex normal in model space (the normal is already normalized).

The output: points (p400... p004) can be found in tescOut[]. The coding is as follows:

p400  = tescOut[0].p0
p310  = tescOut[0].p1
p220  = tescOut[0].p2
p130  = tescOut[0].p3
p040  = tescOut[1].p0
p031  = tescOut[1].p1
p022  = tescOut[1].p2
p013  = tescOut[1].p3
p004  = tescOut[2].p0
p103  = tescOut[2].p1
p202  = tescOut[2].p2
p301  = tescOut[2].p3


Applying this algorithm an a cube with shared normals:

Figure 1: Left: The input cube with 8 vertices. The 12 triangles share these vertices. Hence the normal vectors pointing away from the origin of the cube. Right: The marked points (p400,...,p004) of all 12 triangles.

As this is very confusing, I have only applied the algorithm to one of the triangles:

Figure 2: A triangle from the cube with its points (p400,...,p004).

Initially, it is not visible in this one image: The points that belong to the same edge all lie in a plane that is spanned by the vertex normals of the two corners. This is correct.

But: To me, the scale of the points looks wrong. They seem to be too far away from each other.

does anyone have any ideas on how to fix this?