Max 5 API Reference

jit.vecmath.c

#include "jit.vecmath.h"
#include <math.h> // for finite()
/*
 * jit.vecmath.c
 *
 * Copyright 2001-2005 - Cycling '74
 * Derek Gerstmann - derek@cycling74.com
 *
 * Basic vector math operations for floating point datatypes.
 *
 * All functions were designed to use the first function parameter
 * as the assigned/return value (eg):
 *
 *    // C = A + B
 *    jit_vec3_add(c, a, b);
 *
 * With the notable exception of all "_to_array" & "_to_coord" methods:
 *
 *      // matrix -> array
 *      jit_mat4_to_array(m, a);
 *
 *
 */

// ----------------------------------------------------------------------------

void jit_vec2_fill(t_jit_vec2 *u, float x)
{
    u->x = x;
    u->y = x;
}

void jit_vec2_assign(t_jit_vec2 *u, t_jit_vec2 *v)
{
    if(u == v)
        return;

    u->x = v->x;
    u->y = v->y;
}

void jit_vec2_from_array(t_jit_vec2 *u, float *xy)
{
    u->x = xy[0];
    u->y = xy[1];
}

void jit_vec2_from_coords(t_jit_vec2 *u, float x, float y)
{
    u->x = x;
    u->y = y;
}

void jit_vec2_to_array(t_jit_vec2 *u, float *xy)
{
    xy[0] = u->x;
    xy[1] = u->y;
}

void jit_vec2_to_coords(t_jit_vec2 *u, float *x, float *y)
{
    *x = u->x;
    *y = u->y;
}

long jit_vec2_equal(t_jit_vec2 *u, t_jit_vec2 *v)
{
    return (u->x == v->x && u->y == v->y) ? TRUE : FALSE;
}

long jit_vec2_not_equal(t_jit_vec2 *u, t_jit_vec2 *v)
{
    return (u->x != v->x || u->y != v->y) ? TRUE : FALSE;
}

void jit_vec2_mult( t_jit_vec2 *u,  t_jit_vec2 *a,  t_jit_vec2 *b)
{
    u->x = a->x * b->x;
    u->y = a->y * b->y;
}

void jit_vec2_scale( t_jit_vec2 *u,  t_jit_vec2 *a, float s)
{
    u->x = a->x * s;
    u->y = a->y * s;
}

void jit_vec2_accum_scale( t_jit_vec2 *u,  float s)
{
    u->x *= s;
    u->y *= s;
}

void jit_vec2_div( t_jit_vec2 *u,  t_jit_vec2 *a, float s)
{
    float inv = 0.0f;
    
    if(s != 0.0f)
        inv = 1.0f / s;
    
    u->x = a->x * inv;
    u->y = a->y * inv;
}

void jit_vec2_sub( t_jit_vec2 *u,  t_jit_vec2 *a,  t_jit_vec2 *b)
{
    u->x = a->x - b->x;
    u->y = a->y - b->y;
}

void jit_vec2_accum_sub(t_jit_vec2 *u, t_jit_vec2 *v) 
{
    u->x -= v->x;
    u->y -= v->y;
}

void jit_vec2_add( t_jit_vec2 *u,  t_jit_vec2 *a,  t_jit_vec2 *b)
{
    u->x = a->x + b->x;
    u->y = a->y + b->y;
}

void jit_vec2_accum_add(t_jit_vec2 *u, t_jit_vec2 *v) 
{
    u->x += v->x;
    u->y += v->y;
}

float jit_vec2_get_coord(t_jit_vec2 *u, long i) 
{
    if(i >=0 && i < 2)
        return u->vals[i];
    else
        return 0.0f;
}

float jit_vec2_mag_sqr(t_jit_vec2 *u)
{
    return u->x * u->x + u->y * u->y;
}

float jit_vec2_mag(t_jit_vec2 *u)
{
    float v = u->x * u->x + u->y * u->y;
    if(v != 0.0f)
        return jit_math_sqrt(v);
    else 
        return v;
}

void jit_vec2_negate(t_jit_vec2 *u)
{
    u->x = -u->x;
    u->y = -u->y;
}

void jit_vec2_max(t_jit_vec2 *r,  t_jit_vec2 *a,  t_jit_vec2 *b)
{
    r->x = JIT_MATH_MAX(a->x, b->x);
    r->y = JIT_MATH_MAX(a->y, b->y);
}

void jit_vec2_min(t_jit_vec2 *r,  t_jit_vec2 *a,  t_jit_vec2 *b)
{
    r->x = JIT_MATH_MIN(a->x, b->x);
    r->y = JIT_MATH_MIN(a->y, b->y);
}

void jit_vec2_lerp(t_jit_vec2 *w,  float t,  t_jit_vec2 *u,  t_jit_vec2 *v)
{ 
    w->x = JIT_MATH_LERP(t, u->x, v->x); 
    w->y = JIT_MATH_LERP(t, u->y, v->y); 
}

long jit_vec2_is_valid( t_jit_vec2 *v)
{
    return (!jit_math_is_nan(v->x) && !jit_math_is_nan(v->y) &&
        jit_math_is_finite(v->x) && jit_math_is_finite(v->y) );
}

// ----------------------------------------------------------------------------

void jit_vec3_fill(t_jit_vec3 *u, float x)
{
    u->x = x;
    u->y = x;
    u->z = x;
}

void jit_vec3_assign(t_jit_vec3 *u, t_jit_vec3 *v)
{
    if(u == v)
        return;

    u->x = v->x;
    u->y = v->y;
    u->z = v->z;
}

void jit_vec3_from_array(t_jit_vec3 *u, float *xyz)
{
    u->x = xyz[0];
    u->y = xyz[1];
    u->z = xyz[2];
}

void jit_vec3_from_coords(t_jit_vec3 *u, float x, float y, float z)
{
    u->x = x;
    u->y = y;
    u->z = z;
}

void jit_vec3_to_array(t_jit_vec3 *u, float *xyz)
{
    xyz[0] = u->x;
    xyz[1] = u->y;
    xyz[2] = u->z;
}

void jit_vec3_to_coords(t_jit_vec3 *u, float *x, float *y, float *z)
{
    *x = u->x;
    *y = u->y;
    *z = u->z;
}

long jit_vec3_equal(t_jit_vec3 *u, t_jit_vec3 *v)
{
    return (u->x == v->x && u->y == v->y && u->z == v->z) ? TRUE : FALSE;
}

long jit_vec3_not_equal(t_jit_vec3 *u, t_jit_vec3 *v)
{
    return (u->x != v->x || u->y != v->y || u->z != v->z) ? TRUE : FALSE;
}

void jit_vec3_mult( t_jit_vec3 *u,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    u->x = a->x * b->x;
    u->y = a->y * b->y;
    u->z = a->z * b->z;
}

void jit_vec3_scale( t_jit_vec3 *u,  t_jit_vec3 *a, float s)
{
    u->x = a->x * s;
    u->y = a->y * s;
    u->z = a->z * s;
}

void jit_vec3_accum_scale( t_jit_vec3 *u,  float s)
{
    u->x *= s;
    u->y *= s;
    u->z *= s;
}

void jit_vec3_div( t_jit_vec3 *u,  t_jit_vec3 *a, float s)
{
    float inv = 0.0f;
    
    if(s != 0.0f)
        inv = 1.0f / s;
    
    u->x = a->x * inv;
    u->y = a->y * inv;
    u->z = a->z * inv;
}

void jit_vec3_sub( t_jit_vec3 *u,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    u->x = a->x - b->x;
    u->y = a->y - b->y;
    u->z = a->z - b->z;
}

void jit_vec3_accum_sub(t_jit_vec3 *u, t_jit_vec3 *v) 
{
    u->x -= v->x;
    u->y -= v->y;
    u->z -= v->z;
}

void jit_vec3_add( t_jit_vec3 *u,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    u->x = a->x + b->x;
    u->y = a->y + b->y;
    u->z = a->z + b->z;
}

void jit_vec3_accum_add(t_jit_vec3 *u, t_jit_vec3 *v) 
{
    u->x += v->x;
    u->y += v->y;
    u->z += v->z;
}

void jit_vec3_madd(t_jit_vec3 *u,  t_jit_vec3 *v,  float s)
{
    u->x += v->x*s;
    u->y += v->y*s;
    u->z += v->z*s;
}

float jit_vec3_get_coord(t_jit_vec3 *u, long i) 
{
    if(i <= 0 && i < 3)
        return u->vals[i];
    else
        return 0.0f;
}

float jit_vec3_mag_sqr(t_jit_vec3 *u)
{
    return u->x * u->x + u->y * u->y + u->z * u->z;
}

float jit_vec3_mag(t_jit_vec3 *u)
{
    float v = u->x * u->x + u->y * u->y + u->z * u->z;
    if(v != 0)
        return jit_math_sqrt(v);
    else 
        return v;
}

void jit_vec3_negate(t_jit_vec3 *u)
{
    u->x = -u->x;
    u->y = -u->y;
    u->z = -u->z;
}

void jit_vec3_cross(t_jit_vec3 *u,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    u->x = a->y*b->z - a->z*b->y;
    u->y = a->z*b->x - a->x*b->z;
    u->z = a->x*b->y - a->y*b->x;
}

float jit_vec3_dot(t_jit_vec3 *u,  t_jit_vec3 *v)
{
    return u->x*v->x + u->y*v->y + u->z*v->z;
}

void jit_vec3_reflect(t_jit_vec3 *r,  t_jit_vec3 *n,  t_jit_vec3 *l)
{
    float ndotl = 2.0f * jit_vec3_dot(n, l);
    jit_vec3_scale(r, l, -1.0f);
    jit_vec3_madd(r, n, ndotl);
}

float jit_vec3_normalize( t_jit_vec3 *u )
{
    float s = JIT_MATH_F32_EPS;
    float m = u->x * u->x + u->y * u->y + u->z * u->z;
    
    if (m > JIT_MATH_F32_EPS)
    {
        s = jit_math_sqrt(m);
        s = 1.0f / s;
        u->x *= s;
        u->y *= s;
        u->z *= s;
    }
    else
    {
        u->x = s;
        u->y = s;
        u->z = s;
    }

    return s;

}

void jit_vec3_orthogonalize( t_jit_vec3 *u, t_jit_vec3 *v )
{
    float d;
    t_jit_vec3 p;
    
    jit_vec3_scale(&p, v, jit_vec3_dot( v, u ));
    
    u->x -= p.x;
    u->y -= p.y;
    u->z -= p.y;
}

void jit_vec3_orthonormalize( t_jit_vec3 *u, t_jit_vec3 *v )
{
    jit_vec3_orthogonalize( u, v );
    jit_vec3_normalize( u );       
}

void jit_vec3_max(t_jit_vec3 *r,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    r->x = JIT_MATH_MAX(a->x, b->x);
    r->y = JIT_MATH_MAX(a->y, b->y);
    r->z = JIT_MATH_MAX(a->z, b->z);
}

void jit_vec3_min(t_jit_vec3 *r,  t_jit_vec3 *a,  t_jit_vec3 *b)
{
    r->x = JIT_MATH_MIN(a->x, b->x);
    r->y = JIT_MATH_MIN(a->y, b->y);
    r->z = JIT_MATH_MIN(a->z, b->z);
}

void jit_vec3_lerp(t_jit_vec3 *w,  float t,  t_jit_vec3 *u,  t_jit_vec3 *v)
{ 
    w->x = JIT_MATH_LERP(t, u->x, v->x); 
    w->y = JIT_MATH_LERP(t, u->y, v->y); 
    w->z = JIT_MATH_LERP(t, u->z, v->z); 
}

void jit_vec3_cubemap_normal(t_jit_vec3 *v, int i, int x, int y, int cubesize)
{
    float s, t, sc, tc;
    s = ((float)x + 0.5f) / (float)cubesize;
    t = ((float)y + 0.5f) / (float)cubesize;
    sc = s * 2.0f - 1.0f;
    tc = t * 2.0f - 1.0f;

    switch (i) 
    {
        case 0:
            v->x = 1.0f;
            v->y = -tc;
            v->z = -sc;
            break;
        case 1:
            v->x = -1.0f;
            v->y = -tc;
            v->z = sc;
            break;
        case 2:
            v->x = sc;
            v->y = 1.0f;
            v->z = tc;
            break;
        case 3:
            v->x = sc;
            v->y = -1.0f;
            v->z = -tc;
            break;
        case 4:
            v->x = sc;
            v->y = -tc;
            v->z = 1.0f;
            break;
        case 5:
            v->x = -sc;
            v->y = -tc;
            v->z = -1.0f;
            break;
    }
    jit_vec3_normalize(v);
}

long jit_vec3_is_valid( t_jit_vec3 *v)
{
    return (!jit_math_is_nan(v->x) && !jit_math_is_nan(v->y) && !jit_math_is_nan(v->z) &&
        jit_math_is_finite(v->x) && jit_math_is_finite(v->y) && jit_math_is_finite(v->z));
}

// ----------------------------------------------------------------------------

void jit_vec4_fill(t_jit_vec4 *u, float x)
{
    u->x = x;
    u->y = x;
    u->z = x;
    u->w = x;
}

void jit_vec4_assign(t_jit_vec4 *u, t_jit_vec4 *v)
{
    if(u == v)
        return;
        
    u->x = v->x;
    u->y = v->y;
    u->z = v->z;
    u->w = v->w;
}

void jit_vec4_from_array(t_jit_vec4 *u, float *xyzw)
{
    u->x = xyzw[0];
    u->y = xyzw[1];
    u->z = xyzw[2];
    u->w = xyzw[3];
}

void jit_vec4_from_coords(t_jit_vec4 *u, float x, float y, float z, float w)
{
    u->x = x;
    u->y = y;
    u->z = z;
    u->w = w;
}

void jit_vec4_to_array(t_jit_vec4 *u, float *xyzw)
{
    xyzw[0] = u->x;
    xyzw[1] = u->y;
    xyzw[2] = u->z;
    xyzw[3] = u->w;
}

void jit_vec4_to_coords(t_jit_vec4 *u, float *x, float *y, float *z, float *w)
{
    *x = u->x;
    *y = u->y;
    *z = u->z;
    *w = u->w;
}

long jit_vec4_equal(t_jit_vec4 *u, t_jit_vec4 *v)
{
    return (u->x == v->x && u->y == v->y && u->z == v->z && u->w == v->w) ? TRUE : FALSE;
}

long jit_vec4_not_equal(t_jit_vec4 *u, t_jit_vec4 *v)
{
    return (u->x != v->x || u->y != v->y || u->z != v->z || u->w != v->w) ? TRUE : FALSE;
}

void jit_vec4_mult( t_jit_vec4 *u,  t_jit_vec4 *a,  t_jit_vec4 *b)
{
    u->x = a->x * b->x;
    u->y = a->y * b->y;
    u->z = a->z * b->z;
    u->w = a->w * b->w;
}

void jit_vec4_scale( t_jit_vec4 *u,  t_jit_vec4 *a, float s)
{
    u->x = a->x * s;
    u->y = a->y * s;
    u->z = a->z * s;
    u->w = a->w * s;
}

void jit_vec4_accum_scale( t_jit_vec4 *u,  float s)
{
    u->x *= s;
    u->y *= s;
    u->z *= s;
    u->w *= s;
}

void jit_vec4_div( t_jit_vec4 *u,  t_jit_vec4 *a, float s)
{
    float inv = 0.0f;
    
    if(s != 0.0f)
        inv = 1.0f / s;
    
    u->x = a->x * inv;
    u->y = a->y * inv;
    u->z = a->z * inv;
    u->w = a->z * inv;
}

void jit_vec4_sub( t_jit_vec4 *u,  t_jit_vec4 *a,  t_jit_vec4 *b)
{
    u->x = a->x - b->x;
    u->y = a->y - b->y;
    u->z = a->z - b->z;
    u->w = a->w - b->w;
}

void jit_vec4_accum_sub(t_jit_vec4 *u, t_jit_vec4 *v) 
{
    u->x -= v->x;
    u->y -= v->y;
    u->z -= v->z;
    u->w -= v->w;
}

void jit_vec4_add( t_jit_vec4 *u,  t_jit_vec4 *a,  t_jit_vec4 *b)
{
    u->x = a->x + b->x;
    u->y = a->y + b->y;
    u->z = a->z + b->z;
    u->w = a->w + b->w;
}

void jit_vec4_accum_add(t_jit_vec4 *u, t_jit_vec4 *v) 
{
    u->x += v->x;
    u->y += v->y;
    u->z += v->z;
    u->w += v->w;
}

float jit_vec4_get_coord(t_jit_vec4 *u, long i) 
{
    if(i >= 0 && i < 4)
        return u->vals[i];
    else
        return 0.0f;
}

float jit_vec4_mag_sqr(t_jit_vec4 *u)
{
    return u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
}

float jit_vec4_mag(t_jit_vec4 *u)
{
    float v = u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
    if(v != 0.0f)
        return jit_math_sqrt(v);
    else 
        return v;
}

void jit_vec4_negate(t_jit_vec4 *u)
{
    u->x = -u->x;
    u->y = -u->y;
    u->z = -u->z;
    u->w = -u->w;
}

float jit_vec4_dot( t_jit_vec4 *u,  t_jit_vec4 *v)
{
    return u->x*v->x + u->y*v->y + u->z*v->z + u->w*v->w;
}

float jit_vec4_normalize( t_jit_vec4 *u )
{
    float s = JIT_MATH_F32_EPS;
    float m = u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
    
    if (m > JIT_MATH_F32_EPS)
    {
        s = jit_math_sqrt(m);
        s = 1.0f / s;
        u->x *= s;
        u->y *= s;
        u->z *= s;
        u->w *= s;
    }
    else
    {
        u->x = s;
        u->y = s;
        u->z = s;
        u->w = s;
    }

    return s;
}

void jit_vec4_max(t_jit_vec4 *r,  t_jit_vec4 *a,  t_jit_vec4 *b)
{
    r->x = JIT_MATH_MAX(a->x, b->x);
    r->y = JIT_MATH_MAX(a->y, b->y);
    r->z = JIT_MATH_MAX(a->z, b->z);
    r->w = JIT_MATH_MAX(a->w, b->w);
}

void jit_vec4_min(t_jit_vec4 *r,  t_jit_vec4 *a,  t_jit_vec4 *b)
{
    r->x = JIT_MATH_MIN(a->x, b->x);
    r->y = JIT_MATH_MIN(a->y, b->y);
    r->z = JIT_MATH_MIN(a->z, b->z);
    r->w = JIT_MATH_MIN(a->w, b->w);
}

void jit_vec4_lerp(t_jit_vec4 *w,  float t,  t_jit_vec4 *u,  t_jit_vec4 *v)
{ 
    w->x = JIT_MATH_LERP(t, u->x, v->x); 
    w->y = JIT_MATH_LERP(t, u->y, v->y); 
    w->z = JIT_MATH_LERP(t, u->z, v->z); 
    w->w = JIT_MATH_LERP(t, u->w, v->w); 
}

long jit_vec4_is_valid( t_jit_vec4 *v)
{
    return (!jit_math_is_nan(v->x) && !jit_math_is_nan(v->y) && !jit_math_is_nan(v->z) && !jit_math_is_nan(v->w) &&
        jit_math_is_finite(v->x) && jit_math_is_finite(v->y) && jit_math_is_finite(v->z) && jit_math_is_finite(v->w));
}

// ----------------------------------------------------------------------------

void jit_quat_fill(t_jit_quat *u, float x)
{
    u->x = x;
    u->y = x;
    u->z = x;
    u->w = x;
}

void jit_quat_assign(t_jit_quat *u, t_jit_quat *v)
{
    if(u == v)
        return;
        
    u->x = v->x;
    u->y = v->y;
    u->z = v->z;
    u->w = v->w;
}

void jit_quat_from_array(t_jit_quat *u, float *xyzw)
{
    u->x = xyzw[0];
    u->y = xyzw[1];
    u->z = xyzw[2];
    u->w = xyzw[3];
}

void jit_quat_from_coords(t_jit_quat *u, float x, float y, float z, float w)
{
    u->x = x;
    u->y = y;
    u->z = z;
    u->w = w;
}

void jit_quat_to_array(t_jit_quat *u, float *xyzw)
{
    xyzw[0] = u->x;
    xyzw[1] = u->y;
    xyzw[2] = u->z;
    xyzw[3] = u->w;
}

void jit_quat_to_coords(t_jit_quat *u, float *x, float *y, float *z, float *w)
{
    *x = u->x;
    *y = u->y;
    *z = u->z;
    *w = u->w;
}

long jit_quat_equal(t_jit_quat *u, t_jit_quat *v)
{
    return (u->x == v->x && u->y == v->y && u->z == v->z && u->w == v->w) ? TRUE : FALSE;
}

long jit_quat_not_equal(t_jit_quat *u, t_jit_quat *v)
{
    return (u->x != v->x || u->y != v->y || u->z != v->z || u->w != v->w) ? TRUE : FALSE;
}

void jit_quat_scale( t_jit_quat *u,  t_jit_quat *a, float s)
{
    u->x = a->x * s;
    u->y = a->y * s;
    u->z = a->z * s;
    u->w = a->w * s;
}

void jit_quat_accum_scale( t_jit_quat *u,  float s)
{
    u->x *= s;
    u->y *= s;
    u->z *= s;
    u->w *= s;
}

void jit_quat_div( t_jit_quat *u,  t_jit_quat *a, float s)
{
    float inv = 0.0f;
    
    if(s != 0.0f)
        inv = 1.0f / s;
    
    u->x = a->x * inv;
    u->y = a->y * inv;
    u->z = a->z * inv;
    u->w = a->z * inv;
}

void jit_quat_add(t_jit_quat * r,  t_jit_quat * q1,  t_jit_quat * q2)
{
    t_jit_quat t1, t2;

    t1.x = q1->x * q2->w;
    t1.y = q1->y * q2->w;
    t1.z = q1->z * q2->w;
    t1.w = q1->w;

    t2.x = q2->x * q1->w;
    t2.y = q2->y * q1->w;
    t2.z = q2->z * q1->w;
    t2.w = q2->w;

    r->x = (q2->y * q1->z) - (q2->z * q1->y) + t1.x + t2.x;
    r->y = (q2->z * q1->x) - (q2->x * q1->z) + t1.y + t2.y;
    r->z = (q2->x * q1->y) - (q2->y * q1->x) + t1.z + t2.z;
    r->w = q1->w * q2->w - (q1->x * q2->x + q1->y * q2->y + q1->z * q2->z);
}

void jit_quat_mult( t_jit_quat *r, t_jit_quat * p,  t_jit_quat * q)
{
    r->x = p->w * q->x + p->x * q->w + p->y * q->z - p->z * q->y;
    r->y = p->w * q->y + p->y * q->w + p->z * q->x - p->x * q->z;
    r->z = p->w * q->z + p->z * q->w + p->x * q->y - p->y * q->x;
    r->w = p->w * q->w - p->x * q->x - p->y * q->y - p->z * q->z;
}

float jit_quat_dot( t_jit_quat * q1,  t_jit_quat * q2)
{
    return q1->x*q2->x + q1->y*q2->y + q1->z*q2->z + q1->w*q2->w;
}

float jit_quat_get_coord( t_jit_quat *u, long i) 
{
    if(i >= 0 && i < 4)
        return u->vals[i];
    else
        return 0.0f;
}

float jit_quat_mag_sqr(t_jit_quat *u)
{
    return u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
}

float jit_quat_mag(t_jit_quat *u)
{
    float v = u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
    if(v != 0.0f)
        return jit_math_sqrt(v);
    else 
        return v;
}

void jit_quat_negate(t_jit_quat *u)
{
    u->x = -u->x;
    u->y = -u->y;
    u->z = -u->z;
    u->w = -u->w;
}

float jit_quat_normalize( t_jit_quat *u )
{
    float s = JIT_MATH_F32_EPS;
    float m = u->x * u->x + u->y * u->y + u->z * u->z + u->w * u->w;
    
    if (m > JIT_MATH_F32_EPS)
    {
        s = jit_math_sqrt(m);
        s = 1.0f / s;
        u->x *= s;
        u->y *= s;
        u->z *= s;
        u->w *= s;
    }
    else
    {
        u->x = s;
        u->y = s;
        u->z = s;
        u->w = s;
    }

    return s;
}

void jit_quat_max(t_jit_quat *r,  t_jit_quat *a,  t_jit_quat *b)
{
    r->x = JIT_MATH_MAX(a->x, b->x);
    r->y = JIT_MATH_MAX(a->y, b->y);
    r->z = JIT_MATH_MAX(a->z, b->z);
    r->w = JIT_MATH_MAX(a->w, b->w);
}

void jit_quat_min(t_jit_quat *r,  t_jit_quat *a,  t_jit_quat *b)
{
    r->x = JIT_MATH_MIN(a->x, b->x);
    r->y = JIT_MATH_MIN(a->y, b->y);
    r->z = JIT_MATH_MIN(a->z, b->z);
    r->w = JIT_MATH_MIN(a->w, b->w);
}

void jit_quat_from_mat3( t_jit_quat *q, t_jit_mat3 *m)
{
    static long next[] = { 1, 2, 0 };
    
    long i, j, k;
    float vi, vj, vk;
    float scale;
    
    float trace = m->v00 + m->v11 + m->v22;
    
    if (trace > 0.0f) 
    {
        scale = jit_math_sqrt(trace + 1.0f);
        q->w = 0.5f * scale;
        scale = 0.5f / scale;
        q->x = scale * (m->v21 - m->v12);
        q->y = scale * (m->v02 - m->v20);
        q->z = scale * (m->v10 - m->v01);
    }
    else 
    {
        i = 0;
        if (m->v11 > m->v00)
            i = 1;
            
        vi = jit_mat3_get_coord(m, i, i);
        if (m->v22 > vi)
            i = 2;
        
        j = next[i];
        k = next[j];
        
        vi = jit_mat3_get_coord(m, i, i);
        vj = jit_mat3_get_coord(m, j, j);
        vk = jit_mat3_get_coord(m, k, k);
        
        scale = jit_math_sqrt(vi - vj - vk + 1);

        q->vals[i] = 0.5f * scale;
        scale = 0.5f / scale;

        vi = jit_mat3_get_coord(m, k, j);
        vj = jit_mat3_get_coord(m, j, k);
        q->w = scale * (vi - vj);
        
        vi = jit_mat3_get_coord(m, j, i);
        vj = jit_mat3_get_coord(m, i, j);
        q->vals[j] = scale * (vi + vj);
        
        vi = jit_mat3_get_coord(m, k, i);
        vj = jit_mat3_get_coord(m, i, k);
        q->vals[k] = scale * (vi + vj);
    }
}

void jit_mat3_from_quat(t_jit_mat3 *m, t_jit_quat *q ) 
{
    float x2 = q->x * 2;
    float y2 = q->y * 2;
    float z2 = q->z * 2;
    
    float wx = x2 * q->w;
    float wy = y2 * q->w;
    float wz = z2 * q->w;
    float xx = x2 * q->x;
    float xy = y2 * q->x;
    float xz = z2 * q->x;
    float yy = y2 * q->y;
    float yz = z2 * q->y;
    float zz = z2 * q->z;
    
    m->v00 = 1 - (yy + zz);
    m->v01 = xy - wz;
    m->v02 = xz + wy;
    m->v10 = xy + wz;
    m->v11 = 1 - (xx + zz);
    m->v12 = yz - wx;
    m->v20 = xz - wy;
    m->v21 = yz + wx;
    m->v22 = 1 - (xx + yy);
}

void jit_quat_from_mat4(t_jit_quat * q,  t_jit_mat4 *m)
{
    t_jit_mat3 r;
    jit_mat3_from_mat4(&r, m);
    jit_quat_from_mat3(q, &r);
} 

void jit_quat_from_axisangle(t_jit_quat * q,  t_jit_vec3 *a,  float angle)
{
    t_jit_vec3 t;
    float scale;
    float halfangle;
    
    t.x = a->x;
    t.y = a->y;
    t.z = a->z;
    
    jit_vec3_normalize(&t);
    
    halfangle = angle * 0.5f;
    scale = jit_math_sin(halfangle);
    
    q->x = scale * t.x;
    q->y = scale * t.y;
    q->z = scale * t.z;
    q->w = jit_math_cos(halfangle);
}

void jit_quat_conj(t_jit_quat * p)
{
    p->x = -p->x;
    p->y = -p->y;
    p->z = -p->z;
}

void jit_quat_get_conj(t_jit_quat * p,  t_jit_quat * q)
{
    p->x = -q->x;
    p->y = -q->y;
    p->z = -q->z;
    p->w = q->w;
}

void jit_quat_slerp(t_jit_quat *r, float s,  t_jit_quat * q1,  t_jit_quat * q2)
{
    float cosine;
    float angle;
    float sine;
    float invsine;
    float c1, c2;
    
    cosine = jit_quat_dot(q1, q2);
    if (cosine < -1)
        cosine = -1;
    else if (cosine > 1)
        cosine = 1;
        
    angle = (float)jit_math_acos(cosine);
    if (JIT_MATH_F32_ABS(angle) < JIT_MATH_F32_EPS) 
    {
        r->x = q1->x;
        r->y = q1->y;
        r->z = q1->z;
        r->w = q1->w;
        return;
    }
    
    sine = jit_math_sin(angle);
    invsine = 1.0f / sine;
    
    c1 = jit_math_sin((1.0f - s) * angle) * invsine;
    c2 = jit_math_sin(s * angle) * invsine;
    
    r->x = c1 * q1->x + c2 * q2->x;
    r->y = c1 * q1->y + c2 * q2->y;
    r->z = c1 * q1->z + c2 * q2->z;
    r->w = c1 * q1->w + c2 * q2->w;
}

void jit_quat_trackball(t_jit_quat * q, t_jit_vec2 *p1, t_jit_vec2 *p2, float radius)
{
    t_jit_vec3 d;
    t_jit_vec3 axis;
    t_jit_vec3 v1, v2;

    float t;
    float phi;

    if (p1->x == p2->x && p1->y == p2->y) 
    {
        // zero rotation
        q->x = 0.0f;
        q->y = 0.0f;
        q->z = 0.0f;
        q->w = 1.0f;
        return;
    }

    // get z-coordinates for projection of P1 and P2 to deformed sphere
    jit_vec3_from_coords(&v1, p1->x, p1->y, jit_math_project_to_sphere(radius, p1->x, p1->y));
    jit_vec3_from_coords(&v2, p2->x, p2->y, jit_math_project_to_sphere(radius, p2->x, p2->y));

    //  cross product of P1 and P2
    jit_vec3_cross(&axis, &v1, &v2);

    //  figure out how much to rotate around that axis.
    d.x = v1.x - v2.x;
    d.y = v1.y - v2.y;
    d.z = v1.z - v2.z;
    t = jit_math_sqrt(d.x * d.x + d.y * d.y + d.z * d.z) / (radius);

    // avoid problems with out-of-control values...
    if (t > 1.0f)
        t = 1.0f;
    if (t < -1.0f) 
        t = -1.0f;
    
    phi = 2.0f * (float)jit_math_asin(t);
    jit_quat_from_axisangle(q, &axis, phi);
}

long jit_quat_is_valid( t_jit_quat *v)
{
    return (!jit_math_is_nan(v->x) && !jit_math_is_nan(v->y) && !jit_math_is_nan(v->z) && !jit_math_is_nan(v->w) &&
        jit_math_is_finite(v->x) && jit_math_is_finite(v->y) && jit_math_is_finite(v->z) && jit_math_is_finite(v->w));
}

// ----------------------------------------------------------------------------

void jit_mat3_assign( t_jit_mat3 *m, t_jit_mat3 *a )
{
    // avoid self assignment
    if(m == a)
        return;
        
    m->v00 = a->v00; m->v10 = a->v10; m->v20 = a->v20; 
    m->v01 = a->v01; m->v11 = a->v11; m->v21 = a->v21;
    m->v02 = a->v02; m->v12 = a->v12; m->v22 = a->v22;
}

void jit_mat3_from_array( t_jit_mat3 *m, float *v )
{
    m->v00 = v[0]; m->v10 = v[1]; m->v20 = v[2]; 
    m->v01 = v[3]; m->v11 = v[4]; m->v21 = v[5];
    m->v02 = v[6]; m->v12 = v[7]; m->v22 = v[8];
}

void jit_mat3_from_array3x3( t_jit_mat3 *m, float **v )
{
    m->v00 = v[0][0]; m->v10 = v[1][0]; m->v20 = v[2][0]; 
    m->v01 = v[0][1]; m->v11 = v[1][1]; m->v21 = v[2][1];
    m->v02 = v[0][2]; m->v12 = v[1][2]; m->v22 = v[2][2];
}

void jit_mat3_from_coords( t_jit_mat3 *m,  
    float f0,   float f1,   float f2,  
    float f3,   float f4,   float f5,  
    float f6,   float f7,   float f8 )
{
    m->v00 = f0; m->v10 = f1; m->v20 = f2; 
    m->v01 = f3; m->v11 = f4; m->v21 = f5;
    m->v02 = f6; m->v12 = f7; m->v22 = f8;
}

void jit_mat3_to_array( t_jit_mat3 *m, float *v )
{
    m->v00 = v[0]; m->v10 = v[1]; m->v20 = v[2]; 
    m->v01 = v[3]; m->v11 = v[4]; m->v21 = v[5];
    m->v02 = v[6]; m->v12 = v[7]; m->v22 = v[8];
}

void jit_mat3_to_array3x3( t_jit_mat3 *m, float **v )
{
    v[0][0] = m->v00; v[1][0] = m->v10; v[2][0] = m->v20; 
    v[0][1] = m->v01; v[1][1] = m->v11; v[2][1] = m->v21;
    v[0][2] = m->v02; v[1][2] = m->v12; v[2][2] = m->v22;
}

void jit_mat3_to_coords( t_jit_mat3 *m,  
    float *f0,   float *f1,   float *f2,  
    float *f3,   float *f4,   float *f5,  
    float *f6,   float *f7,   float *f8 )
{
    *f0 = m->v00; *f1 = m->v10; *f2 = m->v20; 
    *f3 = m->v01; *f4 = m->v11; *f5 = m->v21;
    *f6 = m->v02; *f7 = m->v12; *f8 = m->v22;
}

void jit_mat3_get_row( t_jit_vec3 *v, long i, t_jit_mat3 *m )
{
    if( i >= 0 && i < 3 ) 
    {
        v->x = m->vals[i + 0];
        v->y = m->vals[i + 3];
        v->z = m->vals[i + 6];
    }
}

void jit_mat3_set_row( t_jit_mat3 *m, long i,  t_jit_vec3 *v )
{
    if( i >= 0 && i < 3 ) 
    {
        m->vals[i + 0] = v->x;
        m->vals[i + 3] = v->y;
        m->vals[i + 6] = v->z;
    }
}

void jit_mat3_get_col( t_jit_vec3 *v, long i, t_jit_mat3 *m )
{
    if( i >= 0 && i < 3 ) 
    {
        v->x = m->vals[i * 3 + 0];
        v->y = m->vals[i * 3 + 1];
        v->z = m->vals[i * 3 + 2];
    }
}

void jit_mat3_set_col( t_jit_mat3 *m, long i,  t_jit_vec3 *v)
{
    if( i >= 0 && i < 3 ) 
    {
        m->vals[i * 3 + 0] = v->x;
        m->vals[i * 3 + 1] = v->y;
        m->vals[i * 3 + 2] = v->z;
    }
}

float jit_mat3_get_coord( t_jit_mat3 *m, long r, long c )
{
    if(r >= 0 && r < 3 && c >= 0 && c < 3 )
        return m->vals[ c * 3 + r ];
    else
        return 0.0f;
}

void jit_mat3_add(t_jit_mat3 *c,  t_jit_mat3 *a,  t_jit_mat3 *b)
{
    // avoid self assignment
    if (c == a || c == b)
    {
        t_jit_mat3 m;

        m.v00 = a->v00 + b->v00;
        m.v01 = a->v01 + b->v01;
        m.v02 = a->v02 + b->v02;
        m.v10 = a->v10 + b->v10;
        m.v11 = a->v11 + b->v11;
        m.v12 = a->v12 + b->v12;
        m.v20 = a->v20 + b->v20;
        m.v21 = a->v21 + b->v21;
        m.v22 = a->v22 + b->v22;

        jit_mat3_assign(c, &m);
    }
    else
    {
        c->v00 = a->v00 + b->v00;
        c->v01 = a->v01 + b->v01;
        c->v02 = a->v02 + b->v02;
        c->v10 = a->v10 + b->v10;
        c->v11 = a->v11 + b->v11;
        c->v12 = a->v12 + b->v12;
        c->v20 = a->v20 + b->v20;
        c->v21 = a->v21 + b->v21;
        c->v22 = a->v22 + b->v22;
    }
}

void jit_mat3_accum_add(t_jit_mat3 *a,  t_jit_mat3 *b)
{
    a->v00 += b->v00;
    a->v10 += b->v10;
    a->v20 += b->v20;
    a->v01 += b->v01;
    a->v11 += b->v11;
    a->v21 += b->v21;
    a->v02 += b->v02;
    a->v12 += b->v12;
    a->v22 += b->v22;
}

void jit_mat3_accum_scale( t_jit_mat3 *m, float s)
{
    m->v00 *= s; m->v10 *= s; m->v20 *= s; 
    m->v01 *= s; m->v11 *= s; m->v21 *= s;
    m->v02 *= s; m->v12 *= s; m->v22 *= s;
}

void jit_mat3_accum_sub( t_jit_mat3 *m, t_jit_mat3 *a )
{
    m->v00 -= a->v00; m->v10 -= a->v10; m->v20 -= a->v20; 
    m->v01 -= a->v01; m->v11 -= a->v11; m->v21 -= a->v21;
    m->v02 -= a->v02; m->v12 -= a->v12; m->v22 -= a->v22;
}

void jit_mat3_mult(t_jit_mat3 *c, t_jit_mat3 *a, t_jit_mat3 *b)
{
    if (c == a || c == b)
    {
        t_jit_mat3 m;

        m.v00 = a->v00 * b->v00 + a->v01 * b->v10 + a->v02 * b->v20;
        m.v10 = a->v10 * b->v00 + a->v11 * b->v10 + a->v12 * b->v20;
        m.v20 = a->v20 * b->v00 + a->v21 * b->v10 + a->v22 * b->v20;
        m.v01 = a->v00 * b->v01 + a->v01 * b->v11 + a->v02 * b->v21;
        m.v11 = a->v10 * b->v01 + a->v11 * b->v11 + a->v12 * b->v21;
        m.v21 = a->v20 * b->v01 + a->v21 * b->v11 + a->v22 * b->v21;
        m.v02 = a->v00 * b->v02 + a->v01 * b->v12 + a->v02 * b->v22;
        m.v12 = a->v10 * b->v02 + a->v11 * b->v12 + a->v12 * b->v22;
        m.v22 = a->v20 * b->v02 + a->v21 * b->v12 + a->v22 * b->v22;

        jit_mat3_assign(c, &m);
    }
    else
    {
        c->v00 = a->v00 * b->v00 + a->v01 * b->v10 + a->v02 * b->v20;
        c->v10 = a->v10 * b->v00 + a->v11 * b->v10 + a->v12 * b->v20;
        c->v20 = a->v20 * b->v00 + a->v21 * b->v10 + a->v22 * b->v20;
        c->v01 = a->v00 * b->v01 + a->v01 * b->v11 + a->v02 * b->v21;
        c->v11 = a->v10 * b->v01 + a->v11 * b->v11 + a->v12 * b->v21;
        c->v21 = a->v20 * b->v01 + a->v21 * b->v11 + a->v22 * b->v21;
        c->v02 = a->v00 * b->v02 + a->v01 * b->v12 + a->v02 * b->v22;
        c->v12 = a->v10 * b->v02 + a->v11 * b->v12 + a->v12 * b->v22;
        c->v22 = a->v20 * b->v02 + a->v21 * b->v12 + a->v22 * b->v22;
    }
}

void jit_mat3_get_transpose(t_jit_mat3 *b, t_jit_mat3 *a)
{
    b->v00 = a->v00;
    b->v01 = a->v10;
    b->v02 = a->v20;
    b->v10 = a->v01;
    b->v11 = a->v11;
    b->v12 = a->v21;
    b->v20 = a->v02;
    b->v21 = a->v12;
    b->v22 = a->v22;
}

void jit_mat3_transpose(t_jit_mat3 *a)
{
    float t;
    t = a->v01;
    a->v01 = a->v10;
    a->v10 = t;

    t = a->v02;
    a->v02 = a->v20;
    a->v20 = t;

    t = a->v12;
    a->v12 = a->v21;
    a->v21 = t;
}

void jit_mat3_mult_vec3( t_jit_vec3 *r, t_jit_mat3 *m, t_jit_vec3 *v )
{
    r->x = m->v00 * v->x + m->v01 * v->y + m->v02 * v->z;
    r->y = m->v10 * v->x + m->v11 * v->y + m->v12 * v->z;
    r->z = m->v20 * v->x + m->v21 * v->y + m->v22 * v->z;
}

void jit_vec3_mult_mat3( t_jit_vec3 *r, t_jit_vec3 *v, t_jit_mat3 *m )
{
    r->x = m->v00 * v->x + m->v10 * v->y + m->v20 * v->z;
    r->y = m->v01 * v->x + m->v11 * v->y + m->v21 * v->z;
    r->z = m->v02 * v->x + m->v12 * v->y + m->v22 * v->z;
}

void jit_mat3_from_axisangle( t_jit_mat3 *m, t_jit_vec3 *v, float theta ) 
{
    float ct = (float)(jit_math_cos(theta));
    float st = (float)(jit_math_sin(theta));

    float xx = v->x * v->x;
    float yy = v->y * v->y;
    float zz = v->z * v->z;
    float xy = v->x * v->y;
    float xz = v->x * v->z;
    float yz = v->y * v->z;

    m->v00 = xx + ct*(1-xx);
    m->v01 = xy + ct*(-xy) + st*-v->z;
    m->v02 = xz + ct*(-xz) + st*v->y;

    m->v10 = xy + ct*(-xy) + st*v->z;
    m->v11 = yy + ct*(1-yy);
    m->v12 = yz + ct*(-yz) + st*-v->x;

    m->v20 = xz + ct*(-xz) + st*-v->y;
    m->v21 = yz + ct*(-yz) + st*v->x;
    m->v22 = zz + ct*(1-zz);
}

void jit_mat3_from_uv( t_jit_mat3 *m, t_jit_vec3 *u,  t_jit_vec3 *v )
{
    float phi;
    float h, hxy, hxz, hyz;
    float lambda;
    t_jit_vec3 w;

    jit_vec3_cross(&w, u, v);   
    phi = jit_vec3_dot(u, v);
    lambda = jit_vec3_dot(&w, &w);

    if (lambda > JIT_MATH_F32_EPS)
        h = (1.0f - phi) / lambda;
    else
        h = lambda;

    hxy = w.x * w.y * h;
    hxz = w.x * w.z * h;
    hyz = w.y * w.z * h;

    m->v00 = phi + w.x * w.x * h;
    m->v01 = hxy - w.z;
    m->v02 = hxz + w.y;

    m->v10 = hxy + w.z;
    m->v11 = phi + w.y * w.y * h;
    m->v12 = hyz - w.x;

    m->v20 = hxz - w.y;
    m->v21 = hyz + w.x;
    m->v22 = phi + w.z * w.z * h;
}

float jit_mat3_determinant( t_jit_mat3 *m )
{
    return JIT_MATH_DET3x3(
            m->v00, m->v01, m->v02, 
            m->v10, m->v11, m->v12, 
            m->v20, m->v21, m->v22);
}

void jit_mat3_get_inverse(t_jit_mat3 *b,  t_jit_mat3 *a)
{
    float det,oodet;

    b->v00 =  (a->v11 * a->v22 - a->v21 * a->v12);
    b->v10 = -(a->v10 * a->v22 - a->v20 * a->v12);
    b->v20 =  (a->v10 * a->v21 - a->v20 * a->v11);
    b->v01 = -(a->v01 * a->v22 - a->v21 * a->v02);
    b->v11 =  (a->v00 * a->v22 - a->v20 * a->v02);
    b->v21 = -(a->v00 * a->v21 - a->v20 * a->v01);
    b->v02 =  (a->v01 * a->v12 - a->v11 * a->v02);
    b->v12 = -(a->v00 * a->v12 - a->v10 * a->v02);
    b->v22 =  (a->v00 * a->v11 - a->v10 * a->v01);

    det = (a->v00 * b->v00) + (a->v01 * b->v10) + (a->v02 * b->v20);
    
    oodet = 1.0f / det;

    b->v00 *= oodet; b->v01 *= oodet; b->v02 *= oodet;
    b->v10 *= oodet; b->v11 *= oodet; b->v12 *= oodet;
    b->v20 *= oodet; b->v21 *= oodet; b->v22 *= oodet;
}


void jit_mat3_negate( t_jit_mat3 *m )
{
    m->v00 = -m->v00; m->v10 = -m->v10; m->v20 = -m->v20; 
    m->v01 = -m->v01; m->v11 = -m->v11; m->v21 = -m->v21;
    m->v02 = -m->v02; m->v12 = -m->v12; m->v22 = -m->v22;
}

void jit_mat3_tangent_basis(t_jit_mat3 *basis,
    t_jit_vec3 *v0, t_jit_vec3 *v1, t_jit_vec3 *v2,
    t_jit_vec2 *t0, t_jit_vec2 *t1, t_jit_vec2 *t2,  
    t_jit_vec3 *n)
{
    float oonorm;
    t_jit_vec3 cp;
    t_jit_vec3 e0;
    t_jit_vec3 e1;
    
    jit_vec3_from_coords(&e0, v1->x - v0->x, t1->s - t0->s, t1->t - t0->t);
    jit_vec3_from_coords(&e1, v2->x - v0->x, t2->s - t0->s, t2->t - t0->t);

    jit_vec3_cross(&cp, &e0, &e1);
    if ( JIT_MATH_F32_ABS(cp.x) > JIT_MATH_F32_EPS)
    {
        basis->v00 = -cp.y / cp.x;        
        basis->v10 = -cp.z / cp.x;
    }

    e0.x = v1->y - v0->y;
    e1.x = v2->y - v0->y;

    jit_vec3_cross(&cp, &e0, &e1);
    if ( JIT_MATH_F32_ABS(cp.x) > JIT_MATH_F32_EPS)
    {
        basis->v01 = -cp.y / cp.x;        
        basis->v11 = -cp.z / cp.x;
    }

    e0.x = v1->z - v0->z;
    e1.x = v2->z - v0->z;

    jit_vec3_cross(&cp, &e0, &e1);
    if ( JIT_MATH_F32_ABS(cp.x) > JIT_MATH_F32_EPS)
    {
        basis->v02 = -cp.y / cp.x;        
        basis->v12 = -cp.z / cp.x;
    }

    // tangent
    oonorm = 1.0f / jit_math_sqrt(basis->v00 * basis->v00 + basis->v01 * basis->v01 + basis->v02 * basis->v02);
    basis->v00 *= oonorm;
    basis->v01 *= oonorm;
    basis->v02 *= oonorm;

    // binormal
    oonorm = 1.0f / jit_math_sqrt(basis->v10 * basis->v10 + basis->v11 * basis->v11 + basis->v12 * basis->v12);
    basis->v10 *= oonorm;
    basis->v11 *= oonorm;
    basis->v12 *= oonorm;

    // normal via cross product
    basis->v20 = basis->v01*basis->v12 - basis->v02*basis->v11;
    basis->v21 = basis->v02*basis->v10 - basis->v00*basis->v12;
    basis->v22 = basis->v00*basis->v11 - basis->v01*basis->v10;

    oonorm = 1.0f / jit_math_sqrt(basis->v20 * basis->v20 + basis->v21 * basis->v21 + basis->v22 * basis->v22);
    basis->v20 *= oonorm;
    basis->v21 *= oonorm;
    basis->v22 *= oonorm;

    // Gram-Schmidt orthogonalization process
    basis->v10 = basis->v21*basis->v02 - basis->v22*basis->v01;
    basis->v11 = basis->v22*basis->v00 - basis->v20*basis->v02;
    basis->v12 = basis->v20*basis->v01 - basis->v21*basis->v00;

    if (n && (basis->v20 * n->x + basis->v21 * n->y + basis->v22 * n->z < 0.0f))
    {
        basis->v20 = -basis->v20;
        basis->v21 = -basis->v21;
        basis->v22 = -basis->v22;
    }
}

void jit_vec3_mult_mat4_pos(t_jit_vec3 *u,  t_jit_vec3 *v,  t_jit_mat4 *m)
{
    float oow;
    float div = v->x * m->v03 + v->y * m->v13 + v->z * m->v23 + m->v33;
    if (div < JIT_MATH_F32_EPS && div > -JIT_MATH_F32_EPS)
        oow = 1.0f ;
    else
        oow = 1.0f / div;

    u->x = (m->v00 * v->x + m->v10 * v->y + m->v20 * v->z + m->v30) * oow;
    u->y = (m->v01 * v->x + m->v11 * v->y + m->v21 * v->z + m->v31) * oow;
    u->z = (m->v02 * v->x + m->v12 * v->y + m->v22 * v->z + m->v32) * oow;
}

void jit_mat4_mult_vec3_pos(t_jit_vec3 *u,  t_jit_mat4 *m,  t_jit_vec3 *v)
{
    float w;
    float div = v->x * m->v30 + v->y * m->v31 + v->z * m->v32 + m->v33;
    
    if (div < JIT_MATH_F32_EPS && div > -JIT_MATH_F32_EPS)
        w = 1.0f ;
    else
        w = 1.0f / div;
    
    u->x = (m->v00 * v->x + m->v01 * v->y + m->v02 * v->z + m->v03) * w;
    u->y = (m->v10 * v->x + m->v11 * v->y + m->v12 * v->z + m->v13) * w;
    u->z = (m->v20 * v->x + m->v21 * v->y + m->v22 * v->z + m->v23) * w;
}

void jit_vec3_mult_mat4_dir(t_jit_vec3 *u,  t_jit_vec3 *v,  t_jit_mat4 *m)
{
    u->x = m->v00 * v->x + m->v10 * v->y + m->v20 * v->z;
    u->y = m->v01 * v->x + m->v11 * v->y + m->v21 * v->z;
    u->z = m->v02 * v->x + m->v12 * v->y + m->v22 * v->z;
}

void jit_mat4_mult_vec3_dir(t_jit_vec3 *u,  t_jit_mat4 *m,  t_jit_vec3 *v)
{
    u->x = m->v00 * v->x + m->v01 * v->y + m->v02 * v->z;
    u->y = m->v10 * v->x + m->v11 * v->y + m->v12 * v->z;
    u->z = m->v20 * v->x + m->v21 * v->y + m->v22 * v->z;
}

void jit_vec3_mult_mat4(t_jit_vec3 *u,  t_jit_vec3 *v,  t_jit_mat4 *m)
{
    u->x = m->v00 * v->x + m->v10 * v->y + m->v20 * v->z + m->v30;
    u->y = m->v01 * v->x + m->v11 * v->y + m->v21 * v->z + m->v31;
    u->z = m->v02 * v->x + m->v12 * v->y + m->v22 * v->z + m->v32;
}

void jit_mat4_mult_vec3(t_jit_vec3 *u,  t_jit_mat4 *m,  t_jit_vec3 *v)
{
    u->x = m->v00 * v->x + m->v01 * v->y + m->v02 * v->z + m->v03;
    u->y = m->v10 * v->x + m->v11 * v->y + m->v12 * v->z + m->v13;
    u->z = m->v20 * v->x + m->v21 * v->y + m->v22 * v->z + m->v23;
}

// ----------------------------------------------------------------------------

void jit_mat4_assign( t_jit_mat4 *m, t_jit_mat4 *a )
{
    // avoid self assignment
    if(m == a)
        return;
        
    m->v00 = a->v00; m->v10 = a->v10; m->v20 = a->v20; m->v30 = a->v30;
    m->v01 = a->v01; m->v11 = a->v11; m->v21 = a->v21; m->v31 = a->v31;
    m->v02 = a->v02; m->v12 = a->v12; m->v22 = a->v22; m->v32 = a->v32;
    m->v03 = a->v03; m->v13 = a->v13; m->v23 = a->v23; m->v33 = a->v33;
}

void jit_mat4_from_array( t_jit_mat4 *m, float *v )
{
    m->v00 = v[0];  m->v10 = v[1];  m->v20 = v[2];  m->v30 = v[3]; 
    m->v01 = v[4];  m->v11 = v[5];  m->v21 = v[6];  m->v31 = v[7]; 
    m->v02 = v[8];  m->v12 = v[9];  m->v22 = v[10]; m->v32 = v[11]; 
    m->v03 = v[12]; m->v13 = v[13]; m->v23 = v[14]; m->v33 = v[15]; 
}

void jit_mat4_from_array4x4( t_jit_mat4 *m, float **v )
{
    m->v00 = v[0][0];   m->v10 = v[1][0];   m->v20 = v[2][0];   m->v30 = v[3][0]; 
    m->v01 = v[0][1];   m->v11 = v[1][1];   m->v21 = v[2][1];   m->v31 = v[3][1]; 
    m->v02 = v[0][2];   m->v12 = v[1][2];   m->v22 = v[2][2];   m->v32 = v[3][2]; 
    m->v03 = v[0][3];   m->v13 = v[1][3];   m->v23 = v[2][3];   m->v33 = v[3][3]; 
}

void jit_mat4_from_coords( t_jit_mat4 *m, 
    float f0,   float f1,   float f2,   float f3,
    float f4,   float f5,   float f6,   float f7,
    float f8,   float f9,   float f10,  float f11,
    float f12,  float f13,  float f14,  float f15 )
{
    m->v00 = f0;    m->v10 = f1;    m->v20 = f2;    m->v30 = f3; 
    m->v01 = f4;    m->v11 = f5;    m->v21 = f6;    m->v31 = f7; 
    m->v02 = f8;    m->v12 = f9;    m->v22 = f10;   m->v32 = f11; 
    m->v03 = f12;   m->v13 = f13;   m->v23 = f14;   m->v33 = f15; 
}

void jit_mat4_to_array( t_jit_mat4 *m, float *v )
{
    v[0] = m->v00;  v[1] = m->v10;  v[2] = m->v20;  v[3] = m->v30; 
    v[4] = m->v01;  v[5] = m->v11;  v[6] = m->v21;  v[7] = m->v31; 
    v[8] = m->v02;  v[9] = m->v12;  v[10] = m->v22; v[11] = m->v32; 
    v[12] = m->v03; v[13] = m->v13; v[14] = m->v23; v[15] = m->v33; 
}

void jit_mat4_to_array4x4( t_jit_mat4 *m, float **v )
{
    v[0][0] = m->v00;   v[1][0] = m->v10;   v[2][0] = m->v20;   v[3][0] = m->v30; 
    v[0][1] = m->v01;   v[1][1] = m->v11;   v[2][1] = m->v21;   v[3][1] = m->v31; 
    v[0][2] = m->v02;   v[1][2] = m->v12;   v[2][2] = m->v22;   v[3][2] = m->v32; 
    v[0][3] = m->v03;   v[1][3] = m->v13;   v[2][3] = m->v23;   v[3][3] = m->v33; 
}

void jit_mat4_to_coords( t_jit_mat4 *m, 
    float *f0,   float *f1,   float *f2,   float *f3,
    float *f4,   float *f5,   float *f6,   float *f7,
    float *f8,   float *f9,   float *f10,  float *f11,
    float *f12,  float *f13,  float *f14,  float *f15 )
{
    *f0 = m->v00;   *f1 = m->v10;   *f2 = m->v20;   *f3 = m->v30; 
    *f4 = m->v01;   *f5 = m->v11;   *f6 = m->v21;   *f7 = m->v31; 
    *f8 = m->v02;   *f9 = m->v12;   *f10 = m->v22;  *f11 = m->v32; 
    *f12 = m->v03;  *f13 = m->v13;  *f14 = m->v23;  *f15 = m->v33; 
}
 
void jit_mat4_get_row( t_jit_vec4 *v, long i, t_jit_mat4 *m )
{
    if( i >= 0 && i < 4 ) 
    {
        v->x = m->vals[i + 0];
        v->y = m->vals[i + 4];
        v->z = m->vals[i + 8];
        v->w = m->vals[i + 12];
    }
}
 
void jit_mat4_set_row( t_jit_mat4 *m, long i,  t_jit_vec4 *v )
{
    if( i >= 0 && i < 4 ) 
    {
        m->vals[i + 0] = v->x;
        m->vals[i + 4] = v->y;
        m->vals[i + 8] = v->z;
        m->vals[i + 12] = v->w;
    }
}

void jit_mat4_get_col( t_jit_vec4 *v, long i, t_jit_mat4 *m )
{
    if( i >= 0 && i < 4 ) 
    {
        v->x = m->vals[i * 4 + 0];
        v->y = m->vals[i * 4 + 1];
        v->z = m->vals[i * 4 + 2];
        v->w = m->vals[i * 4 + 3];
    }
}

void jit_mat4_set_col( t_jit_mat4 *m, long i,  t_jit_vec4 *v)
{
    if( i >= 0 && i < 4 ) 
    {
        m->vals[i * 4 + 0] = v->x;
        m->vals[i * 4 + 1] = v->y;
        m->vals[i * 4 + 2] = v->z;
        m->vals[i * 4 + 3] = v->w;
    }
}

float jit_mat4_get_coord( t_jit_mat4 *m, long r, long c )
{
    if(c >= 0 && c < 4 && r >= 0 && r < 4)
        return m->vals[ c * 4 + r ];
    else
        return 0.0f;
}

void jit_mat4_mult_vec4(t_jit_vec4 *u,  t_jit_mat4 *m,  t_jit_vec4 *v)
{
    u->x = m->v00 * v->x + m->v01 * v->y + m->v02 * v->z + m->v03 * v->w;
    u->y = m->v10 * v->x + m->v11 * v->y + m->v12 * v->z + m->v13 * v->w;
    u->z = m->v20 * v->x + m->v21 * v->y + m->v22 * v->z + m->v23 * v->w;
    u->w = m->v30 * v->x + m->v31 * v->y + m->v32 * v->z + m->v33 * v->w;
}

void jit_vec4_mult_mat4(t_jit_vec4 *u,  t_jit_vec4 *v,  t_jit_mat4 *m)
{
    u->x = m->v00 * v->x + m->v10 * v->y + m->v20 * v->z + m->v30 * v->w;
    u->y = m->v01 * v->x + m->v11 * v->y + m->v21 * v->z + m->v31 * v->w;
    u->z = m->v02 * v->x + m->v12 * v->y + m->v22 * v->z + m->v32 * v->w;
    u->w = m->v03 * v->x + m->v13 * v->y + m->v23 * v->z + m->v33 * v->w;
}

void jit_mat4_accum_add(t_jit_mat4 *a,  t_jit_mat4 *b)
{
    a->v00 += b->v00;
    a->v10 += b->v10;
    a->v20 += b->v20;
    a->v30 += b->v30;
    a->v01 += b->v01;
    a->v11 += b->v11;
    a->v21 += b->v21;
    a->v31 += b->v31;
    a->v02 += b->v02;
    a->v12 += b->v12;
    a->v22 += b->v22;
    a->v32 += b->v32;
    a->v03 += b->v03;
    a->v13 += b->v13;
    a->v23 += b->v23;
    a->v33 += b->v33;
}

void jit_mat4_add(t_jit_mat4 *c, t_jit_mat4 *a, t_jit_mat4 *b)
{
    if (c == a || c == b)
    {
        t_jit_mat4 m;

        m.v00 = a->v00 + b->v00;
        m.v01 = a->v01 + b->v01;
        m.v02 = a->v02 + b->v02;
        m.v03 = a->v03 + b->v03;
        m.v10 = a->v10 + b->v10;
        m.v11 = a->v11 + b->v11;
        m.v12 = a->v12 + b->v12;
        m.v13 = a->v13 + b->v13;
        m.v20 = a->v20 + b->v20;
        m.v21 = a->v21 + b->v21;
        m.v22 = a->v22 + b->v22;
        m.v23 = a->v23 + b->v23;
        m.v30 = a->v30 + b->v30;
        m.v31 = a->v31 + b->v31;
        m.v32 = a->v32 + b->v32;
        m.v33 = a->v33 + b->v33;
    
        jit_mat4_assign(c, &m);
    }
    else
    {
        c->v00 = a->v00 + b->v00;
        c->v01 = a->v01 + b->v01;
        c->v02 = a->v02 + b->v02;
        c->v03 = a->v03 + b->v03;
        c->v10 = a->v10 + b->v10;
        c->v11 = a->v11 + b->v11;
        c->v12 = a->v12 + b->v12;
        c->v13 = a->v13 + b->v13;
        c->v20 = a->v20 + b->v20;
        c->v21 = a->v21 + b->v21;
        c->v22 = a->v22 + b->v22;
        c->v23 = a->v23 + b->v23;
        c->v30 = a->v30 + b->v30;
        c->v31 = a->v31 + b->v31;
        c->v32 = a->v32 + b->v32;
        c->v33 = a->v33 + b->v33;
    }
}

void jit_mat4_mult(t_jit_mat4 *c,  t_jit_mat4 *a,  t_jit_mat4 *b)
{
    if (c == a || c == b)
    {
        t_jit_mat4 m;

        m.v00 = a->v00 * b->v00 + a->v01 * b->v10 + a->v02 * b->v20 + a->v03 * b->v30;
        m.v10 = a->v10 * b->v00 + a->v11 * b->v10 + a->v12 * b->v20 + a->v13 * b->v30;
        m.v20 = a->v20 * b->v00 + a->v21 * b->v10 + a->v22 * b->v20 + a->v23 * b->v30;
        m.v30 = a->v30 * b->v00 + a->v31 * b->v10 + a->v32 * b->v20 + a->v33 * b->v30;
        m.v01 = a->v00 * b->v01 + a->v01 * b->v11 + a->v02 * b->v21 + a->v03 * b->v31;
        m.v11 = a->v10 * b->v01 + a->v11 * b->v11 + a->v12 * b->v21 + a->v13 * b->v31;
        m.v21 = a->v20 * b->v01 + a->v21 * b->v11 + a->v22 * b->v21 + a->v23 * b->v31;
        m.v31 = a->v30 * b->v01 + a->v31 * b->v11 + a->v32 * b->v21 + a->v33 * b->v31;
        m.v02 = a->v00 * b->v02 + a->v01 * b->v12 + a->v02 * b->v22 + a->v03 * b->v32;
        m.v12 = a->v10 * b->v02 + a->v11 * b->v12 + a->v12 * b->v22 + a->v13 * b->v32;
        m.v22 = a->v20 * b->v02 + a->v21 * b->v12 + a->v22 * b->v22 + a->v23 * b->v32;
        m.v32 = a->v30 * b->v02 + a->v31 * b->v12 + a->v32 * b->v22 + a->v33 * b->v32;
        m.v03 = a->v00 * b->v03 + a->v01 * b->v13 + a->v02 * b->v23 + a->v03 * b->v33;
        m.v13 = a->v10 * b->v03 + a->v11 * b->v13 + a->v12 * b->v23 + a->v13 * b->v33;
        m.v23 = a->v20 * b->v03 + a->v21 * b->v13 + a->v22 * b->v23 + a->v23 * b->v33;
        m.v33 = a->v30 * b->v03 + a->v31 * b->v13 + a->v32 * b->v23 + a->v33 * b->v33;
    
        jit_mat4_assign(c, &m);
    }
    else
    {
        c->v00 = a->v00 * b->v00 + a->v01 * b->v10 + a->v02 * b->v20 + a->v03 * b->v30;
        c->v10 = a->v10 * b->v00 + a->v11 * b->v10 + a->v12 * b->v20 + a->v13 * b->v30;
        c->v20 = a->v20 * b->v00 + a->v21 * b->v10 + a->v22 * b->v20 + a->v23 * b->v30;
        c->v30 = a->v30 * b->v00 + a->v31 * b->v10 + a->v32 * b->v20 + a->v33 * b->v30;
        c->v01 = a->v00 * b->v01 + a->v01 * b->v11 + a->v02 * b->v21 + a->v03 * b->v31;
        c->v11 = a->v10 * b->v01 + a->v11 * b->v11 + a->v12 * b->v21 + a->v13 * b->v31;
        c->v21 = a->v20 * b->v01 + a->v21 * b->v11 + a->v22 * b->v21 + a->v23 * b->v31;
        c->v31 = a->v30 * b->v01 + a->v31 * b->v11 + a->v32 * b->v21 + a->v33 * b->v31;
        c->v02 = a->v00 * b->v02 + a->v01 * b->v12 + a->v02 * b->v22 + a->v03 * b->v32;
        c->v12 = a->v10 * b->v02 + a->v11 * b->v12 + a->v12 * b->v22 + a->v13 * b->v32;
        c->v22 = a->v20 * b->v02 + a->v21 * b->v12 + a->v22 * b->v22 + a->v23 * b->v32;
        c->v32 = a->v30 * b->v02 + a->v31 * b->v12 + a->v32 * b->v22 + a->v33 * b->v32;
        c->v03 = a->v00 * b->v03 + a->v01 * b->v13 + a->v02 * b->v23 + a->v03 * b->v33;
        c->v13 = a->v10 * b->v03 + a->v11 * b->v13 + a->v12 * b->v23 + a->v13 * b->v33;
        c->v23 = a->v20 * b->v03 + a->v21 * b->v13 + a->v22 * b->v23 + a->v23 * b->v33;
        c->v33 = a->v30 * b->v03 + a->v31 * b->v13 + a->v32 * b->v23 + a->v33 * b->v33;
    }
}

void jit_mat4_transpose(t_jit_mat4 *m)
{
    float t;
    t = m->v01;
    m->v01 = m->v10;
    m->v10 = t;

    t = m->v02;
    m->v02 = m->v20;
    m->v20 = t;

    t = m->v03;
    m->v03 = m->v30;
    m->v30 = t;

    t = m->v12;
    m->v12 = m->v21;
    m->v21 = t;

    t = m->v13;
    m->v13 = m->v31;
    m->v31 = t;

    t = m->v23;
    m->v23 = m->v32;
    m->v32 = t;
}

void jit_mat4_get_transpose(t_jit_mat4 *b,  t_jit_mat4 *a)
{
    b->v00 = a->v00;
    b->v01 = a->v10;
    b->v02 = a->v20;
    b->v03 = a->v30;
    b->v10 = a->v01;
    b->v11 = a->v11;
    b->v12 = a->v21;
    b->v13 = a->v31;
    b->v20 = a->v02;
    b->v21 = a->v12;
    b->v22 = a->v22;
    b->v23 = a->v32;
    b->v30 = a->v03;
    b->v31 = a->v13;
    b->v32 = a->v23;
    b->v33 = a->v33;
}

void jit_mat4_get_inverse(t_jit_mat4 *b,  t_jit_mat4 *a)
{
    float det,oodet;

    b->v00 =  JIT_MATH_DET3x3(a->v11, a->v21, a->v31, a->v12, a->v22, a->v32, a->v13, a->v23, a->v33);
    b->v10 = -JIT_MATH_DET3x3(a->v10, a->v20, a->v30, a->v12, a->v22, a->v32, a->v13, a->v23, a->v33);
    b->v20 =  JIT_MATH_DET3x3(a->v10, a->v20, a->v30, a->v11, a->v21, a->v31, a->v13, a->v23, a->v33);
    b->v30 = -JIT_MATH_DET3x3(a->v10, a->v20, a->v30, a->v11, a->v21, a->v31, a->v12, a->v22, a->v32);

    b->v01 = -JIT_MATH_DET3x3(a->v01, a->v21, a->v31, a->v02, a->v22, a->v32, a->v03, a->v23, a->v33);
    b->v11 =  JIT_MATH_DET3x3(a->v00, a->v20, a->v30, a->v02, a->v22, a->v32, a->v03, a->v23, a->v33);
    b->v21 = -JIT_MATH_DET3x3(a->v00, a->v20, a->v30, a->v01, a->v21, a->v31, a->v03, a->v23, a->v33);
    b->v31 =  JIT_MATH_DET3x3(a->v00, a->v20, a->v30, a->v01, a->v21, a->v31, a->v02, a->v22, a->v32);

    b->v02 =  JIT_MATH_DET3x3(a->v01, a->v11, a->v31, a->v02, a->v12, a->v32, a->v03, a->v13, a->v33);
    b->v12 = -JIT_MATH_DET3x3(a->v00, a->v10, a->v30, a->v02, a->v12, a->v32, a->v03, a->v13, a->v33);
    b->v22 =  JIT_MATH_DET3x3(a->v00, a->v10, a->v30, a->v01, a->v11, a->v31, a->v03, a->v13, a->v33);
    b->v32 = -JIT_MATH_DET3x3(a->v00, a->v10, a->v30, a->v01, a->v11, a->v31, a->v02, a->v12, a->v32);

    b->v03 = -JIT_MATH_DET3x3(a->v01, a->v11, a->v21, a->v02, a->v12, a->v22, a->v03, a->v13, a->v23);
    b->v13 =  JIT_MATH_DET3x3(a->v00, a->v10, a->v20, a->v02, a->v12, a->v22, a->v03, a->v13, a->v23);
    b->v23 = -JIT_MATH_DET3x3(a->v00, a->v10, a->v20, a->v01, a->v11, a->v21, a->v03, a->v13, a->v23);
    b->v33 =  JIT_MATH_DET3x3(a->v00, a->v10, a->v20, a->v01, a->v11, a->v21, a->v02, a->v12, a->v22);

    det = (a->v00 * b->v00) + (a->v01 * b->v10) + (a->v02 * b->v20) + (a->v03 * b->v30);
    
    // avoid divide by zero
    if(det < JIT_MATH_F32_EPS)
        det = JIT_MATH_F32_EPS;

    oodet = 1.0f / det;

    b->v00 *= oodet;
    b->v10 *= oodet;
    b->v20 *= oodet;
    b->v30 *= oodet;

    b->v01 *= oodet;
    b->v11 *= oodet;
    b->v21 *= oodet;
    b->v31 *= oodet;

    b->v02 *= oodet;
    b->v12 *= oodet;
    b->v22 *= oodet;
    b->v32 *= oodet;

    b->v03 *= oodet;
    b->v13 *= oodet;
    b->v23 *= oodet;
    b->v33 *= oodet;
}

void jit_mat4_get_inverse_rot_trans(t_jit_mat4 *b,  t_jit_mat4 *a)
{
    b->v00 = a->v00;
    b->v10 = a->v01;
    b->v20 = a->v02;
    b->v30 = a->v30;
    b->v01 = a->v10;
    b->v11 = a->v11;
    b->v21 = a->v12;
    b->v31 = a->v31;
    b->v02 = a->v20;
    b->v12 = a->v21;
    b->v22 = a->v22;
    b->v32 = a->v32;
    b->v03 = - (a->v00 * a->v03 + a->v10 * a->v13 + a->v20 * a->v23);
    b->v13 = - (a->v01 * a->v03 + a->v11 * a->v13 + a->v21 * a->v23);
    b->v23 = - (a->v02 * a->v03 + a->v12 * a->v13 + a->v22 * a->v23);
    b->v33 = a->v33;
}

void jit_mat4_look_at(t_jit_mat4 *m,  t_jit_vec3 *eye,  t_jit_vec3 *center,  t_jit_vec3 *up)
{
    t_jit_vec3 x, y, z;

    // z vector
    z.x = eye->x - center->x;
    z.y = eye->y - center->y;
    z.z = eye->z - center->z;
    jit_vec3_normalize(&z);

    // y vector
    y.x = up->x;
    y.y = up->y;
    y.z = up->z;

    // x vector = u cross z
    jit_vec3_cross(&x, &y, &z);

    // recompute u = z cross x
    jit_vec3_cross(&y, &z, &x);
    
    // renormalize x + y
    jit_vec3_normalize(&x);
    jit_vec3_normalize(&y);

    m->v00 = x.x; m->v01 = x.y; m->v02 = x.z; m->v03 = -x.x * eye->x - x.y * eye->y - x.z*eye->z;
    m->v10 = y.x; m->v11 = y.y; m->v12 = y.z; m->v13 = -y.x * eye->x - y.y * eye->y - y.z*eye->z;
    m->v20 = z.x; m->v21 = z.y; m->v22 = z.z; m->v23 = -z.x * eye->x - z.y * eye->y - z.z*eye->z;
    m->v30 = 0.0f; m->v31 = 0.0f; m->v32 = 0.0f; m->v33 = 1.0f;
}

void jit_mat4_frustum(t_jit_mat4 *m,  
    float left,  float right,  
    float bottom, float top,  
    float n,  float f)
{
    m->v00 = (2.0f*n) / (right-left);
    m->v10 = 0.0;
    m->v20 = 0.0;
    m->v30 = 0.0;

    m->v01 = 0.0;
    m->v11 = (2.0f*n) / (top-bottom);
    m->v21 = 0.0;
    m->v31 = 0.0;

    m->v02 = (right+left) / (right-left);
    m->v12 = (top+bottom) / (top-bottom);
    m->v22 = -(f+n) / (f-n);
    m->v32 = -1.0f;

    m->v03 = 0.0;
    m->v13 = 0.0;
    m->v23 = -(2.0f*f*n) / (f-n);
    m->v33 = 0.0;
}

void jit_mat4_perspective(t_jit_mat4 *m,  float fovy,  float aspect,  float n,  float f)
{
    float xmin, xmax, ymin, ymax;

    ymax = n * jit_math_tan(fovy * JIT_MATH_F32_DEGTORAD * 0.5f);
    ymin = -ymax;

    xmin = ymin * aspect;
    xmax = ymax * aspect;

    jit_mat4_frustum(m, xmin, xmax, ymin, ymax, n, f);
}

void jit_mat4_ortho(t_jit_mat4 *m,  
    float left, float right, 
    float bottom, float top,
    float n, float f)
{
    m->v00 = 2.0f / (right - left);
    m->v01 = 0.0f;
    m->v02 = 0.0f;
    m->v03 = - (right + left) / (right - left);
    m->v10 = 0.0f;
    m->v11 = 2.0f / (top - bottom);
    m->v12 = 0.0f;
    m->v13 = - (top + bottom) / (top - bottom);
    m->v20 = 0.0f;
    m->v21 = 0.0f;
    m->v22 = - 2.0f / (f - n);
    m->v23 = - (f + n) / (f - n);
    m->v30 = 0.0f;
    m->v31 = 0.0f;
    m->v32 = 0.0f;
    m->v33 = 1.0f;
}

void jit_mat4_from_quat( t_jit_mat4 *m, t_jit_quat * q)
{
    t_jit_mat3 r;
    jit_mat3_from_quat(&r, q);
    jit_mat4_from_mat3(m, &r);
}

void jit_mat4_from_axisangle( t_jit_mat4 *m,  t_jit_vec3 *v, float theta) 
{
    float ct = (float)jit_math_cos(theta);
    float st = (float)jit_math_sin(theta);

    float xx = v->x * v->x;
    float yy = v->y * v->y;
    float zz = v->z * v->z;
    float xy = v->x * v->y;
    float xz = v->x * v->z;
    float yz = v->y * v->z;

    m->v00 = xx + ct*(1-xx);
    m->v01 = xy + ct*(-xy) + st*-v->z;
    m->v02 = xz + ct*(-xz) + st*v->y;

    m->v10 = xy + ct*(-xy) + st*v->z;
    m->v11 = yy + ct*(1-yy);
    m->v12 = yz + ct*(-yz) + st*-v->x;

    m->v20 = xz + ct*(-xz) + st*-v->y;
    m->v21 = yz + ct*(-yz) + st*v->x;
    m->v22 = zz + ct*(1-zz);
}

void jit_mat4_from_uv( t_jit_mat4 *m, t_jit_vec3 *u,  t_jit_vec3 *v)
{
    float phi;
    float h, hxy, hxz, hyz;
    float lambda;
    t_jit_vec3 w;

    jit_vec3_cross(&w, u, v);
    phi = jit_vec3_dot(u, v);
    lambda = jit_vec3_dot(&w, &w);
    
    if (lambda > JIT_MATH_F32_EPS)
        h = (1.0f - phi) / lambda;
    else
        h = lambda;
    
    hxy = w.x * w.y * h;
    hxz = w.x * w.z * h;
    hyz = w.y * w.z * h;

    m->v00 = phi + w.x * w.x * h;
    m->v01 = hxy - w.z;
    m->v02 = hxz + w.y;

    m->v10 = hxy + w.z;
    m->v11 = phi + w.y * w.y * h;
    m->v12 = hyz - w.x;

    m->v20 = hxz - w.y;
    m->v21 = hyz + w.x;
    m->v22 = phi + w.z * w.z * h;
}

void jit_mat4_from_mat3( t_jit_mat4 *m, t_jit_mat3 *r)
{
    // copy the 3x3 rotation block
    m->v00 = r->v00; m->v10 = r->v10; m->v20 = r->v20;
    m->v01 = r->v01; m->v11 = r->v11; m->v21 = r->v21;
    m->v02 = r->v02; m->v12 = r->v12; m->v22 = r->v22;
}

void jit_mat3_from_mat4(t_jit_mat3 *r, t_jit_mat4 *m) 
{
    // assign the 3x3 rotation block
    r->v00 = m->v00; r->v10 = m->v10; r->v20 = m->v20;
    r->v01 = m->v01; r->v11 = m->v11; r->v21 = m->v21;
    r->v02 = m->v02; r->v12 = m->v12; r->v22 = m->v22;
}

void jit_mat4_set_scale( t_jit_mat4 *m, t_jit_vec3 *s)
{
    m->v00 = s->x;
    m->v11 = s->y;
    m->v22 = s->z;
}

void jit_mat4_get_scale(t_jit_vec3 *s, t_jit_mat4 *m) 
{
    s->x = m->v00;
    s->y = m->v11;
    s->z = m->v22;
}

void jit_mat4_set_translation(t_jit_mat4 *m, t_jit_vec3 *t)
{
    m->v03 = t->x;
    m->v13 = t->y;
    m->v23 = t->z;
}

void jit_mat4_get_translation(t_jit_vec3 *t, t_jit_mat4 *m) 
{
    t->x = m->v03;
    t->y = m->v13;
    t->z = m->v23;
}

void jit_mat4_negate(t_jit_mat4 *m)
{
    m->v00 = -m->v00;
    m->v10 = -m->v10;
    m->v20 = -m->v20;
    m->v30 = -m->v30;
    m->v01 = -m->v01;
    m->v11 = -m->v11;
    m->v21 = -m->v21;
    m->v31 = -m->v31;
    m->v02 = -m->v02;
    m->v12 = -m->v12;
    m->v22 = -m->v22;
    m->v32 = -m->v32;
    m->v03 = -m->v03;
    m->v13 = -m->v13;
    m->v23 = -m->v23;
    m->v33 = -m->v33;
}

float jit_tri_area( t_jit_vec3 *v1,  t_jit_vec3 *v2,  t_jit_vec3 *v3)
{
    t_jit_vec3 sum;
    t_jit_vec3 cp;
    jit_vec3_cross(&sum, v1, v2);
    jit_vec3_cross(&cp, v2, v3);
    jit_vec3_accum_add(&sum, &cp);
    jit_vec3_cross(&cp, v3, v1);
    jit_vec3_accum_add(&sum, &cp);
    return 0.5f * jit_vec3_mag(&sum); 
}

float jit_tri_perimeter( t_jit_vec3 *v1,  t_jit_vec3 *v2,  t_jit_vec3 *v3)
{
    float p;
    t_jit_vec3 d;
    jit_vec3_sub(&d, v1, v2);
    p = jit_vec3_mag(&d);
    jit_vec3_sub(&d, v2, v3);
    p += jit_vec3_mag(&d);
    jit_vec3_sub(&d, v3, v1);
    p += jit_vec3_mag(&d);
    return p;
}

float jit_tri_find_in_circle(t_jit_vec3 *center,  t_jit_vec3 *v1,  t_jit_vec3 *v2,  t_jit_vec3 *v3)
{
    t_jit_vec3 d;
    float s;
    float perim;
    float area = jit_tri_area(v1, v2, v3);
    
    // if the area is null
    if (area < JIT_MATH_F32_EPS)
    {
        center = v1;
        return 0.0f;
    }

    perim = 1.0f / jit_tri_perimeter(v1, v2, v3);

    jit_vec3_sub(&d, v2, v3);
    s = jit_vec3_mag(&d);
    jit_vec3_scale(center, v1, s);

    jit_vec3_sub(&d, v3, v1);
    s = jit_vec3_mag(&d);
    jit_vec3_madd(center, v2, s);
    
    jit_vec3_sub(&d, v1, v2);
    s = jit_vec3_mag(&d);
    jit_vec3_madd(center, v3, s);

    jit_vec3_accum_scale(center, perim);
    return 2.0f * area * perim;
}

float jit_tri_find_circ_circle( t_jit_vec3 *center,  t_jit_vec3 *v1,  t_jit_vec3 *v2,  t_jit_vec3 *v3)
{
    t_jit_vec3 e0;
    t_jit_vec3 e1;
    float d1, d2, d3;
    float c1, c2, c3, c;

    jit_vec3_sub(&e0, v3, v1);
    jit_vec3_sub(&e1, v2, v1);
    d1 = jit_vec3_dot(&e0, &e1);

    jit_vec3_sub(&e0, v3, v2);
    jit_vec3_sub(&e1, v1, v2);
    d2 = jit_vec3_dot(&e0, &e1);

    jit_vec3_sub(&e0, v1, v3);
    jit_vec3_sub(&e1, v2, v3);
    d3 = jit_vec3_dot(&e0, &e1);

    c1 = d2 * d3;
    c2 = d3 * d1;
    c3 = d1 * d2;
    c = 1.0f / (c1 + c2 + c3);

    jit_vec3_scale(center, v1, c2 + c3);
    jit_vec3_madd(center, v2, c3 + c1);
    jit_vec3_madd(center, v3, c1 + c2);
    
    c *= 0.5f;
    center->x *= c;
    center->y *= c;
    center->z *= c;
 
    return 0.5f * jit_math_sqrt((d1 + d2) * (d2 + d3) * (d3 + d1) * 2 * c);
}

// ----------------------------------------------------------------------------

Copyright © 2008, Cycling '74