Fast Function Approximations lowering.

Makefile

@@ -421,6 +421,7 @@ SOURCE_FILES = \
   AlignLoads.cpp \
   AllocationBoundsInference.cpp \
   ApplySplit.cpp \
+  ApproximationTables.cpp \
   Argument.cpp \
   AssociativeOpsTable.cpp \
   Associativity.cpp \
@@ -479,6 +480,7 @@ SOURCE_FILES = \
   Expr.cpp \
   ExtractTileOperations.cpp \
   FastIntegerDivide.cpp \
+  FastMathFunctions.cpp \
   FindCalls.cpp \
   FindIntrinsics.cpp \
   FlattenNestedRamps.cpp \

src/ApproximationTables.h

@@ -0,0 +1,32 @@
+#ifndef HALIDE_APPROXIMATION_TABLES_H
+#define HALIDE_APPROXIMATION_TABLES_H
+
+#include <vector>
+
+#include "IROperator.h"
+
+namespace Halide {
+namespace Internal {
+
+struct Approximation {
+    ApproximationPrecision::OptimizationObjective objective;
+    struct Metrics {
+        double mse;
+        double mae;
+        double mulpe;
+    } metrics_f32, metrics_f64;
+    std::vector<double> coefficients;
+};
+
+const Approximation *best_atan_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_sin_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_cos_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_tan_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_log_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_exp_approximation(Halide::ApproximationPrecision precision, Type type);
+const Approximation *best_expm1_approximation(Halide::ApproximationPrecision precision, Type type);
+
+}  // namespace Internal
+}  // namespace Halide
+
+#endif

src/CMakeLists.txt

@@ -115,6 +115,7 @@ target_sources(
     ExternFuncArgument.h
     ExtractTileOperations.h
     FastIntegerDivide.h
+    FastMathFunctions.h
     FindCalls.h
     FindIntrinsics.h
     FlattenNestedRamps.h
@@ -220,8 +221,7 @@ target_sources(
     WrapCalls.h
 )
 
-# The sources that go into libHalide. For the sake of IDE support, headers that
-# exist in src/ but are not public should be included here.
+# The sources that go into libHalide.
 target_sources(
     Halide
     PRIVATE
@@ -233,6 +233,7 @@ target_sources(
     AlignLoads.cpp
     AllocationBoundsInference.cpp
     ApplySplit.cpp
+    ApproximationTables.cpp
     Argument.cpp
     AssociativeOpsTable.cpp
     Associativity.cpp
@@ -291,6 +292,7 @@ target_sources(
     Expr.cpp
     ExtractTileOperations.cpp
     FastIntegerDivide.cpp
+    FastMathFunctions.cpp
     FindCalls.cpp
     FindIntrinsics.cpp
     FlattenNestedRamps.cpp

src/CSE.cpp

@@ -33,6 +33,12 @@ bool should_extract(const Expr &e, bool lift_all) {
         return false;
     }
 
+    if (const Call *c = e.as<Call>()) {
+        if (c->type == type_of<ApproximationPrecision *>()) {
+            return false;
+        }
+    }
+
     if (lift_all) {
         return true;
     }

src/CodeGen_Metal_Dev.cpp

@@ -858,6 +858,13 @@ void CodeGen_Metal_Dev::init_module() {
                << "#define acosh_f16 acosh\n"
                << "#define tanh_f16 tanh\n"
                << "#define atanh_f16 atanh\n"
+               << "#define fast_sin_f32 fast::sin \n"
+               << "#define fast_cos_f32 fast::cos \n"
+               << "#define fast_tan_f32 fast::tan \n"
+               << "#define fast_exp_f32 fast::exp \n"
+               << "#define fast_log_f32 fast::log \n"
+               << "#define fast_pow_f32 fast::pow \n"
+               << "#define fast_tanh_f32 fast::tanh \n"
                << "#define fast_inverse_sqrt_f16 rsqrt\n"
                << "constexpr half half_from_bits(unsigned short x) {return as_type<half>(x);}\n"
                << "constexpr half nan_f16() { return half_from_bits(32767); }\n"

src/CodeGen_OpenCL_Dev.cpp

@@ -1148,6 +1148,12 @@ void CodeGen_OpenCL_Dev::init_module() {
                << "#define acosh_f32 acosh \n"
                << "#define tanh_f32 tanh \n"
                << "#define atanh_f32 atanh \n"
+               << "#define fast_sin_f32 native_sin \n"
+               << "#define fast_cos_f32 native_cos \n"
+               << "#define fast_tan_f32 native_tan \n"
+               << "#define fast_exp_f32 native_exp \n"
+               << "#define fast_log_f32 native_log \n"
+               << "#define fast_pow_f32 native_powr \n"
                << "#define fast_inverse_f32 native_recip \n"
                << "#define fast_inverse_sqrt_f32 native_rsqrt \n";
 

src/CodeGen_PTX_Dev.cpp

@@ -572,7 +572,7 @@ string CodeGen_PTX_Dev::mattrs() const {
         return "+ptx70";
     } else if (target.has_feature(Target::CUDACapability70) ||
                target.has_feature(Target::CUDACapability75)) {
-        return "+ptx60";
+        return "+ptx70";
     } else if (target.has_feature(Target::CUDACapability61)) {
         return "+ptx50";
     } else if (target.features_any_of({Target::CUDACapability32,
@@ -728,7 +728,7 @@ vector<char> CodeGen_PTX_Dev::compile_to_src() {
     if (debug::debug_level() >= 2) {
         dump();
     }
-    debug(2) << "Done with CodeGen_PTX_Dev::compile_to_src";
+    debug(2) << "Done with CodeGen_PTX_Dev::compile_to_src\n";
 
     debug(1) << "PTX kernel:\n"
              << outstr.c_str() << "\n";

src/Derivative.cpp

@@ -30,12 +30,20 @@ using FuncKey = Derivative::FuncKey;
 namespace Internal {
 namespace {
 
-bool is_float_extern(const string &op_name,
-                     const string &func_name) {
-    return op_name == (func_name + "_f16") ||
-           op_name == (func_name + "_f32") ||
-           op_name == (func_name + "_f64");
-};
+bool is_math_func(const Call *call,
+                  const string &func_name,
+                  Call::IntrinsicOp intrinsic_op = Call::IntrinsicOp::IntrinsicOpCount) {
+    if (call->is_extern()) {
+        const string &op_name = call->name;
+        return op_name == (func_name + "_f16") ||
+               op_name == (func_name + "_f32") ||
+               op_name == (func_name + "_f64");
+    } else if (call->is_intrinsic() && intrinsic_op != Call::IntrinsicOpCount) {
+        return call->is_intrinsic(intrinsic_op);
+    } else {
+        return false;
+    }
+}
 
 /** Compute derivatives through reverse accumulation
  */
@@ -1058,101 +1066,99 @@ void ReverseAccumulationVisitor::visit(const Select *op) {
 void ReverseAccumulationVisitor::visit(const Call *op) {
     internal_assert(expr_adjoints.find(op) != expr_adjoints.end());
     Expr adjoint = expr_adjoints[op];
-    if (op->is_extern()) {
-        // Math functions
-        if (is_float_extern(op->name, "exp")) {
-            // d/dx exp(x) = exp(x)
-            accumulate(op->args[0], adjoint * exp(op->args[0]));
-        } else if (is_float_extern(op->name, "log")) {
-            // d/dx log(x) = 1 / x
-            accumulate(op->args[0], adjoint / op->args[0]);
-        } else if (is_float_extern(op->name, "sin")) {
-            // d/dx sin(x) = cos(x)
-            accumulate(op->args[0], adjoint * cos(op->args[0]));
-        } else if (is_float_extern(op->name, "asin")) {
-            // d/dx asin(x) = 1 / sqrt(1 - x^2)
-            Expr one = make_one(op->type);
-            accumulate(op->args[0], adjoint / sqrt(one - op->args[0] * op->args[0]));
-        } else if (is_float_extern(op->name, "cos")) {
-            // d/dx cos(x) = -sin(x)
-            accumulate(op->args[0], -adjoint * sin(op->args[0]));
-        } else if (is_float_extern(op->name, "acos")) {
-            // d/dx acos(x) = - 1 / sqrt(1 - x^2)
-            Expr one = make_one(op->type);
-            accumulate(op->args[0], -adjoint / sqrt(one - op->args[0] * op->args[0]));
-        } else if (is_float_extern(op->name, "tan")) {
-            // d/dx tan(x) = 1 / cos(x)^2
-            Expr c = cos(op->args[0]);
-            accumulate(op->args[0], adjoint / (c * c));
-        } else if (is_float_extern(op->name, "atan")) {
-            // d/dx atan(x) = 1 / (1 + x^2)
-            Expr one = make_one(op->type);
-            accumulate(op->args[0], adjoint / (one + op->args[0] * op->args[0]));
-        } else if (is_float_extern(op->name, "atan2")) {
-            Expr x2y2 = op->args[0] * op->args[0] + op->args[1] * op->args[1];
-            // d/dy atan2(y, x) = x / (x^2 + y^2)
-            accumulate(op->args[0], adjoint * (op->args[1] / x2y2));
-            // d/dx atan2(y, x) = -y / (x^2 + y^2)
-            accumulate(op->args[1], adjoint * (-op->args[0] / x2y2));
-        } else if (is_float_extern(op->name, "sinh")) {
-            // d/dx sinh(x) = cosh(x)
-            accumulate(op->args[0], adjoint * cosh(op->args[0]));
-        } else if (is_float_extern(op->name, "asinh")) {
-            // d/dx asin(x) = 1 / sqrt(1 + x^2)
-            Expr one = make_one(op->type);
-            accumulate(op->args[0], adjoint / sqrt(one + op->args[0] * op->args[0]));
-        } else if (is_float_extern(op->name, "cosh")) {
-            // d/dx cosh(x) = sinh(x)
-            accumulate(op->args[0], adjoint * sinh(op->args[0]));
-        } else if (is_float_extern(op->name, "acosh")) {
-            // d/dx acosh(x) = 1 / (sqrt(x - 1) sqrt(x + 1)))
-            Expr one = make_one(op->type);
-            accumulate(op->args[0],
-                       adjoint / (sqrt(op->args[0] - one) * sqrt(op->args[0] + one)));
-        } else if (is_float_extern(op->name, "tanh")) {
-            // d/dx tanh(x) = 1 / cosh(x)^2
-            Expr c = cosh(op->args[0]);
-            accumulate(op->args[0], adjoint / (c * c));
-        } else if (is_float_extern(op->name, "atanh")) {
-            // d/dx atanh(x) = 1 / (1 - x^2)
-            Expr one = make_one(op->type);
-            accumulate(op->args[0], adjoint / (one - op->args[0] * op->args[0]));
-        } else if (is_float_extern(op->name, "ceil")) {
-            // TODO: d/dx = dirac(n) for n in Z ...
-            accumulate(op->args[0], make_zero(op->type));
-        } else if (is_float_extern(op->name, "floor")) {
-            // TODO: d/dx = dirac(n) for n in Z ...
-            accumulate(op->args[0], make_zero(op->type));
-        } else if (is_float_extern(op->name, "round")) {
-            accumulate(op->args[0], make_zero(op->type));
-        } else if (is_float_extern(op->name, "trunc")) {
-            accumulate(op->args[0], make_zero(op->type));
-        } else if (is_float_extern(op->name, "sqrt")) {
-            Expr half = make_const(op->type, 0.5);
-            accumulate(op->args[0], adjoint * (half / sqrt(op->args[0])));
-        } else if (is_float_extern(op->name, "pow")) {
-            Expr one = make_one(op->type);
-            accumulate(op->args[0],
-                       adjoint * op->args[1] * pow(op->args[0], op->args[1] - one));
-            accumulate(op->args[1],
-                       adjoint * pow(op->args[0], op->args[1]) * log(op->args[0]));
-        } else if (is_float_extern(op->name, "fast_inverse")) {
-            // d/dx 1/x = -1/x^2
-            Expr inv_x = fast_inverse(op->args[0]);
-            accumulate(op->args[0], -adjoint * inv_x * inv_x);
-        } else if (is_float_extern(op->name, "fast_inverse_sqrt")) {
-            // d/dx x^(-0.5) = -0.5*x^(-1.5)
-            Expr inv_sqrt_x = fast_inverse_sqrt(op->args[0]);
-            Expr neg_half = make_const(op->type, -0.5);
-            accumulate(op->args[0],
-                       neg_half * adjoint * inv_sqrt_x * inv_sqrt_x * inv_sqrt_x);
-        } else if (op->name == "halide_print") {
-            for (const auto &arg : op->args) {
-                accumulate(arg, make_zero(op->type));
-            }
-        } else {
-            internal_error << "The derivative of " << op->name << " is not implemented.";
+    // Math functions (Can be both intrinsic and extern).
+    if (is_math_func(op, "exp", Call::fast_exp)) {
+        // d/dx exp(x) = exp(x)
+        accumulate(op->args[0], adjoint * exp(op->args[0]));
+    } else if (is_math_func(op, "log", Call::fast_log)) {
+        // d/dx log(x) = 1 / x
+        accumulate(op->args[0], adjoint / op->args[0]);
+    } else if (is_math_func(op, "sin", Call::fast_sin)) {
+        // d/dx sin(x) = cos(x)
+        accumulate(op->args[0], adjoint * cos(op->args[0]));
+    } else if (is_math_func(op, "asin")) {
+        // d/dx asin(x) = 1 / sqrt(1 - x^2)
+        Expr one = make_one(op->type);
+        accumulate(op->args[0], adjoint / sqrt(one - op->args[0] * op->args[0]));
+    } else if (is_math_func(op, "cos", Call::fast_cos)) {
+        // d/dx cos(x) = -sin(x)
+        accumulate(op->args[0], -adjoint * sin(op->args[0]));
+    } else if (is_math_func(op, "acos")) {
+        // d/dx acos(x) = - 1 / sqrt(1 - x^2)
+        Expr one = make_one(op->type);
+        accumulate(op->args[0], -adjoint / sqrt(one - op->args[0] * op->args[0]));
+    } else if (is_math_func(op, "tan", Call::fast_tan)) {
+        // d/dx tan(x) = 1 / cos(x)^2
+        Expr c = cos(op->args[0]);
+        accumulate(op->args[0], adjoint / (c * c));
+    } else if (is_math_func(op, "atan", Call::fast_atan)) {
+        // d/dx atan(x) = 1 / (1 + x^2)
+        Expr one = make_one(op->type);
+        accumulate(op->args[0], adjoint / (one + op->args[0] * op->args[0]));
+    } else if (is_math_func(op, "atan2", Call::fast_atan2)) {
+        Expr x2y2 = op->args[0] * op->args[0] + op->args[1] * op->args[1];
+        // d/dy atan2(y, x) = x / (x^2 + y^2)
+        accumulate(op->args[0], adjoint * (op->args[1] / x2y2));
+        // d/dx atan2(y, x) = -y / (x^2 + y^2)
+        accumulate(op->args[1], adjoint * (-op->args[0] / x2y2));
+    } else if (is_math_func(op, "sinh")) {
+        // d/dx sinh(x) = cosh(x)
+        accumulate(op->args[0], adjoint * cosh(op->args[0]));
+    } else if (is_math_func(op, "asinh")) {
+        // d/dx asin(x) = 1 / sqrt(1 + x^2)
+        Expr one = make_one(op->type);
+        accumulate(op->args[0], adjoint / sqrt(one + op->args[0] * op->args[0]));
+    } else if (is_math_func(op, "cosh")) {
+        // d/dx cosh(x) = sinh(x)
+        accumulate(op->args[0], adjoint * sinh(op->args[0]));
+    } else if (is_math_func(op, "acosh")) {
+        // d/dx acosh(x) = 1 / (sqrt(x - 1) sqrt(x + 1)))
+        Expr one = make_one(op->type);
+        accumulate(op->args[0],
+                   adjoint / (sqrt(op->args[0] - one) * sqrt(op->args[0] + one)));
+    } else if (is_math_func(op, "tanh", Call::fast_tanh)) {
+        // d/dx tanh(x) = 1 / cosh(x)^2
+        Expr c = cosh(op->args[0]);
+        accumulate(op->args[0], adjoint / (c * c));
+    } else if (is_math_func(op, "atanh")) {
+        // d/dx atanh(x) = 1 / (1 - x^2)
+        Expr one = make_one(op->type);
+        accumulate(op->args[0], adjoint / (one - op->args[0] * op->args[0]));
+    } else if (is_math_func(op, "ceil")) {
+        // TODO: d/dx = dirac(n) for n in Z ...
+        accumulate(op->args[0], make_zero(op->type));
+    } else if (is_math_func(op, "floor")) {
+        // TODO: d/dx = dirac(n) for n in Z ...
+        accumulate(op->args[0], make_zero(op->type));
+    } else if (is_math_func(op, "round")) {
+        accumulate(op->args[0], make_zero(op->type));
+    } else if (is_math_func(op, "trunc")) {
+        accumulate(op->args[0], make_zero(op->type));
+    } else if (is_math_func(op, "sqrt")) {
+        Expr half = make_const(op->type, 0.5);
+        accumulate(op->args[0], adjoint * (half / sqrt(op->args[0])));
+    } else if (is_math_func(op, "pow", Call::fast_pow)) {
+        Expr one = make_one(op->type);
+        accumulate(op->args[0],
+                   adjoint * op->args[1] * pow(op->args[0], op->args[1] - one));
+        accumulate(op->args[1],
+                   adjoint * pow(op->args[0], op->args[1]) * log(op->args[0]));
+    } else if (is_math_func(op, "fast_inverse")) {
+        // d/dx 1/x = -1/x^2
+        Expr inv_x = fast_inverse(op->args[0]);
+        accumulate(op->args[0], -adjoint * inv_x * inv_x);
+    } else if (is_math_func(op, "fast_inverse_sqrt")) {
+        // d/dx x^(-0.5) = -0.5*x^(-1.5)
+        Expr inv_sqrt_x = fast_inverse_sqrt(op->args[0]);
+        Expr neg_half = make_const(op->type, -0.5);
+        accumulate(op->args[0],
+                   neg_half * adjoint * inv_sqrt_x * inv_sqrt_x * inv_sqrt_x);
+    } else if (op->is_extern() && op->name == "halide_print") {
+        for (const auto &arg : op->args) {
+            accumulate(arg, make_zero(op->type));
         }
+    } else if (op->is_extern()) {
+        internal_error << "The derivative of " << op->name << " is not implemented.";
     } else if (op->is_intrinsic()) {
         if (op->is_intrinsic(Call::abs)) {
             accumulate(op->args[0],

src/Error.cpp

@@ -106,7 +106,7 @@ namespace Internal {
 
 void unhandled_exception_handler() {
     // Note that we use __cpp_exceptions (rather than HALIDE_WITH_EXCEPTIONS)
-    // to maximize the change of dealing with uncaught exceptions in weird
+    // to maximize the chance of dealing with uncaught exceptions in weird
     // build situations (i.e., exceptions enabled via C++ but HALIDE_WITH_EXCEPTIONS
     // is somehow not set).
 #ifdef __cpp_exceptions