//  (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