text: //==============================================================================
  //         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
  //         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
  //
  //          Distributed under the Boost Software License, Version 1.0.
  //                 See accompanying file LICENSE.txt or copy at
  //                     http://www.boost.org/LICENSE_1_0.txt
  //==============================================================================
  #ifndef NT2_INTEGRATION_FUNCTIONS_QUAD_HPP_INCLUDED
  #define NT2_INTEGRATION_FUNCTIONS_QUAD_HPP_INCLUDED
  
  /*!
   * \file
   * \brief Defines and implements the nt2::quad function
   */
  #include <nt2/integration/interface.hpp>
  
  namespace nt2
  {
    namespace tag
    {
      struct quad_ : ext::unspecified_<quad_>
      {
        typedef ext::unspecified_<quad_> parent;
      };
  
      // definition  of abstol constant for quad method
      BOOST_SIMD_CONSTANT_REGISTER( Quadabstol, double
                                    , 0, 0x3a83126f             //1.0e-3
                                    , 0x3eb0c6f7a0b5ed8dll      //1.0e-6
        );
    }
  
    BOOST_SIMD_CONSTANT_IMPLEMENTATION(tag::Quadabstol, Quadabstol);
  
    // specialization of abstol for quad method
    template<class T, class V> struct integ_params<T, V, tag::quad_>
    : integ_params<T, V, void>
    {
      typedef typename nt2::integ_params<T, V, void>::real_t real_t;
      static real_t        abstol(){return Quadabstol<real_t>(); }
      static size_t    maxintvcnt(){return 512; }
    };
  
    //============================================================================
    /*!
     * Apply quad algorithm to integrate a function over a real interval
     *
     * \param func  Function to optimize
     * \param x    required points in the interval or 2 abscissae a and b
     * \param opt   Options pack related to the tolerance handling
     *
     * \return  a tuple containing the results of the integration, the last error value,
     * the number of required function evaluation and a boolean
     * notifying success of the whole process.
     *
     *    q = quad(fun,a,b) tries to approximate the integral of scalar-valued
     *    function fun from a to b to within a default error of nt2::Quadabstol<real_t>()
     *    where real_t can be float (1.0e-3) or double (1.0e-6) using recursive
     *    adaptive simpson quadrature. fun is a function handle. the function
     *    y=fun(x) should accept a vector argument x and return a vector result
     *    y, the integrand evaluated at each element of x.
     *
     * This routine mimics matlab quad, but has some behavioural differences.
     * including options passing and obtaining results
     *
     * 1 -When matlab subdivides the initial interval [a, b] each subinterval is
     *    given the same tolerance as the initial one which can not insure a total
     *    result with the correct tolerance.
     *    Th nt2 quad function uses a tolerance proportional to the length on such
     *    cases in order the sum of all the partial errors stays above the global
     *    absolute error tolerance.
     *    As this can bring to needlessly bring to too small intervals a theshhold
     *    depending of the maximum of subdivision allowed (maxintvcnt_) is used.
     *
     * 2- a and b can be singular for f,  but contrarily to matlab, there is no
     *    automatic detection and the options singular_a_ and singular_b_ are by
     *    default set to false
     *
     * 3- as in matlab on can give way points through the computation must go
     *    a and b are always added (but never duplicated)
     *    the way points mustn't be singular and in proper order
     *    if the option return_waypoints is set to true the result is a vector of
     *    cumulated integrals from a to wi. Peculiarly result(begin_) is 0 and
     *    result(end_) is the integral from a to b.
     *
     * 4- the abscissae can be complex and the integral a path integral through the
     *    lines following the way points in their given order.
     */
    //============================================================================
  
  
    template<class F, class X> BOOST_FORCEINLINE
    typename details::integration<F, X, tag::quad_>::result_type
    quad(F f, X x)
    {
      return details::integration<F, X, tag::quad_>::call(f, x);
    }
  
    template<class F, class X, class Xpr> BOOST_FORCEINLINE
    typename details::integration<F, X, tag::quad_>::result_type
    quad(F f, X x, nt2::details::option_expr<Xpr> const& opt)
    {
      return details::integration<F, X, tag::quad_>::call(f, x, opt);
    }
  
    template<class F, class A> BOOST_FORCEINLINE
    typename details::integration<F, typename details::h2_t<A>::ab_t, tag::quad_>::result_type
    quad(F f, A a, A b)
    {
      typedef typename details::h2_t<A>::ab_t ab_t;
      return details::integration<F, ab_t, tag::quad_>::call(f, nt2::cath(a, b));
    }
  
    template<class F, class A, class Xpr> BOOST_FORCEINLINE
    typename details::integration<F, typename details::h2_t<A>::ab_t, tag::quad_>::result_type
    quad(F f, A a, A b, nt2::details::option_expr<Xpr> const& opt)
    {
      typedef typename details::h2_t<A>::ab_t ab_t;
      return details::integration<F, ab_t, tag::quad_>::call(f, nt2::cath(a, b), opt);
    }
  
  
  }
  
  #endif
  
  