Max 5 API Reference

jit.noise.simplex.c

/*
* Copyright 2001-2005 - Cycling '74
* Derek Gerstmann - derek@cycling74.com
*
* Coherent simplex based noise functor.
*
* Based on Ken Perlin's replacement for his classic noise algorithm,
* presented at GDC 2001.
*
* "In one dimension, the simplex is the line segment. In two dimensions, 
* the simplex is the convex hull of the equilateral triangle. In three 
* dimensions, the simplex is the convex hull of the tetrahedron. The 
* simplex in four dimensions (the pentatope) is a regular tetrahedron 
* ABCD in which a point E along the fourth dimension through the center 
* of ABCD is chosen so that EA = EB = EC = ED = AB."
*   -- http://mathworld.wolfram.com/Simplex.html
*/

/*************************************************************************/

#include "jit.noise.h"
#include "jit.common.h"
#include "jit.functor.h"
#include "ext_systhread.h"

/*************************************************************************/

// 3d simplex gradient directions
static long jit_functor_noise_simplex_grad3count = 12;
static long jit_functor_noise_simplex_grad3[12][3] =
{
    {0,1,1},{0,1,-1},{0,-1,1},{0,-1,-1},
    {1,0,1},{1,0,-1},{-1,0,1},{-1,0,-1},
    {1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0}
};

// 4d simplex gradient directions
static long jit_functor_noise_simplex_grad4count = 32;
static long jit_functor_noise_simplex_grad4[32][4] =
{
    {0, 1, 1, 1}, {0, 1, 1, -1}, {0, 1, -1, 1}, {0, 1, -1, -1},
    {0, -1, 1, 1}, {0, -1, 1, -1}, {0, -1, -1, 1}, {0, -1, -1, -1},
    {1, 0, 1, 1}, {1, 0, 1, -1}, {1, 0, -1, 1}, {1, 0, -1, -1},
    { -1, 0, 1, 1}, { -1, 0, 1, -1}, { -1, 0, -1, 1}, { -1, 0, -1, -1},
    {1, 1, 0, 1}, {1, 1, 0, -1}, {1, -1, 0, 1}, {1, -1, 0, -1},
    { -1, 1, 0, 1}, { -1, 1, 0, -1}, { -1, -1, 0, 1}, { -1, -1, 0, -1},
    {1, 1, 1, 0}, {1, 1, -1, 0}, {1, -1, 1, 0}, {1, -1, -1, 0},
    { -1, 1, 1, 0}, { -1, 1, -1, 0}, { -1, -1, 1, 0}, { -1, -1, -1, 0}
};

// 4d simplex index table
static long jit_functor_noise_simplex_index[64][4] = 
{
    {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, 
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
    {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, 
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
    {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, 
    {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
    {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, 
    {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 
    {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
    {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 
    {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
    {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, 
    {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}
};

// skewing + unskewing factors
#define F2 ((jit_math_sqrt(3.0)-1.0)/2.0)
#define G2 ((3.0-jit_math_sqrt(3.0))/6.0)

#define F3 (1.0/3.0)
#define G3 (1.0/6.0)

#define F4 ((jit_math_sqrt(5.0)-1.0)/4.0)
#define G4 ((5.0-jit_math_sqrt(5.0))/20.0)

// safe permutation lookup
#define P(i) (x->p[(i)&(x->tablemask)])

/*************************************************************************/

// the noise object
typedef struct _jit_functor_noise_simplex
{
    t_jit_object    ob;

    // pseudo random permutation table
    long    *p;

    // lattice of uniform values [0-1]
    t_jit_functor_combined_dynarray lattice;

    // array of simplex directions
    t_jit_functor_combined_dimvalue *gradients;

    // skew and unskew factors for determining simplex
    t_jit_functor_combined_value    skew;
    t_jit_functor_combined_value    unskew;

    // current dimcount for simplex directions
    long    dimcount;

    // size and mask for tables
    long    tablesize;
    long    tablemask;

    // seed for random number generation
    long    seed;

    // signed mode (0 = unsigned, 1 = signed)
    long    sign;

    // abs mode (0 = normal, 1 = absolute)
    long    abs;

    // lock on our dir init which can be called during parallel evaluation
    long spinwait;

}
t_jit_functor_noise_simplex;

t_class *_jit_functor_noise_simplex_class;
t_symbol *ps_jit_functor_noise_simplex_filter;

/*************************************************************************/

t_jit_object *jit_functor_noise_simplex_new(void);
t_jit_err jit_functor_noise_simplex_free(t_jit_functor_noise_simplex *x);

t_jit_err jit_functor_noise_simplex_abs(t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv);
t_jit_err jit_functor_noise_simplex_sign(t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv);
t_jit_err jit_functor_noise_simplex_seed(t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv);
t_jit_err jit_functor_noise_simplex_tablesize(t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv);
t_jit_err jit_functor_noise_simplex_recalc(t_jit_functor_noise_simplex *x);

double jit_functor_noise_simplex_eval_float64(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals);

double jit_functor_noise_simplex_eval_float64_dim2(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals);

double jit_functor_noise_simplex_eval_float64_dim3(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals);

double jit_functor_noise_simplex_eval_float64_dim4(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals);

float jit_functor_noise_simplex_eval_float32(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals);

float jit_functor_noise_simplex_eval_float32_dim2(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals);

float jit_functor_noise_simplex_eval_float32_dim3(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals);

float jit_functor_noise_simplex_eval_float32_dim4(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals);

long jit_functor_noise_simplex_eval_fixed(
    t_jit_functor_noise_simplex *x, long dimcount, long *vals);

t_jit_err jit_functor_noise_simplex_dir_init(
    t_jit_functor_noise_simplex *x);

/*************************************************************************/

t_jit_err jit_functor_noise_simplex_init(void)
{
    t_jit_object *attr;

    // create class
    _jit_functor_noise_simplex_class = jit_class_new("jit_functor_noise_simplex",
                                       (method)jit_functor_noise_simplex_new, (method)jit_functor_noise_simplex_free,
                                       sizeof(t_jit_functor_noise_simplex), 0L);

    // add attribute methods
    attr = jit_object_new(_jit_sym_jit_attr_offset, "abs", _jit_sym_long, 0,
                          (method)0L, (method)jit_functor_noise_simplex_abs,
                          calcoffset(t_jit_functor_noise_simplex, abs));
    jit_class_addattr(_jit_functor_noise_simplex_class, attr);
    attr = jit_object_new(_jit_sym_jit_attr_offset, "sign", _jit_sym_long, 0,
                          (method)0L, (method)jit_functor_noise_simplex_sign,
                          calcoffset(t_jit_functor_noise_simplex, sign));
    jit_class_addattr(_jit_functor_noise_simplex_class, attr);
    attr = jit_object_new(_jit_sym_jit_attr_offset, "seed", _jit_sym_long, 0,
                          (method)0L, (method)jit_functor_noise_simplex_seed,
                          calcoffset(t_jit_functor_noise_simplex, seed));
    jit_class_addattr(_jit_functor_noise_simplex_class, attr);
    attr = jit_object_new(_jit_sym_jit_attr_offset, "tablesize", _jit_sym_long, 0,
                          (method)0L, (method)jit_functor_noise_simplex_tablesize,
                          calcoffset(t_jit_functor_noise_simplex, tablesize));
    jit_class_addattr(_jit_functor_noise_simplex_class, attr);

    // add evaluation methods
    jit_class_addmethod(_jit_functor_noise_simplex_class,
                        (method)jit_functor_noise_simplex_eval_fixed, "evalfixed", A_CANT, 0L);
    jit_class_addmethod(_jit_functor_noise_simplex_class,
                        (method)jit_functor_noise_simplex_eval_float32, "evalfloat32", A_CANT, 0L);
    jit_class_addmethod(_jit_functor_noise_simplex_class,
                        (method)jit_functor_noise_simplex_eval_float64, "evalfloat64", A_CANT, 0L);

    // register functor *then* class
    jit_functor_setup_class(_jit_functor_noise_simplex_class, "noise", "simplex");
    jit_class_register(_jit_functor_noise_simplex_class);
    return JIT_ERR_NONE;
}

t_jit_object *jit_functor_noise_simplex_new(void)
{
    long i, j;
    t_jit_functor_noise_simplex *x;

    if (x = (t_jit_functor_noise_simplex *)jit_object_alloc(_jit_functor_noise_simplex_class))
    {
        // initialization
        x->p = NULL;
        x->lattice.fixed = NULL;
        x->lattice.float32 = NULL;
        x->lattice.float64 = NULL;

        x->abs = FALSE;
        x->sign = FALSE;

        JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->skew, 0);
        JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->unskew, 0);

        x->seed = JIT_NOISE_DEFAULT_SEED;
        x->tablesize = jit_math_roundup_poweroftwo(JIT_NOISE_DEFAULT_TABLESIZE) * 2;
        x->tablemask = x->tablesize - 1;

        x->dimcount = 0;
        x->gradients = NULL;
        x->spinwait = 0;

        // allocate permutation and uniform tables
        jit_functor_noise_permutations_init(&x->p, x->tablesize, x->seed);
        jit_functor_noise_lattice_init(&x->lattice, x->tablesize, x->seed);
        jit_functor_noise_simplex_recalc(x);
    }

    return (t_jit_object *)x;
}

t_jit_err jit_functor_noise_simplex_free(t_jit_functor_noise_simplex *x)
{
    // free any allocated memory
    jit_functor_combined_dynarray_destroy(&x->lattice);

    if (x->p)
        jit_disposeptr((char*)x->p);
    x->p = NULL;

    if (x->gradients)
        jit_disposeptr((char*)x->gradients);
    x->gradients = NULL;


    return JIT_ERR_NONE;
}

t_jit_err jit_functor_noise_simplex_abs(
    t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv)
{
    long v;

    if (x)
    {
        v = jit_atom_getlong(argv);
        if (x->abs != v)
        {
            x->abs = (v > 0) ? TRUE : FALSE;
            jit_functor_noise_simplex_recalc(x);
        }
        return JIT_ERR_NONE;
    }
    return JIT_ERR_INVALID_PTR;
}

t_jit_err jit_functor_noise_simplex_sign(
    t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv)
{
    long v;

    if (x)
    {
        v = jit_atom_getlong(argv);
        if (x->sign != v)
        {
            x->sign = (v > 0) ? TRUE : FALSE;
            jit_functor_noise_simplex_recalc(x);
        }
        return JIT_ERR_NONE;
    }
    return JIT_ERR_INVALID_PTR;
}

t_jit_err jit_functor_noise_simplex_seed(
    t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv)
{
    long v;

    if (x)
    {
        v = jit_atom_getlong(argv);
        if (x->seed != v)
        {
            x->seed = v;

            // allocate permutation and uniform tables
            jit_functor_noise_permutations_init(&x->p, x->tablesize, x->seed);
            jit_functor_noise_lattice_init(&x->lattice, x->tablesize, x->seed);
            jit_functor_noise_simplex_recalc(x);
        }
        return JIT_ERR_NONE;
    }
    return JIT_ERR_INVALID_PTR;
}

t_jit_err jit_functor_noise_simplex_tablesize(
    t_jit_functor_noise_simplex *x, void *attr, long argc, t_atom *argv)
{
    long v;

    if (x)
    {
        v = jit_atom_getlong(argv);
        if (x->tablesize != v && v > 0)
        {

            x->tablesize = jit_math_roundup_poweroftwo(v) * 2;
            x->tablemask = x->tablesize - 1;

            // allocate permutation and uniform tables
            jit_functor_noise_permutations_init(&x->p, x->tablesize, x->seed);
            jit_functor_noise_lattice_init(&x->lattice, x->tablesize, x->seed);
            jit_functor_noise_simplex_recalc(x);
        }
        return JIT_ERR_NONE;
    }
    return JIT_ERR_INVALID_PTR;
}

t_jit_err jit_functor_noise_simplex_recalc(t_jit_functor_noise_simplex *x)
{
    // calculate intermediary values for efficiency
    return JIT_ERR_NONE;
}

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

float jit_functor_noise_simplex_eval_float32_dim2(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals)
{
    long i, j;
    long i1, j1;
    long ii, jj;
    long gi0, gi1, gi2;

    float s, t;
    float t0, t1, t2;
    float gx, gy;
    float n0, n1, n2;
    float X0, Y0, x0, y0;
    float x1, y1;
    float x2, y2;
    float xin, yin;
    float signal;

    if (!x || !vals || dimcount < 2)
        return JIT_ERR_GENERIC;

    // get input space
    xin = vals[0];
    yin = vals[1];

    // skew input space to determine which simplex we are in
    s = (xin + yin) * x->skew.float32;
    i = JIT_MATH_F32_FLOOR(xin + s);
    j = JIT_MATH_F32_FLOOR(yin + s);
    t = (i + j) * x->unskew.float32;

    // unskew cell origin back to 2d space
    X0 = i - t;
    Y0 = j - t;

    // get distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;

    // determine which simplex we are in.
    if (x0 > y0)
    {
        // lower triangle
        i1 = 1;
        j1 = 0;
    }
    else
    {
        // upper triangle
        i1 = 0;
        j1 = 1;
    }

    // unskew offsets for middle corner
    x1 = x0 - i1 + x->unskew.float32; 
    y1 = y0 - j1 + x->unskew.float32;

    // unskew offsets for last corner
    x2 = x0 - 1.0f + 2.0f * x->unskew.float32; 
    y2 = y0 - 1.0f + 2.0f * x->unskew.float32;

    // get the hashed simplex indices of the three simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;

    // get contribution for first corner
    t0 = 0.5f - (x0 * x0 + y0 * y0);
    if (t0 < 0.0f)
    {
        n0 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj)) % 12;
        gx = x->gradients[gi0].float32[0];
        gy = x->gradients[gi0].float32[1];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0);
    }

    // get contribution for middle corner
    t1 = 0.5f - (x1 * x1 + y1 * y1);
    if (t1 < 0.0f)
    {
        n1 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1)) % 12;
        gx = x->gradients[gi1].float32[0];
        gy = x->gradients[gi1].float32[1];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1);
    }

    // get contribution for last corner
    t2 = 0.5f - (x2 * x2 + y2 * y2);
    if (t2 < 0.0f)
    {
        n2 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + 1 + P(jj + 1)) % 12;
        gx = x->gradients[gi2].float32[0];
        gy = x->gradients[gi2].float32[1];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2);
    }

    // add contributions from each corner, scale to [-1,1].
    signal = 70.0f * (n0 + n1 + n2);
    return signal;
}

float jit_functor_noise_simplex_eval_float32_dim3(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals)
{
    float xin, yin, zin;
    float x0, y0, z0;
    float X0, Y0, Z0;
    float x1, y1, z1;
    float x2, y2, z2;
    float x3, y3, z3;
    float gx, gy, gz;
    float n0, n1, n2, n3;
    float t0, t1, t2, t3, t4;
    float s, t;
    float signal;

    long i, j, k;
    long i1, j1, k1;
    long i2, j2, k2;
    long ii, jj, kk;
    long gi0, gi1, gi2, gi3;

    if (!x || !vals || dimcount < 3)
        return JIT_ERR_GENERIC;

    // get the input space
    xin = vals[0];
    yin = vals[1];
    zin = vals[2];

    // skew the input space to determine which simplex cell we're in
    s = (xin + yin + zin) * x->skew.float32;
    i = JIT_MATH_F32_FLOOR(xin + s);
    j = JIT_MATH_F32_FLOOR(yin + s);
    k = JIT_MATH_F32_FLOOR(zin + s);
    t = (i + j + k) * x->unskew.float32;

    // unskew the cell origin back to input space
    X0 = i - t;
    Y0 = j - t;
    Z0 = k - t;

    // get the distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;
    z0 = zin - Z0;

    // determine which simplex we are in.
    if (x0 >= y0)
    {
        if (y0 >= z0)
        { 
            // X Y Z order
            i1 = 1; j1 = 0; k1 = 0;
            i2 = 1; j2 = 1; k2 = 0;
        }
        else if (x0 >= z0)
        {
            // X Z Y order
            i1 = 1; j1 = 0; k1 = 0;
            i2 = 1; j2 = 0; k2 = 1;
        } 
        else
        {
            // Z X Y order
            i1 = 0; j1 = 0; k1 = 1;
            i2 = 1; j2 = 0; k2 = 1;
        } 
    }
    else
    {   
        if (y0 < z0)
        {
            // Z Y X order
            i1 = 0; j1 = 0; k1 = 1;
            i2 = 0; j2 = 1; k2 = 1;
        }
        else if (x0 < z0)
        {
            // Y Z X order
            i1 = 0; j1 = 1; k1 = 0;
            i2 = 0; j2 = 1; k2 = 1;
        } 
        else
        {
            // Y X Z order
            i1 = 0; j1 = 1; k1 = 0;
            i2 = 1; j2 = 1; k2 = 0;
        } 
    }
    
    // get offsets for second corner
    x1 = x0 - i1 + x->unskew.float32; 
    y1 = y0 - j1 + x->unskew.float32;
    z1 = z0 - k1 + x->unskew.float32;

    // get offsets for third corner
    x2 = x0 - i2 + 2.0f * x->unskew.float32;
    y2 = y0 - j2 + 2.0f * x->unskew.float32;
    z2 = z0 - k2 + 2.0f * x->unskew.float32;

    // get offsets for last corner
    x3 = x0 - 1.0f + 3.0f * x->unskew.float32;
    y3 = y0 - 1.0f + 3.0f * x->unskew.float32;
    z3 = z0 - 1.0f + 3.0f * x->unskew.float32;

    // get the hashed simplex indices of the four simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;
    kk = k & x->tablemask;

    // calculate the contribution from the four corners
    t0 = 0.6f - (x0 * x0 + y0 * y0 + z0 * z0);
    if (t0 < 0.0f)
    {
        n0 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj + P(kk))) % 12;
        gx = x->gradients[gi0].float32[0];
        gy = x->gradients[gi0].float32[1];
        gz = x->gradients[gi0].float32[2];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0 + gz * z0);
    }

    t1 = 0.6f - (x1 * x1 + y1 * y1 + z1 * z1);
    if (t1 < 0.0f)
    {
        n1 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1 + P(kk + k1))) % 12;
        gx = x->gradients[gi1].float32[0];
        gy = x->gradients[gi1].float32[1];
        gz = x->gradients[gi1].float32[2];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1 + gz * z1);
    }

    t2 = 0.6f - (x2 * x2 + y2 * y2 + z2 * z2);
    if (t2 < 0.0f)
    {
        n2 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + i2 + P(jj + j2 + P(kk + k2))) % 12;
        gx = x->gradients[gi2].float32[0];
        gy = x->gradients[gi2].float32[1];
        gz = x->gradients[gi2].float32[2];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2 + gz * z2);
    }

    t3 = 0.6f - (x3 * x3 + y3 * y3 + z3 * z3);
    if (t3 < 0.0f)
    {
        n3 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi3 = P(ii + 1 + P(jj + 1 + P(kk + 1))) % 12;
        gx = x->gradients[gi3].float32[0];
        gy = x->gradients[gi3].float32[1];
        gz = x->gradients[gi3].float32[2];

        // compute contribution
        t3 *= t3;
        n3 = t3 * t3 * (gx * x3 + gy * y3 + gz * z3);
    }

    // add contributions from each corner to get the final noise value, scale to [-1,1]
    signal = 32.0f * (n0 + n1 + n2 + n3);
    return signal;

}

float jit_functor_noise_simplex_eval_float32_dim4(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals)
{
    long i, j, k, l;
    long ii, jj, kk, ll;
    long i1, j1, k1, l1;
    long i2, j2, k2, l2;
    long i3, j3, k3, l3;
    long c, c1, c2, c3, c4, c5, c6;
    long gi0, gi1, gi2, gi3, gi4;

    float s, t;
    float signal;
    float xin, yin, zin, win;
    float X0, Y0, Z0, W0;
    float x0, y0, z0, w0;
    float x1, y1, z1, w1;
    float x2, y2, z2, w2;
    float x3, y3, z3, w3;
    float x4, y4, z4, w4;
    float gx, gy, gz, gw;
    float n0, n1, n2, n3, n4;
    float t0, t1, t2, t3, t4;

    if (!x || !vals || dimcount < 4)
        return JIT_ERR_GENERIC;

    xin = vals[0];
    yin = vals[1];
    zin = vals[2];
    win = vals[3];

    // skew the input space to determine which cell of 24 simplices we're in
    s = (xin + yin + zin + win) * x->skew.float32;
    i = JIT_MATH_F32_FLOOR(xin + s);
    j = JIT_MATH_F32_FLOOR(yin + s);
    k = JIT_MATH_F32_FLOOR(zin + s);
    l = JIT_MATH_F32_FLOOR(win + s);
    t = (i + j + k + l) * x->unskew.float32;

    // unskew the cell origin back to (x,y,z,w) space
    X0 = i - t; 
    Y0 = j - t;
    Z0 = k - t;
    W0 = l - t;

    // get the distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;
    z0 = zin - Z0;
    w0 = win - W0;

    // To find out which of the 24 possible simplices we're in, we need to
    // determine the magnitude ordering of x0, y0, z0 and w0.
    // The method below is a good way of finding the ordering of x,y,z,w and
    // then find the correct traversal order for the simplex we�re in.
    // First, six pair-wise comparisons are performed between each possible pair
    // of the four coordinates, and the results are used to add up binary bits
    // for an integer index.
    c1 = (x0 > y0) ? 32 : 0;
    c2 = (x0 > z0) ? 16 : 0;
    c3 = (y0 > z0) ? 8 : 0;
    c4 = (x0 > w0) ? 4 : 0;
    c5 = (y0 > w0) ? 2 : 0;
    c6 = (z0 > w0) ? 1 : 0;

    // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
    c = c1 + c2 + c3 + c4 + c5 + c6;

    // the number 3 in the "simplex" array is at the position of the largest coordinate.
    i1 = jit_functor_noise_simplex_index[c][0] >= 3 ? 1 : 0;
    j1 = jit_functor_noise_simplex_index[c][1] >= 3 ? 1 : 0;
    k1 = jit_functor_noise_simplex_index[c][2] >= 3 ? 1 : 0;
    l1 = jit_functor_noise_simplex_index[c][3] >= 3 ? 1 : 0;

    // the number 2 in the "simplex" array is at the second largest coordinate.
    i2 = jit_functor_noise_simplex_index[c][0] >= 2 ? 1 : 0;
    j2 = jit_functor_noise_simplex_index[c][1] >= 2 ? 1 : 0;
    k2 = jit_functor_noise_simplex_index[c][2] >= 2 ? 1 : 0;
    l2 = jit_functor_noise_simplex_index[c][3] >= 2 ? 1 : 0;

    // the number 1 in the "simplex" array is at the second smallest coordinate.
    i3 = jit_functor_noise_simplex_index[c][0] >= 1 ? 1 : 0;
    j3 = jit_functor_noise_simplex_index[c][1] >= 1 ? 1 : 0;
    k3 = jit_functor_noise_simplex_index[c][2] >= 1 ? 1 : 0;
    l3 = jit_functor_noise_simplex_index[c][3] >= 1 ? 1 : 0;

    // skip the fifth corner, it has all coordinate offsets = 1, so no need to look that up.

    // get offsets for second corner
    x1 = x0 - i1 + x->unskew.float32; 
    y1 = y0 - j1 + x->unskew.float32;
    z1 = z0 - k1 + x->unskew.float32;
    w1 = w0 - l1 + x->unskew.float32;

    // get offsets for third corner
    x2 = x0 - i2 + 2.0f * x->unskew.float32; 
    y2 = y0 - j2 + 2.0f * x->unskew.float32;
    z2 = z0 - k2 + 2.0f * x->unskew.float32;
    w2 = w0 - l2 + 2.0f * x->unskew.float32;

    // get offsets for fourth corner
    x3 = x0 - i3 + 3.0f * x->unskew.float32;
    y3 = y0 - j3 + 3.0f * x->unskew.float32;
    z3 = z0 - k3 + 3.0f * x->unskew.float32;
    w3 = w0 - l3 + 3.0f * x->unskew.float32;

    // get offsets for last corner
    x4 = x0 - 1.0f + 4.0f * x->unskew.float32;
    y4 = y0 - 1.0f + 4.0f * x->unskew.float32;
    z4 = z0 - 1.0f + 4.0f * x->unskew.float32;
    w4 = w0 - 1.0f + 4.0f * x->unskew.float32;

    // get the hashed simplex indices of the five simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;
    kk = k & x->tablemask;
    ll = l & x->tablemask;

    // calculate the contribution from the five corners
    t0 = 0.6f - (x0 * x0 + y0 * y0 + z0 * z0 + w0 * w0);
    if (t0 < 0.0f)
    {
        n0 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj + P(kk + P(ll)))) % 32;
        gx = x->gradients[gi0].float32[0];
        gy = x->gradients[gi0].float32[1];
        gz = x->gradients[gi0].float32[2];
        gw = x->gradients[gi0].float32[3];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0 + gz * z0 + gw * w0);
    }

    t1 = 0.6 - (x1 * x1 + y1 * y1 + z1 * z1 + w1 * w1);
    if (t1 < 0.0f)
    {
        n1 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1 + P(kk + k1 + P(ll + l1)))) % 32;
        gx = x->gradients[gi1].float32[0];
        gy = x->gradients[gi1].float32[1];
        gz = x->gradients[gi1].float32[2];
        gw = x->gradients[gi1].float32[3];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1 + gz * z1 + gw * w1);
    }

    t2 = 0.6 - (x2 * x2 + y2 * y2 + z2 * z2 + w2 * w2);
    if (t2 < 0.0f)
    {
        n2 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + i2 + P(jj + j2 + P(kk + k2 + P(ll + l2)))) % 32;
        gx = x->gradients[gi2].float32[0];
        gy = x->gradients[gi2].float32[1];
        gz = x->gradients[gi2].float32[2];
        gw = x->gradients[gi2].float32[3];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2 + gz * z2 + gw * w2);
    }

    t3 = 0.6f - (x3 * x3 + y3 * y3 + z3 * z3 + w3 * w3);
    if (t3 < 0.0f)
    {
        n3 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi3 = P(ii + i3 + P(jj + j3 + P(kk + k3 + P(ll + l3)))) % 32;
        gx = x->gradients[gi3].float32[0];
        gy = x->gradients[gi3].float32[1];
        gz = x->gradients[gi3].float32[2];
        gw = x->gradients[gi3].float32[3];

        // compute contribution
        t3 *= t3;
        n3 = t3 * t3 * (gx * x3 + gy * y3 + gz * z3 + gw * w3);
    }

    t4 = 0.6f - (x4 * x4 + y4 * y4 + z4 * z4 + w4 * w4);
    if (t4 < 0.0f)
    {
        n4 = 0.0f;
    }
    else
    {
        // get gradient value via permutation table
        gi4 = P(ii + 1 + P(jj + 1 + P(kk + 1 + P(ll + 1)))) % 32;
        gx = x->gradients[gi4].float32[0];
        gy = x->gradients[gi4].float32[1];
        gz = x->gradients[gi4].float32[2];
        gw = x->gradients[gi4].float32[3];

        // compute contribution
        t4 *= t4;
        n4 = t4 * t4 * (gx * x4 + gy * y4 + gz * z4 + gw * w4);
    }

    // sum up and scale the result to cover the range [-1,1]
    signal = 27.0f * (n0 + n1 + n2 + n3 + n4);
    return signal;
}

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

double jit_functor_noise_simplex_eval_float64_dim2(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals)
{
    long i, j;
    long i1, j1;
    long ii, jj;
    long gi0, gi1, gi2;

    double s, t;
    double t0, t1, t2;
    double gx, gy;
    double n0, n1, n2;
    double X0, Y0, x0, y0;
    double x1, y1;
    double x2, y2;
    double xin, yin;
    double signal;

    if (!x || !vals || dimcount < 2)
        return JIT_ERR_GENERIC;

    // get input space
    xin = vals[0];
    yin = vals[1];

    // skew input space to determine which simplex we are in
    s = (xin + yin) * x->skew.float64;
    i = JIT_MATH_F64_FLOOR(xin + s);
    j = JIT_MATH_F64_FLOOR(yin + s);
    t = (i + j) * x->unskew.float64;

    // unskew cell origin back to 2d space
    X0 = i - t;
    Y0 = j - t;

    // get distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;

    // determine which simplex we are in.
    if (x0 > y0)
    {
        // lower triangle
        i1 = 1;
        j1 = 0;
    }
    else
    {
        // upper triangle
        i1 = 0;
        j1 = 1;
    }

    // unskew offsets for middle corner
    x1 = x0 - i1 + x->unskew.float64; 
    y1 = y0 - j1 + x->unskew.float64;

    // unskew offsets for last corner
    x2 = x0 - 1.0 + 2.0 * x->unskew.float64; 
    y2 = y0 - 1.0 + 2.0 * x->unskew.float64;

    // get the hashed simplex indices of the three simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;

    // get contribution for first corner
    t0 = 0.5f - (x0 * x0 + y0 * y0);
    if (t0 < 0.0)
    {
        n0 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj)) % 12;
        gx = x->gradients[gi0].float64[0];
        gy = x->gradients[gi0].float64[1];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0);
    }

    // get contribution for middle corner
    t1 = 0.5f - (x1 * x1 + y1 * y1);
    if (t1 < 0.0)
    {
        n1 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1)) % 12;
        gx = x->gradients[gi1].float64[0];
        gy = x->gradients[gi1].float64[1];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1);
    }

    // get contribution for last corner
    t2 = 0.5f - (x2 * x2 + y2 * y2);
    if (t2 < 0.0)
    {
        n2 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + 1 + P(jj + 1)) % 12;
        gx = x->gradients[gi2].float64[0];
        gy = x->gradients[gi2].float64[1];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2);
    }

    // add contributions from each corner, scale to [-1,1].
    signal = 70.0 * (n0 + n1 + n2);

    return signal;
}

double jit_functor_noise_simplex_eval_float64_dim3(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals)
{
    double xin, yin, zin;
    double x0, y0, z0;
    double X0, Y0, Z0;
    double x1, y1, z1;
    double x2, y2, z2;
    double x3, y3, z3;
    double gx, gy, gz;
    double n0, n1, n2, n3;
    double t0, t1, t2, t3, t4;
    double s, t;
    double signal;

    long i, j, k;
    long i1, j1, k1;
    long i2, j2, k2;
    long ii, jj, kk;
    long gi0, gi1, gi2, gi3;

    if (!x || !vals || dimcount < 3)
        return JIT_ERR_GENERIC;

    // get the input space
    xin = vals[0];
    yin = vals[1];
    zin = vals[2];

    // skew the input space to determine which simplex cell we're in
    s = (xin + yin + zin) * x->skew.float64;
    i = JIT_MATH_F64_FLOOR(xin + s);
    j = JIT_MATH_F64_FLOOR(yin + s);
    k = JIT_MATH_F64_FLOOR(zin + s);
    t = (i + j + k) * x->unskew.float64;

    // unskew the cell origin back to input space
    X0 = i - t;
    Y0 = j - t;
    Z0 = k - t;

    // get the distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;
    z0 = zin - Z0;

    // determine which simplex we are in.
    if (x0 >= y0)
    {
        if (y0 >= z0)
        { 
            // X Y Z order
            i1 = 1; j1 = 0; k1 = 0;
            i2 = 1; j2 = 1; k2 = 0;
        }
        else if (x0 >= z0)
        {
            // X Z Y order
            i1 = 1; j1 = 0; k1 = 0;
            i2 = 1; j2 = 0; k2 = 1;
        } 
        else
        {
            // Z X Y order
            i1 = 0; j1 = 0; k1 = 1;
            i2 = 1; j2 = 0; k2 = 1;
        } 
    }
    else
    {   
        if (y0 < z0)
        {
            // Z Y X order
            i1 = 0; j1 = 0; k1 = 1;
            i2 = 0; j2 = 1; k2 = 1;
        }
        else if (x0 < z0)
        {
            // Y Z X order
            i1 = 0; j1 = 1; k1 = 0;
            i2 = 0; j2 = 1; k2 = 1;
        } 
        else
        {
            // Y X Z order
            i1 = 0; j1 = 1; k1 = 0;
            i2 = 1; j2 = 1; k2 = 0;
        } 
    }
    
    // get offsets for second corner
    x1 = x0 - i1 + x->unskew.float64; 
    y1 = y0 - j1 + x->unskew.float64;
    z1 = z0 - k1 + x->unskew.float64;

    // get offsets for third corner
    x2 = x0 - i2 + 2.0 * x->unskew.float64;
    y2 = y0 - j2 + 2.0 * x->unskew.float64;
    z2 = z0 - k2 + 2.0 * x->unskew.float64;

    // get offsets for last corner
    x3 = x0 - 1.0 + 3.0 * x->unskew.float64;
    y3 = y0 - 1.0 + 3.0 * x->unskew.float64;
    z3 = z0 - 1.0 + 3.0 * x->unskew.float64;

    // get the hashed simplex indices of the four simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;
    kk = k & x->tablemask;

    // calculate the contribution from the four corners
    t0 = 0.6 - (x0 * x0 + y0 * y0 + z0 * z0);
    if (t0 < 0)
    {
        n0 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj + P(kk))) % 12;
        gx = x->gradients[gi0].float64[0];
        gy = x->gradients[gi0].float64[1];
        gz = x->gradients[gi0].float64[2];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0 + gz * z0);
    }

    t1 = 0.6 - (x1 * x1 + y1 * y1 + z1 * z1);
    if (t1 < 0)
    {
        n1 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1 + P(kk + k1))) % 12;
        gx = x->gradients[gi1].float64[0];
        gy = x->gradients[gi1].float64[1];
        gz = x->gradients[gi1].float64[2];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1 + gz * z1);
    }

    t2 = 0.6 - (x2 * x2 + y2 * y2 + z2 * z2);
    if (t2 < 0)
    {
        n2 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + i2 + P(jj + j2 + P(kk + k2))) % 12;
        gx = x->gradients[gi2].float64[0];
        gy = x->gradients[gi2].float64[1];
        gz = x->gradients[gi2].float64[2];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2 + gz * z2);
    }

    t3 = 0.6 - (x3 * x3 + y3 * y3 + z3 * z3);
    if (t3 < 0)
    {
        n3 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi3 = P(ii + 1 + P(jj + 1 + P(kk + 1))) % 12;
        gx = x->gradients[gi3].float64[0];
        gy = x->gradients[gi3].float64[1];
        gz = x->gradients[gi3].float64[2];

        // compute contribution
        t3 *= t3;
        n3 = t3 * t3 * (gx * x3 + gy * y3 + gz * z3);
    }

    // add contributions from each corner to get the final noise value, scale to [-1,1]
    signal = 32.0 * (n0 + n1 + n2 + n3);

    return signal;

}

double jit_functor_noise_simplex_eval_float64_dim4(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals)
{
    long i, j, k, l;
    long ii, jj, kk, ll;
    long i1, j1, k1, l1;
    long i2, j2, k2, l2;
    long i3, j3, k3, l3;
    long c, c1, c2, c3, c4, c5, c6;
    long gi0, gi1, gi2, gi3, gi4;

    double s, t;
    double signal;
    double xin, yin, zin, win;
    double X0, Y0, Z0, W0;
    double x0, y0, z0, w0;
    double x1, y1, z1, w1;
    double x2, y2, z2, w2;
    double x3, y3, z3, w3;
    double x4, y4, z4, w4;
    double gx, gy, gz, gw;
    double n0, n1, n2, n3, n4;
    double t0, t1, t2, t3, t4;

    if (!x || !vals || dimcount < 4)
        return JIT_ERR_GENERIC;

    xin = vals[0];
    yin = vals[1];
    zin = vals[2];
    win = vals[3];

    // skew the input space to determine which cell of 24 simplices we're in
    s = (xin + yin + zin + win) * x->skew.float64;
    i = JIT_MATH_F64_FLOOR(xin + s);
    j = JIT_MATH_F64_FLOOR(yin + s);
    k = JIT_MATH_F64_FLOOR(zin + s);
    l = JIT_MATH_F64_FLOOR(win + s);
    t = (i + j + k + l) * x->unskew.float64;

    // unskew the cell origin back to (x,y,z,w) space
    X0 = i - t; 
    Y0 = j - t;
    Z0 = k - t;
    W0 = l - t;

    // get the distances from the cell origin
    x0 = xin - X0;
    y0 = yin - Y0;
    z0 = zin - Z0;
    w0 = win - W0;

    // To find out which of the 24 possible simplices we're in, we need to
    // determine the magnitude ordering of x0, y0, z0 and w0.
    // The method below is a good way of finding the ordering of x,y,z,w and
    // then find the correct traversal order for the simplex we�re in.
    // First, six pair-wise comparisons are performed between each possible pair
    // of the four coordinates, and the results are used to add up binary bits
    // for an integer index.
    c1 = (x0 > y0) ? 32 : 0;
    c2 = (x0 > z0) ? 16 : 0;
    c3 = (y0 > z0) ? 8 : 0;
    c4 = (x0 > w0) ? 4 : 0;
    c5 = (y0 > w0) ? 2 : 0;
    c6 = (z0 > w0) ? 1 : 0;

    // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
    c = c1 + c2 + c3 + c4 + c5 + c6;

    // the number 3 in the "simplex" array is at the position of the largest coordinate.
    i1 = jit_functor_noise_simplex_index[c][0] >= 3 ? 1 : 0;
    j1 = jit_functor_noise_simplex_index[c][1] >= 3 ? 1 : 0;
    k1 = jit_functor_noise_simplex_index[c][2] >= 3 ? 1 : 0;
    l1 = jit_functor_noise_simplex_index[c][3] >= 3 ? 1 : 0;

    // the number 2 in the "simplex" array is at the second largest coordinate.
    i2 = jit_functor_noise_simplex_index[c][0] >= 2 ? 1 : 0;
    j2 = jit_functor_noise_simplex_index[c][1] >= 2 ? 1 : 0;
    k2 = jit_functor_noise_simplex_index[c][2] >= 2 ? 1 : 0;
    l2 = jit_functor_noise_simplex_index[c][3] >= 2 ? 1 : 0;

    // the number 1 in the "simplex" array is at the second smallest coordinate.
    i3 = jit_functor_noise_simplex_index[c][0] >= 1 ? 1 : 0;
    j3 = jit_functor_noise_simplex_index[c][1] >= 1 ? 1 : 0;
    k3 = jit_functor_noise_simplex_index[c][2] >= 1 ? 1 : 0;
    l3 = jit_functor_noise_simplex_index[c][3] >= 1 ? 1 : 0;

    // skip the fifth corner, it has all coordinate offsets = 1, so no need to look that up.

    // get offsets for second corner
    x1 = x0 - i1 + x->unskew.float64; 
    y1 = y0 - j1 + x->unskew.float64;
    z1 = z0 - k1 + x->unskew.float64;
    w1 = w0 - l1 + x->unskew.float64;

    // get offsets for third corner
    x2 = x0 - i2 + 2.0 * x->unskew.float64; 
    y2 = y0 - j2 + 2.0 * x->unskew.float64;
    z2 = z0 - k2 + 2.0 * x->unskew.float64;
    w2 = w0 - l2 + 2.0 * x->unskew.float64;

    // get offsets for fourth corner
    x3 = x0 - i3 + 3.0 * x->unskew.float64;
    y3 = y0 - j3 + 3.0 * x->unskew.float64;
    z3 = z0 - k3 + 3.0 * x->unskew.float64;
    w3 = w0 - l3 + 3.0 * x->unskew.float64;

    // get offsets for last corner
    x4 = x0 - 1.0 + 4.0 * x->unskew.float64;
    y4 = y0 - 1.0 + 4.0 * x->unskew.float64;
    z4 = z0 - 1.0 + 4.0 * x->unskew.float64;
    w4 = w0 - 1.0 + 4.0 * x->unskew.float64;

    // get the hashed simplex indices of the five simplex corners
    ii = i & x->tablemask;
    jj = j & x->tablemask;
    kk = k & x->tablemask;
    ll = l & x->tablemask;

    // calculate the contribution from the five corners
    t0 = 0.6 - (x0 * x0 + y0 * y0 + z0 * z0 + w0 * w0);
    if (t0 < 0)
    {
        n0 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi0 = P(ii + P(jj + P(kk + P(ll)))) % 32;
        gx = x->gradients[gi0].float64[0];
        gy = x->gradients[gi0].float64[1];
        gz = x->gradients[gi0].float64[2];
        gw = x->gradients[gi0].float64[3];

        // compute contribution
        t0 *= t0;
        n0 = t0 * t0 * (gx * x0 + gy * y0 + gz * z0 + gw * w0);
    }

    t1 = 0.6 - (x1 * x1 + y1 * y1 + z1 * z1 + w1 * w1);
    if (t1 < 0)
    {
        n1 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi1 = P(ii + i1 + P(jj + j1 + P(kk + k1 + P(ll + l1)))) % 32;
        gx = x->gradients[gi1].float64[0];
        gy = x->gradients[gi1].float64[1];
        gz = x->gradients[gi1].float64[2];
        gw = x->gradients[gi1].float64[3];

        // compute contribution
        t1 *= t1;
        n1 = t1 * t1 * (gx * x1 + gy * y1 + gz * z1 + gw * w1);
    }

    t2 = 0.6 - (x2 * x2 + y2 * y2 + z2 * z2 + w2 * w2);
    if (t2 < 0)
    {
        n2 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi2 = P(ii + i2 + P(jj + j2 + P(kk + k2 + P(ll + l2)))) % 32;
        gx = x->gradients[gi2].float64[0];
        gy = x->gradients[gi2].float64[1];
        gz = x->gradients[gi2].float64[2];
        gw = x->gradients[gi2].float64[3];

        // compute contribution
        t2 *= t2;
        n2 = t2 * t2 * (gx * x2 + gy * y2 + gz * z2 + gw * w2);
    }

    t3 = 0.6 - (x3 * x3 + y3 * y3 + z3 * z3 + w3 * w3);
    if (t3 < 0)
    {
        n3 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi3 = P(ii + i3 + P(jj + j3 + P(kk + k3 + P(ll + l3)))) % 32;
        gx = x->gradients[gi3].float64[0];
        gy = x->gradients[gi3].float64[1];
        gz = x->gradients[gi3].float64[2];
        gw = x->gradients[gi3].float64[3];

        // compute contribution
        t3 *= t3;
        n3 = t3 * t3 * (gx * x3 + gy * y3 + gz * z3 + gw * w3);
    }

    t4 = 0.6 - (x4 * x4 + y4 * y4 + z4 * z4 + w4 * w4);
    if (t4 < 0)
    {
        n4 = 0.0;
    }
    else
    {
        // get gradient value via permutation table
        gi4 = P(ii + 1 + P(jj + 1 + P(kk + 1 + P(ll + 1)))) % 32;
        gx = x->gradients[gi4].float64[0];
        gy = x->gradients[gi4].float64[1];
        gz = x->gradients[gi4].float64[2];
        gw = x->gradients[gi4].float64[3];

        // compute contribution
        t4 *= t4;
        n4 = t4 * t4 * (gx * x4 + gy * y4 + gz * z4 + gw * w4);
    }

    // sum up and scale the result to cover the range [-1,1]
    signal = 27.0 * (n0 + n1 + n2 + n3 + n4);

    // remap
    return signal;

}

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

t_jit_err jit_functor_noise_simplex_dir_init(
    t_jit_functor_noise_simplex *x)
{
    long i, j;
    long gradcount = 0;
    t_jit_err err = JIT_ERR_NONE;

    if (x->gradients)
        jit_disposeptr((char*)x->gradients);
    x->gradients = NULL;

    jit_rand_setseed(x->seed);
    switch (x->dimcount)
    {
    case 1:
    case 2:
        {
            gradcount = jit_functor_noise_simplex_grad3count;
            x->gradients = (t_jit_functor_combined_dimvalue*)jit_newptr(gradcount * sizeof(t_jit_functor_combined_dimvalue));
            if (!x->gradients)
            {
                err = JIT_ERR_INVALID_PTR;
                goto out;
            }

            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->skew, F2);
            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->unskew, G2);

            // use 3d gradients
            for (i = 0; i < gradcount; i++)
            {
                for (j = 0; j < x->dimcount; j++)
                    x->gradients[i].float64[j] = jit_functor_noise_simplex_grad3[i][j] * 4.0f - 1.0f;
            }
            break;
        }
    case 3:
        {
            gradcount = jit_functor_noise_simplex_grad3count;
            x->gradients = (t_jit_functor_combined_dimvalue*)jit_newptr(gradcount * sizeof(t_jit_functor_combined_dimvalue));
            if (!x->gradients)
            {
                err = JIT_ERR_INVALID_PTR;
                goto out;
            }

            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->skew, F3);
            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->unskew, G3);

            // use 3d gradients
            for (i = 0; i < gradcount; i++)
            {
                for (j = 0; j < x->dimcount; j++)
                    x->gradients[i].float64[j] = jit_functor_noise_simplex_grad3[i][j] * 4.0f - 1.0f;
            }
            break;
        }
    case 4:
        {
            gradcount = jit_functor_noise_simplex_grad4count;
            x->gradients = (t_jit_functor_combined_dimvalue*)jit_newptr(gradcount * sizeof(t_jit_functor_combined_dimvalue));
            if (!x->gradients)
            {
                err = JIT_ERR_INVALID_PTR;
                goto out;
            }

            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->skew, F4);
            JIT_FUNCTOR_COMBINED_VALUE_SETALL(x->unskew, G4);

            // use 4d gradients
            for (i = 0; i < gradcount; i++)
            {
                for (j = 0; j < x->dimcount; j++)
                    x->gradients[i].float64[j] = jit_functor_noise_simplex_grad4[i][j] * 4.0f - 1.0f;
            }
            break;
        }
    default:
        {
            // TODO!  n-dim simplex gradient dir ?
            jit_object_error((t_object *)x,"jit.noise.simplex: dimcount too high -- max dimcount is 4");
        }
    };

    // copy and convert
    for (i = 0; i < gradcount; i++)
    {
        for (j = 0; j < x->dimcount; j++)
        {
            x->gradients[i].float32[j] = (float)x->gradients[i].float64[j];
            x->gradients[i].fixed[j] = DoubleToFixed(x->gradients[i].float64[j]);
        }
    }

out:
    return err;
}

long jit_functor_noise_simplex_eval_fixed(
    t_jit_functor_noise_simplex *x, long dimcount, long *vals)
{
    long i;
    float floatvals[JIT_MATRIX_MAX_DIMCOUNT];
    
    for(i = 0; i < dimcount; i++)
        floatvals[i] = FixedToFloat(vals[i]);

    return FloatToFixed(jit_functor_noise_simplex_eval_float32(x, dimcount, floatvals));
}

float jit_functor_noise_simplex_eval_float32(
    t_jit_functor_noise_simplex *x, long dimcount, float *vals)
{
    long i;
    float v2[2];
    float signal = 0;
    
    // safety
    if (!x->p)
        jit_functor_noise_permutations_init(&x->p, x->tablesize, x->seed);

    if (!x->lattice.float32)
        jit_functor_noise_lattice_init(&x->lattice, x->tablesize, x->seed);

    if (!x->gradients || dimcount != x->dimcount)
    {
        x->dimcount = dimcount;
        jit_functor_noise_simplex_dir_init(x);
    }

    switch(dimcount)
    {
    case 1:
        v2[0] = vals[0];
        v2[1] = 0;
        signal = jit_functor_noise_simplex_eval_float32_dim2(x, 2, v2);
        break;
    case 2:
        signal = jit_functor_noise_simplex_eval_float32_dim2(x, dimcount, vals);
        break;
    case 3:
        signal = jit_functor_noise_simplex_eval_float32_dim3(x, dimcount, vals);
        break;
    case 4:
        signal = jit_functor_noise_simplex_eval_float32_dim4(x, dimcount, vals);
        break;
    default:
        // iterate downward 
        for (i = 0; i < dimcount;i++) 
        {
            signal = jit_functor_noise_simplex_eval_float32_dim4(x, dimcount-1, vals);
        }
    };

    if (!x->sign)
        signal = 0.5f * signal + 0.5f;

    if (x->abs)
        signal = JIT_MATH_F32_ABS(signal);

    return signal;
}

double jit_functor_noise_simplex_eval_float64(
    t_jit_functor_noise_simplex *x, long dimcount, double *vals)
{
    long i;
    double v2[2];
    double signal = 0;
    
    // safety
    if (!x->p)
        jit_functor_noise_permutations_init(&x->p, x->tablesize, x->seed);

    if (!x->lattice.float64)
        jit_functor_noise_lattice_init(&x->lattice, x->tablesize, x->seed);

    if (!x->gradients || dimcount != x->dimcount)
    {
        x->dimcount = dimcount;
        jit_functor_noise_simplex_dir_init(x);
    }

    switch(dimcount)
    {
    case 1:
        v2[0] = vals[0];
        v2[1] = 0;
        signal = jit_functor_noise_simplex_eval_float64_dim2(x, 2, v2);
        break;
    case 2:
        signal = jit_functor_noise_simplex_eval_float64_dim2(x, dimcount, vals);
        break;
    case 3:
        signal = jit_functor_noise_simplex_eval_float64_dim3(x, dimcount, vals);
        break;
    case 4:
        signal = jit_functor_noise_simplex_eval_float64_dim4(x, dimcount, vals);
        break;
    default:
        // iterate downward 
        for (i = 0; i < dimcount;i++) 
        {
            signal = jit_functor_noise_simplex_eval_float64_dim4(x, dimcount-1, vals);
        }
    };

    if (!x->sign)
        signal = 0.5 * signal + 0.5;

    if (x->abs)
        signal = JIT_MATH_F64_ABS(signal);

    return signal;
}

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

Copyright © 2008, Cycling '74