integration.hpp

Go to the documentation of this file.

/*
 * $Id: integration.hpp 293 2012-12-17 20:27:36Z jdl3 $
 * Copyright (C) 2012, 2019, 2020 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
#ifndef CCGSL_INTEGRATION_HPP
#define CCGSL_INTEGRATION_HPP
 
#include<new>
#include<algorithm>
#include<gsl/gsl_integration.h>
#include"exception.hpp"
#include"function_scl.hpp"
 
namespace gsl {
  namespace integration {
    class workspace {
    public:
      workspace(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit workspace( size_t const n ){
    ccgsl_pointer = gsl_integration_workspace_alloc( n );
    // just plausibly we could allocate workspace but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_workspace_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit workspace( gsl_integration_workspace* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      workspace( workspace const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      workspace& operator=( workspace const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_workspace_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~workspace(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_workspace_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      workspace( workspace&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      workspace& operator=( workspace&& v ){
    workspace( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( workspace const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( workspace const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( workspace const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( workspace const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( workspace const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( workspace const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( workspace& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_workspace* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_workspace* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
    class fixed_workspace {
    public:
      fixed_workspace(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit fixed_workspace( const gsl_integration_fixed_type * type, const size_t n,
                const double a, const double b, const double alpha,
                const double beta ){
    ccgsl_pointer = gsl_integration_fixed_alloc( type, n, a, b, alpha, beta );
    // just plausibly we could allocate workspace but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_fixed_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit fixed_workspace( gsl_integration_fixed_workspace* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      fixed_workspace( fixed_workspace const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      fixed_workspace& operator=( fixed_workspace const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_fixed_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~fixed_workspace(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_fixed_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      fixed_workspace( fixed_workspace&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      fixed_workspace& operator=( fixed_workspace&& v ){
    fixed_workspace( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( fixed_workspace const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( fixed_workspace const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( fixed_workspace const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( fixed_workspace const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( fixed_workspace const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( fixed_workspace const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( fixed_workspace& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_fixed_workspace* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_fixed_workspace* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
    class qaws_table {
    public:
      qaws_table(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit qaws_table( double const alpha, double const beta, double const mu,
               double const nu ){
    ccgsl_pointer = gsl_integration_qaws_table_alloc( alpha, beta, mu, nu );
    // just plausibly we could allocate qaws_table but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_qaws_table_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit qaws_table( gsl_integration_qaws_table* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      qaws_table( qaws_table const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      qaws_table& operator=( qaws_table const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_qaws_table_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~qaws_table(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_qaws_table_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      qaws_table( qaws_table&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      qaws_table& operator=( qaws_table&& v ){
    qaws_table( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( qaws_table const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( qaws_table const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( qaws_table const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( qaws_table const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( qaws_table const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( qaws_table const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( qaws_table& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_qaws_table* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_qaws_table* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
    inline int qaws_table_set( qaws_table& t, double const alpha, double const beta,
                   int const mu, int const nu ){
      return gsl_integration_qaws_table_set( t.get(), alpha, beta, mu, nu ); }
 
    enum qawo_enum { 
      SINE   = GSL_INTEG_SINE,   
      COSINE = GSL_INTEG_COSINE  
    };
 
    class qawo_table {
    public:
      qawo_table(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit qawo_table( double const omega, double const L,
               qawo_enum const sine, size_t const n ){
    ccgsl_pointer
      = gsl_integration_qawo_table_alloc( omega, L,
                          static_cast<gsl_integration_qawo_enum>( sine ), n );
    // just plausibly we could allocate qawo_table but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_qawo_table_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit qawo_table( gsl_integration_qawo_table* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      qawo_table( qawo_table const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      qawo_table& operator=( qawo_table const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_qawo_table_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~qawo_table(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_qawo_table_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      qawo_table( qawo_table&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      qawo_table& operator=( qawo_table&& v ){
    qawo_table( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( qawo_table const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( qawo_table const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( qawo_table const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( qawo_table const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( qawo_table const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( qawo_table const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( qawo_table& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_qawo_table* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_qawo_table* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
    inline int qawo_table_set( qawo_table& t, double const omega, double const L,
                   qawo_enum const sine ){
      return gsl_integration_qawo_table_set( t.get(), omega, L,
                         static_cast<gsl_integration_qawo_enum>( sine ) ); }
    
    inline int qawo_table_set_length( qawo_table& t, double L ){
      return gsl_integration_qawo_table_set_length( t.get(), L ); }
 
    inline void qk15( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk15( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk15( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk15( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
    inline void qk21( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk21( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk21( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk21( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
    inline void qk31( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk31( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk31( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk31( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
    inline void qk41( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk41( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk41( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk41( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
    inline void qk51( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk51( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk51( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk51( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
    inline void qk61( function_scl const& f, double a, double b, double& result, double& abserr,
              double& resabs, double& resasc ){
      gsl_integration_qk61( &f, a, b, &result, &abserr, &resabs, &resasc ); }
 
#ifndef DOXYGEN_SKIP 
    inline void qk61( gsl_function const* f, double a, double b, double* result, double* abserr,
              double* resabs, double* resasc ){
      gsl_integration_qk61( f, a, b, result, abserr, resabs, resasc ); }
#endif /* DOXYGEN_SKIP */
 
#ifndef DOXYGEN_SKIP 
    inline void qcheb( function_scl& f, double a, double b, double* cheb12, double* cheb24 ){
      gsl_integration_qcheb( &f, a, b, cheb12, cheb24 ); }
    
    inline void qcheb( gsl_function* f, double a, double b, double* cheb12, double* cheb24 ){
      gsl_integration_qcheb( f, a, b, cheb12, cheb24 ); }
#endif /* DOXYGEN_SKIP */
 
    template<typename T>
    inline void qcheb( function_scl& f, double a, double b, T& cheb12, T& cheb24 ){
#ifndef GSL_RANGE_CHECK_OFF
      // Check that vectors are long enough
      if( cheb12.size() < 13 ){
    gsl_error( "expected cheb12 of length 13 (or more)",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return;
      }
      if( cheb24.size() < 25 ){
    gsl_error( "expected cheb24 of length 25 (or more)",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return;
      }
#endif
      gsl_integration_qcheb( &f, a, b, cheb12.data(), cheb24.data() );
    }
 
#ifndef DOXYGEN_SKIP 
    template<typename T>
    inline void qcheb( gsl_function* f, double a, double b, T& cheb12, T& cheb24 ){
#ifndef GSL_RANGE_CHECK_OFF
      // Check that vectors are long enough
      if( cheb12.size() < 13 ){
    gsl_error( "expected cheb12 of length 13 (or more)",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return;
      }
      if( cheb24.size() < 25 ){
    gsl_error( "expected cheb24 of length 25 (or more)",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return;
      }
#endif
      gsl_integration_qcheb( f, a, b, cheb12.data(), cheb24.data() );
    }
#endif /* DOXYGEN_SKIP */
 
    inline void qk( int const n, double const xgk[], double const wg[], double const wgk[],
            double fv1[], double fv2[], function_scl const& f, double a, double b,
            double& result, double& abserr, double& resabs, double& resasc ){
      gsl_integration_qk( n, xgk, wg, wgk, fv1, fv2, &f, a, b, &result, &abserr, &resabs, &resasc );
    }    
 
#ifndef DOXYGEN_SKIP 
    inline void qk( int const n, double const xgk[], double const wg[], double const wgk[],
            double fv1[], double fv2[], gsl_function const* f, double a, double b,
            double* result, double* abserr, double* resabs, double* resasc ){
      gsl_integration_qk( n, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc );
    }    
#endif /* DOXYGEN_SKIP */
 
    template<typename T>
    inline void qk( T const& xgk, T const& wg, T const& wgk,
            T& fv1, T& fv2, function_scl const& f, double a, double b,
            double& result, double& abserr, double& resabs, double& resasc ){
#ifndef GSL_RANGE_CHECK_OFF
      size_t const n = xgk.size();
      // Check that vectors are same length
      if( n == 0 ){ gsl_error( "expected xgk of nonzero size",
                   __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( wg.size() != n ){ gsl_error( "size mismatch: xgk and wg",
                       __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( wgk.size() != n ){ gsl_error( "size mismatch: xgk and wgk",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( fv1.size() != n ){ gsl_error( "size mismatch: xgk and fv1",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( fv2.size() != n ){ gsl_error( "size mismatch: xgk and fv2",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
#endif
      gsl_integration_qk( n, xgk.data(), wg.data(), wgk.data(), fv1.data(), fv2.data(),
              &f, a, b, &result, &abserr, &resabs, &resasc );
    }
 
#ifndef DOXYGEN_SKIP 
    template<typename T>
    inline void qk( T const& xgk, T const& wg, T const& wgk,
            T& fv1, T& fv2, gsl_function const* f, double a, double b,
            double* result, double* abserr, double* resabs, double* resasc ){
#ifndef GSL_RANGE_CHECK_OFF
      size_t const n = xgk.size();
      // Check that vectors are same length
      if( n == 0 ){ gsl_error( "expected xgk of nonzero size",
                   __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( wg.size() != n ){ gsl_error( "size mismatch: xgk and wg",
                       __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( wgk.size() != n ){ gsl_error( "size mismatch: xgk and wgk",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( fv1.size() != n ){ gsl_error( "size mismatch: xgk and fv1",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
      if( fv2.size() != n ){ gsl_error( "size mismatch: xgk and fv2",
                    __FILE__, __LINE__, exception::GSL_EFAILED ); return; }
#endif
      gsl_integration_qk( n, xgk.data(), wg.data(), wgk.data(), fv1.data(), fv2.data(),
              f, a, b, result, abserr, resabs, resasc );
    }
#endif /* DOXYGEN_SKIP */
 
    inline int qng( function_scl const& f, double a, double b, double epsabs, double epsrel,
            double& result, double& abserr, size_t& neval ){
      return gsl_integration_qng( &f, a, b, epsabs, epsrel, &result, &abserr, &neval ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qng( gsl_function const* f, double a, double b, double epsabs, double epsrel,
            double* result, double* abserr, size_t* neval ){
      return gsl_integration_qng( f, a, b, epsabs, epsrel, result, abserr, neval ); }
#endif /* DOXYGEN_SKIP */
    
    enum qag_key { 
      GAUSS15   = 1,   
      GAUSS21   = 2,   
      GAUSS31   = 3,   
      GAUSS41   = 4,   
      GAUSS51   = 5,   
      GAUSS61   = 6    
    };
 
    inline int qag( function_scl& f, double const a, double const b, double const epsabs,
            double const epsrel, size_t const limit, qag_key const key, workspace& workspace,
            double& result, double& abserr ){
#ifndef GSL_RANGE_CHECK_OFF
      if( limit > (workspace.get() == 0 ? 0 : workspace.get()->limit) ){
    gsl_error( "limit must not exceed size of workspace",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return GSL_EFAILED; }
#endif
      return gsl_integration_qag( &f, a, b, epsabs, epsrel, limit, static_cast<int>( key ),
                  workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qag( gsl_function* f, double const a, double const b, double const epsabs,
            double const epsrel, size_t const limit, qag_key const key, workspace& workspace,
            double* result, double* abserr ){
#ifndef GSL_RANGE_CHECK_OFF
      if( limit > (workspace.get() == 0 ? 0 : workspace.get()->limit) ){
    gsl_error( "limit must not exceed size of workspace",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return GSL_EFAILED; }
#endif
      return gsl_integration_qag( f, a, b, epsabs, epsrel, limit, static_cast<int>( key ),
                  workspace.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
    inline int qagi( function_scl& f, double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qagi( &f, epsabs, epsrel, limit, workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qagi( gsl_function* f, double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qagi( f, epsabs, epsrel, limit, workspace.get(), result, abserr ); }
 
#endif /* DOXYGEN_SKIP */
    inline int qagiu( function_scl& f, double const a, double const epsabs, double const epsrel,
              size_t const limit, workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qagiu( &f, a, epsabs, epsrel, limit, workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qagiu( gsl_function* f, double const a, double const epsabs, double const epsrel,
              size_t const limit, workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qagiu( f, a, epsabs, epsrel, limit, workspace.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qagil( function_scl& f, double const b, double const epsabs, double const epsrel,
              size_t const limit, workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qagil( &f, b, epsabs, epsrel, limit, workspace.get(),
                    &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qagil( gsl_function* f, double const b, double const epsabs, double const epsrel,
              size_t const limit, workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qagil( f, b, epsabs, epsrel, limit, workspace.get(),
                    result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qags( function_scl& f, double const a, double const b, double const epsabs,
             double const epsrel, size_t const limit, workspace& workspace, double& result,
             double& abserr ){
      return gsl_integration_qags( &f, a, b, epsabs, epsrel, limit, workspace.get(),
                   &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qags( gsl_function* f, double const a, double const b, double const epsabs,
             double const epsrel, size_t const limit, workspace& workspace, double* result,
             double* abserr ){
      return gsl_integration_qags( f, a, b, epsabs, epsrel, limit, workspace.get(),
                   result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qagp( function_scl& f, double const* pts, size_t const npts, double const epsabs,
             double const epsrel,
             size_t const limit, workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qagp( &f, const_cast<double*>( pts ), npts, epsabs, epsrel, limit,
                   workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qagp( gsl_function* f, double const* pts, size_t const npts, double const epsabs,
             double const epsrel,
             size_t const limit, workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qagp( f, const_cast<double*>( pts ), npts, epsabs, epsrel, limit,
                   workspace.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    template<typename T>
    inline int qagp( function_scl& f, T const& pts, double const epsabs, double const epsrel,
             size_t const limit, workspace& workspace, double& result, double& abserr ){
      size_t const npts = pts.size();
#ifndef GSL_RANGE_CHECK_OFF
      // Check that vectors are long enough
      if( pts.size() < 2 ){
    gsl_error( "expected pts of length 2 or greater",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return GSL_EFAILED;
      }
#endif
      return gsl_integration_qagp( &f, const_cast<double*>( pts.data() ), npts, epsabs, epsrel,
                   limit, workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    template<typename T>
    inline int qagp( gsl_function* f, T const& pts, double const epsabs, double const epsrel,
             size_t const limit, workspace& workspace, double* result, double* abserr ){
      size_t const npts = pts.size();
#ifndef GSL_RANGE_CHECK_OFF
      // Check that vectors are long enough
      if( pts.size() < 2 ){
    gsl_error( "expected pts of length 2 or greater",
           __FILE__, __LINE__, exception::GSL_EFAILED );
    return GSL_EFAILED;
      }
#endif
      return gsl_integration_qagp( f, const_cast<double*>( pts.data() ), npts, epsabs, epsrel,
                   limit, workspace.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qawc( function_scl& f, double const a, double const b, double const c,
             double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qawc( &f, a, b, c, epsabs, epsrel, limit, workspace.get(),
                   &result, &abserr ); }
                         
#ifndef DOXYGEN_SKIP 
    inline int qawc( gsl_function* f, double const a, double const b, double const c,
             double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qawc( f, a, b, c, epsabs, epsrel, limit, workspace.get(),
                   result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qaws( function_scl& f, double const a, double const b, qaws_table& t,
             double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double& result, double& abserr ){
      return gsl_integration_qaws( &f, a, b, t.get(), epsabs, epsrel, limit,
                   workspace.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qaws( gsl_function* f, double const a, double const b, qaws_table& t,
             double const epsabs, double const epsrel, size_t const limit,
             workspace& workspace, double* result, double* abserr ){
      return gsl_integration_qaws( f, a, b, t.get(), epsabs, epsrel, limit,
                   workspace.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qawo( function_scl& f, double const a, double const epsabs, double const epsrel,
             size_t const limit, workspace& workspace, qawo_table& wf,
             double& result, double& abserr ){
      return gsl_integration_qawo( &f, a, epsabs, epsrel, limit, workspace.get(), wf.get(),
                   &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qawo( gsl_function* f, double const a, double const epsabs, double const epsrel,
             size_t const limit, workspace& workspace, qawo_table& wf,
             double* result, double* abserr ){
      return gsl_integration_qawo( f, a, epsabs, epsrel, limit, workspace.get(), wf.get(),
                   result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    inline int qawf( function_scl& f, double const a, double const epsabs, size_t const limit,
             integration::workspace& workspace, integration::workspace& cycle_workspace,
             qawo_table& wf, double& result, double& abserr ){
      return gsl_integration_qawf( &f, a, epsabs, limit, workspace.get(), cycle_workspace.get(),
                   wf.get(), &result, &abserr ); }
 
#ifndef DOXYGEN_SKIP 
    inline int qawf( gsl_function* f, double const a, double const epsabs, size_t const limit,
             integration::workspace& workspace, integration::workspace& cycle_workspace,
             qawo_table& wf, double* result, double* abserr ){
      return gsl_integration_qawf( f, a, epsabs, limit, workspace.get(), cycle_workspace.get(),
                   wf.get(), result, abserr ); }
#endif /* DOXYGEN_SKIP */
 
    class glfixed_table {
    public:
      glfixed_table(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit glfixed_table( size_t const n ){
    ccgsl_pointer = gsl_integration_glfixed_table_alloc( n );
    // just plausibly we could allocate glfixed_table but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_glfixed_table_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit glfixed_table( gsl_integration_glfixed_table* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      glfixed_table( glfixed_table const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      glfixed_table& operator=( glfixed_table const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_glfixed_table_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~glfixed_table(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_glfixed_table_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      glfixed_table( glfixed_table&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      glfixed_table& operator=( glfixed_table&& v ){
    glfixed_table( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( glfixed_table const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( glfixed_table const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( glfixed_table const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( glfixed_table const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( glfixed_table const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( glfixed_table const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( glfixed_table& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_glfixed_table* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_glfixed_table* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
    inline double glfixed( function_scl const& f, double const a, double const b,
               glfixed_table const& t ){
      return gsl_integration_glfixed( &f, a, b, t.get() ); }
    
#ifndef DOXYGEN_SKIP 
    inline double glfixed( gsl_function const* f, double const a, double const b,
               glfixed_table const& t ){
      return gsl_integration_glfixed( f, a, b, t.get() ); }
#endif /* DOXYGEN_SKIP */
 
#ifndef DOXYGEN_SKIP 
    inline int glfixed_point( double const a, double const b, size_t const i, double* xi, double* wi,
                  glfixed_table const& t ){
      return gsl_integration_glfixed_point( a, b, i, xi, wi, t.get() ); } 
#endif /* DOXYGEN_SKIP */
    inline int glfixed_point( double const a, double const b, size_t const i,
                  double& xi, double& wi, glfixed_table const& t ){
      return gsl_integration_glfixed_point( a, b, i, &xi, &wi, t.get() ); } 
 
    class cquad_workspace {
    public:
      cquad_workspace(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit cquad_workspace( size_t const n ){
    ccgsl_pointer = gsl_integration_cquad_workspace_alloc( n );
    // just plausibly we could allocate cquad_workspace but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_cquad_workspace_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit cquad_workspace( gsl_integration_cquad_workspace* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      cquad_workspace( cquad_workspace const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      cquad_workspace& operator=( cquad_workspace const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_cquad_workspace_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~cquad_workspace(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_cquad_workspace_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      cquad_workspace( cquad_workspace&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      cquad_workspace& operator=( cquad_workspace&& v ){
    cquad_workspace( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( cquad_workspace const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( cquad_workspace const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( cquad_workspace const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( cquad_workspace const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( cquad_workspace const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( cquad_workspace const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( cquad_workspace& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_cquad_workspace* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_cquad_workspace* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
 
    inline int cquad( function_scl const& f, double const a, double const b, double const epsabs,
              double const epsrel, cquad_workspace& ws,
              double& result, double& abserr, size_t& nevals ){
      return gsl_integration_cquad( &f, a, b, epsabs, epsrel, ws.get(), &result, &abserr, &nevals ); }
    
#ifndef DOXYGEN_SKIP 
    inline int cquad( gsl_function const* f, double const a, double const b, double const epsabs,
              double const epsrel, cquad_workspace& ws,
              double* result, double* abserr, size_t* nevals ){
      return gsl_integration_cquad( f, a, b, epsabs, epsrel, ws.get(), result, abserr, nevals ); }
#endif /* DOXYGEN_SKIP */
 
    class romberg_workspace {
    public:
      romberg_workspace(){
    ccgsl_pointer = 0;
    count = 0; // initially nullptr will do
      }
      // Refines random access container
      // Refines assignable
      explicit romberg_workspace( size_t const n ){
    ccgsl_pointer = gsl_integration_romberg_alloc( n );
    // just plausibly we could allocate romberg_workspace but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_integration_romberg_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit romberg_workspace( gsl_integration_romberg_workspace* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      romberg_workspace( romberg_workspace const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      romberg_workspace& operator=( romberg_workspace const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_integration_romberg_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~romberg_workspace(){
    if( count == 0 or --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_integration_romberg_free( ccgsl_pointer );
      delete count;
    }
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      romberg_workspace( romberg_workspace&& v ) : ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
    std::swap( count, v.count );
    v.ccgsl_pointer = nullptr;
      }
      romberg_workspace& operator=( romberg_workspace&& v ){
    romberg_workspace( std::move( v ) ).swap( *this );
    return *this;
      }
#endif
      // Refines equality comparable
      // == operator
      bool operator==( romberg_workspace const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( romberg_workspace const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( romberg_workspace const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( romberg_workspace const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( romberg_workspace const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( romberg_workspace const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( romberg_workspace& v ){
    std::swap( ccgsl_pointer, v.ccgsl_pointer );
    std::swap( count, v.count );
      }
    private:
      gsl_integration_romberg_workspace* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_integration_romberg_workspace* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
    inline int romberg( function_scl const& f, double const a, double const b,
            double const epsabs, double const epsrel, double& result,
            size_t& neval, romberg_workspace& w ){
      return gsl_integration_romberg( &f, a, b, epsabs, epsrel, &result, &neval, w.get() ); }
#ifndef DOXYGEN_SKIP 
    inline int romberg( gsl_function const* f, double const a, double const b,
            double const epsabs, double const epsrel, double* result,
            size_t* neval, romberg_workspace& w ){
      return gsl_integration_romberg( f, a, b, epsabs, epsrel, result, neval, w.get() ); }
#endif /* DOXYGEN_SKIP */
    inline size_t fixed_n( fixed_workspace const& w ){
      return gsl_integration_fixed_n( w.get() ); }
    inline double* fixed_nodes( fixed_workspace const& w ){
      return gsl_integration_fixed_nodes( w.get() ); }
    inline double* fixed_weights( fixed_workspace const& w ){
      return gsl_integration_fixed_weights( w.get() ); }
 
    inline int fixed( function_scl const& func, double* result, fixed_workspace& w ){
      return gsl_integration_fixed( &func, result, w.get() ); }
#ifndef DOXYGEN_SKIP 
    inline int fixed( gsl_function const* func, double* result, fixed_workspace& w ){
      return gsl_integration_fixed( func, result, w.get() ); }
#endif /* DOXYGEN_SKIP */
  }
}
 
#endif