diff --git a/src/Makefile.sources b/src/Makefile.sources
index 24ce798..70b80be 100644
--- a/src/Makefile.sources
+++ b/src/Makefile.sources
@@ -353,6 +353,7 @@
 	hb-subset-cff2.hh \
 	hb-subset-input.cc \
 	hb-subset-input.hh \
+	hb-subset-instancer-solver.cc \
 	hb-subset-accelerator.hh \
 	hb-subset-plan.cc \
 	hb-subset-plan.hh \
diff --git a/src/harfbuzz-subset.cc b/src/harfbuzz-subset.cc
index d485233..196e32f 100644
--- a/src/harfbuzz-subset.cc
+++ b/src/harfbuzz-subset.cc
@@ -53,6 +53,7 @@
 #include "hb-subset-cff1.cc"
 #include "hb-subset-cff2.cc"
 #include "hb-subset-input.cc"
+#include "hb-subset-instancer-solver.cc"
 #include "hb-subset-plan.cc"
 #include "hb-subset-repacker.cc"
 #include "hb-subset.cc"
diff --git a/src/hb-subset-instancer-solver.cc b/src/hb-subset-instancer-solver.cc
new file mode 100644
index 0000000..89360d0
--- /dev/null
+++ b/src/hb-subset-instancer-solver.cc
@@ -0,0 +1,405 @@
+/*
+ * Copyright © 2023  Behdad Esfahbod
+ *
+ *  This is part of HarfBuzz, a text shaping library.
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and its documentation for any purpose, provided that the
+ * above copyright notice and the following two paragraphs appear in
+ * all copies of this software.
+ *
+ * IN NO EVENT SHALL THE COPYRIGHT HOLDER BE LIABLE TO ANY PARTY FOR
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
+ * ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN
+ * IF THE COPYRIGHT HOLDER HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
+ * DAMAGE.
+ *
+ * THE COPYRIGHT HOLDER SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING,
+ * BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+ * FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS
+ * ON AN "AS IS" BASIS, AND THE COPYRIGHT HOLDER HAS NO OBLIGATION TO
+ * PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+ */
+
+#include "hb.hh"
+
+/* This file is a straight port of the following:
+ *
+ * https://github.com/fonttools/fonttools/blob/f73220816264fc383b8a75f2146e8d69e455d398/Lib/fontTools/varLib/instancer/solver.py
+ *
+ * See XXX markers for unfinished parts.
+ *
+ * Where that file returns None for a triple, we return Triple{}.
+ * This should be safe.
+ */
+
+constexpr static float EPSILON = 1.f / (1 << 14);
+constexpr static float MAX_F2DOT14 = float (0x7FFF) / (1 << 14);
+
+struct Triple {
+
+  Triple () :
+    minimum (0.f), middle (0.f), maximum (0.f) {}
+
+  Triple (float minimum_, float middle_, float maximum_) :
+    minimum (minimum_), middle (middle_), maximum (maximum_) {}
+
+  float minimum;
+  float middle;
+  float maximum;
+};
+
+static inline Triple _reverse_negate(const Triple &v)
+{ return {-v.maximum, -v.middle, -v.minimum}; }
+
+using result_item_t = hb_pair_t<float, Triple>;
+using result_t = hb_vector_t<result_item_t>;
+
+static inline result_t
+_solve (Triple tent, Triple axisLimit, bool negative = false)
+{
+  float axisMin = axisLimit.minimum;
+  float axisDef = axisLimit.middle;
+  float axisMax = axisLimit.maximum;
+  float lower = tent.minimum;
+  float peak  = tent.middle;
+  float upper = tent.maximum;
+
+  // Mirror the problem such that axisDef <= peak
+  if (axisDef > peak)
+  {
+    result_t vec = _solve (_reverse_negate (tent),
+			   _reverse_negate (axisLimit),
+			   !negative);
+
+    for (auto &p : vec)
+      p = hb_pair (p.first, _reverse_negate (p.second));
+
+    return vec;
+  }
+  // axisDef <= peak
+
+  /* case 1: The whole deltaset falls outside the new limit; we can drop it
+   *
+   *                                          peak
+   *  1.........................................o..........
+   *                                           / \
+   *                                          /   \
+   *                                         /     \
+   *                                        /       \
+   *  0---|-----------|----------|-------- o         o----1
+   *    axisMin     axisDef    axisMax   lower     upper
+   */
+  if (axisMax <= lower && axisMax < peak)
+      return result_t{};  // No overlap
+
+  /* case 2: Only the peak and outermost bound fall outside the new limit;
+   * we keep the deltaset, update peak and outermost bound and and scale deltas
+   * by the scalar value for the restricted axis at the new limit, and solve
+   * recursively.
+   *
+   *                                  |peak
+   *  1...............................|.o..........
+   *                                  |/ \
+   *                                  /   \
+   *                                 /|    \
+   *                                / |     \
+   *  0--------------------------- o  |      o----1
+   *                           lower  |      upper
+   *                                  |
+   *                                axisMax
+   *
+   * Convert to:
+   *
+   *  1............................................
+   *                                  |
+   *                                  o peak
+   *                                 /|
+   *                                /x|
+   *  0--------------------------- o  o upper ----1
+   *                           lower  |
+   *                                  |
+   *                                axisMax
+   */
+  if (axisMax < peak)
+  {
+    float mult = 1.f; //XXX supportScalar({"tag": axisMax}, {"tag": tent})
+    tent = Triple{lower, axisMax, axisMax};
+
+    result_t vec = _solve (tent, axisLimit);
+
+    for (auto &p : vec)
+      p = hb_pair (p.first * mult, p.second);
+
+    return vec;
+  }
+
+  // lower <= axisDef <= peak <= axisMax
+
+  float gain = 1.0f; // XXX supportScalar({"tag": axisDef}, {"tag": tent})
+  result_t out {hb_pair (gain, Triple{})};
+
+  // First, the positive side
+
+  // outGain is the scalar of axisMax at the tent.
+  float outGain = 1.f; // XXX supportScalar({"tag": axisMax}, {"tag": tent})
+
+  /* Case 3a: Gain is more than outGain. The tent down-slope crosses
+   * the axis into negative. We have to split it into multiples.
+   *
+   *                      | peak  |
+   *  1...................|.o.....|..............
+   *                      |/x\_   |
+   *  gain................+....+_.|..............
+   *                     /|    |y\|
+   *  ................../.|....|..+_......outGain
+   *                   /  |    |  | \
+   *  0---|-----------o   |    |  |  o----------1
+   *    axisMin    lower  |    |  |   upper
+   *                      |    |  |
+   *                axisDef    |  axisMax
+   *                           |
+   *                      crossing
+   */
+  if (gain > outGain)
+  {
+    // Crossing point on the axis.
+    float crossing = peak + ((1 - gain) * (upper - peak) / (1 - outGain));
+
+    Triple loc{peak, peak, crossing};
+    float scalar = 1.f;
+
+    // The part before the crossing point.
+    out.push (hb_pair (scalar - gain, loc));
+
+    /* The part after the crossing point may use one or two tents,
+     * depending on whether upper is before axisMax or not, in one
+     * case we need to keep it down to eternity.
+     *
+     * Case 3a1, similar to case 1neg; just one tent needed, as in
+     * the drawing above.
+     */
+    if (upper >= axisMax)
+    {
+      Triple loc {crossing, axisMax, axisMax};
+      float scalar = 1.f; // XXX supportScalar({"tag": axisMax}, {"tag": tent})
+
+      out.push (hb_pair (scalar - gain, loc));
+    }
+
+    /* Case 3a2: Similar to case 2neg; two tents needed, to keep
+     * down to eternity.
+     *
+     *                      | peak             |
+     *  1...................|.o................|...
+     *                      |/ \_              |
+     *  gain................+....+_............|...
+     *                     /|    | \xxxxxxxxxxy|
+     *                    / |    |  \_xxxxxyyyy|
+     *                   /  |    |    \xxyyyyyy|
+     *  0---|-----------o   |    |     o-------|--1
+     *    axisMin    lower  |    |      upper  |
+     *                      |    |             |
+     *                axisDef    |             axisMax
+     *                           |
+     *                      crossing
+     */
+    else
+    {
+      // A tent's peak cannot fall on axis default. Nudge it.
+      if (upper == axisDef)
+	upper += EPSILON;
+
+      // Downslope.
+      Triple loc1 {crossing, upper, axisMax};
+      float scalar1 = 0.f;
+
+      // Eternity justify.
+      Triple loc2 {upper, axisMax, axisMax};
+      float scalar2 = 1.f; // supportScalar({"tag": axisMax}, {"tag": tent})
+
+      out.push (hb_pair (scalar1 - gain, loc1));
+      out.push (hb_pair (scalar2 - gain, loc2));
+    }
+  }
+
+  /* Case 3: Outermost limit still fits within F2Dot14 bounds;
+   * we keep deltas as is and only scale the axes bounds. Deltas beyond -1.0
+   * or +1.0 will never be applied as implementations must clamp to that range.
+   *
+   * A second tent is needed for cases when gain is positive, though we add it
+   * unconditionally and it will be dropped because scalar ends up 0.
+   *
+   * TODO: See if we can just move upper closer to adjust the slope, instead of
+   * second tent.
+   *
+   *            |           peak |
+   *  1.........|............o...|..................
+   *            |           /x\  |
+   *            |          /xxx\ |
+   *            |         /xxxxx\|
+   *            |        /xxxxxxx+
+   *            |       /xxxxxxxx|\
+   *  0---|-----|------oxxxxxxxxx|xo---------------1
+   *    axisMin |  lower         | upper
+   *            |                |
+   *          axisDef          axisMax
+   */
+  else if (axisDef + (axisMax - axisDef) * 2 >= upper)
+  {
+    if (!negative && axisDef + (axisMax - axisDef) * MAX_F2DOT14 < upper)
+    {
+      // we clamp +2.0 to the max F2Dot14 (~1.99994) for convenience
+      upper = axisDef + (axisMax - axisDef) * MAX_F2DOT14;
+      assert (peak < upper);
+    }
+
+    // Special-case if peak is at axisMax.
+    if (axisMax == peak)
+	upper = peak;
+
+    Triple loc1 {hb_max (axisDef, lower), peak, upper};
+    float scalar1 = 1.f;
+
+    Triple loc2 {peak, upper, upper};
+    float scalar2 = 0.f;
+
+    // Don't add a dirac delta!
+    if (axisDef < upper)
+	out.push (hb_pair (scalar1 - gain, loc1));
+    if (peak < upper)
+	out.push (hb_pair (scalar2 - gain, loc2));
+  }
+
+  /* Case 4: New limit doesn't fit; we need to chop into two tents,
+   * because the shape of a triangle with part of one side cut off
+   * cannot be represented as a triangle itself.
+   *
+   *            |   peak |
+   *  1.........|......o.|...................
+   *            |     /x\|
+   *            |    |xxy|\_
+   *            |   /xxxy|  \_
+   *            |  |xxxxy|    \_
+   *            |  /xxxxy|      \_
+   *  0---|-----|-oxxxxxx|        o----------1
+   *    axisMin | lower  |        upper
+   *            |        |
+   *          axisDef  axisMax
+   */
+  else
+  {
+    Triple loc1 {hb_max (axisDef, lower), peak, axisMax};
+    float scalar1 = 1.f;
+
+    Triple loc2 {peak, axisMax, axisMax};
+    float scalar2 = 1.f; // XXX supportScalar({"tag": axisMax}, {"tag": tent})
+
+    out.push (hb_pair (scalar1 - gain, loc1));
+    // Don't add a dirac delta!
+    if (peak < axisMax)
+      out.push (hb_pair (scalar2 - gain, loc2));
+  }
+
+  /* Now, the negative side
+   *
+   * Case 1neg: Lower extends beyond axisMin: we chop. Simple.
+   *
+   *                     |   |peak
+   *  1..................|...|.o.................
+   *                     |   |/ \
+   *  gain...............|...+...\...............
+   *                     |x_/|    \
+   *                     |/  |     \
+   *                   _/|   |      \
+   *  0---------------o  |   |       o----------1
+   *              lower  |   |       upper
+   *                     |   |
+   *               axisMin   axisDef
+   */
+  if (lower <= axisMin)
+  {
+    Triple loc {axisMin, axisMin, axisDef};
+    float scalar = 1.f; // XXX supportScalar({"tag": axisMin}, {"tag": tent})
+
+    out.push (hb_pair (scalar - gain, loc));
+  }
+
+  /* Case 2neg: Lower is betwen axisMin and axisDef: we add two
+   * tents to keep it down all the way to eternity.
+   *
+   *      |               |peak
+   *  1...|...............|.o.................
+   *      |               |/ \
+   *  gain|...............+...\...............
+   *      |yxxxxxxxxxxxxx/|    \
+   *      |yyyyyyxxxxxxx/ |     \
+   *      |yyyyyyyyyyyx/  |      \
+   *  0---|-----------o   |       o----------1
+   *    axisMin    lower  |       upper
+   *                      |
+   *                    axisDef
+   */
+  else
+  {
+    // A tent's peak cannot fall on axis default. Nudge it.
+    if (lower == axisDef)
+      lower -= EPSILON;
+
+    // Downslope.
+    Triple loc1 {axisMin, lower, axisDef};
+    float scalar1 = 0.f;
+
+    // Eternity justify.
+    Triple loc2 {axisMin, axisMin, lower};
+    float scalar2 = 0.f;
+
+    out.push (hb_pair (scalar1 - gain, loc1));
+    out.push (hb_pair (scalar2 - gain, loc2));
+  }
+
+  return out;
+}
+
+/* Given a tuple (lower,peak,upper) "tent" and new axis limits
+ * (axisMin,axisDefault,axisMax), solves how to represent the tent
+ * under the new axis configuration.  All values are in normalized
+ * -1,0,+1 coordinate system. Tent values can be outside this range.
+ *
+ * Return value: a list of tuples. Each tuple is of the form
+ * (scalar,tent), where scalar is a multipler to multiply any
+ * delta-sets by, and tent is a new tent for that output delta-set.
+ * If tent value is Triple{}, that is a special deltaset that should
+ * be always-enabled (called "gain").
+ */
+HB_INTERNAL result_t rebase_tent (Triple tent, Triple axisLimit);
+
+result_t
+rebase_tent (Triple tent, Triple axisLimit)
+{
+  float axisMin = axisLimit.minimum;
+  float axisDef = axisLimit.middle;
+  float axisMax = axisLimit.maximum;
+  float lower = tent.minimum;
+  float peak  = tent.middle;
+  float upper = tent.maximum;
+
+  assert (-1.f <= axisMin && axisMin <= axisDef && axisDef <= axisMax && axisMax <= +1.f);
+  assert (-2.f <= lower && lower <= peak && peak <= upper && upper <= +2.f);
+  assert (peak != 0.f);
+
+  result_t sols = _solve (tent, axisLimit);
+
+#if 0
+  // XXX
+  n = lambda v: normalizeValue(v, axisLimit, extrapolate=True)
+  sols = [
+      (scalar, (n(v[0]), n(v[1]), n(v[2])) if v is not None else None)
+      for scalar, v in sols
+      if scalar
+  ]
+#endif
+
+  return sols;
+}
diff --git a/src/meson.build b/src/meson.build
index b678b5b..27d9a31 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -353,6 +353,7 @@
   'hb-subset-cff2.hh',
   'hb-subset-input.cc',
   'hb-subset-input.hh',
+  'hb-subset-instancer-solver.cc',
   'hb-subset-plan.cc',
   'hb-subset-plan.hh',
   'hb-subset-repacker.cc',