diff --git a/src/core/ulib_math.cpp b/src/core/ulib_math.cpp index 29fe61f7f..23e206531 100644 --- a/src/core/ulib_math.cpp +++ b/src/core/ulib_math.cpp @@ -157,6 +157,52 @@ DLL_QUERY libmath_query( Chuck_DL_Query * QUERY ) QUERY->add_arg( QUERY, "float", "y" ); QUERY->doc_func( QUERY, "Compute the euclidean distance sqrt(x*x+y*y)." ); + // Fast versions of trig functions!!! + // Warning: not always accurate + + // sin + QUERY->add_sfun( QUERY, fast_sin_impl, "float", "ssin" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of sine (measured in radians). Uses a Padé approximant, and is recommended to only use inputs between -pi and +pi." ); + + // cos + QUERY->add_sfun( QUERY, fast_cos_impl, "float", "scos" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of cosine (measured in radians). Uses a Padé approximant, and is recommended to only use inputs between -pi and +pi." ); + + // tan + QUERY->add_sfun( QUERY, fast_tan_impl, "float", "stan" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of tangent (measured in radians). Uses a Padé approximant, and is recommended to only use inputs between -pi/2 and +pi/2." ); + + // sinh + QUERY->add_sfun( QUERY, fast_sinh_impl, "float", "ssinh" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of hyperbolic sine. Uses a Padé approximant, and is recommended to only use inputs between -5 and +5." ); + + // cosh + QUERY->add_sfun( QUERY, fast_cosh_impl, "float", "scosh" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of the hyperbolic cosine. Uses a Padé approximant, and is recommended to only use inputs between -5 and +5." ); + + // tanh + QUERY->add_sfun( QUERY, fast_tanh_impl, "float", "stanh" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of the hyperbolic tangent. Uses a Padé approximant, and is recommended to only use inputs between -5 and +5." ); + + // exp + QUERY->add_sfun( QUERY, fast_exp_impl, "float", "sexp" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "Speedy approximation of e^x, the base-e exponential. Uses a Padé approximant, and is recommended to only use inputs between -6 and +4." ); + + // fast inverse-square root + QUERY->add_sfun( QUERY, fast_inv_sqrt, "float", "sinsqrt" ); + QUERY->add_arg( QUERY, "float", "x" ); + QUERY->doc_func( QUERY, "A speedy (fast) approximation of the inverse square root. Popularized by Quake" ); + + + // End of fast implementations of trig functions + // pow QUERY->add_sfun( QUERY, pow_impl, "float", "pow" ); QUERY->add_arg( QUERY, "float", "x" ); @@ -602,6 +648,104 @@ CK_DLL_SFUN( tanh_impl ) RETURN->v_float = ::tanh( GET_CK_FLOAT(ARGS) ); } +// Fast versions of trig functions!! + +// sin +CK_DLL_SFUN( fast_sin_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = -x * (-11511339840 + x2 * (1640635920 + x2 * (-52785432 + x2 * 479249))); + t_CKFLOAT denominator = 11511339840 + x2 * (277920720 + x2 * (3177720 + x2 * 18361)); + RETURN->v_float = numerator / denominator; +} + +// cos +CK_DLL_SFUN( fast_cos_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = -(-39251520 + x2 * (18471600 + x2 * (-1075032 + 14615 * x2))); + t_CKFLOAT denominator = 39251520 + x2 * (1154160 + x2 * (16632 + x2 * 127)); + RETURN->v_float = numerator / denominator; +} + +// tan +CK_DLL_SFUN( fast_tan_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = x * (-135135 + x2 * (17325 + x2 * (-378 + x2))); + t_CKFLOAT denominator = -135135 + x2 * (62370 + x2 * (-3150 + 28 * x2)); + RETURN->v_float = numerator / denominator; +} + + +// sinh +CK_DLL_SFUN( fast_sinh_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = -x * (11511339840 + x2 * (1640635920 + x2 * (52785432 + x2 * 479249))); + t_CKFLOAT denominator = -11511339840 + x2 * (277920720 + x2 * (-3177720 + x2 * 18361)); + RETURN->v_float = numerator / denominator; +} + +// cosh +CK_DLL_SFUN( fast_cosh_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = -(39251520 + x2 * (18471600 + x2 * (1075032 + 14615 * x2))); + t_CKFLOAT denominator = -39251520 + x2 * (1154160 + x2 * (-16632 + 127 * x2)); + RETURN->v_float = numerator / denominator; +} + +// tanh +CK_DLL_SFUN( fast_tanh_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT x2 = x * x; + t_CKFLOAT numerator = x * (135135 + x2 * (17325 + x2 * (378 + x2))); + t_CKFLOAT denominator = 135135 + x2 * (62370 + x2 * (3150 + 28 * x2)); + RETURN->v_float = numerator / denominator; +} + +// exp +CK_DLL_SFUN( fast_exp_impl ) +{ + t_CKFLOAT x = GET_CK_FLOAT(ARGS); + t_CKFLOAT numerator = 1680 + x * (840 + x * (180 + x * (20 + x))); + t_CKFLOAT denominator = 1680 + x *(-840 + x * (180 + x * (-20 + x))); + RETURN->v_float = numerator / denominator; +} + + +// fast inverse square root (of quake fame) +CK_DLL_SFUN( fast_inv_sqrt ) +{ + t_CKFLOAT num = GET_CK_FLOAT(ARGS); + + int i; + float x, y; + x = num * 0.5; + y = num; + i = *(int*)&y; // copy bits from y into i so we can bit shift + + i = 0x5f3759df - (i >> 1); // original, has .175% error + i = 0x5F375A86 - (i >> 1); // Chris Lomont's updated magic number, .09% error after 1 newton iteration + + y = *(float*)&i; // copy back into float + + // newton iteration + y = (y * (1.5 - (x * y * y))); + y = (y * (1.5 - (x * y * y))); // optional, improves accuracy + RETURN->v_float = num * y; +} + + +// end of fast versions of trig functions + // hypot CK_DLL_SFUN( hypot_impl ) { diff --git a/src/core/ulib_math.h b/src/core/ulib_math.h index 4a00d79b3..6893d2bf7 100644 --- a/src/core/ulib_math.h +++ b/src/core/ulib_math.h @@ -52,6 +52,18 @@ CK_DLL_SFUN( cosh_impl ); CK_DLL_SFUN( tanh_impl ); CK_DLL_SFUN( hypot_impl ); +// fast approximations of the above trig functions +CK_DLL_SFUN( fast_sin_impl ); +CK_DLL_SFUN( fast_cos_impl ); +CK_DLL_SFUN( fast_tan_impl ); +CK_DLL_SFUN( fast_cot_impl ); +CK_DLL_SFUN( fast_sinh_impl ); +CK_DLL_SFUN( fast_cosh_impl ); +CK_DLL_SFUN( fast_tanh_impl ); +CK_DLL_SFUN( fast_exp_impl ); +CK_DLL_SFUN( fast_inv_sqrt ); + +// Normal not-fast not-approximations CK_DLL_SFUN( pow_impl ); CK_DLL_SFUN( sqrt_impl ); CK_DLL_SFUN( exp_impl ); diff --git a/src/test/01-Basic/276-speedy-trig.ck b/src/test/01-Basic/276-speedy-trig.ck new file mode 100644 index 000000000..c50dbdc9f --- /dev/null +++ b/src/test/01-Basic/276-speedy-trig.ck @@ -0,0 +1,84 @@ +0.01 => float inc; + + +-1 * (Math.PI / 2) => float i; +Math.PI / 2 => float dest; + +while (i < dest) { + testEquality(i, Math.sin(i), Math.ssin(i), "sin"); + inc +=> i; +} + +-1 * (Math.PI / 2) => i; +Math.PI / 2 => dest; + +while (i < dest) { + testEquality(i, Math.cos(i), Math.scos(i), "cos"); + inc +=> i; +} + +-1 * (Math.PI / 2) + 0.2 => i; +Math.PI / 2 - 0.1 => dest; + +while (i < dest) { + testEquality(i, Math.tan(i), Math.stan(i), "tan"); + inc +=> i; +} + +-1 * (Math.PI / 2) => i; +Math.PI / 2 => dest; + +while (i < dest) { + testEquality(i, Math.sinh(i), Math.ssinh(i), "sinh"); + inc +=> i; +} + + +-1 * (Math.PI / 2) => i; +Math.PI / 2 => dest; + +while (i < dest) { + testEquality(i, Math.cosh(i), Math.scosh(i), "cosh"); + inc +=> i; +} + +-1 * (Math.PI / 2) => i; +Math.PI / 2 => dest; + +while (i < dest) { + testEquality(i, Math.tanh(i), Math.stanh(i), "tanh"); + inc +=> i; +} + +-1 * (Math.PI / 2) => i; +Math.PI / 2 - 0.35 => dest; + +while (i < dest) { + testEquality(i, Math.exp(i), Math.sexp(i), "exp"); + inc +=> i; +} + +<<< "success" >>>; + + +fun int testEquality(float input, float x, float y, string fname) { + if (!fequals(x,y)) { + cherr <= "ERROR: speed approximation of " <= fname <= " are too different at input " + <= input <= ". Expected " + <= x <= ", but got " <= y <= IO.newline(); + return false; + } + return true; +} + +// custom float approximate equality test +fun int fequals(float x, float y) { + 0.000001 => float epsilon; // our error bound + + Math.fabs(x) => x; + Math.fabs(y) => y; + + x-y => Math.fabs => float diff; + + return diff < epsilon; +} \ No newline at end of file