xLADIV implementation using the complex division algorithm described in http://www.open-std.org/jtc1/sc22/wg11/docs/n435.pdf. This is the same algorithm used by LLVM's libc++.
PiperOrigin-RevId: 381083249
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);
+}
