/************************************************************************/
/*                                                                      */
/*               Copyright 1998-2002 by Ullrich Koethe                  */
/*       Cognitive Systems Group, University of Hamburg, Germany        */
/*                                                                      */
/*    This file is part of the VIGRA computer vision library.           */
/*    ( Version 1.3.2, Jan 27 2005 )                                    */
/*    You may use, modify, and distribute this software according       */
/*    to the terms stated in the LICENSE file included in               */
/*    the VIGRA distribution.                                           */
/*                                                                      */
/*    The VIGRA Website is                                              */
/*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
/*    Please direct questions, bug reports, and contributions to        */
/*        koethe@informatik.uni-hamburg.de                              */
/*                                                                      */
/*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
/*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
/*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
/*                                                                      */
/************************************************************************/

#ifndef VIGRA_MATHUTIL_HXX
#define VIGRA_MATHUTIL_HXX

#include <cmath>
#include <cstdlib>
#include "vigra/config.hxx"
#include "vigra/numerictraits.hxx"

/*! \page MathConstants Mathematical Constants

    <TT>M_PI, M_SQRT2</TT>

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"

    Since <TT>M_PI</TT> and <TT>M_SQRT2</TT> are not officially standardized,
    we provide definitions here for those compilers that don't support them.

    \code
    #ifndef M_PI
    #    define M_PI     3.14159265358979323846
    #endif

    #ifndef M_SQRT2
    #    define M_SQRT2  1.41421356237309504880
    #endif
    \endcode
*/
#ifndef M_PI
#    define M_PI     3.14159265358979323846
#endif

#ifndef M_SQRT2
#    define M_SQRT2  1.41421356237309504880
#endif

namespace vigra {

#ifndef __sun 

/** \addtogroup MathFunctions Mathematical Functions

    Useful mathematical functions and functors.
*/
//@{
/*! The error function.

    With the exception of Solaris (where <tt>erf()</tt> is provided as an extension of the 
    C math library), VIGRA implements <tt>erf()</tt> as an approximation of the error 
    function
    
    \f[
        \mbox{erf}(x) = \int_0^x e^{-x^2} dx
    \f]
    
    according to the formula given in Press et al. "Numerical Recipes".

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    Namespace: vigra
*/
template <class T>
double erf(T x)
{
    double t = 1.0/(1.0+0.5*VIGRA_CSTD::fabs(x));
    double ans = t*VIGRA_CSTD::exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+
                                    t*(0.09678418+t*(-0.18628806+t*(0.27886807+
                                    t*(-1.13520398+t*(1.48851587+t*(-0.82215223+
                                    t*0.17087277)))))))));
    if (x >= 0.0)
        return 1.0 - ans;
    else
        return ans - 1.0;
}

#else

using VIGRA_CSTD::erf;

#endif

// import functions into namespace vigra which VIGRA is going to overload
using std::abs;  
using VIGRA_CSTD::pow;  
using VIGRA_CSTD::floor;  
using VIGRA_CSTD::ceil;  

/*! The square function.

    sq(x) is needed so often that it makes sense to define it as a function.

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    Namespace: vigra
*/
template <class T>
inline 
typename NumericTraits<T>::Promote sq(T t)
{
    return t*t;
}

#ifdef VIGRA_NO_HYPOT
/*! Compute the Euclidean distance.

    The  hypot()  function  returns  the  sqrt(a*a  +  b*b).
    It is implemented in a way that minimizes round-off error.

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    Namespace: vigra
*/
template <class T>
T hypot(T a, T b) 
{ 
    T absa = abs(a), absb = abs(b);
    if (absa > absb) 
        return absa * VIGRA_CSTD::sqrt(1.0 + sq(absb/absa)); 
    else 
        return absb == NumericTraits<T>::zero()
                   ? NumericTraits<T>::zero()
                   : absb * VIGRA_CSTD::sqrt(1.0 + sq(absa/absb)); 
}

#else

using ::hypot;

#endif

/*! The sign function.

    Returns 1, 0, or -1 depending on the signm of \a t.

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    Namespace: vigra
*/
template <class T>
T sign(T t) 
{ 
    return t > NumericTraits<T>::zero()
               ? NumericTraits<T>::one()
               : t < NumericTraits<T>::zero()
                    ? -NumericTraits<T>::one()
                    : NumericTraits<T>::zero();
}

/*! The binary sign function.

    Transfers the sign of \a t2 to \a t1.

    <b>\#include</b> "<a href="mathutil_8hxx-source.html">vigra/mathutil.hxx</a>"<br>
    Namespace: vigra
*/
template <class T1, class T2>
T1 sign(T1 t1, T2 t2) 
{ 
    return t2 >= NumericTraits<T2>::zero()
               ? abs(t1)
               : -abs(t1);
}

//@}

} // namespace vigra


#endif /* VIGRA_MATHUTIL_HXX */
