diff --git a/src/index.d.ts b/src/index.d.ts
index b64e05f..c32e8f0 100644
--- a/src/index.d.ts
+++ b/src/index.d.ts
@@ -4,5 +4,5 @@ export = bezier;
 declare function bezier(x1: number, y1: number, x2: number, y2: number): bezier.EasingFunction;
 
 declare namespace bezier {
-    export interface EasingFunction { (input: number): number }
+    export type EasingFunction = (input: number) => number
 }
diff --git a/src/index.js b/src/index.js
index aa859a4..54e573c 100644
--- a/src/index.js
+++ b/src/index.js
@@ -4,101 +4,65 @@
  * by Gaëtan Renaudeau 2014 - 2015 – MIT License
  */
 
-// These values are established by empiricism with tests (tradeoff: performance VS precision)
-var NEWTON_ITERATIONS = 4;
-var NEWTON_MIN_SLOPE = 0.001;
-var SUBDIVISION_PRECISION = 0.0000001;
-var SUBDIVISION_MAX_ITERATIONS = 10;
-
-var kSplineTableSize = 11;
-var kSampleStepSize = 1.0 / (kSplineTableSize - 1.0);
-
-var float32ArraySupported = typeof Float32Array === 'function';
-
-function A (aA1, aA2) { return 1.0 - 3.0 * aA2 + 3.0 * aA1; }
-function B (aA1, aA2) { return 3.0 * aA2 - 6.0 * aA1; }
-function C (aA1)      { return 3.0 * aA1; }
-
-// Returns x(t) given t, x1, and x2, or y(t) given t, y1, and y2.
-function calcBezier (aT, aA1, aA2) { return ((A(aA1, aA2) * aT + B(aA1, aA2)) * aT + C(aA1)) * aT; }
-
-// Returns dx/dt given t, x1, and x2, or dy/dt given t, y1, and y2.
-function getSlope (aT, aA1, aA2) { return 3.0 * A(aA1, aA2) * aT * aT + 2.0 * B(aA1, aA2) * aT + C(aA1); }
-
-function binarySubdivide (aX, aA, aB, mX1, mX2) {
-  var currentX, currentT, i = 0;
-  do {
-    currentT = aA + (aB - aA) / 2.0;
-    currentX = calcBezier(currentT, mX1, mX2) - aX;
-    if (currentX > 0.0) {
-      aB = currentT;
-    } else {
-      aA = currentT;
-    }
-  } while (Math.abs(currentX) > SUBDIVISION_PRECISION && ++i < SUBDIVISION_MAX_ITERATIONS);
-  return currentT;
-}
-
-function newtonRaphsonIterate (aX, aGuessT, mX1, mX2) {
- for (var i = 0; i < NEWTON_ITERATIONS; ++i) {
-   var currentSlope = getSlope(aGuessT, mX1, mX2);
-   if (currentSlope === 0.0) {
-     return aGuessT;
-   }
-   var currentX = calcBezier(aGuessT, mX1, mX2) - aX;
-   aGuessT -= currentX / currentSlope;
- }
- return aGuessT;
-}
-
-function LinearEasing (x) {
-  return x;
+function LinearEasing(x) {
+	return x;
 }
 
-module.exports = function bezier (mX1, mY1, mX2, mY2) {
-  if (!(0 <= mX1 && mX1 <= 1 && 0 <= mX2 && mX2 <= 1)) {
-    throw new Error('bezier x values must be in [0, 1] range');
-  }
-
-  if (mX1 === mY1 && mX2 === mY2) {
-    return LinearEasing;
-  }
-
-  // Precompute samples table
-  var sampleValues = float32ArraySupported ? new Float32Array(kSplineTableSize) : new Array(kSplineTableSize);
-  for (var i = 0; i < kSplineTableSize; ++i) {
-    sampleValues[i] = calcBezier(i * kSampleStepSize, mX1, mX2);
-  }
-
-  function getTForX (aX) {
-    var intervalStart = 0.0;
-    var currentSample = 1;
-    var lastSample = kSplineTableSize - 1;
-
-    for (; currentSample !== lastSample && sampleValues[currentSample] <= aX; ++currentSample) {
-      intervalStart += kSampleStepSize;
-    }
-    --currentSample;
-
-    // Interpolate to provide an initial guess for t
-    var dist = (aX - sampleValues[currentSample]) / (sampleValues[currentSample + 1] - sampleValues[currentSample]);
-    var guessForT = intervalStart + dist * kSampleStepSize;
-
-    var initialSlope = getSlope(guessForT, mX1, mX2);
-    if (initialSlope >= NEWTON_MIN_SLOPE) {
-      return newtonRaphsonIterate(aX, guessForT, mX1, mX2);
-    } else if (initialSlope === 0.0) {
-      return guessForT;
-    } else {
-      return binarySubdivide(aX, intervalStart, intervalStart + kSampleStepSize, mX1, mX2);
-    }
-  }
-
-  return function BezierEasing (x) {
-    // Because JavaScript number are imprecise, we should guarantee the extremes are right.
-    if (x === 0 || x === 1) {
-      return x;
-    }
-    return calcBezier(getTForX(x), mY1, mY2);
-  };
+const { cbrt, sqrt, PI: π } = Math;
+
+const x2t = (x, a, b, c, d) => {
+	const q = a + b * x;
+	const s = q ** 2 + c;
+	if (s > 0) {
+		const root = sqrt(s);
+		return cbrt(q + root) + cbrt(q - root) - d;
+	}
+	const l = cbrt(sqrt(q * q - s));
+	const angle = q ? Math.atan(sqrt(-s) / q) : -π / 2;
+	let φ;
+	if (b < 0) {
+		φ = (q > 0 ? 2 * π : π) - angle;
+	} else if (d < 0) {
+		φ = (q > 0 ? 2 * π : -3 * π) + angle;
+	} else {
+		φ = (q > 0 ? 0 : π) + angle;
+	}
+	return 2 * l * Math.cos(φ / 3) - d;
+};
+const Y = (t, ay, by, cy) => ((ay * t + 3 * by) * t + cy) * t;
+
+module.exports = function bezier(mX1, mY1, mX2, mY2) {
+	if (!(0 <= mX1 && mX1 <= 1 && 0 <= mX2 && mX2 <= 1)) {
+		throw new Error("Bezier x values must be in [0, 1] range");
+	}
+
+	if (mX1 === mY1 && mX2 === mY2) {
+		return LinearEasing;
+	}
+
+	const a = 6 * (3 * mX1 - 3 * mX2 + 1);
+	const b = 6 * (mX2 - 2 * mX1);
+	const c = 3 * mX1;
+
+	const a2 = a * a;
+	const b2 = b * b;
+
+	const d = b / a;
+	const e = (3 * b * c) / a2 - (b2 * b) / (a2 * a);
+	const w1 = (2 * c) / a - b2 / a2;
+	const w = w1 * w1 * w1;
+	const o = 3 / a;
+
+	const ay = 3 * mY1 - 3 * mY2 + 1;
+	const by = mY2 - 2 * mY1;
+	const cy = 3 * mY1;
+
+	const X2T = a ? x2t : LinearEasing;
+
+	return (x) => {
+		if (x === 0 || x === 1) {
+			return x;
+		}
+		return Y(X2T(x, e, o, w, d), ay, by, cy);
+	};
 };