Speedy trig functions

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 )
 {

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 );

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;
+}