Boost C++ Libraries

boost/math/special_functions/hypot.hpp

//  (C) Copyright John Maddock 2005.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_HYPOT_INCLUDED
#define BOOST_MATH_HYPOT_INCLUDED

#include <cmath>
#include <boost/limits.hpp>
#include <algorithm> // swap

#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
#  include <boost/static_assert.hpp>
#else
#  include <boost/assert.hpp>
#endif

#ifdef BOOST_NO_STDC_NAMESPACE
namespace std{ using ::sqrt; using ::fabs; }
#endif


namespace boost{ namespace math{

template <class T>
T hypot(T x, T y)
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
   BOOST_STATIC_ASSERT(::std::numeric_limits<T>::is_specialized);
#else
   BOOST_ASSERT(std::numeric_limits<T>::is_specialized);
#endif

   //
   // normalize x and y, so that both are positive and x >= y:
   //
   x = (std::fabs)(x);
   y = (std::fabs)(y);

   // special case, see C99 Annex F:
   if(std::numeric_limits<T>::has_infinity
      && ((x == std::numeric_limits<T>::infinity())
      || (y == std::numeric_limits<T>::infinity())))
      return std::numeric_limits<T>::infinity();

   if(y > x) 
      (std::swap)(x, y);
   //
   // figure out overflow and underflow limits:
   //
   T safe_upper = (std::sqrt)((std::numeric_limits<T>::max)()) / 2;
   T safe_lower = (std::sqrt)((std::numeric_limits<T>::min)());
   static const T one = 1;
   //
   // Now handle special cases:
   //
   if(x >= safe_upper)
   {
      if(y <= one)
      {
         // y is neligible:
         return x;
      }
      return (std::sqrt)(x) * (std::sqrt)(y) * (std::sqrt)(x/y + y/x);
   }
   else if(y <= safe_lower)
   {
      if((x >= one) || (y == 0))
      {
         // y is negligible:
         return x;
      }
      return (std::sqrt)(x) * (std::sqrt)(y) * (std::sqrt)(x/y + y/x);
   }
   //
   // If we get here then x^2+y^2 will not overflow or underflow:
   //
   return (std::sqrt)(x*x + y*y);
}

} } // namespaces

#endif // BOOST_MATH_HYPOT_INCLUDED