/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */
/*
 * hypot() function for DJGPP.
 *
 * hypot() computes sqrt(x^2 + y^2).  The problem with the obvious
 * naive implementation is that it might fail for very large or
 * very small arguments.  For instance, for large x or y the result
 * might overflow even if the value of the function should not,
 * because squaring a large number might trigger an overflow.  For
 * very small numbers, their square might underflow and will be
 * silently replaced by zero; this won't cause an exception, but might
 * have an adverse effect on the accuracy of the result.
 *
 * This implementation tries to avoid the above pitfals, without
 * inflicting too much of a performance hit.
 *
 */

/// #include <float.h>
#include <errno.h>
#include <math.h>

/* Approximate square roots of DBL_MAX and DBL_MIN.  Numbers
   between these two shouldn't neither overflow nor underflow
   when squared.  */
#define __SQRT_DBL_MAX 1.3e+154
#define __SQRT_DBL_MIN 2.3e-162

double hypot(double x, double y)
{
    double abig = fabs(x), asmall = fabs(y);
    double ratio;

    /* Make abig = max(|x|, |y|), asmall = min(|x|, |y|).  */
    if (abig < asmall) {
        double temp = abig;

        abig = asmall;
        asmall = temp;
    }

    /* Trivial case.  */
    if (asmall == 0.)
        return abig;

    /* Scale the numbers as much as possible by using its ratio.
       For example, if both ABIG and ASMALL are VERY small, then
       X^2 + Y^2 might be VERY inaccurate due to loss of
       significant digits.  Dividing ASMALL by ABIG scales them
       to a certain degree, so that accuracy is better.  */

    if ((ratio = asmall / abig) > __SQRT_DBL_MIN && abig < __SQRT_DBL_MAX)
        return abig * sqrt(1.0 + ratio * ratio);
    else {
        /* Slower but safer algorithm due to Moler and Morrison.  Never
           produces any intermediate result greater than roughly the
           larger of X and Y.  Should converge to machine-precision
           accuracy in 3 iterations.  */

        double r = ratio * ratio, t, s, p = abig, q = asmall;

        do {
            t = 4. + r;
            if (t == 4.)
                break;
            s = r / t;
            p += 2. * s * p;
            q *= s;
            r = (q / p) * (q / p);
        } while (1);

        return p;
    }
}

#ifdef TEST

#include <stdio.h>

int main(void)
{
    printf("hypot(3, 4) =\t\t\t %25.17e\n", hypot(3., 4.));
    printf("hypot(3*10^150, 4*10^150) =\t %25.17g\n", hypot(3.e+150, 4.e+150));
    printf("hypot(3*10^306, 4*10^306) =\t %25.17g\n", hypot(3.e+306, 4.e+306));
    printf("hypot(3*10^-320, 4*10^-320) =\t %25.17g\n",
        hypot(3.e-320, 4.e-320));
    printf("hypot(0.7*DBL_MAX, 0.7*DBL_MAX) =%25.17g\n",
        hypot(0.7 * DBL_MAX, 0.7 * DBL_MAX));
    printf("hypot(DBL_MAX, 1.0) =\t\t %25.17g\n", hypot(DBL_MAX, 1.0));
    printf("hypot(1.0, DBL_MAX) =\t\t %25.17g\n", hypot(1.0, DBL_MAX));
    printf("hypot(0.0, DBL_MAX) =\t\t %25.17g\n", hypot(0.0, DBL_MAX));

    return 0;
}

#endif