Add C++ implementations for xLAPY2 and xLAPY3. They are straight translations of the corresponding Fortran code, even though xLAPY2 can just call std:hypot.

The code currently lives in eigen3/lapack/experimental. The plan is to finish development here, and once its done upstream it replacing all the f2c'ed code.

PiperOrigin-RevId: 380897335

lapack/experimental/README.md

diff --git a/lapack/experimental/README.md b/lapack/experimental/README.md
new file mode 100644
index 0000000..86aaace
--- /dev/null
+++ b/lapack/experimental/README.md
@@ -0,0 +1,3 @@
+This directory contains C++11 replacment for the f2c-ed functions in the lapack
+directory. This is a temporary directory, which is being used to develop and
+stage these changes, once this code is upstreamed, it will be deleted.

lapack/experimental/xlapy2.cpp

diff --git a/lapack/experimental/xlapy2.cpp b/lapack/experimental/xlapy2.cpp
new file mode 100644
index 0000000..5992f9e
--- /dev/null
+++ b/lapack/experimental/xlapy2.cpp
@@ -0,0 +1,24 @@
+#include <algorithm>
+#include <cmath>
+
+extern "C" {
+float slapy2_(float* x, float* y);
+double dlapy2_(double* x, double* y);
+}
+
+// Compute sqrt(x*x + y*y) taking care not to cause unnecessary overflow.
+template <typename T>
+T xlapy2(const T x, const T y) {
+  const T xabs = std::abs(x);
+  const T yabs = std::abs(y);
+  const T w = std::fmax(xabs, yabs);
+  const T z = std::fmin(xabs, yabs);
+  if (z > T(0)) {
+    const T t = z / w;
+    return w * std::sqrt(T(1) + t * t);
+  }
+  return w;
+}
+
+float slapy2_(float* x, float* y) { return xlapy2(*x, *y); }
+double dlapy2_(double* x, double* y) { return xlapy2(*x, *y); }

lapack/experimental/xlapy3.cpp

diff --git a/lapack/experimental/xlapy3.cpp b/lapack/experimental/xlapy3.cpp
new file mode 100644
index 0000000..9c935fa
--- /dev/null
+++ b/lapack/experimental/xlapy3.cpp
@@ -0,0 +1,30 @@
+#include <algorithm>
+#include <cmath>
+
+extern "C" {
+float slapy3_(float* x, float* y, float* z);
+double dlapy3_(double* x, double* y, double* z);
+}
+
+// Compute sqrt(x*x + y*y + z*z) taking care not to cause unnecessary overflow.
+template <typename T>
+T xlapy3(const T x, const T y, const T z) {
+  const T xabs = std::abs(x);
+  const T yabs = std::abs(y);
+  const T zabs = std::abs(z);
+  const T w = std::fmax(std::fmax(xabs, yabs), zabs);
+
+  if (w == T(0)) {
+    // W can be zero for max(0,nan,0) adding all three entries together will
+    // make sure NaN will not disappear.
+    return (xabs + yabs + zabs);
+  }
+
+  const T tx = xabs / w;
+  const T ty = yabs / w;
+  const T tz = zabs / w;
+  return w * std::sqrt(tx * tx + ty * ty + tz * tz);
+}
+
+float slapy3_(float* x, float* y, float* z) { return xlapy3(*x, *y, *z); }
+double dlapy3_(double* x, double* y, double* z) { return xlapy3(*x, *y, *z); }