C implementation of Havel-Herout

typedef struct {
    vec3 n0; float d0;
    vec3 n1; float d1;
    vec3 n2; float d2;
} isect_hh_data;
void
isect_hh_pre(vec3 v0, vec3 v1, vec3 v2, isect_hh_data *D) {
    vec3 e1 = v3_sub(v1, v0);
    vec3 e2 = v3_sub(v2, v0);
    D->n0 = v3_cross(e1, e2);
    D->d0 = v3_dot(D->n0, v0);

    float inv_denom = 1 / v3_dot(D->n0, D->n0);

    D->n1 = v3_scale(v3_cross(e2, D->n0), inv_denom);
    D->d1 = -v3_dot(D->n1, v0);

    D->n2 = v3_scale(v3_cross(D->n0, e1), inv_denom);
    D->d2 = -v3_dot(D->n2, v0);
}
inline bool
isect_hh(vec3 o, vec3 d, float *t, vec2 *uv, isect_hh_data *D) {
    float det = v3_dot(D->n0, d);
    float dett = D->d0 - v3_dot(o, D->n0);
    vec3 wr = v3_add(v3_scale(o, det), v3_scale(d, dett));
    uv->x = v3_dot(wr, D->n1) + det * D->d1;
    uv->y = v3_dot(wr, D->n2) + det * D->d2;
    float tmpdet0 = det - uv->x - uv->y;
    int pdet0 = ((int_or_float)tmpdet0).i;
    int pdetu = ((int_or_float)uv->x).i;
    int pdetv = ((int_or_float)uv->y).i;
    pdet0 = pdet0 ^ pdetu;
    pdet0 = pdet0 | (pdetu ^ pdetv);
    if (pdet0 & 0x80000000)
        return false;
    float rdet = 1 / det;
    uv->x *= rdet;
    uv->y *= rdet;
    *t = dett * rdet;
    return *t >= ISECT_NEAR && *t <= ISECT_FAR;
}