#include "cminpack.h"
#include <math.h>
#include "cminpackP.h"

/*
  About the values for rdwarf and rgiant.

  The original values, both in signe-precision FORTRAN source code and in double-precision code were:
#define rdwarf 3.834e-20
#define rgiant 1.304e19
  See for example:
    http://www.netlib.org/slatec/src/denorm.f
    http://www.netlib.org/slatec/src/enorm.f
  However, rdwarf is smaller than sqrt(FLT_MIN) = 1.0842021724855044e-19, so that rdwarf**2 will
  underflow. This contradicts the constraints expressed in the comments below.

  We changed these constants to be sqrt(dpmpar(2))*0.9 and sqrt(dpmpar(3))*0.9, as proposed by the
  implementation found in MPFIT http://cow.physics.wisc.edu/~craigm/idl/fitting.html
*/

#define double_dwarf (1.4916681462400413e-154*0.9)
#define double_giant (1.3407807929942596e+154*0.9)
#define float_dwarf (1.0842021724855044e-19f*0.9f)
#define float_giant (1.8446743523953730e+19f*0.9f)
#define half_dwarf (2.4414062505039999e-4f*0.9f)
#define half_giant (255.93749236874225497222f*0.9f)

#define dwarf(type) _dwarf(type)
#define _dwarf(type) type ## _dwarf
#define giant(type) _giant(type)
#define _giant(type) type ## _giant

#define rdwarf dwarf(real)
#define rgiant giant(real)

__cminpack_attr__
real __cminpack_func__(enorm)(int n, const real *x)
{
#ifdef USE_CBLAS
    return cblas_dnrm2(n, x, 1);
#else /* !USE_CBLAS */
    /* System generated locals */
    real ret_val, d1;

    /* Local variables */
    int i;
    real s1, s2, s3, xabs, x1max, x3max, agiant, floatn;

/*     ********** */

/*     function enorm */

/*     given an n-vector x, this function calculates the */
/*     euclidean norm of x. */

/*     the euclidean norm is computed by accumulating the sum of */
/*     squares in three different sums. the sums of squares for the */
/*     small and large components are scaled so that no overflows */
/*     occur. non-destructive underflows are permitted. underflows */
/*     and overflows do not occur in the computation of the unscaled */
/*     sum of squares for the intermediate components. */
/*     the definitions of small, intermediate and large components */
/*     depend on two constants, rdwarf and rgiant. the main */
/*     restrictions on these constants are that rdwarf**2 not */
/*     underflow and rgiant**2 not overflow. the constants */
/*     given here are suitable for every known computer. */

/*     the function statement is */

/*       double precision function enorm(n,x) */

/*     where */

/*       n is a positive integer input variable. */

/*       x is an input array of length n. */

/*     subprograms called */

/*       fortran-supplied ... dabs,dsqrt */

/*     argonne national laboratory. minpack project. march 1980. */
/*     burton s. garbow, kenneth e. hillstrom, jorge j. more */

/*     ********** */

    s1 = 0.;
    s2 = 0.;
    s3 = 0.;
    x1max = 0.;
    x3max = 0.;
    floatn = (real) (n);
    agiant = rgiant / floatn;
    for (i = 0; i < n; ++i) {
	xabs = fabs(x[i]);
	if (xabs <= rdwarf || xabs >= agiant) {
            if (xabs > rdwarf) {

/*              sum for large components. */

                if (xabs > x1max) {
                    /* Computing 2nd power */
                    d1 = x1max / xabs;
                    s1 = 1. + s1 * (d1 * d1);
                    x1max = xabs;
                } else {
                    /* Computing 2nd power */
                    d1 = xabs / x1max;
                    s1 += d1 * d1;
                }
            } else {

/*              sum for small components. */

                if (xabs > x3max) {
                    /* Computing 2nd power */
                    d1 = x3max / xabs;
                    s3 = 1. + s3 * (d1 * d1);
                    x3max = xabs;
                } else {
                    if (xabs != 0.) {
                        /* Computing 2nd power */
                        d1 = xabs / x3max;
                        s3 += d1 * d1;
                    }
                }
            }
	} else {

/*           sum for intermediate components. */

            /* Computing 2nd power */
            s2 += xabs * xabs;
        }
    }

/*     calculation of norm. */

    if (s1 != 0.) {
        ret_val = x1max * sqrt(s1 + (s2 / x1max) / x1max);
    } else {
        if (s2 != 0.) {
            if (s2 >= x3max) {
                ret_val = sqrt(s2 * (1. + (x3max / s2) * (x3max * s3)));
            } else {
                ret_val = sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
            }
        } else {
            ret_val = x3max * sqrt(s3);
        }
    }
    return ret_val;

/*     last card of function enorm. */
#endif /* !USE_CBLAS */
} /* enorm_ */