lapack/experimental/xladiv.cpp

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