do some optimizations

src/sdf/context.cpp

@@ -82,9 +82,7 @@ Operand Context::addInstruction(Instruction inst) {
   return result;
 }
 
-// bypass the cache because we don't expect to have more common subexpressions
-// after optimizations
-Operand Context::addInstructionNoCache(Instruction inst) {
+std::optional<Operand> Context::trySimplify(Instruction inst) {
   // constant choice
   auto op = inst.op;
   auto &operands = inst.operands;
@@ -161,11 +159,68 @@ Operand Context::addInstructionNoCache(Instruction inst) {
     return addConstant(result);
   }
 
+  // simple simplifications
+  if (op == OpCode::ADD) {
+    // add is commutative, so if there is a constant, it must be on the left
+    // 0 + x => x
+    if (operands[0].isConst() && constants[operands[0].toConstIndex()] == 0.0)
+      return operands[1];
+  }
+  if (op == OpCode::SUB) {
+    // x - 0 => x
+    if (operands[1].isConst() && constants[operands[1].toConstIndex()] == 0.0)
+      return operands[0];
+  }
+  if (op == OpCode::MUL) {
+    // mul is commutative, so if there is a constant, it must be on the left
+    // 0 * x => 0
+    if (operands[0].isConst() && constants[operands[0].toConstIndex()] == 0.0)
+      return operands[0];
+    // 1 * x => x
+    if (operands[0].isConst() && constants[operands[0].toConstIndex()] == 1.0)
+      return operands[1];
+  }
+  if (op == OpCode::DIV) {
+    if (operands[1].isConst() && constants[operands[1].toConstIndex()] == 1.0)
+      return operands[0];
+  }
+
+  return {};
+}
+
+Instruction Context::strengthReduction(Instruction inst) {
+  // strength reduction: reduce instructions to simpler variants
+  // not very helpful for point evaluation in a vm because instruction decoding
+  // is the most time consuming part.
+  // This can be useful if we want to do JIT, interval evaluation or bulk
+  // evaluation.
+  if (inst.op == OpCode::MUL && inst.operands[1].isConst() &&
+      constants[inst.operands[1].toConstIndex()] == 2.0) {
+    // x * 2 => x + x
+    return {OpCode::ADD, {inst.operands[0], inst.operands[1], Operand::none()}};
+  }
+  if (inst.op == OpCode::DIV && inst.operands[1].isConst()) {
+    // x / c => x * (1/c)
+    return {OpCode::MUL,
+            {inst.operands[0],
+             addConstant(1.0 / constants[inst.operands[1].toConstIndex()]),
+             Operand::none()}};
+  }
+  return inst;
+}
+
+// bypass the cache because we don't expect to have more common subexpressions
+// after optimizations
+Operand Context::addInstructionNoCache(Instruction inst) {
+  auto simplified = trySimplify(inst);
+  if (simplified.has_value()) return simplified.value();
+  inst = strengthReduction(inst);
+
   size_t i = instructions.size();
-  instructions.push_back({op, operands});
+  instructions.push_back(inst);
   opUses.emplace_back();
   // update uses
-  for (auto operand : operands) {
+  for (auto operand : inst.operands) {
     auto target = getUses(operand);
     if (target == nullptr) continue;
     // avoid duplicates

src/sdf/context.h

@@ -13,6 +13,7 @@
 // limitations under the License.
 #pragma once
 
+#include <optional>
 #include <utility>
 #include <vector>
 
@@ -104,6 +105,8 @@ class Context {
   std::vector<UsesVector> opUses;
   unordered_map<Instruction, Operand> cache;
 
+  std::optional<Operand> trySimplify(Instruction);
+  Instruction strengthReduction(Instruction);
   Operand addInstructionNoCache(Instruction);
 
   UsesVector* getUses(Operand operand) {

src/sdf/interval.h

@@ -27,6 +27,9 @@ struct Interval {
   Domain lower;
   Domain upper;
 
+  static constexpr Domain zero = static_cast<Domain>(0);
+  static constexpr Domain one = static_cast<Domain>(1);
+
   Interval()
       : lower(-std::numeric_limits<Domain>::infinity()),
         upper(std::numeric_limits<Domain>::infinity()) {}
@@ -42,14 +45,22 @@ struct Interval {
 
   Interval operator-(const Interval &other) const { return *this + (-other); }
 
+  Interval operator*(Domain d) const {
+    if (d > zero) return {lower * d, upper * d};
+    return {upper * d, lower * d};
+  }
+
   Interval operator*(const Interval &other) const {
+    if (is_const()) return other * lower;
+    if (other.is_const()) return *this * other.lower;
+
     Domain a1b1 = lower * other.lower;
     Domain a2b2 = upper * other.upper;
     // we can write more "fast paths", but at some point it will become slower
     // than just going the general path...
-    if (lower >= 0.0 && other.lower >= 0.0)
+    if (lower >= zero && other.lower >= zero)
       return {a1b1, a2b2};
-    else if (upper <= 0.0 && other.upper <= 0.0)
+    else if (upper <= zero && other.upper <= zero)
       return {a2b2, a1b1};
 
     Domain a1b2 = lower * other.upper;
@@ -58,53 +69,47 @@ struct Interval {
             std::max(std::max(a1b1, a1b2), std::max(a2b1, a2b2))};
   }
 
-  Interval operator*(double d) const {
-    if (d > 0) return {lower * d, upper * d};
-    return {upper * d, lower * d};
-  }
-
   Interval operator/(const Interval &other) const {
     if (other.is_const()) return *this / other.lower;
-    constexpr Domain zero = static_cast<Domain>(0);
     constexpr Domain infty = std::numeric_limits<Domain>::infinity();
     Interval reci;
     if (other.lower >= zero || other.upper <= zero) {
-      reci.lower = other.upper == zero ? -infty : (1 / other.upper);
-      reci.upper = other.lower == zero ? infty : (1 / other.lower);
+      reci.lower = other.upper == zero ? -infty : (one / other.upper);
+      reci.upper = other.lower == zero ? infty : (one / other.lower);
     } else {
       reci.lower = -infty;
       reci.upper = infty;
     }
     return *this * reci;
   }
 
-  Interval operator/(double d) const {
-    if (d > 0) return {lower / d, upper / d};
+  Interval operator/(Domain d) const {
+    if (d > zero) return {lower / d, upper / d};
     return {upper / d, lower / d};
   }
 
   constexpr bool is_const() const { return lower == upper; }
 
   Interval operator==(const Interval &other) const {
     if (is_const() && other.is_const() && lower == other.lower)
-      return constant(1);  // must be equal
+      return constant(one);  // must be equal
     if (lower > other.upper || upper < other.lower)
-      return constant(0);  // disjoint, cannot possibly be equal
-    return {0, 1};
+      return constant(zero);  // disjoint, cannot possibly be equal
+    return {zero, one};
   }
 
-  constexpr bool operator==(double d) const { return is_const() && lower == d; }
+  constexpr bool operator==(Domain d) const { return is_const() && lower == d; }
 
   Interval operator>(const Interval &other) const {
-    if (lower > other.upper) return constant(1);
-    if (upper < other.lower) return constant(0);
-    return {0, 1};
+    if (lower > other.upper) return constant(one);
+    if (upper < other.lower) return constant(zero);
+    return {zero, one};
   }
 
   Interval operator<(const Interval &other) const {
-    if (upper < other.lower) return constant(1);
-    if (lower > other.upper) return constant(0);
-    return {0, 1};
+    if (upper < other.lower) return constant(one);
+    if (lower > other.upper) return constant(zero);
+    return {zero, one};
   }
 
   Interval min(const Interval &other) const {
@@ -132,70 +137,76 @@ struct Interval {
   }
 
   Interval abs() const {
-    if (lower >= 0) return *this;
-    if (upper <= 0) return {-upper, -lower};
-    return {0.0, std::max(-lower, upper)};
+    if (lower >= zero) return *this;
+    if (upper <= zero) return {-upper, -lower};
+    return {zero, std::max(-lower, upper)};
   }
 
-  Interval mod(double m) const {
+  Interval mod(Domain m) const {
     // FIXME: cannot deal with negative m right now...
     Domain diff = std::fmod(lower, m);
-    if (diff < 0) diff += m;
+    if (diff < zero) diff += m;
     Domain cycle_min = lower - diff;
     // may be disjoint intervals, but we don't deal with that...
-    if (upper - cycle_min >= m) return {0.0, m};
+    if (upper - cycle_min >= m) return {zero, m};
     return {diff, upper - cycle_min};
   }
 
   Interval logical_and(const Interval &other) const {
-    return {lower == 0.0 || other.lower == 0.0 ? 0.0 : 1.0,
-            upper == 1.0 && other.upper == 1.0 ? 1.0 : 0.0};
+    return {lower == zero || other.lower == zero ? zero : one,
+            upper == one && other.upper == one ? one : zero};
   }
 
   Interval logical_or(const Interval &other) const {
-    return {lower == 0.0 && other.lower == 0.0 ? 0.0 : 1.0,
-            upper == 1.0 || other.upper == 1.0 ? 1.0 : 0.0};
+    return {lower == zero && other.lower == zero ? zero : one,
+            upper == one || other.upper == one ? one : zero};
   }
 
   Interval sin() const {
     if (is_const()) return constant(std::sin(lower));
     // largely similar to cos
-    int64_t min_pis = static_cast<int64_t>(std::floor((lower - kHalfPi) / kPi));
-    int64_t max_pis = static_cast<int64_t>(std::floor((upper - kHalfPi) / kPi));
+    int64_t min_pis = static_cast<int64_t>(std::floor(
+        (lower - static_cast<Domain>(kHalfPi)) / static_cast<Domain>(kPi)));
+    int64_t max_pis = static_cast<int64_t>(std::floor(
+        (upper - static_cast<Domain>(kHalfPi)) / static_cast<Domain>(kPi)));
 
     bool not_cross_pos_1 =
         (min_pis % 2 == 0) ? max_pis - min_pis <= 1 : max_pis == min_pis;
     bool not_cross_neg_1 =
         (min_pis % 2 == 0) ? max_pis == min_pis : max_pis - min_pis <= 1;
 
     Domain new_min =
-        not_cross_neg_1 ? std::min(std::sin(lower), std::sin(upper)) : -1.0;
+        not_cross_neg_1 ? std::min(std::sin(lower), std::sin(upper)) : -one;
     Domain new_max =
-        not_cross_pos_1 ? std::max(std::sin(lower), std::sin(upper)) : 1.0;
+        not_cross_pos_1 ? std::max(std::sin(lower), std::sin(upper)) : one;
     return {new_min, new_max};
   }
 
   Interval cos() const {
     if (is_const()) return constant(std::cos(lower));
-    int64_t min_pis = static_cast<int64_t>(std::floor(lower / kPi));
-    int64_t max_pis = static_cast<int64_t>(std::floor(upper / kPi));
+    int64_t min_pis =
+        static_cast<int64_t>(std::floor(lower / static_cast<Domain>(kPi)));
+    int64_t max_pis =
+        static_cast<int64_t>(std::floor(upper / static_cast<Domain>(kPi)));
 
     bool not_cross_pos_1 =
         (min_pis % 2 == 0) ? max_pis - min_pis <= 1 : max_pis == min_pis;
     bool not_cross_neg_1 =
         (min_pis % 2 == 0) ? max_pis == min_pis : max_pis - min_pis <= 1;
 
     Domain new_min =
-        not_cross_neg_1 ? std::min(std::cos(lower), std::cos(upper)) : -1.0;
+        not_cross_neg_1 ? std::min(std::cos(lower), std::cos(upper)) : -one;
     Domain new_max =
-        not_cross_pos_1 ? std::max(std::cos(lower), std::cos(upper)) : 1.0;
+        not_cross_pos_1 ? std::max(std::cos(lower), std::cos(upper)) : one;
     return {new_min, new_max};
   }
 
   Interval tan() const {
     if (is_const()) return constant(std::tan(lower));
-    int64_t min_pis = static_cast<int64_t>(std::floor((lower + kHalfPi) / kPi));
-    int64_t max_pis = static_cast<int64_t>(std::floor((upper + kHalfPi) / kPi));
+    int64_t min_pis = static_cast<int64_t>(std::floor(
+        (lower + static_cast<Domain>(kHalfPi)) / static_cast<Domain>(kPi)));
+    int64_t max_pis = static_cast<int64_t>(std::floor(
+        (upper + static_cast<Domain>(kHalfPi)) / static_cast<Domain>(kPi)));
     if (min_pis != max_pis)
       return {-std::numeric_limits<Domain>::infinity(),
               std::numeric_limits<Domain>::infinity()};