lapack/experimental/xlapy2.cpp

#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); }