dmikoss · October 30, 2020 07:58 · Aug 8, 2020 · Mar 29, 2020 · Mar 28, 2020 · Jan 20, 2020
diff --git a/gjk.c b/gjk.c
@@ -59,7 +59,7 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
     struct gjk_vertex *vert = &s->v[s->cnt];
     f3cpy(vert->a, sup->a);
     f3cpy(vert->b, sup->b);
-    f3cpy(vert->p, sup->d);
+    f3sub(vert->p, sup->b, sup->a);
     vert->aid = sup->aid;
     vert->bid = sup->bid;
     s->bc[s->cnt++] = 1.0f;

diff --git a/gjk.h b/gjk.h
@@ -7,7 +7,6 @@ struct gjk_support {
     int aid, bid;   /* in */
     float a[3];     /* in */
     float b[3];     /* in */
-    float d[3];     /* in */
     float da[3];    /* out */
     float db[3];    /* out */
 };

diff --git a/xample_simple.c b/xample_simple.c
@@ -32,14 +32,12 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     struct gjk_support s = {0};
     f3cpy(s.a, verts);
     f3cpy(s.b, ca);
-    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
     while (gjk(&gsx, &s)) {
         s.aid = polyhedron_support(s.a, s.da, verts, cnt);
         s.bid = line_support(s.b, s.db, ca, cb);
-        f3sub(s.d, s.b, s.a);
     }
     /* check distance between closest points */
     struct gjk_result res;

diff --git a/xample_transform.c b/xample_transform.c
@@ -38,7 +38,6 @@ polyhedron_intersect_polyhedron(
     struct gjk_support s = {0};
     transform(s.a, averts, arot, apos);
     transform(s.b, bverts, brot, bpos);
-    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
@@ -53,8 +52,7 @@ polyhedron_intersect_polyhedron(
 
         /* calculate distance vector on transformed points */
         transform(s.a, sa, arot, apos);
-        transform(s.b, sb, brot, bpos);       
-        f3sub(s.d, s.b, s.a);
+        transform(s.b, sb, brot, bpos);
     }
     return gsx.hit;
 }
diff --git a/xample_xdebug.c b/xample_xdebug.c
@@ -25,12 +25,10 @@ polyhedron_intersect_sphere_debug(
     struct gjk_support s = {0};
     f3cpy(s.a, verts);
     f3cpy(s.b, sc);
-    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     while (gjk(&gsx[n++], &s)) {
         s.aid = polyhedron_support(s.a, s.da, verts, cnt);
-        f3sub(s.d, s.b, s.a);
         gsx[n] = gsx[n-1];
     }
     return n;

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -38,7 +38,6 @@ polyhedron_intersect_polyhedron(
     struct gjk_support s = {0};
     transform(s.a, averts, arot, apos);
     transform(s.b, bverts, brot, bpos);
-    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
@@ -58,7 +57,6 @@ polyhedron_intersect_polyhedron(
         /* calculate distance vector */
         if (f3dot(s.da, ta) > 0) f3add(s.a, s.a, ta);
         if (f3dot(s.db, tb) > 0) f3add(s.b, s.b, tb);
-        f3sub(s.d, s.b, s.a);
     }
     return gsx.hit;
 }
diff --git a/gjk.c b/gjk.c
@@ -345,38 +345,17 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
         float v_abcd = f3box(c, a, d) * volume;
         float w_abcd = f3box(d, a, b) * volume;
         float x_abcd = f3box(b, a, c) * volume;
-        {
-            /* check minkowski distance to origin */
-            float pa[3]; f3mul(a, s->v[0].p, denom * u_abcd);
-            float pb[3]; f3mul(b, s->v[1].p, denom * v_abcd);
-            float pc[3]; f3mul(c, s->v[2].p, denom * w_abcd);
-            float pd[3]; f3mul(d, s->v[3].p, denom * x_abcd);
-
-            float pnt[3];
-            f3add(pnt, pa, pb);
-            f3add(pnt, pnt, pc);
-            f3add(pnt, pnt, pd);
-
-            if (f3dot(pnt,pnt) < GJK_EPSILON * GJK_EPSILON) {
-                s->bc[0] = u_abcd;
-                s->bc[1] = v_abcd;
-                s->bc[2] = w_abcd;
-                s->bc[3] = x_abcd;
-                s->cnt = 4;
-                return 1;
-            }
-        }
 
         /* check faces for closest point */
-        if (x_abcd <= 0.0f && u_abc > 0.0f && v_abc > 0.0f && w_abc > 0.0f) {
+        if (x_abcd < 0.0f && u_abc > 0.0f && v_abc > 0.0f && w_abc > 0.0f) {
             /* region ABC */
             s->bc[0] = u_abc;
             s->bc[1] = v_abc;
             s->bc[2] = w_abc;
             s->cnt = 3;
             break;
         }
-        if (u_abcd <= 0.0f && u_cbd > 0.0f && v_cbd > 0.0f && w_cbd > 0.0f) {
+        if (u_abcd < 0.0f && u_cbd > 0.0f && v_cbd > 0.0f && w_cbd > 0.0f) {
             /* region CBD */
             s->v[0] = s->v[2];
             s->v[2] = s->v[3];
@@ -386,7 +365,7 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
             s->cnt = 3;
             break;
         }
-        if (v_abcd <= 0.0f && u_acd > 0.0f && v_acd > 0.0f && w_acd > 0.0f) {
+        if (v_abcd < 0.0f && u_acd > 0.0f && v_acd > 0.0f && w_acd > 0.0f) {
             /* region ACD */
             s->v[1] = s->v[2];
             s->v[2] = s->v[3];
@@ -396,7 +375,7 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
             s->cnt = 3;
             break;
         }
-        if (w_abcd <= 0.0f && u_adb > 0.0f && v_adb > 0.0f && w_adb > 0.0f) {
+        if (w_abcd < 0.0f && u_adb > 0.0f && v_adb > 0.0f && w_adb > 0.0f) {
             /* region ADB */
             s->v[2] = s->v[1];
             s->v[1] = s->v[3];
@@ -407,7 +386,7 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
             break;
         }
         /* region ABCD */
-        assert(u_abcd > 0.0f && v_abcd > 0.0f && w_abcd > 0.0f && x_abcd > 0.0f);
+        assert(u_abcd >= 0.0f && v_abcd >= 0.0f && w_abcd >= 0.0f && x_abcd >= 0.0f);
         s->bc[0] = u_abcd;
         s->bc[1] = v_abcd;
         s->bc[2] = w_abcd;

diff --git a/gjk.c b/gjk.c
@@ -345,6 +345,27 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
         float v_abcd = f3box(c, a, d) * volume;
         float w_abcd = f3box(d, a, b) * volume;
         float x_abcd = f3box(b, a, c) * volume;
+        {
+            /* check minkowski distance to origin */
+            float pa[3]; f3mul(a, s->v[0].p, denom * u_abcd);
+            float pb[3]; f3mul(b, s->v[1].p, denom * v_abcd);
+            float pc[3]; f3mul(c, s->v[2].p, denom * w_abcd);
+            float pd[3]; f3mul(d, s->v[3].p, denom * x_abcd);
+
+            float pnt[3];
+            f3add(pnt, pa, pb);
+            f3add(pnt, pnt, pc);
+            f3add(pnt, pnt, pd);
+
+            if (f3dot(pnt,pnt) < GJK_EPSILON * GJK_EPSILON) {
+                s->bc[0] = u_abcd;
+                s->bc[1] = v_abcd;
+                s->bc[2] = w_abcd;
+                s->bc[3] = x_abcd;
+                s->cnt = 4;
+                return 1;
+            }
+        }
 
         /* check faces for closest point */
         if (x_abcd <= 0.0f && u_abc > 0.0f && v_abc > 0.0f && w_abc > 0.0f) {
@@ -421,17 +442,6 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
 
         f3add(pnt, a, b);
         f3add(pnt, pnt, c);
-    } break;
-    case 4: {
-        /* ----- Tetrahedron ----- */
-        float a[3]; f3mul(a, s->v[0].p, denom * s->bc[0]);
-        float b[3]; f3mul(b, s->v[1].p, denom * s->bc[1]);
-        float c[3]; f3mul(c, s->v[2].p, denom * s->bc[2]);
-        float d[3]; f3mul(d, s->v[3].p, denom * s->bc[3]);
-
-        f3add(pnt, a, b);
-        f3add(pnt, pnt, c);
-        f3add(pnt, pnt, d);
     } break;}
 
     float d2 = f3dot(pnt, pnt);

diff --git a/_GJK.md b/_GJK.md
@@ -123,7 +123,7 @@ debug drawing.
 References
 --------------
 
-- [Erin Cattos presentation on 2D GJK](http://box2d.org/files/GDC2010/GDC2010_Catto_Erin_GJK.pdf)
+- [Erin Cattos presentation on 2D GJK](https://box2d.org/files/ErinCatto_GJK_GDC2010.pdf)
 - [Randy Gauls 3D GJK implementation]( https://bitbucket.org/rgaul/lm/src/806fe4427db32f98e514c646b51f99ce1cfc5c14/src/lm/collision/query/lmDistance.cpp?at=default&fileviewer=file-view-default)
 - [Casey Muratoris video lecture](https://mollyrocket.com/849)
 - Christer Ericson Real Time Collision Detection

diff --git a/xample_transform.c b/xample_transform.c
@@ -48,9 +48,8 @@ polyhedron_intersect_polyhedron(
         float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed directions  */
-        float sa[3], sb[3];
-        s.aid = polyhedron_support(sa, da, averts, acnt);
-        s.bid = polyhedron_support(sb, db, bverts, bcnt);
+        float sa[3]; s.aid = polyhedron_support(sa, da, averts, acnt);
+        float sb[3]; s.bid = polyhedron_support(sb, db, bverts, bcnt);
 
         /* calculate distance vector on transformed points */
         transform(s.a, sa, arot, apos);

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -48,9 +48,8 @@ polyhedron_intersect_polyhedron(
         float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed into model space directions  */
-        float sa[3], sb[3];
-        s.aid = polyhedron_support(sa, da, averts, acnt);
-        s.bid = polyhedron_support(sb, db, bverts, bcnt);
+        float sa[3]; s.aid = polyhedron_support(sa, da, averts, acnt);
+        float sb[3]; s.bid = polyhedron_support(sb, db, bverts, bcnt);
 
         /* transform points to world space */
         transform(s.a, sa, arot, apos);

diff --git a/xample_transform.c b/xample_transform.c
@@ -48,12 +48,13 @@ polyhedron_intersect_polyhedron(
         float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed directions  */
-        s.aid = polyhedron_support(s.a, da, averts, acnt);
-        s.bid = polyhedron_support(s.b, db, bverts, bcnt);
+        float sa[3], sb[3];
+        s.aid = polyhedron_support(sa, da, averts, acnt);
+        s.bid = polyhedron_support(sb, db, bverts, bcnt);
 
         /* calculate distance vector on transformed points */
-        transform(s.a, s.a, arot, apos);
-        transform(s.b, s.b, brot, bpos);       
+        transform(s.a, sa, arot, apos);
+        transform(s.b, sb, brot, bpos);       
         f3sub(s.d, s.b, s.a);
     }
     return gsx.hit;

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -48,12 +48,13 @@ polyhedron_intersect_polyhedron(
         float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed into model space directions  */
-        s.aid = polyhedron_support(s.a, da, averts, acnt);
-        s.bid = polyhedron_support(s.b, db, bverts, bcnt);
+        float sa[3], sb[3];
+        s.aid = polyhedron_support(sa, da, averts, acnt);
+        s.bid = polyhedron_support(sb, db, bverts, bcnt);
 
         /* transform points to world space */
-        transform(s.a, s.a, arot, apos);
-        transform(s.b, s.b, brot, bpos);  
+        transform(s.a, sa, arot, apos);
+        transform(s.b, sb, brot, bpos);  
 
         /* calculate distance vector */
         if (f3dot(s.da, ta) > 0) f3add(s.a, s.a, ta);

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -56,8 +56,8 @@ polyhedron_intersect_polyhedron(
         transform(s.b, s.b, brot, bpos);  
 
         /* calculate distance vector */
-        if (f3dot(n, ta) > 0) f3add(s.a, s.a, ta);
-        if (f3dot(d, tb) > 0) f3add(s.b, s.b, tb);
+        if (f3dot(s.da, ta) > 0) f3add(s.a, s.a, ta);
+        if (f3dot(s.db, tb) > 0) f3add(s.b, s.b, tb);
         f3sub(s.d, s.b, s.a);
     }
     return gsx.hit;

diff --git a/gjk.c b/gjk.c
@@ -464,7 +464,7 @@ gjk(struct gjk_simplex *s, struct gjk_support *sup)
     }}
     if (f3dot(d,d) < GJK_EPSILON * GJK_EPSILON)
         return 0;
-    
+
     f3mul(sup->da, d, -1.0f);
     f3cpy(sup->db, d);
     return 1;

diff --git a/gjk.h b/gjk.h
@@ -7,7 +7,7 @@ struct gjk_support {
     int aid, bid;   /* in */
     float a[3];     /* in */
     float b[3];     /* in */
-    float d;        /* in */
+    float d[3];     /* in */
     float da[3];    /* out */
     float db[3];    /* out */
 };
@@ -33,4 +33,4 @@ extern int gjk(struct gjk_simplex *s, struct gjk_support *sup);
 extern void gjk_analyze(struct gjk_result *res, const struct gjk_simplex *s);
 extern void gjk_quad(struct gjk_result *res, float a_radius, float b_radius);
 
-#endif
+#endif
diff --git a/gjk.c b/gjk.c
@@ -36,12 +36,11 @@ inv_sqrt(float n)
     return conv.f;
 }
 extern int
-gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
+gjk(struct gjk_simplex *s, struct gjk_support *sup)
 {
     assert(s);
-    assert(dv);
     assert(sup);
-    if (!s || !sup || !dv) return 0;
+    if (!s || !sup) return 0;
     if (s->max_iter > 0 && s->iter >= s->max_iter)
         return 0;
 
@@ -60,7 +59,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     struct gjk_vertex *vert = &s->v[s->cnt];
     f3cpy(vert->a, sup->a);
     f3cpy(vert->b, sup->b);
-    f3cpy(vert->p, dv);
+    f3cpy(vert->p, sup->d);
     vert->aid = sup->aid;
     vert->bid = sup->bid;
     s->bc[s->cnt++] = 1.0f;
@@ -440,30 +439,34 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     s->D = d2;
 
     /* VII.) New search direction */
+    float d[3] = {0};
     switch (s->cnt) {
     default: assert(0); break;
     case 1: {
         /* --------- Point -------- */
-        f3mul(dv, s->v[0].p, -1);
+        f3mul(d, s->v[0].p, -1);
     } break;
     case 2: {
         /* ------ Line segment ---- */
         float ba[3]; f3sub(ba, s->v[1].p, s->v[0].p);
         float b0[3]; f3mul(b0, s->v[1].p, -1);
         float t[3];  f3cross(t, ba, b0);
-        f3cross(dv, t, ba);
+        f3cross(d, t, ba);
     } break;
     case 3: {
         /* ------- Triangle ------- */
         float ab[3]; f3sub(ab, s->v[1].p, s->v[0].p);
         float ac[3]; f3sub(ac, s->v[2].p, s->v[0].p);
         float n[3];  f3cross(n, ab, ac);
         if (f3dot(n, s->v[0].p) <= 0.0f)
-            f3cpy(dv, n);
-        else f3mul(dv, n, -1);
+            f3cpy(d, n);
+        else f3mul(d, n, -1);
     }}
-    if (f3dot(dv,dv) < GJK_EPSILON * GJK_EPSILON)
+    if (f3dot(d,d) < GJK_EPSILON * GJK_EPSILON)
         return 0;
+
+    f3mul(sup->da, d, -1.0f);
+    f3cpy(sup->db, d);
     return 1;
 }
 extern void

diff --git a/gjk.h b/gjk.h
@@ -4,9 +4,12 @@
 #define GJK_MAX_ITERATIONS 20
 
 struct gjk_support {
-    int aid, bid;
-    float a[3];
-    float b[3];
+    int aid, bid;   /* in */
+    float a[3];     /* in */
+    float b[3];     /* in */
+    float d;        /* in */
+    float da[3];    /* out */
+    float db[3];    /* out */
 };
 struct gjk_simplex {
     int max_iter, iter;
@@ -26,7 +29,7 @@ struct gjk_result {
     float distance_squared;
     int iterations;
 };
-extern int gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv);
+extern int gjk(struct gjk_simplex *s, struct gjk_support *sup);
 extern void gjk_analyze(struct gjk_result *res, const struct gjk_simplex *s);
 extern void gjk_quad(struct gjk_result *res, float a_radius, float b_radius);
 

diff --git a/xample_simple.c b/xample_simple.c
@@ -29,19 +29,17 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     const float *ca, const float *cb, float cr)
 {
     /* initial guess */
-    float d[3] = {0};
     struct gjk_support s = {0};
     f3cpy(s.a, verts);
     f3cpy(s.b, ca);
-    f3sub(d, s.b, s.a);
+    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
-    while (gjk(&gsx, &s, d)) {
-        float n[3]; f3mul(n, d, -1);
-        s.aid = polyhedron_support(s.a, n, verts, cnt);
-        s.bid = line_support(s.b, d, ca, cb);
-        f3sub(d, s.b, s.a);
+    while (gjk(&gsx, &s)) {
+        s.aid = polyhedron_support(s.a, s.da, verts, cnt);
+        s.bid = line_support(s.b, s.db, ca, cb);
+        f3sub(s.d, s.b, s.a);
     }
     /* check distance between closest points */
     struct gjk_result res;

diff --git a/xample_transform.c b/xample_transform.c
@@ -35,19 +35,17 @@ polyhedron_intersect_polyhedron(
     const float *bverts, int bcnt, const float *bpos, const float *brot)
 {
     /* initial guess */
-    float d[3] = {0};
     struct gjk_support s = {0};
     transform(s.a, averts, arot, apos);
     transform(s.b, bverts, brot, bpos);
-    f3sub(d, s.b, s.a);
+    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
-    while (gjk(&gsx, &s, d)) {
+    while (gjk(&gsx, &s)) {
         /* transform direction */
-        float n[3]; f3mul(n, d, -1);
-        float da[3]; transformI(da, n, arot, apos);
-        float db[3]; transformI(db, d, brot, bpos);
+        float da[3]; transformI(da, s.da, arot, apos);
+        float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed directions  */
         s.aid = polyhedron_support(s.a, da, averts, acnt);
@@ -56,7 +54,7 @@ polyhedron_intersect_polyhedron(
         /* calculate distance vector on transformed points */
         transform(s.a, s.a, arot, apos);
         transform(s.b, s.b, brot, bpos);       
-        f3sub(d, s.b, s.a);
+        f3sub(s.d, s.b, s.a);
     }
     return gsx.hit;
 }
diff --git a/xample_xdebug.c b/xample_xdebug.c
@@ -22,18 +22,15 @@ polyhedron_intersect_sphere_debug(
 {
     /* initial guess */
     int n = 0;
-    float d[3] = {0};
     struct gjk_support s = {0};
-
     f3cpy(s.a, verts);
     f3cpy(s.b, sc);
-    f3sub(d, s.b, s.a);
+    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
-    while (gjk(&gsx[n++], &s, d)) {
-        float nd[3]; f3mul(nd, d, -1);
-        s.aid = polyhedron_support(s.a, nd, verts, cnt);
-        f3sub(d, s.b, s.a);
+    while (gjk(&gsx[n++], &s)) {
+        s.aid = polyhedron_support(s.a, s.da, verts, cnt);
+        f3sub(s.d, s.b, s.a);
         gsx[n] = gsx[n-1];
     }
     return n;

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -35,19 +35,17 @@ polyhedron_intersect_polyhedron(
     const float *bverts, int bcnt, const float *bpos, const float *brot, const float *tb)
 {
     /* initial guess */
-    float d[3] = {0};
     struct gjk_support s = {0};
     transform(s.a, averts, arot, apos);
     transform(s.b, bverts, brot, bpos);
-    f3sub(d, s.b, s.a);
+    f3sub(s.d, s.b, s.a);
 
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
-    while (gjk(&gsx, &s, d)) {
+    while (gjk(&gsx, &s)) {
         /* transform direction */
-        float n[3]; f3mul(n, d, -1);
-        float da[3]; transformI(da, n, arot, apos);
-        float db[3]; transformI(db, d, brot, bpos);
+        float da[3]; transformI(da, s.da, arot, apos);
+        float db[3]; transformI(db, s.db, brot, bpos);
 
         /* run support function on tranformed into model space directions  */
         s.aid = polyhedron_support(s.a, da, averts, acnt);
@@ -60,7 +58,7 @@ polyhedron_intersect_polyhedron(
         /* calculate distance vector */
         if (f3dot(n, ta) > 0) f3add(s.a, s.a, ta);
         if (f3dot(d, tb) > 0) f3add(s.b, s.b, tb);
-        f3sub(d, s.b, s.a);
+        f3sub(s.d, s.b, s.a);
     }
     return gsx.hit;
 }
diff --git a/xample_transform.c b/xample_transform.c
@@ -3,17 +3,17 @@
 static void
 transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    r3[0] = v3[0] * r33[0*3+0] + v3[0] * r33[0*3+1] + v3[0] * r33[0*3+2] + t3[0];
-    r3[1] = v3[1] * r33[1*3+0] + v3[1] * r33[1*3+1] + v3[1] * r33[1*3+2] + t3[1];
-    r3[2] = v3[2] * r33[2*3+0] + v3[2] * r33[2*3+1] + v3[2] * r33[2*3+2] + t3[2];
+    r3[0] = v3[0] * r33[0*3+0] + v3[1] * r33[1*3+0] + v3[2] * r33[2*3+0] + t3[0];
+    r3[1] = v3[0] * r33[0*3+1] + v3[1] * r33[1*3+1] + v3[2] * r33[2*3+1] + t3[1];
+    r3[2] = v3[0] * r33[0*3+2] + v3[1] * r33[1*3+2] + v3[2] * r33[2*3+2] + t3[2];
 }
 static void
-transformI(float *r, const float *v, const float *r33, const float *t3)
+transformI(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    const float v[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
-    r3[0] = v[0] * r33[0*3+0] + v[0] * r33[1*3+0] + v[0] * r33[2*3+0];
-    r3[1] = v[1] * r33[0*3+1] + v[1] * r33[1*3+1] + v[1] * r33[2*3+1];
-    r3[2] = v[2] * r33[0*3+2] + v[2] * r33[1*3+2] + v[2] * r33[2*3+2];
+    const float p[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
+    r3[0] = p[0] * r33[0*3+0] + p[1] * r33[0*3+1] + p[2] * r33[0*3+2];
+    r3[1] = p[0] * r33[1*3+0] + p[1] * r33[1*3+1] + p[2] * r33[1*3+2];
+    r3[2] = p[0] * r33[2*3+0] + p[1] * r33[2*3+1] + p[2] * r33[2*3+2];
 }
 static int
 polyhedron_support(float *support, const float *d,

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -3,17 +3,17 @@
 static void
 transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    r3[0] = v3[0] * r33[0*3+0] + v3[0] * r33[0*3+1] + v3[0] * r33[0*3+2] + t3[0];
-    r3[1] = v3[1] * r33[1*3+0] + v3[1] * r33[1*3+1] + v3[1] * r33[1*3+2] + t3[1];
-    r3[2] = v3[2] * r33[2*3+0] + v3[2] * r33[2*3+1] + v3[2] * r33[2*3+2] + t3[2];
+    r3[0] = v3[0] * r33[0*3+0] + v3[1] * r33[1*3+0] + v3[2] * r33[2*3+0] + t3[0];
+    r3[1] = v3[0] * r33[0*3+1] + v3[1] * r33[1*3+1] + v3[2] * r33[2*3+1] + t3[1];
+    r3[2] = v3[0] * r33[0*3+2] + v3[1] * r33[1*3+2] + v3[2] * r33[2*3+2] + t3[2];
 }
 static void
-transformI(float *r, const float *v, const float *r33, const float *t3)
+transformI(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    const float v[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
-    r3[0] = v[0] * r33[0*3+0] + v[0] * r33[1*3+0] + v[0] * r33[2*3+0];
-    r3[1] = v[1] * r33[0*3+1] + v[1] * r33[1*3+1] + v[1] * r33[2*3+1];
-    r3[2] = v[2] * r33[0*3+2] + v[2] * r33[1*3+2] + v[2] * r33[2*3+2];
+    const float p[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
+    r3[0] = p[0] * r33[0*3+0] + p[1] * r33[0*3+1] + p[2] * r33[0*3+2];
+    r3[1] = p[0] * r33[1*3+0] + p[1] * r33[1*3+1] + p[2] * r33[1*3+2];
+    r3[2] = p[0] * r33[2*3+0] + p[1] * r33[2*3+1] + p[2] * r33[2*3+2];
 }
 static int
 polyhedron_support(float *support, const float *d,

diff --git a/_GJK.md b/_GJK.md
@@ -86,7 +86,7 @@ Examples
     translation `t` in context of our current direction `d`. Only if both point in the same direction we need to consider
     the translated object (Real time Collision Detection Page.: 409).
 
-Support
+Support Functions
 -----------
 
 - Support functions are the core of the algorithm and main abstraction over convex polytopes

diff --git a/xample_transform.c b/xample_transform.c
@@ -3,24 +3,17 @@
 static void
 transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    for (int i = 0; i < 3; ++i) {
-        const float v = v3[i];
-        const float t = t3[i];
-        const float r0 = v * r33[i*3+0];
-       	const float r1 = v * r33[i*3+1];
-        const float r2 = v * r33[i*3+2];
-        r3[i] = r0 + r1 + r2 + t;
-    }
+    r3[0] = v3[0] * r33[0*3+0] + v3[0] * r33[0*3+1] + v3[0] * r33[0*3+2] + t3[0];
+    r3[1] = v3[1] * r33[1*3+0] + v3[1] * r33[1*3+1] + v3[1] * r33[1*3+2] + t3[1];
+    r3[2] = v3[2] * r33[2*3+0] + v3[2] * r33[2*3+1] + v3[2] * r33[2*3+2] + t3[2];
 }
 static void
-transformT(float *r, const float *v, const float *r33, const float *t3)
+transformI(float *r, const float *v, const float *r33, const float *t3)
 {
-    for (int i = 0; i < 3; ++i) {
-        float p = v[i] - t3[i];
-        r[i]  = p * r33[0*3+i];
-        r[i] += p * r33[1*3+i];
-        r[i] += p * r33[2*3+i];
-    }
+    const float v[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
+    r3[0] = v[0] * r33[0*3+0] + v[0] * r33[1*3+0] + v[0] * r33[2*3+0];
+    r3[1] = v[1] * r33[0*3+1] + v[1] * r33[1*3+1] + v[1] * r33[2*3+1];
+    r3[2] = v[2] * r33[0*3+2] + v[2] * r33[1*3+2] + v[2] * r33[2*3+2];
 }
 static int
 polyhedron_support(float *support, const float *d,
@@ -53,8 +46,8 @@ polyhedron_intersect_polyhedron(
     while (gjk(&gsx, &s, d)) {
         /* transform direction */
         float n[3]; f3mul(n, d, -1);
-        float da[3]; transformT(da, n, arot, apos);
-        float db[3]; transformT(db, d, brot, bpos);
+        float da[3]; transformI(da, n, arot, apos);
+        float db[3]; transformI(db, d, brot, bpos);
 
         /* run support function on tranformed directions  */
         s.aid = polyhedron_support(s.a, da, averts, acnt);

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -3,24 +3,17 @@
 static void
 transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
-    for (int i = 0; i < 3; ++i) {
-        const float v = v3[i];
-        const float t = t3[i];
-        const float r0 = v * r33[i*3+0];
-       	const float r1 = v * r33[i*3+1];
-        const float r2 = v * r33[i*3+2];
-        r3[i] = r0 + r1 + r2 + t;
-    }
+    r3[0] = v3[0] * r33[0*3+0] + v3[0] * r33[0*3+1] + v3[0] * r33[0*3+2] + t3[0];
+    r3[1] = v3[1] * r33[1*3+0] + v3[1] * r33[1*3+1] + v3[1] * r33[1*3+2] + t3[1];
+    r3[2] = v3[2] * r33[2*3+0] + v3[2] * r33[2*3+1] + v3[2] * r33[2*3+2] + t3[2];
 }
 static void
-transformT(float *r, const float *v, const float *r33, const float *t3)
+transformI(float *r, const float *v, const float *r33, const float *t3)
 {
-    for (int i = 0; i < 3; ++i) {
-        float p = v[i] - t3[i];
-        r[i]  = p * r33[0*3+i];
-        r[i] += p * r33[1*3+i];
-        r[i] += p * r33[2*3+i];
-    }
+    const float v[3] = {v3[0] - t3[0], v3[1] - t3[1], v3[2] - t3[2] };
+    r3[0] = v[0] * r33[0*3+0] + v[0] * r33[1*3+0] + v[0] * r33[2*3+0];
+    r3[1] = v[1] * r33[0*3+1] + v[1] * r33[1*3+1] + v[1] * r33[2*3+1];
+    r3[2] = v[2] * r33[0*3+2] + v[2] * r33[1*3+2] + v[2] * r33[2*3+2];
 }
 static int
 polyhedron_support(float *support, const float *d,
@@ -53,8 +46,8 @@ polyhedron_intersect_polyhedron(
     while (gjk(&gsx, &s, d)) {
         /* transform direction */
         float n[3]; f3mul(n, d, -1);
-        float da[3]; transformT(da, n, arot, apos);
-        float db[3]; transformT(db, d, brot, bpos);
+        float da[3]; transformI(da, n, arot, apos);
+        float db[3]; transformI(db, d, brot, bpos);
 
         /* run support function on tranformed into model space directions  */
         s.aid = polyhedron_support(s.a, da, averts, acnt);

diff --git a/xample_transform.c b/xample_transform.c
@@ -1,22 +1,18 @@
 #include "gjk.h"
 
 static void
-transform(float *r, const float *v, const float *r33, const float *t3)
+transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
     for (int i = 0; i < 3; ++i) {
-        r[i]  = v[i] * r33[i*3+0];
-        r[i] += v[i] * r33[i*3+1];
-        r[i] += v[i] * r33[i*3+2];
-        r[i] += t3[i];
+        const float v = v3[i];
+        const float t = t3[i];
+        const float r0 = v * r33[i*3+0];
+       	const float r1 = v * r33[i*3+1];
+        const float r2 = v * r33[i*3+2];
+        r3[i] = r0 + r1 + r2 + t;
     }
 }
 static void
-transformS(float *v, const float *r33, const float *t3)
-{
-    float tmp[3]; f3cpy(tmp, v);
-    transform(v, tmp, r33, t3);
-}
-static void
 transformT(float *r, const float *v, const float *r33, const float *t3)
 {
     for (int i = 0; i < 3; ++i) {
@@ -48,11 +44,8 @@ polyhedron_intersect_polyhedron(
     /* initial guess */
     float d[3] = {0};
     struct gjk_support s = {0};
-    f3cpy(s.a, averts);
-    f3cpy(s.b, bverts);
-
-    transformS(s.a, arot, apos);
-    transformS(s.b, brot, bpos);
+    transform(s.a, averts, arot, apos);
+    transform(s.b, bverts, brot, bpos);
     f3sub(d, s.b, s.a);
 
     /* run gjk algorithm */
@@ -68,8 +61,8 @@ polyhedron_intersect_polyhedron(
         s.bid = polyhedron_support(s.b, db, bverts, bcnt);
 
         /* calculate distance vector on transformed points */
-        transformS(s.a, arot, apos);
-        transformS(s.b, brot, bpos);       
+        transform(s.a, s.a, arot, apos);
+        transform(s.b, s.b, brot, bpos);       
         f3sub(d, s.b, s.a);
     }
     return gsx.hit;

diff --git a/xample_xmov.c b/xample_xmov.c
@@ -1,22 +1,18 @@
 #include "gjk.h"
 
 static void
-transform(float *r, const float *v, const float *r33, const float *t3)
+transform(float *r3, const float *v3, const float *r33, const float *t3)
 {
     for (int i = 0; i < 3; ++i) {
-        r[i]  = v[i] * r33[i*3+0];
-        r[i] += v[i] * r33[i*3+1];
-        r[i] += v[i] * r33[i*3+2];
-        r[i] += t3[i];
+        const float v = v3[i];
+        const float t = t3[i];
+        const float r0 = v * r33[i*3+0];
+       	const float r1 = v * r33[i*3+1];
+        const float r2 = v * r33[i*3+2];
+        r3[i] = r0 + r1 + r2 + t;
     }
 }
 static void
-transformS(float *v, const float *r33, const float *t3)
-{
-    float tmp[3]; f3cpy(tmp, v);
-    transform(v, tmp, r33, t3);
-}
-static void
 transformT(float *r, const float *v, const float *r33, const float *t3)
 {
     for (int i = 0; i < 3; ++i) {
@@ -48,11 +44,8 @@ polyhedron_intersect_polyhedron(
     /* initial guess */
     float d[3] = {0};
     struct gjk_support s = {0};
-    f3cpy(s.a, averts);
-    f3cpy(s.b, bverts);
-
-    transformS(s.a, arot, apos);
-    transformS(s.b, brot, bpos);
+    transform(s.a, averts, arot, apos);
+    transform(s.b, bverts, brot, bpos);
     f3sub(d, s.b, s.a);
 
     /* run gjk algorithm */
@@ -68,8 +61,8 @@ polyhedron_intersect_polyhedron(
         s.bid = polyhedron_support(s.b, db, bverts, bcnt);
 
         /* transform points to world space */
-        transformS(s.a, arot, apos);
-        transformS(s.b, brot, bpos);  
+        transform(s.a, s.a, arot, apos);
+        transform(s.b, s.b, brot, bpos);  
 
         /* calculate distance vector */
         if (f3dot(n, ta) > 0) f3add(s.a, s.a, ta);

diff --git a/xample_simple.c b/xample_simple.c
@@ -29,16 +29,15 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     const float *ca, const float *cb, float cr)
 {
     /* initial guess */
+    float d[3] = {0};
     struct gjk_support s = {0};
     f3cpy(s.a, verts);
     f3cpy(s.b, ca);
-
+    f3sub(d, s.b, s.a);
+
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
-    f3sub(gsx.d, s.b, s.a);
-
-    float d[3] = {0};
-    while (gjk(&d, &gsx, &s)) {
+    while (gjk(&gsx, &s, d)) {
         float n[3]; f3mul(n, d, -1);
         s.aid = polyhedron_support(s.a, n, verts, cnt);
         s.bid = line_support(s.b, d, ca, cb);
@@ -49,4 +48,4 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     gjk_analyze(&res, &gsx);
     gjk_quad(&res, 0, cr);
     return res.hit;
-}
+}
diff --git a/xample_simple.c b/xample_simple.c
@@ -29,15 +29,16 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     const float *ca, const float *cb, float cr)
 {
     /* initial guess */
-    float d[3] = {0};
     struct gjk_support s = {0};
     f3cpy(s.a, verts);
     f3cpy(s.b, ca);
-    f3sub(d, s.b, s.a);
-
+
     /* run gjk algorithm */
     struct gjk_simplex gsx = {0};
-    while (gjk(&gsx, &s, d)) {
+    f3sub(gsx.d, s.b, s.a);
+
+    float d[3] = {0};
+    while (gjk(&d, &gsx, &s)) {
         float n[3]; f3mul(n, d, -1);
         s.aid = polyhedron_support(s.a, n, verts, cnt);
         s.bid = line_support(s.b, d, ca, cb);

diff --git a/gjk.c b/gjk.c
@@ -6,23 +6,20 @@
 
 static const float f3z[3];
 #define fop(r,e,a,p,b,i,s) (r) e ((a) p (b)) i (s)
-#define f3op(r,e,a,p,b,i,s)\
-    do{fop((r)[0],e,(a)[0],p,(b)[0],i,s),\
+#define f3op(r,e,a,p,b,i,s) do {\
+    fop((r)[0],e,(a)[0],p,(b)[0],i,s),\
     fop((r)[1],e,(a)[1],p,(b)[1],i,s),\
     fop((r)[2],e,(a)[2],p,(b)[2],i,s);}while(0)
 #define f3cpy(d,s) (d)[0]=(s)[0],(d)[1]=(s)[1],(d)[2]=(s)[2]
 #define f3add(d,a,b) f3op(d,=,a,+,b,+,0)
 #define f3sub(d,a,b) f3op(d,=,a,-,b,+,0)
 #define f3mul(d,a,s) f3op(d,=,a,+,f3z,*,s)
 #define f3dot(a,b) ((a)[0]*(b)[0]+(a)[1]*(b)[1]+(a)[2]*(b)[2])
+#define f3cross(d,a,b) do {\
+    (d)[0] = ((a)[1]*(b)[2]) - ((a)[2]*(b)[1]),\
+    (d)[1] = ((a)[2]*(b)[0]) - ((a)[0]*(b)[2]),\
+    (d)[2] = ((a)[0]*(b)[1]) - ((a)[1]*(b)[0]);}while(0)
 
-static inline void
-f3cross(float *r, const float *a, const float *b)
-{
-    r[0] = (a[1]*b[2]) - (a[2]*b[1]);
-    r[1] = (a[2]*b[0]) - (a[0]*b[2]);
-    r[2] = (a[0]*b[1]) - (a[1]*b[0]);
-}
 static inline float
 f3box(const float *a, const float *b, const float *c)
 {

diff --git a/gjk.c b/gjk.c
@@ -26,7 +26,8 @@ f3cross(float *r, const float *a, const float *b)
 static inline float
 f3box(const float *a, const float *b, const float *c)
 {
-    float n[3]; f3cross(n, a, b);
+    float n[3];
+    f3cross(n, a, b);
     return f3dot(n, c);
 }
 static float

diff --git a/xample_simple.c b/xample_simple.c
@@ -46,6 +46,6 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     /* check distance between closest points */
     struct gjk_result res;
     gjk_analyze(&res, &gsx);
-    gjk_quad(res, 0, cr);
+    gjk_quad(&res, 0, cr);
     return res.hit;
 }
diff --git a/xample_xdebug.c b/xample_xdebug.c
@@ -49,7 +49,7 @@ debug_draw_polyhedron_intersect_sphere(int key_frame,
 
     struct gjk_result res;
     gjk_analyze(&gsx[key_frame]);
-    gjk_quad(res, 0, sr);
+    gjk_quad(&res, 0, sr);
 
     glBox(res.p0, 0.05f, 0.05f, 0.05f);
     glBox(res.p1, 0.05f, 0.05f, 0.05f);

diff --git a/xample_simple.c b/xample_simple.c
@@ -44,7 +44,8 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
         f3sub(d, s.b, s.a);
     }
     /* check distance between closest points */
-    struct gjk_result res = gjk_analyze(&gsx);
+    struct gjk_result res;
+    gjk_analyze(&res, &gsx);
     gjk_quad(res, 0, cr);
     return res.hit;
 }
diff --git a/xample_xdebug.c b/xample_xdebug.c
@@ -47,8 +47,10 @@ debug_draw_polyhedron_intersect_sphere(int key_frame,
     int n = polyhedron_intersect_sphere_debug(gsx, verts, cnt, sc, sr);
     key_frame = clamp(0, key_frame, n-1);
 
-    struct gjk_result res = gjk_analyze(&gsx[key_frame]);
+    struct gjk_result res;
+    gjk_analyze(&gsx[key_frame]);
     gjk_quad(res, 0, sr);
+
     glBox(res.p0, 0.05f, 0.05f, 0.05f);
     glBox(res.p1, 0.05f, 0.05f, 0.05f);
     glLine(res.p0, res.p1);

diff --git a/gjk.c b/gjk.c
@@ -7,9 +7,9 @@
 static const float f3z[3];
 #define fop(r,e,a,p,b,i,s) (r) e ((a) p (b)) i (s)
 #define f3op(r,e,a,p,b,i,s)\
-    {fop((r)[0],e,(a)[0],p,(b)[0],i,s),\
+    do{fop((r)[0],e,(a)[0],p,(b)[0],i,s),\
     fop((r)[1],e,(a)[1],p,(b)[1],i,s),\
-    fop((r)[2],e,(a)[2],p,(b)[2],i,s);}
+    fop((r)[2],e,(a)[2],p,(b)[2],i,s);}while(0)
 #define f3cpy(d,s) (d)[0]=(s)[0],(d)[1]=(s)[1],(d)[2]=(s)[2]
 #define f3add(d,a,b) f3op(d,=,a,+,b,+,0)
 #define f3sub(d,a,b) f3op(d,=,a,-,b,+,0)
@@ -26,8 +26,7 @@ f3cross(float *r, const float *a, const float *b)
 static inline float
 f3box(const float *a, const float *b, const float *c)
 {
-    float n[3];
-    f3cross(n, a, b);
+    float n[3]; f3cross(n, a, b);
     return f3dot(n, c);
 }
 static float
@@ -49,7 +48,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         return 0;
 
     /* I.) Initialize */
-    if (!s->cnt) {
+    if (s->cnt == 0) {
         s->D = GJK_FLT_MAX;
         s->max_iter = !s->max_iter ? GJK_MAX_ITERATIONS: s->max_iter;
     }
@@ -73,14 +72,12 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     case 1: break;
     case 2: {
         /* -------------------- Line ----------------------- */
-        float a[3], b[3];
-        f3cpy(a, s->v[0].p);
-        f3cpy(b, s->v[1].p);
+        float a[3]; f3cpy(a, s->v[0].p);
+        float b[3]; f3cpy(b, s->v[1].p);
 
         /* compute barycentric coordinates */
-        float ab[3], ba[3];
-        f3sub(ab, a, b);
-        f3sub(ba, b, a);
+        float ab[3]; f3sub(ab, a, b);
+        float ba[3]; f3sub(ba, b, a);
 
         float u = f3dot(b, ba);
         float v = f3dot(a, ab);
@@ -104,18 +101,16 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     } break;
     case 3: {
         /* -------------------- Triangle ----------------------- */
-        float a[3], b[3], c[3];
-        f3cpy(a, s->v[0].p);
-        f3cpy(b, s->v[1].p);
-        f3cpy(c, s->v[2].p);
-
-        float ab[3], ba[3], bc[3], cb[3], ca[3], ac[3];
-        f3sub(ab, a, b);
-        f3sub(ba, b, a);
-        f3sub(bc, b, c);
-        f3sub(cb, c, b);
-        f3sub(ca, c, a);
-        f3sub(ac, a, c);
+        float a[3]; f3cpy(a, s->v[0].p);
+        float b[3]; f3cpy(b, s->v[1].p);
+        float c[3]; f3cpy(c, s->v[2].p);
+
+        float ab[3]; f3sub(ab, a, b);
+        float ba[3]; f3sub(ba, b, a);
+        float bc[3]; f3sub(bc, b, c);
+        float cb[3]; f3sub(cb, c, b);
+        float ca[3]; f3sub(ca, c, a);
+        float ac[3]; f3sub(ac, a, c);
 
         /* compute barycentric coordinates */
         float u_ab = f3dot(b, ba);
@@ -148,11 +143,10 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             break;
         }
         /* calculate fractional area */
-        float n[3], n1[3], n2[3], n3[3];
-        f3cross(n, ba, ca);
-        f3cross(n1, b, c);
-        f3cross(n2, c, a);
-        f3cross(n3, a, b);
+        float n[3]; f3cross(n, ba, ca);
+        float n1[3]; f3cross(n1, b, c);
+        float n2[3]; f3cross(n2, c, a);
+        float n3[3]; f3cross(n3, a, b);
 
         float u_abc = f3dot(n1, n);
         float v_abc = f3dot(n2, n);
@@ -192,28 +186,24 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     } break;
     case 4: {
         /* -------------------- Tetrahedron ----------------------- */
-        float a[3], b[3], c[3], d[3];
-        f3cpy(a, s->v[0].p);
-        f3cpy(b, s->v[1].p);
-        f3cpy(c, s->v[2].p);
-        f3cpy(d, s->v[3].p);
-
-        float ab[3], ba[3], bc[3], cb[3], ca[3], ac[3];
-        float db[3], bd[3], dc[3], cd[3], ad[3], da[3];
-
-        f3sub(ab, a, b);
-        f3sub(ba, b, a);
-        f3sub(bc, b, c);
-        f3sub(cb, c, b);
-        f3sub(ca, c, a);
-        f3sub(ac, a, c);
-
-        f3sub(db, d, b);
-        f3sub(bd, b, d);
-        f3sub(dc, d, c);
-        f3sub(cd, c, d);
-        f3sub(da, d, a);
-        f3sub(ad, a, d);
+        float a[3]; f3cpy(a, s->v[0].p);
+        float b[3]; f3cpy(b, s->v[1].p);
+        float c[3]; f3cpy(c, s->v[2].p);
+        float d[3]; f3cpy(d, s->v[3].p);
+
+        float ab[3]; f3sub(ab, a, b);
+        float ba[3]; f3sub(ba, b, a);
+        float bc[3]; f3sub(bc, b, c);
+        float cb[3]; f3sub(cb, c, b);
+        float ca[3]; f3sub(ca, c, a);
+        float ac[3]; f3sub(ac, a, c);
+
+        float db[3]; f3sub(db, d, b);
+        float bd[3]; f3sub(bd, b, d);
+        float dc[3]; f3sub(dc, d, c);
+        float cd[3]; f3sub(cd, c, d);
+        float da[3]; f3sub(da, d, a);
+        float ad[3]; f3sub(ad, a, d);
 
         /* compute barycentric coordinates */
         float u_ab = f3dot(b, ba);
@@ -263,11 +253,10 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             break;
         }
         /* calculate fractional area */
-        float n[3], n1[3], n2[3], n3[3];
-        f3cross(n, da, ba);
-        f3cross(n1, d, b);
-        f3cross(n2, b, a);
-        f3cross(n3, a, d);
+        float n[3]; f3cross(n, da, ba);
+        float n1[3]; f3cross(n1, d, b);
+        float n2[3]; f3cross(n2, b, a);
+        float n3[3]; f3cross(n3, a, d);
 
         float u_adb = f3dot(n1, n);
         float v_adb = f3dot(n2, n);
@@ -414,7 +403,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         return 0;
     }
     /* VI.) Ensure closing in on origin to prevent multi-step cycling */
-    {float pnt[3], denom = 0;
+    float pnt[3], denom = 0;
     for (int i = 0; i < s->cnt; ++i)
         denom += s->bc[i];
     denom = 1.0f / denom;
@@ -423,28 +412,25 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     case 1: f3cpy(pnt, s->v[0].p); break;
     case 2: {
         /* --------- Line -------- */
-        float a[3], b[3];
-        f3mul(a, s->v[0].p, denom * s->bc[0]);
-        f3mul(b, s->v[1].p, denom * s->bc[1]);
+        float a[3]; f3mul(a, s->v[0].p, denom * s->bc[0]);
+        float b[3]; f3mul(b, s->v[1].p, denom * s->bc[1]);
         f3add(pnt, a, b);
     } break;
     case 3: {
         /* ------- Triangle ------ */
-        float a[3], b[3], c[3];
-        f3mul(a, s->v[0].p, denom * s->bc[0]);
-        f3mul(b, s->v[1].p, denom * s->bc[1]);
-        f3mul(c, s->v[2].p, denom * s->bc[2]);
+        float a[3]; f3mul(a, s->v[0].p, denom * s->bc[0]);
+        float b[3]; f3mul(b, s->v[1].p, denom * s->bc[1]);
+        float c[3]; f3mul(c, s->v[2].p, denom * s->bc[2]);
 
         f3add(pnt, a, b);
         f3add(pnt, pnt, c);
     } break;
     case 4: {
         /* ----- Tetrahedron ----- */
-        float a[3], b[3], c[3], d[3];
-        f3mul(a, s->v[0].p, denom * s->bc[0]);
-        f3mul(b, s->v[1].p, denom * s->bc[1]);
-        f3mul(c, s->v[2].p, denom * s->bc[2]);
-        f3mul(d, s->v[3].p, denom * s->bc[3]);
+        float a[3]; f3mul(a, s->v[0].p, denom * s->bc[0]);
+        float b[3]; f3mul(b, s->v[1].p, denom * s->bc[1]);
+        float c[3]; f3mul(c, s->v[2].p, denom * s->bc[2]);
+        float d[3]; f3mul(d, s->v[3].p, denom * s->bc[3]);
 
         f3add(pnt, a, b);
         f3add(pnt, pnt, c);
@@ -453,7 +439,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
 
     float d2 = f3dot(pnt, pnt);
     if (d2 >= s->D) return 0;
-    s->D = d2;}
+    s->D = d2;
 
     /* VII.) New search direction */
     switch (s->cnt) {
@@ -464,18 +450,16 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     } break;
     case 2: {
         /* ------ Line segment ---- */
-        float b0[3], ba[3], t[3];
-        f3sub(ba, s->v[1].p, s->v[0].p);
-        f3mul(b0, s->v[1].p, -1);
-        f3cross(t, ba, b0);
+        float ba[3]; f3sub(ba, s->v[1].p, s->v[0].p);
+        float b0[3]; f3mul(b0, s->v[1].p, -1);
+        float t[3];  f3cross(t, ba, b0);
         f3cross(dv, t, ba);
     } break;
     case 3: {
         /* ------- Triangle ------- */
-        float n[3], ab[3], ac[3];
-        f3sub(ab, s->v[1].p, s->v[0].p);
-        f3sub(ac, s->v[2].p, s->v[0].p);
-        f3cross(n, ab, ac);
+        float ab[3]; f3sub(ab, s->v[1].p, s->v[0].p);
+        float ac[3]; f3sub(ac, s->v[2].p, s->v[0].p);
+        float n[3];  f3cross(n, ab, ac);
         if (f3dot(n, s->v[0].p) <= 0.0f)
             f3cpy(dv, n);
         else f3mul(dv, n, -1);
@@ -484,12 +468,11 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         return 0;
     return 1;
 }
-extern struct gjk_result
-gjk_analyze(const struct gjk_simplex *s)
+extern void
+gjk_analyze(struct gjk_result *res, const struct gjk_simplex *s)
 {
-    struct gjk_result r = {0};
-    r.iterations = s->iter;
-    r.hit = s->hit;
+    res->iterations = s->iter;
+    res->hit = s->hit;
 
     /* calculate normalization denominator */
     float denom = 0;
@@ -502,92 +485,88 @@ gjk_analyze(const struct gjk_simplex *s)
     default: assert(0); break;
     case 1: {
         /* Point */
-        f3cpy(r.p0, s->v[0].a);
-        f3cpy(r.p1, s->v[0].b);
+        f3cpy(res->p0, s->v[0].a);
+        f3cpy(res->p1, s->v[0].b);
     } break;
     case 2: {
         /* Line */
         float as = denom * s->bc[0];
         float bs = denom * s->bc[1];
 
-        float a[3], b[3], c[3], d[3];
-        f3mul(a, s->v[0].a, as);
-        f3mul(b, s->v[1].a, bs);
-        f3mul(c, s->v[0].b, as);
-        f3mul(d, s->v[1].b, bs);
+        float a[3]; f3mul(a, s->v[0].a, as);
+        float b[3]; f3mul(b, s->v[1].a, bs);
+        float c[3]; f3mul(c, s->v[0].b, as);
+        float d[3]; f3mul(d, s->v[1].b, bs);
 
-        f3add(r.p0, a, b);
-        f3add(r.p1, c, d);
+        f3add(res->p0, a, b);
+        f3add(res->p1, c, d);
     } break;
     case 3: {
         /* Triangle */
         float as = denom * s->bc[0];
         float bs = denom * s->bc[1];
         float cs = denom * s->bc[2];
 
-        float a[3], b[3], c[3];
-        f3mul(a, s->v[0].a, as);
-        f3mul(b, s->v[1].a, bs);
-        f3mul(c, s->v[2].a, cs);
+        float a[3]; f3mul(a, s->v[0].a, as);
+        float b[3]; f3mul(b, s->v[1].a, bs);
+        float c[3]; f3mul(c, s->v[2].a, cs);
 
-        float d[3], e[3], f[3];
-        f3mul(d, s->v[0].b, as);
-        f3mul(e, s->v[1].b, bs);
-        f3mul(f, s->v[2].b, cs);
+        float d[3]; f3mul(d, s->v[0].b, as);
+        float e[3]; f3mul(e, s->v[1].b, bs);
+        float f[3]; f3mul(f, s->v[2].b, cs);
 
-        f3add(r.p0, a, b);
-        f3add(r.p0, r.p0, c);
+        f3add(res->p0, a, b);
+        f3add(res->p0, res->p0, c);
 
-        f3add(r.p1, d, e);
-        f3add(r.p1, r.p1, f);
+        f3add(res->p1, d, e);
+        f3add(res->p1, res->p1, f);
     } break;
     case 4: {
         /* Tetrahedron */
-        float a[3], b[3], c[3], d[3];
-        f3mul(a, s->v[0].a, denom * s->bc[0]);
-        f3mul(b, s->v[1].a, denom * s->bc[1]);
-        f3mul(c, s->v[2].a, denom * s->bc[2]);
-        f3mul(d, s->v[3].a, denom * s->bc[3]);
-
-        f3add(r.p0, a, b);
-        f3add(r.p0, r.p0, c);
-        f3add(r.p0, r.p0, d);
-        f3cpy(r.p1, r.p0);
+        float a[3]; f3mul(a, s->v[0].a, denom * s->bc[0]);
+        float b[3]; f3mul(b, s->v[1].a, denom * s->bc[1]);
+        float c[3]; f3mul(c, s->v[2].a, denom * s->bc[2]);
+        float d[3]; f3mul(d, s->v[3].a, denom * s->bc[3]);
+
+        f3add(res->p0, a, b);
+        f3add(res->p0, res->p0, c);
+        f3add(res->p0, res->p0, d);
+        f3cpy(res->p1, res->p0);
     } break;}
 
-    if (!r.hit) {
+    if (!res->hit) {
         /* compute distance */
-        float d[3]; f3sub(d, r.p1, r.p0);
-        r.distance_squared = f3dot(d, d);
-    } else r.distance_squared = 0;
-    return r;
+        float d[3]; f3sub(d, res->p1, res->p0);
+        res->distance_squared = f3dot(d, d);
+    } else res->distance_squared = 0;
 }
 extern void
-gjk_quad(struct gjk_result *r, float ar, float br)
+gjk_quad(struct gjk_result *res, float a_radius, float b_radius)
 {
-    float rad = ar + br;
-    if (r->distance_squared > GJK_EPSILON &&
-        r->distance_squared > rad*rad) {
-        r->distance_squared -= rad*rad;
+    float radius = a_radius + b_radius;
+    float radius_squared = radius * radius;
+    if (res->distance_squared > GJK_EPSILON &&
+        res->distance_squared > radius_squared) {
+        res->distance_squared -= radius_squared;
 
         /* calculate normal */
-        float n[3]; f3sub(n, r->p1, r->p0);
+        float n[3]; f3sub(n, res->p1, res->p0);
         float l2 = f3dot(n, n);
         if (l2 != 0.0f) {
             float il = inv_sqrt(l2);
             f3mul(n,n,il);
         }
-        float da[3]; f3mul(da, n, ar);
-        float db[3]; f3mul(db, n, br);
+        float da[3]; f3mul(da, n, a_radius);
+        float db[3]; f3mul(db, n, b_radius);
 
         /* calculate new collision points */
-        f3add(r->p0, r->p0, da);
-        f3sub(r->p1, r->p1, db);
+        f3add(res->p0, res->p0, da);
+        f3sub(res->p1, res->p1, db);
     } else {
-        float p[3]; f3add(p, r->p0, r->p1);
-        f3mul(r->p0, p, 0.5f);
-        f3cpy(r->p1, r->p0);
-        r->distance_squared = 0;
-        r->hit = 1;
+        float p[3]; f3add(p, res->p0, res->p1);
+        f3mul(res->p0, p, 0.5f);
+        f3cpy(res->p1, res->p0);
+        res->distance_squared = 0;
+        res->hit = 1;
     }
 }
diff --git a/gjk.h b/gjk.h
@@ -27,7 +27,7 @@ struct gjk_result {
     int iterations;
 };
 extern int gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv);
-extern struct gjk_result gjk_analyze(const struct gjk_simplex *s);
-extern void gjk_quad(struct gjk_result *r, float ar, float br);
+extern void gjk_analyze(struct gjk_result *res, const struct gjk_simplex *s);
+extern void gjk_quad(struct gjk_result *res, float a_radius, float b_radius);
 
 #endif
diff --git a/gjk.c b/gjk.c
@@ -33,7 +33,7 @@ f3box(const float *a, const float *b, const float *c)
 static float
 inv_sqrt(float n)
 {
-    union {unsigned u; float f;} conv = {.f = n};
+    union {unsigned u; float f;} conv; conv.f = n;
     conv.u = 0x5f375A84 - (conv.u >> 1);
     conv.f = conv.f * (1.5f - (n * 0.5f * conv.f * conv.f));
     return conv.f;

diff --git a/gjk.c b/gjk.c
@@ -31,15 +31,11 @@ f3box(const float *a, const float *b, const float *c)
     return f3dot(n, c);
 }
 static float
-inv_sqrt(float i)
+inv_sqrt(float n)
 {
-    float x2;
-    const float threehalfs = 1.5f;
-    union {unsigned i; float f;} conv = {0};
-    conv.f = i;
-    x2 = i * 0.5f;
-    conv.i = 0x5f375A84 - (conv.i >> 1);
-    conv.f = conv.f * (threehalfs - (x2 * conv.f * conv.f));
+    union {unsigned u; float f;} conv = {.f = n};
+    conv.u = 0x5f375A84 - (conv.u >> 1);
+    conv.f = conv.f * (1.5f - (n * 0.5f * conv.f * conv.f));
     return conv.f;
 }
 extern int

diff --git a/gjk.c b/gjk.c
@@ -42,14 +42,6 @@ inv_sqrt(float i)
     conv.f = conv.f * (threehalfs - (x2 * conv.f * conv.f));
     return conv.f;
 }
-static inline void
-f3unit(float *r, const float *v)
-{
-    float il, l2 = f3dot(v, v);
-    if (l2 == 0.0f) return;
-    il = inv_sqrt(l2);
-    f3mul(r,v,il);
-}
 extern int
 gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
 {
@@ -583,7 +575,12 @@ gjk_quad(struct gjk_result *r, float ar, float br)
         r->distance_squared -= rad*rad;
 
         /* calculate normal */
-        float n[3]; f3sub(n, r->p1, r->p0); f3unit(n,n);
+        float n[3]; f3sub(n, r->p1, r->p0);
+        float l2 = f3dot(n, n);
+        if (l2 != 0.0f) {
+            float il = inv_sqrt(l2);
+            f3mul(n,n,il);
+        }
         float da[3]; f3mul(da, n, ar);
         float db[3]; f3mul(db, n, br);
 

diff --git a/gjk.c b/gjk.c
@@ -30,6 +30,26 @@ f3box(const float *a, const float *b, const float *c)
     f3cross(n, a, b);
     return f3dot(n, c);
 }
+static float
+inv_sqrt(float i)
+{
+    float x2;
+    const float threehalfs = 1.5f;
+    union {unsigned i; float f;} conv = {0};
+    conv.f = i;
+    x2 = i * 0.5f;
+    conv.i = 0x5f375A84 - (conv.i >> 1);
+    conv.f = conv.f * (threehalfs - (x2 * conv.f * conv.f));
+    return conv.f;
+}
+static inline void
+f3unit(float *r, const float *v)
+{
+    float il, l2 = f3dot(v, v);
+    if (l2 == 0.0f) return;
+    il = inv_sqrt(l2);
+    f3mul(r,v,il);
+}
 extern int
 gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
 {
@@ -563,7 +583,7 @@ gjk_quad(struct gjk_result *r, float ar, float br)
         r->distance_squared -= rad*rad;
 
         /* calculate normal */
-        float n[3]; f3sub(n, r->p1, r->p0); f3norm(n);
+        float n[3]; f3sub(n, r->p1, r->p0); f3unit(n,n);
         float da[3]; f3mul(da, n, ar);
         float db[3]; f3mul(db, n, br);
 

diff --git a/gjk.h b/gjk.h
@@ -9,7 +9,7 @@ struct gjk_support {
     float b[3];
 };
 struct gjk_simplex {
-    int iter, max_iter;
+    int max_iter, iter;
     int hit, cnt;
     struct gjk_vertex {
         float a[3];

diff --git a/gjk.h b/gjk.h
@@ -8,16 +8,15 @@ struct gjk_support {
     float a[3];
     float b[3];
 };
-struct gjk_vertex {
-    float a[3];
-    float b[3];
-    float p[3];
-    int aid, bid;
-};
 struct gjk_simplex {
     int iter, max_iter;
     int hit, cnt;
-    struct gjk_vertex v[4];
+    struct gjk_vertex {
+        float a[3];
+        float b[3];
+        float p[3];
+        int aid, bid;
+    } v[4];
     float bc[4], D;
 };
 struct gjk_result {

diff --git a/gjk.c b/gjk.c
@@ -41,27 +41,27 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         return 0;
 
     /* I.) Initialize */
-    if (!s->vcnt) {
+    if (!s->cnt) {
         s->D = GJK_FLT_MAX;
         s->max_iter = !s->max_iter ? GJK_MAX_ITERATIONS: s->max_iter;
     }
     /* II.) Check for duplications */
-    for (int i = 0; i < s->vcnt; ++i) {
+    for (int i = 0; i < s->cnt; ++i) {
         if (sup->aid != s->v[i].aid) continue;
         if (sup->bid != s->v[i].bid) continue;
         return 0;
     }
     /* III.) Add vertex into simplex */
-    struct gjk_vertex *vert = &s->v[s->vcnt];
+    struct gjk_vertex *vert = &s->v[s->cnt];
     f3cpy(vert->a, sup->a);
     f3cpy(vert->b, sup->b);
     f3cpy(vert->p, dv);
     vert->aid = sup->aid;
     vert->bid = sup->bid;
-    s->bc[s->vcnt++] = 1.0f;
+    s->bc[s->cnt++] = 1.0f;
 
     /* IV.) Find closest simplex point */
-    switch (s->vcnt) {
+    switch (s->cnt) {
     case 1: break;
     case 2: {
         /* -------------------- Line ----------------------- */
@@ -79,20 +79,20 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         if (v <= 0.0f) {
             /* region A */
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u <= 0.0f) {
             /* region B */
             s->v[0] = s->v[1];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         /* region AB */
         s->bc[0] = u;
         s->bc[1] = v;
-        s->vcnt = 2;
+        s->cnt = 2;
     } break;
     case 3: {
         /* -------------------- Triangle ----------------------- */
@@ -122,21 +122,21 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         if (v_ab <= 0.0f && u_ca <= 0.0f) {
             /* region A */
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u_ab <= 0.0f && v_bc <= 0.0f) {
             /* region B */
             s->v[0] = s->v[1];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u_bc <= 0.0f && v_ca <= 0.0f) {
             /* region C */
             s->v[0] = s->v[2];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         /* calculate fractional area */
@@ -154,7 +154,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             /* region AB */
             s->bc[0] = u_ab;
             s->bc[1] = v_ab;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (u_bc > 0.0f && v_bc > 0.0f && u_abc <= 0.0f) {
@@ -163,7 +163,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[1] = s->v[2];
             s->bc[0] = u_bc;
             s->bc[1] = v_bc;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (u_ca > 0.0f && v_ca > 0.0f && v_abc <= 0.0f) {
@@ -172,15 +172,15 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[0] = s->v[2];
             s->bc[0] = u_ca;
             s->bc[1] = v_ca;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         /* region ABC */
         assert(u_abc > 0.0f && v_abc > 0.0f && w_abc > 0.0f);
         s->bc[0] = u_abc;
         s->bc[1] = v_abc;
         s->bc[2] = w_abc;
-        s->vcnt = 3;
+        s->cnt = 3;
     } break;
     case 4: {
         /* -------------------- Tetrahedron ----------------------- */
@@ -230,28 +230,28 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         if (v_ab <= 0.0f && u_ca <= 0.0f && v_ad <= 0.0f) {
             /* region A */
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u_ab <= 0.0f && v_bc <= 0.0f && v_bd <= 0.0f) {
             /* region B */
             s->v[0] = s->v[1];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u_bc <= 0.0f && v_ca <= 0.0f && u_dc <= 0.0f) {
             /* region C */
             s->v[0] = s->v[2];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         if (u_bd <= 0.0f && v_dc <= 0.0f && u_ad <= 0.0f) {
             /* region D */
             s->v[0] = s->v[3];
             s->bc[0] = 1.0f;
-            s->vcnt = 1;
+            s->cnt = 1;
             break;
         }
         /* calculate fractional area */
@@ -297,7 +297,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             /* region AB */
             s->bc[0] = u_ab;
             s->bc[1] = v_ab;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (u_abc <= 0.0f && w_cbd <= 0.0f && u_bc > 0.0f && v_bc > 0.0f) {
@@ -306,7 +306,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[1] = s->v[2];
             s->bc[0] = u_bc;
             s->bc[1] = v_bc;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (v_abc <= 0.0f && w_acd <= 0.0f && u_ca > 0.0f && v_ca > 0.0f) {
@@ -315,7 +315,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[0] = s->v[2];
             s->bc[0] = u_ca;
             s->bc[1] = v_ca;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (v_cbd <= 0.0f && u_acd <= 0.0f && u_dc > 0.0f && v_dc > 0.0f) {
@@ -324,15 +324,15 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[1] = s->v[2];
             s->bc[0] = u_dc;
             s->bc[1] = v_dc;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (v_acd <= 0.0f && w_adb <= 0.0f && u_ad > 0.0f && v_ad > 0.0f) {
             /* region AD */
             s->v[1] = s->v[3];
             s->bc[0] = u_ad;
             s->bc[1] = v_ad;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         if (u_cbd <= 0.0f && u_adb <= 0.0f && u_bd > 0.0f && v_bd > 0.0f) {
@@ -341,7 +341,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->v[1] = s->v[3];
             s->bc[0] = u_bd;
             s->bc[1] = v_bd;
-            s->vcnt = 2;
+            s->cnt = 2;
             break;
         }
         /* calculate fractional volume (volume can be negative!) */
@@ -358,7 +358,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->bc[0] = u_abc;
             s->bc[1] = v_abc;
             s->bc[2] = w_abc;
-            s->vcnt = 3;
+            s->cnt = 3;
             break;
         }
         if (u_abcd <= 0.0f && u_cbd > 0.0f && v_cbd > 0.0f && w_cbd > 0.0f) {
@@ -368,7 +368,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->bc[0] = u_cbd;
             s->bc[1] = v_cbd;
             s->bc[2] = w_cbd;
-            s->vcnt = 3;
+            s->cnt = 3;
             break;
         }
         if (v_abcd <= 0.0f && u_acd > 0.0f && v_acd > 0.0f && w_acd > 0.0f) {
@@ -378,7 +378,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->bc[0] = u_acd;
             s->bc[1] = v_acd;
             s->bc[2] = w_acd;
-            s->vcnt = 3;
+            s->cnt = 3;
             break;
         }
         if (w_abcd <= 0.0f && u_adb > 0.0f && v_adb > 0.0f && w_adb > 0.0f) {
@@ -388,7 +388,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             s->bc[0] = u_adb;
             s->bc[1] = v_adb;
             s->bc[2] = w_adb;
-            s->vcnt = 3;
+            s->cnt = 3;
             break;
         }
         /* region ABCD */
@@ -397,21 +397,21 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
         s->bc[1] = v_abcd;
         s->bc[2] = w_abcd;
         s->bc[3] = x_abcd;
-        s->vcnt = 4;
+        s->cnt = 4;
     } break;}
 
     /* V.) Check if origin is enclosed by tetrahedron */
-    if (s->vcnt == 4) {
+    if (s->cnt == 4) {
         s->hit = 1;
         return 0;
     }
     /* VI.) Ensure closing in on origin to prevent multi-step cycling */
     {float pnt[3], denom = 0;
-    for (int i = 0; i < s->vcnt; ++i)
+    for (int i = 0; i < s->cnt; ++i)
         denom += s->bc[i];
     denom = 1.0f / denom;
 
-    switch (s->vcnt) {
+    switch (s->cnt) {
     case 1: f3cpy(pnt, s->v[0].p); break;
     case 2: {
         /* --------- Line -------- */
@@ -448,7 +448,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     s->D = d2;}
 
     /* VII.) New search direction */
-    switch (s->vcnt) {
+    switch (s->cnt) {
     default: assert(0); break;
     case 1: {
         /* --------- Point -------- */
@@ -485,12 +485,12 @@ gjk_analyze(const struct gjk_simplex *s)
 
     /* calculate normalization denominator */
     float denom = 0;
-    for (int i = 0; i < s->vcnt; ++i)
+    for (int i = 0; i < s->cnt; ++i)
         denom += s->bc[i];
     denom = 1.0f / denom;
 
     /* compute closest points */
-    switch (s->vcnt) {
+    switch (s->cnt) {
     default: assert(0); break;
     case 1: {
         /* Point */

diff --git a/gjk.h b/gjk.h
@@ -16,7 +16,7 @@ struct gjk_vertex {
 };
 struct gjk_simplex {
     int iter, max_iter;
-    int hit, vcnt;
+    int hit, cnt;
     struct gjk_vertex v[4];
     float bc[4], D;
 };

diff --git a/gjk.c b/gjk.c
@@ -42,14 +42,13 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
 
     /* I.) Initialize */
     if (!s->vcnt) {
-        s->scnt = 0;
         s->D = GJK_FLT_MAX;
         s->max_iter = !s->max_iter ? GJK_MAX_ITERATIONS: s->max_iter;
     }
     /* II.) Check for duplications */
-    for (int i = 0; i < s->scnt; ++i) {
-        if (sup->aid != s->saveA[i]) continue;
-        if (sup->bid != s->saveB[i]) continue;
+    for (int i = 0; i < s->vcnt; ++i) {
+        if (sup->aid != s->v[i].aid) continue;
+        if (sup->bid != s->v[i].bid) continue;
         return 0;
     }
     /* III.) Add vertex into simplex */
@@ -66,12 +65,14 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     case 1: break;
     case 2: {
         /* -------------------- Line ----------------------- */
-        float a[3]; f3cpy(a, s->v[0].p);
-        float b[3]; f3cpy(b, s->v[1].p);
+        float a[3], b[3];
+        f3cpy(a, s->v[0].p);
+        f3cpy(b, s->v[1].p);
 
         /* compute barycentric coordinates */
-        float ab[3]; f3sub(ab, a, b);
-        float ba[3]; f3sub(ba, b, a);
+        float ab[3], ba[3];
+        f3sub(ab, a, b);
+        f3sub(ba, b, a);
 
         float u = f3dot(b, ba);
         float v = f3dot(a, ab);
@@ -95,16 +96,18 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     } break;
     case 3: {
         /* -------------------- Triangle ----------------------- */
-        float a[3]; f3cpy(a, s->v[0].p);
-        float b[3]; f3cpy(b, s->v[1].p);
-        float c[3]; f3cpy(c, s->v[2].p);
-
-        float ab[3]; f3sub(ab, a, b);
-        float ba[3]; f3sub(ba, b, a);
-        float bc[3]; f3sub(bc, b, c);
-        float cb[3]; f3sub(cb, c, b);
-        float ca[3]; f3sub(ca, c, a);
-        float ac[3]; f3sub(ac, a, c);
+        float a[3], b[3], c[3];
+        f3cpy(a, s->v[0].p);
+        f3cpy(b, s->v[1].p);
+        f3cpy(c, s->v[2].p);
+
+        float ab[3], ba[3], bc[3], cb[3], ca[3], ac[3];
+        f3sub(ab, a, b);
+        f3sub(ba, b, a);
+        f3sub(bc, b, c);
+        f3sub(cb, c, b);
+        f3sub(ca, c, a);
+        f3sub(ac, a, c);
 
         /* compute barycentric coordinates */
         float u_ab = f3dot(b, ba);
@@ -137,10 +140,11 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
             break;
         }
         /* calculate fractional area */
-        float n[3]; f3cross(n, ba, ca);
-        float n1[3]; f3cross(n1, b, c);
-        float n2[3]; f3cross(n2, c, a);
-        float n3[3]; f3cross(n3, a, b);
+        float n[3], n1[3], n2[3], n3[3];
+        f3cross(n, ba, ca);
+        f3cross(n1, b, c);
+        f3cross(n2, c, a);
+        f3cross(n3, a, b);
 
         float u_abc = f3dot(n1, n);
         float v_abc = f3dot(n2, n);
@@ -180,24 +184,28 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     } break;
     case 4: {
         /* -------------------- Tetrahedron ----------------------- */
-        float a[3]; f3cpy(a, s->v[0].p);
-        float b[3]; f3cpy(b, s->v[1].p);
-        float c[3]; f3cpy(c, s->v[2].p);
-        float d[3]; f3cpy(d, s->v[3].p);
-
-        float ab[3]; f3sub(ab, a, b);
-        float ba[3]; f3sub(ba, b, a);
-        float bc[3]; f3sub(bc, b, c);
-        float cb[3]; f3sub(cb, c, b);
-        float ca[3]; f3sub(ca, c, a);
-        float ac[3]; f3sub(ac, a, c);
-
-        float db[3]; f3sub(db, d, b);
-        float bd[3]; f3sub(bd, b, d);
-        float dc[3]; f3sub(dc, d, c);
-        float cd[3]; f3sub(cd, c, d);
-        float da[3]; f3sub(da, d, a);
-        float ad[3]; f3sub(ad, a, d);
+        float a[3], b[3], c[3], d[3];
+        f3cpy(a, s->v[0].p);
+        f3cpy(b, s->v[1].p);
+        f3cpy(c, s->v[2].p);
+        f3cpy(d, s->v[3].p);
+
+        float ab[3], ba[3], bc[3], cb[3], ca[3], ac[3];
+        float db[3], bd[3], dc[3], cd[3], ad[3], da[3];
+
+        f3sub(ab, a, b);
+        f3sub(ba, b, a);
+        f3sub(bc, b, c);
+        f3sub(cb, c, b);
+        f3sub(ca, c, a);
+        f3sub(ac, a, c);
+
+        f3sub(db, d, b);
+        f3sub(bd, b, d);
+        f3sub(dc, d, c);
+        f3sub(cd, c, d);
+        f3sub(da, d, a);
+        f3sub(ad, a, d);
 
         /* compute barycentric coordinates */
         float u_ab = f3dot(b, ba);
@@ -466,13 +474,7 @@ gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv)
     }}
     if (f3dot(dv,dv) < GJK_EPSILON * GJK_EPSILON)
         return 0;
-
-    /* VIII.) Save ids for next duplicate check */
-    s->scnt = s->vcnt; s->iter++;
-    for (int i = 0; i < s->scnt; ++i) {
-        s->saveA[i] = s->v[i].aid;
-        s->saveB[i] = s->v[i].bid;
-    } return 1;
+    return 1;
 }
 extern struct gjk_result
 gjk_analyze(const struct gjk_simplex *s)
@@ -575,4 +577,4 @@ gjk_quad(struct gjk_result *r, float ar, float br)
         r->distance_squared = 0;
         r->hit = 1;
     }
-}
+}
diff --git a/gjk.h b/gjk.h
@@ -15,11 +15,9 @@ struct gjk_vertex {
     int aid, bid;
 };
 struct gjk_simplex {
-    int hit;
     int iter, max_iter;
-    int vcnt, scnt;
-    int saveA[5], saveB[5];
-    struct gjk_vertex v[5];
+    int hit, vcnt;
+    struct gjk_vertex v[4];
     float bc[4], D;
 };
 struct gjk_result {
@@ -29,8 +27,8 @@ struct gjk_result {
     float distance_squared;
     int iterations;
 };
-int gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv);
-struct gjk_result gjk_analyze(const struct gjk_simplex *s);
-void gjk_quad(struct gjk_result *r, float ar, float br);
+extern int gjk(struct gjk_simplex *s, const struct gjk_support *sup, float *dv);
+extern struct gjk_result gjk_analyze(const struct gjk_simplex *s);
+extern void gjk_quad(struct gjk_result *r, float ar, float br);
 
-#endif
+#endif
diff --git a/xample_simple.c b/xample_simple.c
@@ -45,5 +45,6 @@ polyhedron_intersect_capsule(const float *verts, int cnt,
     }
     /* check distance between closest points */
     struct gjk_result res = gjk_analyze(&gsx);
-    return res.distance_squared <= cr*cr;
+    gjk_quad(res, 0, cr);
+    return res.hit;
 }
diff --git a/gjk.c b/gjk.c
@@ -561,7 +561,7 @@ gjk_quad(struct gjk_result *r, float ar, float br)
         r->distance_squared -= rad*rad;
 
         /* calculate normal */
-        float n[3]; f3sub(n, r->p1, r->p0);
+        float n[3]; f3sub(n, r->p1, r->p0); f3norm(n);
         float da[3]; f3mul(da, n, ar);
         float db[3]; f3mul(db, n, br);
 

diff --git a/xample_xdebug.c b/xample_xdebug.c
@@ -48,6 +48,7 @@ debug_draw_polyhedron_intersect_sphere(int key_frame,
     key_frame = clamp(0, key_frame, n-1);
 
     struct gjk_result res = gjk_analyze(&gsx[key_frame]);
+    gjk_quad(res, 0, sr);
     glBox(res.p0, 0.05f, 0.05f, 0.05f);
     glBox(res.p1, 0.05f, 0.05f, 0.05f);
     glLine(res.p0, res.p1);