diff --git a/lapack/experimental/xladiv.cpp b/lapack/experimental/xladiv.cpp
new file mode 100644
index 0000000..381950f
--- /dev/null
+++ b/lapack/experimental/xladiv.cpp
@@ -0,0 +1,64 @@
+#include <cmath>
+#include <complex>
+
+extern "C" {
+void sladiv_(float* a, float* b, float* c, float* d, float* p, float* q);
+void dladiv_(double* a, double* b, double* c, double* d, double* p, double* q);
+std::complex<float> cladiv_(std::complex<float>* x, std::complex<float>* y);
+std::complex<double> zladiv_(std::complex<double>* x, std::complex<double>* y);
+}
+
+// This implementation is from pg. 21
+// http://www.open-std.org/jtc1/sc22/wg11/docs/n435.pdf and LLVM's libc++.
+template <typename T>
+void complex_divide(T a, T b, T c, T d, T& x, T& y) {
+  int ilogbw = 0;
+  const T logbw = std::logb(std::fmax(std::abs(c), std::abs(d)));
+  if (std::isfinite(logbw)) {
+    ilogbw = static_cast<int>(logbw);
+    c = scalbn(c, -ilogbw);
+    d = scalbn(d, -ilogbw);
+  }
+  const T denom = c * c + d * d;
+
+  x = scalbn((a * c + b * d) / denom, -ilogbw);
+  y = scalbn((b * c - a * d) / denom, -ilogbw);
+  if (std::isnan(x) && std::isnan(y)) {
+    if ((denom == T(0)) && (!std::isnan(a) || !std::isnan(b))) {
+      x = std::copysign(std::numeric_limits<T>::infinity(), c) * a;
+      y = std::copysign(std::numeric_limits<T>::infinity(), c) * b;
+    } else if ((std::isinf(a) || std::isinf(b)) && std::isfinite(c) &&
+               std::isfinite(d)) {
+      a = std::copysign(std::isinf(a) ? T(1) : T(0), a);
+      b = std::copysign(std::isinf(b) ? T(1) : T(0), b);
+      x = std::numeric_limits<T>::infinity() * (a * c + b * d);
+      y = std::numeric_limits<T>::infinity() * (b * c - a * d);
+    } else if (std::isinf(logbw) && logbw > T(0) && std::isfinite(a) &&
+               std::isfinite(b)) {
+      c = std::copysign(std::isinf(c) ? T(1) : T(0), c);
+      d = std::copysign(std::isinf(d) ? T(1) : T(0), d);
+      x = T(0) * (a * c + b * d);
+      y = T(0) * (b * c - a * d);
+    }
+  }
+}
+
+void sladiv_(float* a, float* b, float* c, float* d, float* p, float* q) {
+  complex_divide(*a, *b, *c, *d, *p, *q);
+}
+
+void dladiv_(double* a, double* b, double* c, double* d, double* p, double* q) {
+  complex_divide(*a, *b, *c, *d, *p, *q);
+}
+
+std::complex<float> cladiv_(std::complex<float>* x, std::complex<float>* y) {
+  float p, q;
+  complex_divide(x->real(), x->imag(), y->real(), y->imag(), p, q);
+  return std::complex<float>(p, q);
+}
+
+std::complex<double> zladiv_(std::complex<double>* x, std::complex<double>* y) {
+  double p, q;
+  complex_divide(x->real(), x->imag(), y->real(), y->imag(), p, q);
+  return std::complex<double>(p, q);
+}