ABACUS/src/MATRIX/pythag.cc

#include "ABACUS.h"
using namespace std;

DP ABACUS::pythag(DP a, DP b)
{
  DP absa, absb;

  absa = fabs(a);
  absb = fabs(b);
  if (absa > absb) return absa * sqrt(1.0 + (absb * absb/(absa * absa)));
  else return (absb == 0.0 ? 0.0 : absb * sqrt(1.0 + (absa * absa/(absb * absb))));
}