Max 5 API Reference

jit.math.c

/* 
    jit.math.c
    Copyright 2001-2005 - Cycling '74
    Joshua Kit Clayton jkc@cycling74.com
*/

#include "jit.common.h"
#include <math.h>
#include <float.h>

#ifdef WIN_VERSION
double acosh(double x);
double asinh(double x);
double atanh(double x);
double expm1(double x);
double exp2(double x);
double log2(double x);
double round(double x);
double trunc(double x);
#endif //WIN_VERSION

/**
 * @defgroup mathmod Math Module
    @ingroup jitter
 */


/**
 * Calculates the cosine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_cos(double x)
{
    return cos(x);
}

/**
 * Calculates the sine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_sin(double x)
{
    return sin(x);
}

/**
 * Calculates the tangent.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_tan(double x)
{
    return tan(x);
}

/**
 * Calculates the arccosine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_acos(double x)
{
    return acos(x);
}

/**
 * Calculates the arcsine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_asin(double x)
{
    return asin(x);
}

/**
 * Calculates the arctangent.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_atan(double x)
{
    return atan(x);
}

/**
 * Calculates the four quadrant arctangent.
 *
 * @ingroup mathmod 
 *
 * @param   y       input
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_atan2(double y, double x)
{
    return atan2(y,x);
}

/**
 * Calculates the hyperbolic cosine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_cosh(double x)
{
    return cosh(x);
}

/**
 * Calculates the hyperbolic sine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_sinh(double x)
{
    return sinh(x);
}

/**
 * Calculates the hyperbolic tangent.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_tanh(double x)
{
    return tanh(x);
}

/**
 * Calculates the hyperbolic arccosine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_acosh(double x)
{
    return acosh(x);
}

/**
 * Calculates the hyperbolic arcsine.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_asinh(double x)
{
    return asinh(x);
}

/**
 * Calculates the hyperbolic arctangent.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_atanh(double x)
{
    return atanh(x);
}

/**
 * Calculates the exponent.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_exp(double x)
{
    return exp(x);
}

/**
 * Calculates the exponent minus 1.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_expm1(double x)
{
    return expm1(x);
}

/**
 * Calculates the exponent base 2.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_exp2(double x)
{
    return exp2(x);
}

/**
 * Calculates the logarithm.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_log(double x)
{
    return log(x);
}

/**
 * Calculates the logarithm base 2.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_log2(double x)
{
    return log2(x);
}

/**
 * Calculates the logarithm base 10.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_log10(double x)
{
    return log10(x);
}

/**
 * Calculates the hypotenuse.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @param   y       input
 * @return  output 
 *
 */
double jit_math_hypot(double x, double y)
{
    return hypot(x,y);
}

/**
 * Calculates x raised to the y power.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @param   y       input
 * @return  output 
 *
 */
double jit_math_pow(double x, double y)
{
    return pow(x,y);
}

/**
 * Calculates the square root.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_sqrt(double x)
{
    return sqrt(x);
}

/**
 * Calculates the ceiling.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_ceil(double x)
{
    return ceil(x);
}

/**
 * Calculates the floor.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_floor(double x)
{
    return floor(x);
}

/**
 * Rounds the input.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_round(double x)
{
    return round(x);
}

/**
 * Truncates the input.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_trunc(double x)
{
    return trunc(x);
}

/**
 * Calculates the floating point x modulo y.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @param   y       input
 * @return  output 
 *
 */
double jit_math_fmod(double x, double y)
{
    double d;

    if (y) {
        if (y<0) y = -y; // modulus needs to be absolute value
        //don't do float->int->float if can be done w/ one add or sub       
        if (x>=y) {
            if (x>=(y*2.)) {
                d = x / y;
                d = d - (long) d;
                x = d * y;
            } else {
                x -= y;
            }
        } else if (x<=(-y)) {
            if (x<=(-y*2.)) {
                d = x / y;
                d = d - (long) d;
                x = d * y;
            }  else {
                x += y;
            }
        }
    } else x = 0.; //don't divide by zero

    return x;
}


/**
 * Calculates the fold of x between lo and hi.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @param   lo      lower bound
 * @param   hi      upper bound
 * @return  output 
 *
 */
double jit_math_fold(double x, double lo, double hi)
{
    long di;
    double m,d,tmp;
    
    if (lo > hi) {
        tmp = lo;
        lo  = hi;
        hi  = tmp;
    }
    if (lo)
        x -= lo;
    m = hi-lo;
    if (m) {            
        if (x < 0.) 
            x = -x;
        if (x>m) {
            if (x>(m*2.)) {
                d = x / m;
                di = (long) d;
                d = d - (double) di;
                if (di%2) {
                    if (d < 0) {
                        d = -1. - d;
                    } else {                
                        d = 1. - d;
                    }
                }
                x = d * m;
                if (x < 0.) 
                    x = m+x;
            } else {
                x = m-(x-m);
            }
        }           
    } else x = 0.; //don't divide by zero
    
    return x + lo; 
}

/**
 * Calculates the wrap of x between lo and hi.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @param   lo      lower bound
 * @param   hi      upper bound
 * @return  output 
 *
 */
double jit_math_wrap(double x, double lo, double hi)
{
    double m,d,tmp;
    long di;

    if (lo > hi) {
        tmp = lo;
        lo  = hi;
        hi  = tmp;
    }
    if (lo)
        x -= lo;
    m = hi-lo;
    if (m) {        
        if (x>m) {
            if (x>(m*2.)) {
                d = x / m;
                di = (long) d;
                d = d - (double) di;
                x = d * m;
            } else {
                x -= m;
            }
        } else if (x<0.) {
            if (x<(-m)) {
                d = x / m;
                di = (long) d;
                d = d - (double) di;
                x = d * m;
                x += m;
            } else {
                x += m;
            }
        }
    } else x = 0.; //don't divide by zero
    
    return x + lo; 
}


/*
 *  Derek Gerstmann - derek@cycling74.com
 *
 *  The following Bessel functions were taken from the free NUMLIBC
 *  library since CodeWarrior doesn't define them:
 *
 *  http://www.math.ntnu.no/num/nnm/Program/Numlibc/#header
 *
 *  -[dg]
 */
 
/**
 * Calcuates the j1_0 Bessel function.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_j1_0(double x)
{
    static int J1_MAX = 16; 

    static double a[ 16] = {  1.29671754121052984,
                             -1.19180116054121687,
                              1.28799409885767762,
                             -0.66144393413454325,
                              0.17770911723972828,
                             -0.02917552480615421,
                              0.00324027018268386,
                             -0.00026044438934858,
                              0.00001588701923993,
                             -0.00000076175878054,
                              0.00000002949707007,
                             -0.00000000094242130,
                              0.00000000002528124,
                             -0.00000000000057774,
                              0.00000000000001139,
                             -0.00000000000000020 };

    int k;
    double  arg, t, value, b0, b1, b2;

    arg = .125 * x;
    t = 2. * (2.*arg*arg - 1.);

    b2 = 0.;
    b1 = 0.;
    b0 = a[J1_MAX-1];

    for (k = J1_MAX-2; k >= 0; k--) {
            b2 = b1;
            b1 = b0;
            b0 = t * b1 - b2 + a[k];
        }

    value = .5 *  (b0 - b2);

    return(arg * value);

}

/**
 * Calcuates the p1 Bessel function.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_p1( double x)
{
    static int P1_MAX = 13;

    static double a[ 13] = { 2.00180608172002740,
                             0.00089898983308594,
                            -0.00000398728430049,
                             0.00000006177633961,
                            -0.00000000187189075,
                             0.00000000008816899,
                            -0.00000000000570486,
                             0.00000000000046992,
                            -0.00000000000004684,
                             0.00000000000000545,
                            -0.00000000000000072,
                             0.00000000000000011,
                            -0.00000000000000002 };

    int k;
    double  arg, t, value, b0, b1, b2;

    arg = 8. / x; 
    t = 2. * (2.*arg*arg - 1.);

    b2 = 0.;
    b1 = 0.;
    b0 = a[P1_MAX-1];

    for (k = P1_MAX-2; k >= 0; k--) {
            b2 = b1;
            b1 = b0;
            b0 = t * b1 - b2 + a[k];
        }

    value = .5 *  (b0 - b2);

    return(value);
}

/**
 * Calcuates the q1 Bessel function.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_q1(double x)
{
    static int Q1_MAX = 13;

    static double a[ 13] = { 0.09355557413907065,
                            -0.00009627723549157,
                             0.00000091386152580,
                            -0.00000002095978138,
                             0.00000000082291933,
                            -0.00000000004683664,
                             0.00000000000351522,
                            -0.00000000000032643,
                             0.00000000000003597,
                            -0.00000000000000456,
                             0.00000000000000065,
                            -0.00000000000000010,
                             0.00000000000000002 };

    int k;
    double  arg, t, value, b0, b1, b2;

    arg = 8. / x;
    t = 2. * (2.*arg*arg - 1.);

    b2 = 0.;
    b1 = 0.;
    b0 = a[Q1_MAX-1];

    for (k = Q1_MAX-2; k >= 0; k--) {
            b2 = b1;
            b1 = b0;
            b0 = t * b1 - b2 + a[k];
        }

    value = .5 *  (b0 - b2);

    return(arg * value);
}

/**
 * Calcuates the j1 Bessel function.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
double jit_math_j1(double x)
{
    static double cof1 =  .79788456080286536;  /* sqrt(2/pi) */
    static double cof2 = 2.35619449019234493;  /* 3*pi/4 */

    double value, ksi1;
    int flag;

    flag = 1;
    if (x >= -8. && x <= 8.)
        value = jit_math_j1_0(x);
    else if (x > 8.) {
        ksi1 = x - cof2;
        value = jit_math_p1(x) * jit_math_cos(ksi1) - jit_math_q1(x) * jit_math_sin(ksi1);
        value *= cof1/sqrt(x);
    }
    else {
        flag = 0;
        value = 0.;
    }
    return(value);

}

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

/**
 * Rounds up to the nearest power of two.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
unsigned long jit_math_roundup_poweroftwo(unsigned long x) 
{
    x--;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    return x+1;
}

/**
 * Checks if input is finite.
 *
 * @ingroup mathmod 
 *
 * @param   v       input
 * @return  1 if finite. Otherwise, 0. 
 *
 */
long jit_math_is_finite(float v)
{
#ifdef WIN_VERSION
    return _finite(v);
#else
    return !isinf(v);
#endif
}

/**
 * Checks if input is not a number (NaN).
 *
 * @ingroup mathmod 
 *
 * @param   v       input
 * @return  1 if not a number. Otherwise, 0. 
 *
 */
long jit_math_is_nan(float v)
{
#ifdef WIN_VERSION
    return _isnan(v);
#else
    return isnan(v);
#endif
}

/**
 * Checks if input is both finite and a number.
 *
 * @ingroup mathmod 
 *
 * @param   v       input
 * @return  1 if vaild. Otherwise, 0. 
 *
 */
long jit_math_is_valid(float v)
{
    return (!jit_math_is_nan(v) && jit_math_is_finite(v));
}

/**
 * Checks if input is a power of two.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  1 if finite. Otherwise, 0. 
 *
 */
long jit_math_is_poweroftwo(long x) 
{
    return (x & (x - 1)) == 0;
}

float jit_math_project_to_sphere(float r, float x, float y)
{
    float d, t, z;

    d = jit_math_sqrt(x*x + y*y);
    if (d < r * 0.70710678118654752440) 
    {
        // inside sphere
        z = jit_math_sqrt(r*r - d*d);
    } else {
        // on hyperbola
        t = r / (float)1.41421356237309504880;
        z = t*t / d;
    }
    return z;
}

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

/*
 *  Derek Gerstmann - derek@cycling74.com
 *
 *  The following square root table methods were adapted from examples
 *  from NVidia's public domain fast math utility functions:
 *
 *  http://www.nvidia.com/object/fast_math_routines.html
 *
 *  -[dg]
 */
 
static unsigned int g_jit_math_fast_sqrt_table[0x10000]; 
static long g_jit_math_fast_sqrt_table_init = FALSE;

typedef union _jit_math_fast_sqrt_union
{
  float f;
  unsigned int i;
} t_jit_math_fast_sqrt_union;

void jit_math_fast_sqrt_table_init()
{
  unsigned int i;
  t_jit_math_fast_sqrt_union s;
  
  for (i = 0; i <= 0x7FFF; i++)
  {
        // build a float with the bit pattern i as mantissa and an exponent of 0, stored as 127
        s.i = (i << 8) | (0x7F << 23);
        s.f = (float)jit_math_sqrt(s.f);
    
        // take the square root then strip the first 7 bits of the mantissa into the table
        g_jit_math_fast_sqrt_table[i + 0x8000] = (s.i & 0x7FFFFF);
    
        // repeat the process, this time with an exponent of 1, stored as 128
        s.i = (i << 8) | (0x80 << 23);
        s.f = (float)jit_math_sqrt(s.f);
    
        g_jit_math_fast_sqrt_table[i] = (s.i & 0x7FFFFF);
    }

    g_jit_math_fast_sqrt_table_init = TRUE;
}

/**
 * Calculates the square root by fast approximation.
 *
 * @ingroup mathmod 
 *
 * @param   n       input
 * @return  output 
 *
 */
float jit_math_fast_sqrt(float n)
{
    if(g_jit_math_fast_sqrt_table_init == FALSE)
        jit_math_fast_sqrt_table_init();

    if (JIT_MATH_F32_BITS(n) == 0)
        return 0.0f;
  
    JIT_MATH_F32_BITS(n) = g_jit_math_fast_sqrt_table[(JIT_MATH_F32_BITS(n) >> 8) & 0xFFFF] | ((((JIT_MATH_F32_BITS(n) - 0x3F800000) >> 1) + 0x3F800000) & 0x7F800000); 
    return n;
}

/**
 * Calculates the inverse square root by fast approximation.
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_invsqrt(float x)
{
    float half = 0.5f*x;
    int i  = *(int*)&x;
    i = 0x5f3759df - (i >> 1);
    x = *(float*)&i;
    x = x*(1.5f - half*x*x);
    return x;
}

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

/*
 *  Derek Gerstmann - derek@cycling74.com
 *
 *  The following polynomial approximations for the various trig
 *  functions use constants derived by David Eberly of 
 *  Magic Software in his free Wild Magic Library:
 *
 *  http://www.magic-software.com/Math.html
 *
 *  -[dg]
 */

/**
 * Calculates the sine by fast approximation.
 *
 * Absolute error of 1.7e-04 for [0, PI/2]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_sin(float x)
{
    float x2 = x*x;
    float rv = (float)7.61e-03;
    rv *= x2;
    rv -= (float)1.6605e-01;
    rv *= x2;
    rv += (float)1.0;
    rv *= x;
    return rv;
}

/**
 * Calculates the cosine by fast approximation.
 *
 * Absolute error of 1.2e-03 for [0, PI/2]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_cos(float x)
{
    float x2 = x*x;
    float rv = (float)3.705e-02;
    rv *= x2;
    rv -= (float)4.967e-01;
    rv *= x2;
    rv += (float)1.0;
    return rv;
}

/**
 * Calculates the tangent by fast approximation.
 *
 * Absolute error of 1.9e-00 for [0, PI/4]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_tan(float x)
{
    float x2 = x*x;
    float rv = (float)9.5168091e-03;
    rv *= x2;
    rv += (float)2.900525e-03;
    rv *= x2;
    rv += (float)2.45650893e-02;
    rv *= x2;
    rv += (float)5.33740603e-02;
    rv *= x2;
    rv += (float)1.333923995e-01;
    rv *= x2;
    rv += (float)3.333314036e-01;
    rv *= x2;
    rv += (float)1.0;
    rv *= x;
    return rv;
}

/**
 * Calculates the arcsine by fast approximation.
 *
 * Absolute error of 6.8e-05 for [0, 1]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_asin(float x)
{
    float root = jit_math_sqrt(((float)1.0)-x);
    float rv = -(float)0.0187293;
    rv *= x;
    rv += (float)0.0742610;
    rv *= x;
    rv -= (float)0.2121144;
    rv *= x;
    rv += (float)1.5707288;
    rv = JIT_MATH_F32_HALF_PI - root*rv;
    return rv;
}

/**
 * Calculates the arccosine by fast approximation.
 *
 * Absolute error of 6.8e-05 for [0, 1]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_acos(float x)
{
    float root = jit_math_sqrt(((float)1.0)-x);
    float rv = -(float)0.0187293;
    rv *= x;
    rv += (float)0.0742610;
    rv *= x;
    rv -= (float)0.2121144;
    rv *= x;
    rv += (float)1.5707288;
    rv *= root;
    return rv;
}

/**
 * Calculates the arctangent by fast approximation.
 *
 * Absolute error of 1.43-08 for [-1, 1]
 *
 * @ingroup mathmod 
 *
 * @param   x       input
 * @return  output 
 *
 */
float jit_math_fast_atan(float x)
{
    float x2 = x*x;
    float rv = (float)0.0028662257;
    rv *= x2;
    rv -= (float)0.0161657367;
    rv *= x2;
    rv += (float)0.0429096138;
    rv *= x2;
    rv -= (float)0.0752896400;
    rv *= x2;
    rv += (float)0.1065626393;
    rv *= x2;
    rv -= (float)0.1420889944;
    rv *= x2;
    rv += (float)0.1999355085;
    rv *= x2;
    rv -= (float)0.3333314528;
    rv *= x2;
    rv += (float)1.0;
    rv *= x;
    return rv;
}

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

#ifdef WIN_VERSION
double acosh(double x)
{
    return log(x + sqrt(x * x - 1));
}

double asinh(double x)
{
    return log(x + sqrt(x * x + 1));
}

double atanh(double x)
{
    return log((1 + x) / (1 - x)) / 2;
}

double expm1(double x)
{
    return exp(x)-1;
}

double exp2(double x)
{
    return pow(2,x);
}

double log2(double x)
{
    return log(x)/log(2);
}

double round(double x)
{
    long c=x;
    if (x<0)
        return ((x-c)<=-0.5)?c-1:c;
    else 
        return ((x-c)>=0.5)?c+1:c;
}

double trunc(double x)
{
    long c=x;
    return c;
}

#endif //WIN_VERSION

Copyright © 2008, Cycling '74