diff --git a/patched/any.hpp b/patched/any.hpp new file mode 100644 index 0000000..7be5231 --- /dev/null +++ b/patched/any.hpp @@ -0,0 +1,239 @@ +// See http://www.boost.org/libs/any for Documentation. + +#ifndef BOOST_ANY_INCLUDED +#define BOOST_ANY_INCLUDED + +// what: variant type boost::any +// who: contributed by Kevlin Henney, +// with features contributed and bugs found by +// Ed Brey, Mark Rodgers, Peter Dimov, and James Curran +// when: July 2001 +// where: tested with BCC 5.5, MSVC 6.0, and g++ 2.95 + +#include +#include + +#include "boost/config.hpp" +#include +#include +#include +#include +#include + +namespace boost +{ + class any + { + public: // structors + + any() + : content(0) + { + } + + template + any(const ValueType & value) + : content(new holder(value)) + { + } + + any(const any & other) + : content(other.content ? other.content->clone() : 0) + { + } + + ~any() + { + delete content; + } + + public: // modifiers + + any & swap(any & rhs) + { + std::swap(content, rhs.content); + return *this; + } + + template + any & operator=(const ValueType & rhs) + { + any(rhs).swap(*this); + return *this; + } + + any & operator=(any rhs) + { + rhs.swap(*this); + return *this; + } + + public: // queries + + bool empty() const + { + return !content; + } + + type_index type() const + { + return content ? content->type() : type_id(); + } + +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS + private: // types +#else + public: // types (public so any_cast can be non-friend) +#endif + + class placeholder + { + public: // structors + + virtual ~placeholder() + { + } + + public: // queries + + virtual const type_index & type() const = 0; + + virtual placeholder * clone() const = 0; + + }; + + template + class holder : public placeholder + { + public: // structors + + holder(const ValueType & value) + : held(value) + { + } + + public: // queries + + virtual const type_index & type() const + { + return type_id(); + } + + virtual placeholder * clone() const + { + return new holder(held); + } + + public: // representation + + ValueType held; + + private: // intentionally left unimplemented + holder & operator=(const holder &); + }; + +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS + + private: // representation + + template + friend ValueType * any_cast(any *); + + template + friend ValueType * unsafe_any_cast(any *); + +#else + + public: // representation (public so any_cast can be non-friend) + +#endif + + placeholder * content; + + }; + + class bad_any_cast : public std::bad_cast + { + public: + virtual const char * what() const throw() + { + return "boost::bad_any_cast: " + "failed conversion using boost::any_cast"; + } + }; + + template + ValueType * any_cast(any * operand) + { + return operand && + operand->type() == type_id() + ? &static_cast *>(operand->content)->held + : 0; + } + + template + inline const ValueType * any_cast(const any * operand) + { + return any_cast(const_cast(operand)); + } + + template + ValueType any_cast(any & operand) + { + typedef BOOST_DEDUCED_TYPENAME remove_reference::type nonref; + +#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION + // If 'nonref' is still reference type, it means the user has not + // specialized 'remove_reference'. + + // Please use BOOST_BROKEN_COMPILER_TYPE_TRAITS_SPECIALIZATION macro + // to generate specialization of remove_reference for your class + // See type traits library documentation for details + BOOST_STATIC_ASSERT(!is_reference::value); +#endif + + nonref * result = any_cast(&operand); + if(!result) + boost::throw_exception(bad_any_cast()); + return *result; + } + + template + inline ValueType any_cast(const any & operand) + { + typedef BOOST_DEDUCED_TYPENAME remove_reference::type nonref; + +#ifdef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION + // The comment in the above version of 'any_cast' explains when this + // assert is fired and what to do. + BOOST_STATIC_ASSERT(!is_reference::value); +#endif + + return any_cast(const_cast(operand)); + } + + // Note: The "unsafe" versions of any_cast are not part of the + // public interface and may be removed at any time. They are + // required where we know what type is stored in the any and can't + // use type_id<>() comparison, e.g., when our types may travel across + // different shared libraries. + template + inline ValueType * unsafe_any_cast(any * operand) + { + return &static_cast *>(operand->content)->held; + } + + template + inline const ValueType * unsafe_any_cast(const any * operand) + { + return unsafe_any_cast(const_cast(operand)); + } +} + +// Copyright Kevlin Henney, 2000, 2001, 2002. All rights reserved. +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#endif diff --git a/patched/detail/REMOVE_sp_typeinfo.hpp b/patched/detail/REMOVE_sp_typeinfo.hpp new file mode 100644 index 0000000..50da7aa --- /dev/null +++ b/patched/detail/REMOVE_sp_typeinfo.hpp @@ -0,0 +1,135 @@ +#ifndef BOOST_DETAIL_SP_TYPEINFO_HPP_INCLUDED +#define BOOST_DETAIL_SP_TYPEINFO_HPP_INCLUDED + +// MS compatible compilers support #pragma once + +#if defined(_MSC_VER) && (_MSC_VER >= 1020) +# pragma once +#endif + +// detail/sp_typeinfo.hpp +// +// Copyright 2007 Peter Dimov +// +// Distributed under the Boost Software License, Version 1.0. +// See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#include + +#if defined( BOOST_NO_TYPEID ) + +#include +#include + +namespace boost +{ + +namespace detail +{ + +class sp_typeinfo +{ +private: + + sp_typeinfo( sp_typeinfo const& ); + sp_typeinfo& operator=( sp_typeinfo const& ); + + char const * name_; + +public: + + explicit sp_typeinfo( char const * name ): name_( name ) + { + } + + bool operator==( sp_typeinfo const& rhs ) const + { + return this == &rhs; + } + + bool operator!=( sp_typeinfo const& rhs ) const + { + return this != &rhs; + } + + bool before( sp_typeinfo const& rhs ) const + { + return std::less< sp_typeinfo const* >()( this, &rhs ); + } + + char const* name() const + { + return name_; + } +}; + +template struct sp_typeid_ +{ + static sp_typeinfo ti_; + + static char const * name() + { + return BOOST_CURRENT_FUNCTION; + } +}; + +#if defined(__SUNPRO_CC) +// see #4199, the Sun Studio compiler gets confused about static initialization +// constructor arguments. But an assignment works just fine. +template sp_typeinfo sp_typeid_< T >::ti_ = sp_typeid_< T >::name(); +#else +template sp_typeinfo sp_typeid_< T >::ti_(sp_typeid_< T >::name()); +#endif + +template struct sp_typeid_< T & >: sp_typeid_< T > +{ +}; + +template struct sp_typeid_< T const >: sp_typeid_< T > +{ +}; + +template struct sp_typeid_< T volatile >: sp_typeid_< T > +{ +}; + +template struct sp_typeid_< T const volatile >: sp_typeid_< T > +{ +}; + +} // namespace detail + +} // namespace boost + +#define BOOST_SP_TYPEID(T) (boost::detail::sp_typeid_::ti_) + +#else + +#include + +namespace boost +{ + +namespace detail +{ + +#if defined( BOOST_NO_STD_TYPEINFO ) + +typedef ::type_info sp_typeinfo; + +#else + +typedef std::type_info sp_typeinfo; + +#endif + +} // namespace detail + +} // namespace boost + +#define BOOST_SP_TYPEID(T) typeid(T) + +#endif + +#endif // #ifndef BOOST_DETAIL_SP_TYPEINFO_HPP_INCLUDED diff --git a/patched/function/function_base.hpp b/patched/function/function_base.hpp new file mode 100644 index 0000000..a95792b --- /dev/null +++ b/patched/function/function_base.hpp @@ -0,0 +1,874 @@ +// Boost.Function library + +// Copyright Douglas Gregor 2001-2006 +// Copyright Emil Dotchevski 2007 +// Use, modification and distribution is subject to the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// For more information, see http://www.boost.org + +#ifndef BOOST_FUNCTION_BASE_HEADER +#define BOOST_FUNCTION_BASE_HEADER + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifndef BOOST_NO_SFINAE +# include "boost/utility/enable_if.hpp" +#else +# include "boost/mpl/bool.hpp" +#endif +#include +#include + +#if defined(BOOST_MSVC) +# pragma warning( push ) +# pragma warning( disable : 4793 ) // complaint about native code generation +# pragma warning( disable : 4127 ) // "conditional expression is constant" +#endif + +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 || defined(__ICL) && __ICL <= 600 || defined(__MWERKS__) && __MWERKS__ < 0x2406 && !defined(BOOST_STRICT_CONFIG) +# define BOOST_FUNCTION_TARGET_FIX(x) x +#else +# define BOOST_FUNCTION_TARGET_FIX(x) +#endif // not MSVC + +#if !BOOST_WORKAROUND(__BORLANDC__, < 0x5A0) +# define BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor,Type) \ + typename ::boost::enable_if_c<(::boost::type_traits::ice_not< \ + (::boost::is_integral::value)>::value), \ + Type>::type +#else +// BCC doesn't recognize this depends on a template argument and complains +// about the use of 'typename' +# define BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor,Type) \ + ::boost::enable_if_c<(::boost::type_traits::ice_not< \ + (::boost::is_integral::value)>::value), \ + Type>::type +#endif + +namespace boost { + namespace detail { + namespace function { + class X; + + /** + * A buffer used to store small function objects in + * boost::function. It is a union containing function pointers, + * object pointers, and a structure that resembles a bound + * member function pointer. + */ + union function_buffer + { + // For pointers to function objects + mutable void* obj_ptr; + + // For pointers to std::type_info objects + struct type_t { + // (get_functor_type_tag, check_functor_type_tag). + type_index type; + + // Whether the type is const-qualified. + bool const_qualified; + // Whether the type is volatile-qualified. + bool volatile_qualified; + } type; + + // For function pointers of all kinds + mutable void (*func_ptr)(); + + // For bound member pointers + struct bound_memfunc_ptr_t { + void (X::*memfunc_ptr)(int); + void* obj_ptr; + } bound_memfunc_ptr; + + // For references to function objects. We explicitly keep + // track of the cv-qualifiers on the object referenced. + struct obj_ref_t { + mutable void* obj_ptr; + bool is_const_qualified; + bool is_volatile_qualified; + } obj_ref; + + // To relax aliasing constraints + mutable char data; + }; + + /** + * The unusable class is a placeholder for unused function arguments + * It is also completely unusable except that it constructable from + * anything. This helps compilers without partial specialization to + * handle Boost.Function objects returning void. + */ + struct unusable + { + unusable() {} + template unusable(const T&) {} + }; + + /* Determine the return type. This supports compilers that do not support + * void returns or partial specialization by silently changing the return + * type to "unusable". + */ + template struct function_return_type { typedef T type; }; + + template<> + struct function_return_type + { + typedef unusable type; + }; + + // The operation type to perform on the given functor/function pointer + enum functor_manager_operation_type { + clone_functor_tag, + move_functor_tag, + destroy_functor_tag, + check_functor_type_tag, + get_functor_type_tag + }; + + // Tags used to decide between different types of functions + struct function_ptr_tag {}; + struct function_obj_tag {}; + struct member_ptr_tag {}; + struct function_obj_ref_tag {}; + + template + class get_function_tag + { + typedef typename mpl::if_c<(is_pointer::value), + function_ptr_tag, + function_obj_tag>::type ptr_or_obj_tag; + + typedef typename mpl::if_c<(is_member_pointer::value), + member_ptr_tag, + ptr_or_obj_tag>::type ptr_or_obj_or_mem_tag; + + typedef typename mpl::if_c<(is_reference_wrapper::value), + function_obj_ref_tag, + ptr_or_obj_or_mem_tag>::type or_ref_tag; + + public: + typedef or_ref_tag type; + }; + + // The trivial manager does nothing but return the same pointer (if we + // are cloning) or return the null pointer (if we are deleting). + template + struct reference_manager + { + static inline void + manage(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op) + { + switch (op) { + case clone_functor_tag: + out_buffer.obj_ref = in_buffer.obj_ref; + return; + + case move_functor_tag: + out_buffer.obj_ref = in_buffer.obj_ref; + in_buffer.obj_ref.obj_ptr = 0; + return; + + case destroy_functor_tag: + out_buffer.obj_ref.obj_ptr = 0; + return; + + case check_functor_type_tag: + // Check whether we have the same type. We can add + // cv-qualifiers, but we can't take them away. + if (out_buffer.type.type == type_id()) + && (!in_buffer.obj_ref.is_const_qualified + || out_buffer.type.const_qualified) + && (!in_buffer.obj_ref.is_volatile_qualified + || out_buffer.type.volatile_qualified)) + out_buffer.obj_ptr = in_buffer.obj_ref.obj_ptr; + else + out_buffer.obj_ptr = 0; + return; + + case get_functor_type_tag: + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = in_buffer.obj_ref.is_const_qualified; + out_buffer.type.volatile_qualified = in_buffer.obj_ref.is_volatile_qualified; + return; + } + } + }; + + /** + * Determine if boost::function can use the small-object + * optimization with the function object type F. + */ + template + struct function_allows_small_object_optimization + { + BOOST_STATIC_CONSTANT + (bool, + value = ((sizeof(F) <= sizeof(function_buffer) && + (alignment_of::value + % alignment_of::value == 0)))); + }; + + template + struct functor_wrapper: public F, public A + { + functor_wrapper( F f, A a ): + F(f), + A(a) + { + } + + functor_wrapper(const functor_wrapper& f) : + F(static_cast(f)), + A(static_cast(f)) + { + } + }; + + /** + * The functor_manager class contains a static function "manage" which + * can clone or destroy the given function/function object pointer. + */ + template + struct functor_manager_common + { + typedef Functor functor_type; + + // Function pointers + static inline void + manage_ptr(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op) + { + if (op == clone_functor_tag) + out_buffer.func_ptr = in_buffer.func_ptr; + else if (op == move_functor_tag) { + out_buffer.func_ptr = in_buffer.func_ptr; + in_buffer.func_ptr = 0; + } else if (op == destroy_functor_tag) + out_buffer.func_ptr = 0; + else if (op == check_functor_type_tag) { + if (out_buffer.type.type == type_id()) + out_buffer.obj_ptr = &in_buffer.func_ptr; + else + out_buffer.obj_ptr = 0; + } else /* op == get_functor_type_tag */ { + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + } + } + + // Function objects that fit in the small-object buffer. + static inline void + manage_small(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op) + { + if (op == clone_functor_tag || op == move_functor_tag) { + const functor_type* in_functor = + reinterpret_cast(&in_buffer.data); + new (reinterpret_cast(&out_buffer.data)) functor_type(*in_functor); + + if (op == move_functor_tag) { + functor_type* f = reinterpret_cast(&in_buffer.data); + (void)f; // suppress warning about the value of f not being used (MSVC) + f->~Functor(); + } + } else if (op == destroy_functor_tag) { + // Some compilers (Borland, vc6, ...) are unhappy with ~functor_type. + functor_type* f = reinterpret_cast(&out_buffer.data); + (void)f; // suppress warning about the value of f not being used (MSVC) + f->~Functor(); + } else if (op == check_functor_type_tag) { + if (out_buffer.type.type == type_id()) + out_buffer.obj_ptr = &in_buffer.data; + else + out_buffer.obj_ptr = 0; + } else /* op == get_functor_type_tag */ { + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + } + } + }; + + template + struct functor_manager + { + private: + typedef Functor functor_type; + + // Function pointers + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, function_ptr_tag) + { + functor_manager_common::manage_ptr(in_buffer,out_buffer,op); + } + + // Function objects that fit in the small-object buffer. + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, mpl::true_) + { + functor_manager_common::manage_small(in_buffer,out_buffer,op); + } + + // Function objects that require heap allocation + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, mpl::false_) + { + if (op == clone_functor_tag) { + // Clone the functor + // GCC 2.95.3 gets the CV qualifiers wrong here, so we + // can't do the static_cast that we should do. + // jewillco: Changing this to static_cast because GCC 2.95.3 is + // obsolete. + const functor_type* f = + static_cast(in_buffer.obj_ptr); + functor_type* new_f = new functor_type(*f); + out_buffer.obj_ptr = new_f; + } else if (op == move_functor_tag) { + out_buffer.obj_ptr = in_buffer.obj_ptr; + in_buffer.obj_ptr = 0; + } else if (op == destroy_functor_tag) { + /* Cast from the void pointer to the functor pointer type */ + functor_type* f = + static_cast(out_buffer.obj_ptr); + delete f; + out_buffer.obj_ptr = 0; + } else if (op == check_functor_type_tag) { + if (out_buffer.type.type == type_id()) + out_buffer.obj_ptr = in_buffer.obj_ptr; + else + out_buffer.obj_ptr = 0; + } else /* op == get_functor_type_tag */ { + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + } + } + + // For function objects, we determine whether the function + // object can use the small-object optimization buffer or + // whether we need to allocate it on the heap. + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, function_obj_tag) + { + manager(in_buffer, out_buffer, op, + mpl::bool_<(function_allows_small_object_optimization::value)>()); + } + + // For member pointers, we use the small-object optimization buffer. + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, member_ptr_tag) + { + manager(in_buffer, out_buffer, op, mpl::true_()); + } + + public: + /* Dispatch to an appropriate manager based on whether we have a + function pointer or a function object pointer. */ + static inline void + manage(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op) + { + typedef typename get_function_tag::type tag_type; + switch (op) { + case get_functor_type_tag: + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + return; + + default: + manager(in_buffer, out_buffer, op, tag_type()); + return; + } + } + }; + + template + struct functor_manager_a + { + private: + typedef Functor functor_type; + + // Function pointers + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, function_ptr_tag) + { + functor_manager_common::manage_ptr(in_buffer,out_buffer,op); + } + + // Function objects that fit in the small-object buffer. + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, mpl::true_) + { + functor_manager_common::manage_small(in_buffer,out_buffer,op); + } + + // Function objects that require heap allocation + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, mpl::false_) + { + typedef functor_wrapper functor_wrapper_type; + typedef typename Allocator::template rebind::other + wrapper_allocator_type; + typedef typename wrapper_allocator_type::pointer wrapper_allocator_pointer_type; + + if (op == clone_functor_tag) { + // Clone the functor + // GCC 2.95.3 gets the CV qualifiers wrong here, so we + // can't do the static_cast that we should do. + const functor_wrapper_type* f = + static_cast(in_buffer.obj_ptr); + wrapper_allocator_type wrapper_allocator(static_cast(*f)); + wrapper_allocator_pointer_type copy = wrapper_allocator.allocate(1); + wrapper_allocator.construct(copy, *f); + + // Get back to the original pointer type + functor_wrapper_type* new_f = static_cast(copy); + out_buffer.obj_ptr = new_f; + } else if (op == move_functor_tag) { + out_buffer.obj_ptr = in_buffer.obj_ptr; + in_buffer.obj_ptr = 0; + } else if (op == destroy_functor_tag) { + /* Cast from the void pointer to the functor_wrapper_type */ + functor_wrapper_type* victim = + static_cast(in_buffer.obj_ptr); + wrapper_allocator_type wrapper_allocator(static_cast(*victim)); + wrapper_allocator.destroy(victim); + wrapper_allocator.deallocate(victim,1); + out_buffer.obj_ptr = 0; + } else if (op == check_functor_type_tag) { + if (out_buffer.type.type == type_id()) + out_buffer.obj_ptr = in_buffer.obj_ptr; + else + out_buffer.obj_ptr = 0; + } else /* op == get_functor_type_tag */ { + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + } + } + + // For function objects, we determine whether the function + // object can use the small-object optimization buffer or + // whether we need to allocate it on the heap. + static inline void + manager(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op, function_obj_tag) + { + manager(in_buffer, out_buffer, op, + mpl::bool_<(function_allows_small_object_optimization::value)>()); + } + + public: + /* Dispatch to an appropriate manager based on whether we have a + function pointer or a function object pointer. */ + static inline void + manage(const function_buffer& in_buffer, function_buffer& out_buffer, + functor_manager_operation_type op) + { + typedef typename get_function_tag::type tag_type; + switch (op) { + case get_functor_type_tag: + out_buffer.type.type = type_id(); + out_buffer.type.const_qualified = false; + out_buffer.type.volatile_qualified = false; + return; + + default: + manager(in_buffer, out_buffer, op, tag_type()); + return; + } + } + }; + + // A type that is only used for comparisons against zero + struct useless_clear_type {}; + +#ifdef BOOST_NO_SFINAE + // These routines perform comparisons between a Boost.Function + // object and an arbitrary function object (when the last + // parameter is mpl::bool_) or against zero (when the + // last parameter is mpl::bool_). They are only necessary + // for compilers that don't support SFINAE. + template + bool + compare_equal(const Function& f, const Functor&, int, mpl::bool_) + { return f.empty(); } + + template + bool + compare_not_equal(const Function& f, const Functor&, int, + mpl::bool_) + { return !f.empty(); } + + template + bool + compare_equal(const Function& f, const Functor& g, long, + mpl::bool_) + { + if (const Functor* fp = f.template target()) + return function_equal(*fp, g); + else return false; + } + + template + bool + compare_equal(const Function& f, const reference_wrapper& g, + int, mpl::bool_) + { + if (const Functor* fp = f.template target()) + return fp == g.get_pointer(); + else return false; + } + + template + bool + compare_not_equal(const Function& f, const Functor& g, long, + mpl::bool_) + { + if (const Functor* fp = f.template target()) + return !function_equal(*fp, g); + else return true; + } + + template + bool + compare_not_equal(const Function& f, + const reference_wrapper& g, int, + mpl::bool_) + { + if (const Functor* fp = f.template target()) + return fp != g.get_pointer(); + else return true; + } +#endif // BOOST_NO_SFINAE + + /** + * Stores the "manager" portion of the vtable for a + * boost::function object. + */ + struct vtable_base + { + void (*manager)(const function_buffer& in_buffer, + function_buffer& out_buffer, + functor_manager_operation_type op); + }; + } // end namespace function + } // end namespace detail + +/** + * The function_base class contains the basic elements needed for the + * function1, function2, function3, etc. classes. It is common to all + * functions (and as such can be used to tell if we have one of the + * functionN objects). + */ +class function_base +{ +public: + function_base() : vtable(0) { } + + /** Determine if the function is empty (i.e., has no target). */ + bool empty() const { return !vtable; } + + /** Retrieve the type of the stored function object, or type_id() + if this is empty. */ + type_index target_type() const + { + if (!vtable) return type_id(); + + detail::function::function_buffer type; + get_vtable()->manager(functor, type, detail::function::get_functor_type_tag); + return *type.type.type; + } + + template + Functor* target() + { + if (!vtable) return 0; + + detail::function::function_buffer type_result; + type_result.type.type = type_id(); + type_result.type.const_qualified = is_const::value; + type_result.type.volatile_qualified = is_volatile::value; + get_vtable()->manager(functor, type_result, + detail::function::check_functor_type_tag); + return static_cast(type_result.obj_ptr); + } + + template +#if defined(BOOST_MSVC) && BOOST_WORKAROUND(BOOST_MSVC, < 1300) + const Functor* target( Functor * = 0 ) const +#else + const Functor* target() const +#endif + { + if (!vtable) return 0; + + detail::function::function_buffer type_result; + type_result.type.type = type_id(); + type_result.type.const_qualified = true; + type_result.type.volatile_qualified = is_volatile::value; + get_vtable()->manager(functor, type_result, + detail::function::check_functor_type_tag); + // GCC 2.95.3 gets the CV qualifiers wrong here, so we + // can't do the static_cast that we should do. + return static_cast(type_result.obj_ptr); + } + + template + bool contains(const F& f) const + { +#if defined(BOOST_MSVC) && BOOST_WORKAROUND(BOOST_MSVC, < 1300) + if (const F* fp = this->target( (F*)0 )) +#else + if (const F* fp = this->template target()) +#endif + { + return function_equal(*fp, f); + } else { + return false; + } + } + +#if defined(__GNUC__) && __GNUC__ == 3 && __GNUC_MINOR__ <= 3 + // GCC 3.3 and newer cannot copy with the global operator==, due to + // problems with instantiation of function return types before it + // has been verified that the argument types match up. + template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator==(Functor g) const + { + if (const Functor* fp = target()) + return function_equal(*fp, g); + else return false; + } + + template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator!=(Functor g) const + { + if (const Functor* fp = target()) + return !function_equal(*fp, g); + else return true; + } +#endif + +public: // should be protected, but GCC 2.95.3 will fail to allow access + detail::function::vtable_base* get_vtable() const { + return reinterpret_cast( + reinterpret_cast(vtable) & ~static_cast(0x01)); + } + + bool has_trivial_copy_and_destroy() const { + return reinterpret_cast(vtable) & 0x01; + } + + detail::function::vtable_base* vtable; + mutable detail::function::function_buffer functor; +}; + +/** + * The bad_function_call exception class is thrown when a boost::function + * object is invoked + */ +class bad_function_call : public std::runtime_error +{ +public: + bad_function_call() : std::runtime_error("call to empty boost::function") {} +}; + +#ifndef BOOST_NO_SFINAE +inline bool operator==(const function_base& f, + detail::function::useless_clear_type*) +{ + return f.empty(); +} + +inline bool operator!=(const function_base& f, + detail::function::useless_clear_type*) +{ + return !f.empty(); +} + +inline bool operator==(detail::function::useless_clear_type*, + const function_base& f) +{ + return f.empty(); +} + +inline bool operator!=(detail::function::useless_clear_type*, + const function_base& f) +{ + return !f.empty(); +} +#endif + +#ifdef BOOST_NO_SFINAE +// Comparisons between boost::function objects and arbitrary function objects +template + inline bool operator==(const function_base& f, Functor g) + { + typedef mpl::bool_<(is_integral::value)> integral; + return detail::function::compare_equal(f, g, 0, integral()); + } + +template + inline bool operator==(Functor g, const function_base& f) + { + typedef mpl::bool_<(is_integral::value)> integral; + return detail::function::compare_equal(f, g, 0, integral()); + } + +template + inline bool operator!=(const function_base& f, Functor g) + { + typedef mpl::bool_<(is_integral::value)> integral; + return detail::function::compare_not_equal(f, g, 0, integral()); + } + +template + inline bool operator!=(Functor g, const function_base& f) + { + typedef mpl::bool_<(is_integral::value)> integral; + return detail::function::compare_not_equal(f, g, 0, integral()); + } +#else + +# if !(defined(__GNUC__) && __GNUC__ == 3 && __GNUC_MINOR__ <= 3) +// Comparisons between boost::function objects and arbitrary function +// objects. GCC 3.3 and before has an obnoxious bug that prevents this +// from working. +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator==(const function_base& f, Functor g) + { + if (const Functor* fp = f.template target()) + return function_equal(*fp, g); + else return false; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator==(Functor g, const function_base& f) + { + if (const Functor* fp = f.template target()) + return function_equal(g, *fp); + else return false; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator!=(const function_base& f, Functor g) + { + if (const Functor* fp = f.template target()) + return !function_equal(*fp, g); + else return true; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator!=(Functor g, const function_base& f) + { + if (const Functor* fp = f.template target()) + return !function_equal(g, *fp); + else return true; + } +# endif + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator==(const function_base& f, reference_wrapper g) + { + if (const Functor* fp = f.template target()) + return fp == g.get_pointer(); + else return false; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator==(reference_wrapper g, const function_base& f) + { + if (const Functor* fp = f.template target()) + return g.get_pointer() == fp; + else return false; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator!=(const function_base& f, reference_wrapper g) + { + if (const Functor* fp = f.template target()) + return fp != g.get_pointer(); + else return true; + } + +template + BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL(Functor, bool) + operator!=(reference_wrapper g, const function_base& f) + { + if (const Functor* fp = f.template target()) + return g.get_pointer() != fp; + else return true; + } + +#endif // Compiler supporting SFINAE + +namespace detail { + namespace function { + inline bool has_empty_target(const function_base* f) + { + return f->empty(); + } + +#if BOOST_WORKAROUND(BOOST_MSVC, <= 1310) + inline bool has_empty_target(const void*) + { + return false; + } +#else + inline bool has_empty_target(...) + { + return false; + } +#endif + } // end namespace function +} // end namespace detail +} // end namespace boost + +#undef BOOST_FUNCTION_ENABLE_IF_NOT_INTEGRAL + +#if defined(BOOST_MSVC) +# pragma warning( pop ) +#endif + +#endif // BOOST_FUNCTION_BASE_HEADER diff --git a/patched/graph/graphml.hpp b/patched/graph/graphml.hpp new file mode 100644 index 0000000..05e8d60 --- /dev/null +++ b/patched/graph/graphml.hpp @@ -0,0 +1,352 @@ +// Copyright (C) 2006 Tiago de Paula Peixoto +// Copyright (C) 2004 The Trustees of Indiana University. +// +// Use, modification and distribution is subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// +// Authors: Douglas Gregor +// Andrew Lumsdaine +// Tiago de Paula Peixoto + +#ifndef BOOST_GRAPH_GRAPHML_HPP +#define BOOST_GRAPH_GRAPHML_HPP + +#include +#include +#include +#include +#include +#include // for exceptions +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost +{ + +///////////////////////////////////////////////////////////////////////////// +// Graph reader exceptions +///////////////////////////////////////////////////////////////////////////// +struct parse_error: public graph_exception +{ + parse_error(const std::string& err) {error = err; statement = "parse error: " + error;} + virtual ~parse_error() throw() {} + virtual const char* what() const throw() {return statement.c_str();} + std::string statement; + std::string error; +}; + + +class mutate_graph +{ +public: + virtual ~mutate_graph() {} + virtual bool is_directed() const = 0; + + virtual boost::any do_add_vertex() = 0; + virtual std::pair do_add_edge(boost::any source, boost::any target) = 0; + + virtual void + set_graph_property(const std::string& name, const std::string& value, const std::string& value_type) = 0; + + virtual void + set_vertex_property(const std::string& name, boost::any vertex, const std::string& value, const std::string& value_type) = 0; + + virtual void + set_edge_property(const std::string& name, boost::any edge, const std::string& value, const std::string& value_type) = 0; +}; + +template +class mutate_graph_impl : public mutate_graph +{ + typedef typename graph_traits::vertex_descriptor vertex_descriptor; + typedef typename graph_traits::edge_descriptor edge_descriptor; + + public: + mutate_graph_impl(MutableGraph& g, dynamic_properties& dp) + : m_g(g), m_dp(dp) { } + + bool is_directed() const + { + return is_convertible::directed_category, + directed_tag>::value; + } + + virtual any do_add_vertex() + { + return any(add_vertex(m_g)); + } + + virtual std::pair do_add_edge(any source, any target) + { + std::pair retval = add_edge(any_cast(source), + any_cast(target), m_g); + return std::make_pair(any(retval.first), retval.second); + } + + virtual void + set_graph_property(const std::string& name, const std::string& value, const std::string& value_type) + { + bool type_found = false; + try + { + mpl::for_each(put_property + (name, m_dp, m_g, value, value_type, m_type_names, type_found)); + } + catch (bad_lexical_cast) + { + BOOST_THROW_EXCEPTION( + parse_error("invalid value \"" + value + "\" for key " + + name + " of type " + value_type)); + } + if (!type_found) + { + BOOST_THROW_EXCEPTION( + parse_error("unrecognized type \"" + value_type + + "\" for key " + name)); + } + + } + + virtual void + set_vertex_property(const std::string& name, any vertex, const std::string& value, const std::string& value_type) + { + bool type_found = false; + try + { + mpl::for_each(put_property + (name, m_dp, any_cast(vertex), + value, value_type, m_type_names, type_found)); + } + catch (bad_lexical_cast) + { + BOOST_THROW_EXCEPTION( + parse_error("invalid value \"" + value + "\" for key " + + name + " of type " + value_type)); + } + if (!type_found) + { + BOOST_THROW_EXCEPTION( + parse_error("unrecognized type \"" + value_type + + "\" for key " + name)); + } + + } + + virtual void + set_edge_property(const std::string& name, any edge, const std::string& value, const std::string& value_type) + { + bool type_found = false; + try + { + mpl::for_each(put_property + (name, m_dp, any_cast(edge), + value, value_type, m_type_names, type_found)); + } + catch (bad_lexical_cast) + { + BOOST_THROW_EXCEPTION( + parse_error("invalid value \"" + value + "\" for key " + + name + " of type " + value_type)); + } + if (!type_found) + { + BOOST_THROW_EXCEPTION( + parse_error("unrecognized type \"" + value_type + + "\" for key " + name)); + } + } + + template + class put_property + { + public: + put_property(const std::string& name, dynamic_properties& dp, const Key& key, + const std::string& value, const std::string& value_type, + const char** type_names, bool& type_found) + : m_name(name), m_dp(dp), m_key(key), m_value(value), + m_value_type(value_type), m_type_names(type_names), + m_type_found(type_found) {} + template + void operator()(Value) + { + if (m_value_type == m_type_names[mpl::find::type::pos::value]) + { + put(m_name, m_dp, m_key, lexical_cast(m_value)); + m_type_found = true; + } + } + private: + const std::string& m_name; + dynamic_properties& m_dp; + const Key& m_key; + const std::string& m_value; + const std::string& m_value_type; + const char** m_type_names; + bool& m_type_found; + }; + +protected: + MutableGraph& m_g; + dynamic_properties& m_dp; + typedef mpl::vector value_types; + static const char* m_type_names[]; +}; + +template +const char* mutate_graph_impl::m_type_names[] = {"boolean", "int", "long", "float", "double", "string"}; + +void BOOST_GRAPH_DECL +read_graphml(std::istream& in, mutate_graph& g); + +template +void +read_graphml(std::istream& in, MutableGraph& g, dynamic_properties& dp) +{ + mutate_graph_impl mg(g,dp); + read_graphml(in, mg); +} + +template +class get_type_name +{ +public: + get_type_name(const type_index& type, const char** type_names, std::string& type_name) + : m_type(type), m_type_names(type_names), m_type_name(type_name) {} + template + void operator()(Type) + { + if (type_id() == m_type) + m_type_name = m_type_names[mpl::find::type::pos::value]; + } +private: + const type_index m_type; + const char** m_type_names; + std::string &m_type_name; +}; + + +template +void +write_graphml(std::ostream& out, const Graph& g, VertexIndexMap vertex_index, + const dynamic_properties& dp, bool ordered_vertices=false) +{ + typedef typename graph_traits::directed_category directed_category; + typedef typename graph_traits::edge_descriptor edge_descriptor; + typedef typename graph_traits::vertex_descriptor vertex_descriptor; + + using boost::property_tree::xml_parser::encode_char_entities; + + BOOST_STATIC_CONSTANT(bool, + graph_is_directed = + (is_convertible::value)); + + out << "\n" + << "\n"; + + typedef mpl::vector value_types; + const char* type_names[] = {"boolean", "int", "int", "int", "int", "long", "long", "long", "long", "float", "double", "double", "string"}; + std::map graph_key_ids; + std::map vertex_key_ids; + std::map edge_key_ids; + int key_count = 0; + + // Output keys + for (dynamic_properties::const_iterator i = dp.begin(); i != dp.end(); ++i) + { + std::string key_id = "key" + lexical_cast(key_count++); + if (i->second->key() == type_id()) + graph_key_ids[i->first] = key_id; + else if (i->second->key() == type_id()) + vertex_key_ids[i->first] = key_id; + else if (i->second->key() == type_id()) + edge_key_ids[i->first] = key_id; + else + continue; + std::string type_name = "string"; + mpl::for_each(get_type_name(i->second->value(), type_names, type_name)); + out << " second->key() == type_id() ? "graph" : (i->second->key() == type_id() ? "node" : "edge")) << "\"" + << " attr.name=\"" << i->first << "\"" + << " attr.type=\"" << type_name << "\"" + << " />\n"; + } + + out << " \n"; + + // Output graph data + for (dynamic_properties::const_iterator i = dp.begin(); i != dp.end(); ++i) + { + if (i->second->key() == type_id()) + { + // The const_cast here is just to get typeid correct for property + // map key; the graph should not be mutated using it. + out << " first] << "\">" + << encode_char_entities(i->second->get_string(const_cast(&g))) << "\n"; + } + } + + typedef typename graph_traits::vertex_iterator vertex_iterator; + vertex_iterator v, v_end; + for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) + { + out << " \n"; + // Output data + for (dynamic_properties::const_iterator i = dp.begin(); i != dp.end(); ++i) + { + if (i->second->key() == type_id()) + { + out << " first] << "\">" + << encode_char_entities(i->second->get_string(*v)) << "\n"; + } + } + out << " \n"; + } + + typedef typename graph_traits::edge_iterator edge_iterator; + edge_iterator e, e_end; + typename graph_traits::edges_size_type edge_count = 0; + for (boost::tie(e, e_end) = edges(g); e != e_end; ++e) + { + out << " \n"; + + // Output data + for (dynamic_properties::const_iterator i = dp.begin(); i != dp.end(); ++i) + { + if (i->second->key() == type_id()) + { + out << " first] << "\">" + << encode_char_entities(i->second->get_string(*e)) << "\n"; + } + } + out << " \n"; + } + + out << " \n" + << "\n"; +} + + +template +void +write_graphml(std::ostream& out, const Graph& g, const dynamic_properties& dp, + bool ordered_vertices=false) +{ + write_graphml(out, g, get(vertex_index, g), dp, ordered_vertices); +} + +} // boost namespace + +#endif // BOOST_GRAPH_GRAPHML_HPP diff --git a/patched/graph/graphviz.hpp b/patched/graph/graphviz.hpp new file mode 100644 index 0000000..aff74a9 --- /dev/null +++ b/patched/graph/graphviz.hpp @@ -0,0 +1,858 @@ +//======================================================================= +// Copyright 2001 University of Notre Dame. +// Copyright 2003 Jeremy Siek +// Authors: Lie-Quan Lee, Jeremy Siek, and Douglas Gregor +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= +#ifndef BOOST_GRAPHVIZ_HPP +#define BOOST_GRAPHVIZ_HPP + +#include +#include +#include +#include +#include +#include // for FILE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + + template + struct graphviz_io_traits { + static std::string name() { + return "digraph"; + } + static std::string delimiter() { + return "->"; + } }; + + template <> + struct graphviz_io_traits { + static std::string name() { + return "graph"; + } + static std::string delimiter() { + return "--"; + } + }; + + struct default_writer { + void operator()(std::ostream&) const { + } + template + void operator()(std::ostream&, const VorE&) const { + } + }; + + template + inline std::string escape_dot_string(const T& obj) { + using namespace boost::xpressive; + static sregex valid_unquoted_id = (((alpha | '_') >> *_w) | (!as_xpr('-') >> (('.' >> *_d) | (+_d >> !('.' >> *_d))))); + std::string s(boost::lexical_cast(obj)); + if (regex_match(s, valid_unquoted_id)) { + return s; + } else { + boost::algorithm::replace_all(s, "\"", "\\\""); + return "\"" + s + "\""; + } + } + + template + class label_writer { + public: + label_writer(Name _name) : name(_name) {} + template + void operator()(std::ostream& out, const VertexOrEdge& v) const { + out << "[label=" << escape_dot_string(get(name, v)) << "]"; + } + private: + Name name; + }; + template + inline label_writer + make_label_writer(Name n) { + return label_writer(n); + } + + enum edge_attribute_t { edge_attribute = 1111 }; + enum vertex_attribute_t { vertex_attribute = 2222 }; + enum graph_graph_attribute_t { graph_graph_attribute = 3333 }; + enum graph_vertex_attribute_t { graph_vertex_attribute = 4444 }; + enum graph_edge_attribute_t { graph_edge_attribute = 5555 }; + + BOOST_INSTALL_PROPERTY(edge, attribute); + BOOST_INSTALL_PROPERTY(vertex, attribute); + BOOST_INSTALL_PROPERTY(graph, graph_attribute); + BOOST_INSTALL_PROPERTY(graph, vertex_attribute); + BOOST_INSTALL_PROPERTY(graph, edge_attribute); + + + template + inline void write_attributes(const Attribute& attr, std::ostream& out) { + typename Attribute::const_iterator i, iend; + i = attr.begin(); + iend = attr.end(); + + while ( i != iend ) { + out << i->first << "=" << escape_dot_string(i->second); + ++i; + if ( i != iend ) + out << ", "; + } + } + + template + inline void write_all_attributes(Attributes attributes, + const std::string& name, + std::ostream& out) + { + typename Attributes::const_iterator i = attributes.begin(), + end = attributes.end(); + if (i != end) { + out << name << " [\n"; + write_attributes(attributes, out); + out << "];\n"; + } + } + + inline void write_all_attributes(detail::error_property_not_found, + const std::string&, + std::ostream&) + { + // Do nothing - no attributes exist + } + + + + + template + struct graph_attributes_writer + { + graph_attributes_writer(GraphGraphAttributes gg, + GraphNodeAttributes gn, + GraphEdgeAttributes ge) + : g_attributes(gg), n_attributes(gn), e_attributes(ge) { } + + void operator()(std::ostream& out) const { + write_all_attributes(g_attributes, "graph", out); + write_all_attributes(n_attributes, "node", out); + write_all_attributes(e_attributes, "edge", out); + } + GraphGraphAttributes g_attributes; + GraphNodeAttributes n_attributes; + GraphEdgeAttributes e_attributes; + }; + + template + graph_attributes_writer + make_graph_attributes_writer(const GAttrMap& g_attr, const NAttrMap& n_attr, + const EAttrMap& e_attr) { + return graph_attributes_writer + (g_attr, n_attr, e_attr); + } + + + template + graph_attributes_writer + ::type, + typename graph_property::type, + typename graph_property::type> + make_graph_attributes_writer(const Graph& g) + { + typedef typename graph_property::type + GAttrMap; + typedef typename graph_property::type + NAttrMap; + typedef typename graph_property::type + EAttrMap; + GAttrMap gam = get_property(g, graph_graph_attribute); + NAttrMap nam = get_property(g, graph_vertex_attribute); + EAttrMap eam = get_property(g, graph_edge_attribute); + graph_attributes_writer writer(gam, nam, eam); + return writer; + } + + template + struct attributes_writer { + attributes_writer(AttributeMap attr) + : attributes(attr) { } + + template + void operator()(std::ostream& out, const VorE& e) const { + this->write_attribute(out, attributes[e]); + } + + private: + template + void write_attribute(std::ostream& out, + const AttributeSequence& seq) const + { + if (!seq.empty()) { + out << "["; + write_attributes(seq, out); + out << "]"; + } + } + + void write_attribute(std::ostream&, + detail::error_property_not_found) const + { + } + AttributeMap attributes; + }; + + template + attributes_writer + ::const_type> + make_edge_attributes_writer(const Graph& g) + { + typedef typename property_map::const_type + EdgeAttributeMap; + return attributes_writer(get(edge_attribute, g)); + } + + template + attributes_writer + ::const_type> + make_vertex_attributes_writer(const Graph& g) + { + typedef typename property_map::const_type + VertexAttributeMap; + return attributes_writer(get(vertex_attribute, g)); + } + + template + inline void + write_graphviz + (std::ostream& out, const Graph& g, + VertexPropertiesWriter vpw, + EdgePropertiesWriter epw, + GraphPropertiesWriter gpw, + VertexID vertex_id + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + BOOST_CONCEPT_ASSERT((EdgeListGraphConcept)); + + typedef typename graph_traits::directed_category cat_type; + typedef graphviz_io_traits Traits; + std::string name = "G"; + out << Traits::name() << " " << escape_dot_string(name) << " {" << std::endl; + + gpw(out); //print graph properties + + typename graph_traits::vertex_iterator i, end; + + for(boost::tie(i,end) = vertices(g); i != end; ++i) { + out << escape_dot_string(get(vertex_id, *i)); + vpw(out, *i); //print vertex attributes + out << ";" << std::endl; + } + typename graph_traits::edge_iterator ei, edge_end; + for(boost::tie(ei, edge_end) = edges(g); ei != edge_end; ++ei) { + out << escape_dot_string(get(vertex_id, source(*ei, g))) << Traits::delimiter() << escape_dot_string(get(vertex_id, target(*ei, g))) << " "; + epw(out, *ei); //print edge attributes + out << ";" << std::endl; + } + out << "}" << std::endl; + } + + template + inline void + write_graphviz(std::ostream& out, const Graph& g, + VertexPropertiesWriter vpw, + EdgePropertiesWriter epw, + GraphPropertiesWriter gpw + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { write_graphviz(out, g, vpw, epw, gpw, get(vertex_index, g)); } + +#if !defined(BOOST_MSVC) || BOOST_MSVC > 1300 + // ambiguous overload problem with VC++ + template + inline void + write_graphviz(std::ostream& out, const Graph& g + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + default_writer dw; + default_writer gw; + write_graphviz(out, g, dw, dw, gw); + } +#endif + + template + inline void + write_graphviz(std::ostream& out, const Graph& g, VertexWriter vw + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + default_writer dw; + default_writer gw; + write_graphviz(out, g, vw, dw, gw); + } + + template + inline void + write_graphviz(std::ostream& out, const Graph& g, + VertexWriter vw, EdgeWriter ew + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + default_writer gw; + write_graphviz(out, g, vw, ew, gw); + } + + namespace detail { + + template + void write_graphviz_subgraph (std::ostream& out, + const subgraph& g, + RandomAccessIterator vertex_marker, + RandomAccessIterator edge_marker, + VertexID vertex_id) + { + typedef subgraph Graph; + typedef typename graph_traits::vertex_descriptor Vertex; + typedef typename graph_traits::directed_category cat_type; + typedef graphviz_io_traits Traits; + + typedef typename graph_property::type NameType; + const NameType& g_name = get_property(g, graph_name); + + if ( g.is_root() ) + out << Traits::name() ; + else + out << "subgraph"; + + out << " " << escape_dot_string(g_name) << " {" << std::endl; + + typename Graph::const_children_iterator i_child, j_child; + + //print graph/node/edge attributes +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 + typedef typename graph_property::type + GAttrMap; + typedef typename graph_property::type + NAttrMap; + typedef typename graph_property::type + EAttrMap; + GAttrMap gam = get_property(g, graph_graph_attribute); + NAttrMap nam = get_property(g, graph_vertex_attribute); + EAttrMap eam = get_property(g, graph_edge_attribute); + graph_attributes_writer writer(gam, nam, eam); + writer(out); +#else + make_graph_attributes_writer(g)(out); +#endif + + //print subgraph + for ( boost::tie(i_child,j_child) = g.children(); + i_child != j_child; ++i_child ) + write_graphviz_subgraph(out, *i_child, vertex_marker, edge_marker, + vertex_id); + + // Print out vertices and edges not in the subgraphs. + + typename graph_traits::vertex_iterator i, end; + typename graph_traits::edge_iterator ei, edge_end; + + for(boost::tie(i,end) = vertices(g); i != end; ++i) { + Vertex v = g.local_to_global(*i); + int pos = get(vertex_id, v); + if ( vertex_marker[pos] ) { + vertex_marker[pos] = false; + out << escape_dot_string(pos); +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 + typedef typename property_map::const_type + VertexAttributeMap; + attributes_writer vawriter(get(vertex_attribute, + g.root())); + vawriter(out, v); +#else + make_vertex_attributes_writer(g.root())(out, v); +#endif + out << ";" << std::endl; + } + } + + for (boost::tie(ei, edge_end) = edges(g); ei != edge_end; ++ei) { + Vertex u = g.local_to_global(source(*ei,g)), + v = g.local_to_global(target(*ei, g)); + int pos = get(get(edge_index, g.root()), g.local_to_global(*ei)); + if ( edge_marker[pos] ) { + edge_marker[pos] = false; + out << escape_dot_string(get(vertex_id, u)) << " " << Traits::delimiter() + << " " << escape_dot_string(get(vertex_id, v)); +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 + typedef typename property_map::const_type + EdgeAttributeMap; + attributes_writer eawriter(get(edge_attribute, g)); + eawriter(out, *ei); +#else + make_edge_attributes_writer(g)(out, *ei); //print edge properties +#endif + out << ";" << std::endl; + } + } + out << "}" << std::endl; + } + } // namespace detail + + // requires graph_name graph property + template + void write_graphviz(std::ostream& out, const subgraph& g) { + std::vector edge_marker(num_edges(g), true); + std::vector vertex_marker(num_vertices(g), true); + + detail::write_graphviz_subgraph(out, g, + vertex_marker.begin(), + edge_marker.begin(), + get(vertex_index, g)); + } + + template + void write_graphviz(const std::string& filename, const subgraph& g) { + std::ofstream out(filename.c_str()); + std::vector edge_marker(num_edges(g), true); + std::vector vertex_marker(num_vertices(g), true); + + detail::write_graphviz_subgraph(out, g, + vertex_marker.begin(), + edge_marker.begin(), + get(vertex_index, g)); + } + + template + void write_graphviz(std::ostream& out, const subgraph& g, + VertexID vertex_id) + { + std::vector edge_marker(num_edges(g), true); + std::vector vertex_marker(num_vertices(g), true); + + detail::write_graphviz_subgraph(out, g, + vertex_marker.begin(), + edge_marker.begin(), + vertex_id); + } + + template + void write_graphviz(const std::string& filename, const subgraph& g, + VertexID vertex_id) + { + std::ofstream out(filename.c_str()); + std::vector edge_marker(num_edges(g), true); + std::vector vertex_marker(num_vertices(g), true); + + detail::write_graphviz_subgraph(out, g, + vertex_marker.begin(), + edge_marker.begin(), + vertex_id); + } + +#if 0 + // This interface has not worked for a long time + typedef std::map GraphvizAttrList; + + typedef property + GraphvizVertexProperty; + + typedef property > + GraphvizEdgeProperty; + + typedef property > > > + GraphvizGraphProperty; + + typedef subgraph > + GraphvizDigraph; + + typedef subgraph > + GraphvizGraph; + + // These four require linking the BGL-Graphviz library: libbgl-viz.a + // from the /src directory. + // Library has not existed for a while + extern void read_graphviz(const std::string& file, GraphvizDigraph& g); + extern void read_graphviz(FILE* file, GraphvizDigraph& g); + + extern void read_graphviz(const std::string& file, GraphvizGraph& g); + extern void read_graphviz(FILE* file, GraphvizGraph& g); +#endif + + class dynamic_properties_writer + { + public: + dynamic_properties_writer(const dynamic_properties& dp) : dp(&dp) { } + + template + void operator()(std::ostream& out, Descriptor key) const + { + bool first = true; + for (dynamic_properties::const_iterator i = dp->begin(); + i != dp->end(); ++i) { + if (type_id() == i->second->key()) { + if (first) out << " ["; + else out << ", "; + first = false; + + out << i->first << "=" << escape_dot_string(i->second->get_string(key)); + } + } + + if (!first) out << "]"; + } + + private: + const dynamic_properties* dp; + }; + + class dynamic_vertex_properties_writer + { + public: + dynamic_vertex_properties_writer(const dynamic_properties& dp, + const std::string& node_id) + : dp(&dp), node_id(&node_id) { } + + template + void operator()(std::ostream& out, Descriptor key) const + { + bool first = true; + for (dynamic_properties::const_iterator i = dp->begin(); + i != dp->end(); ++i) { + if (type_id(key) == i->second->key() + && i->first != *node_id) { + if (first) out << " ["; + else out << ", "; + first = false; + + out << i->first << "=" << escape_dot_string(i->second->get_string(key)); + } + } + + if (!first) out << "]"; + } + + private: + const dynamic_properties* dp; + const std::string* node_id; + }; + + namespace graph { namespace detail { + + template + struct node_id_property_map + { + typedef std::string value_type; + typedef value_type reference; + typedef Vertex key_type; + typedef readable_property_map_tag category; + + node_id_property_map() {} + + node_id_property_map(const dynamic_properties& dp, + const std::string& node_id) + : dp(&dp), node_id(&node_id) { } + + const dynamic_properties* dp; + const std::string* node_id; + }; + + template + inline std::string + get(node_id_property_map pm, + typename node_id_property_map::key_type v) + { return get(*pm.node_id, *pm.dp, v); } + + } } // end namespace graph::detail + + template + inline void + write_graphviz_dp(std::ostream& out, const Graph& g, + const dynamic_properties& dp, + const std::string& node_id = "node_id" + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + typedef typename graph_traits::vertex_descriptor Vertex; + write_graphviz_dp(out, g, dp, node_id, + graph::detail::node_id_property_map(dp, node_id)); + } + + template + void + write_graphviz_dp(std::ostream& out, const Graph& g, + const dynamic_properties& dp, const std::string& node_id, + VertexID id + BOOST_GRAPH_ENABLE_IF_MODELS_PARM(Graph,vertex_list_graph_tag)) + { + write_graphviz + (out, g, + /*vertex_writer=*/dynamic_vertex_properties_writer(dp, node_id), + /*edge_writer=*/dynamic_properties_writer(dp), + /*graph_writer=*/default_writer(), + id); + } + +///////////////////////////////////////////////////////////////////////////// +// Graph reader exceptions +///////////////////////////////////////////////////////////////////////////// +struct graph_exception : public std::exception { + virtual ~graph_exception() throw() {} + virtual const char* what() const throw() = 0; +}; + +struct bad_parallel_edge : public graph_exception { + std::string from; + std::string to; + mutable std::string statement; + bad_parallel_edge(const std::string& i, const std::string& j) : + from(i), to(j) {} + + virtual ~bad_parallel_edge() throw() {} + const char* what() const throw() { + if(statement.empty()) + statement = + std::string("Failed to add parallel edge: (") + + from + "," + to + ")\n"; + + return statement.c_str(); + } +}; + +struct directed_graph_error : public graph_exception { + virtual ~directed_graph_error() throw() {} + virtual const char* what() const throw() { + return + "read_graphviz: " + "Tried to read a directed graph into an undirected graph."; + } +}; + +struct undirected_graph_error : public graph_exception { + virtual ~undirected_graph_error() throw() {} + virtual const char* what() const throw() { + return + "read_graphviz: " + "Tried to read an undirected graph into a directed graph."; + } +}; + +struct bad_graphviz_syntax: public graph_exception { + std::string errmsg; + bad_graphviz_syntax(const std::string& errmsg) + : errmsg(errmsg) {} + const char* what() const throw () {return errmsg.c_str();} + ~bad_graphviz_syntax() throw () {}; +}; + +namespace detail { namespace graph { + +typedef std::string id_t; +typedef id_t node_t; + +// edges are not uniquely determined by adjacent nodes +class edge_t { + int idx_; + explicit edge_t(int i) : idx_(i) {} +public: + static edge_t new_edge() { + static int idx = 0; + return edge_t(idx++); + }; + + bool operator==(const edge_t& rhs) const { + return idx_ == rhs.idx_; + } + bool operator<(const edge_t& rhs) const { + return idx_ < rhs.idx_; + } +}; + +class mutate_graph +{ + public: + virtual ~mutate_graph() {} + virtual bool is_directed() const = 0; + virtual void do_add_vertex(const node_t& node) = 0; + + virtual void + do_add_edge(const edge_t& edge, const node_t& source, const node_t& target) + = 0; + + virtual void + set_node_property(const id_t& key, const node_t& node, const id_t& value) = 0; + + virtual void + set_edge_property(const id_t& key, const edge_t& edge, const id_t& value) = 0; + + virtual void // RG: need new second parameter to support BGL subgraphs + set_graph_property(const id_t& key, const id_t& value) = 0; +}; + +template +class mutate_graph_impl : public mutate_graph +{ + typedef typename graph_traits::vertex_descriptor bgl_vertex_t; + typedef typename graph_traits::edge_descriptor bgl_edge_t; + + public: + mutate_graph_impl(MutableGraph& graph, dynamic_properties& dp, + std::string node_id_prop) + : graph_(graph), dp_(dp), node_id_prop_(node_id_prop) { } + + ~mutate_graph_impl() {} + + bool is_directed() const + { + return + boost::is_convertible< + typename boost::graph_traits::directed_category, + boost::directed_tag>::value; + } + + virtual void do_add_vertex(const node_t& node) + { + // Add the node to the graph. + bgl_vertex_t v = add_vertex(graph_); + + // Set up a mapping from name to BGL vertex. + bgl_nodes.insert(std::make_pair(node, v)); + + // node_id_prop_ allows the caller to see the real id names for nodes. + put(node_id_prop_, dp_, v, node); + } + + void + do_add_edge(const edge_t& edge, const node_t& source, const node_t& target) + { + std::pair result = + add_edge(bgl_nodes[source], bgl_nodes[target], graph_); + + if(!result.second) { + // In the case of no parallel edges allowed + boost::throw_exception(bad_parallel_edge(source, target)); + } else { + bgl_edges.insert(std::make_pair(edge, result.first)); + } + } + + void + set_node_property(const id_t& key, const node_t& node, const id_t& value) + { + put(key, dp_, bgl_nodes[node], value); + } + + void + set_edge_property(const id_t& key, const edge_t& edge, const id_t& value) + { + put(key, dp_, bgl_edges[edge], value); + } + + void + set_graph_property(const id_t& key, const id_t& value) + { + /* RG: pointer to graph prevents copying */ + put(key, dp_, &graph_, value); + } + + + protected: + MutableGraph& graph_; + dynamic_properties& dp_; + std::string node_id_prop_; + std::map bgl_nodes; + std::map bgl_edges; +}; + +} } } // end namespace boost::detail::graph + +#ifdef BOOST_GRAPH_USE_SPIRIT_PARSER +# ifndef BOOST_GRAPH_READ_GRAPHVIZ_ITERATORS +# define BOOST_GRAPH_READ_GRAPHVIZ_ITERATORS +# endif +# include +#else // New default parser +# include +#endif // BOOST_GRAPH_USE_SPIRIT_PARSER + +namespace boost { + +// Parse the passed string as a GraphViz dot file. +template +bool read_graphviz(const std::string& data, + MutableGraph& graph, + dynamic_properties& dp, + std::string const& node_id = "node_id") { +#ifdef BOOST_GRAPH_USE_SPIRIT_PARSER + return read_graphviz_spirit(data.begin(), data.end(), graph, dp, node_id); +#else // Non-Spirit parser + return read_graphviz_new(data,graph,dp,node_id); +#endif +} + +// Parse the passed iterator range as a GraphViz dot file. +template +bool read_graphviz(InputIterator user_first, + InputIterator user_last, + MutableGraph& graph, + dynamic_properties& dp, + std::string const& node_id = "node_id") { +#ifdef BOOST_GRAPH_USE_SPIRIT_PARSER + typedef InputIterator is_t; + typedef boost::spirit::classic::multi_pass iterator_t; + + iterator_t first(boost::spirit::classic::make_multi_pass(user_first)); + iterator_t last(boost::spirit::classic::make_multi_pass(user_last)); + + return read_graphviz_spirit(first, last, graph, dp, node_id); +#else // Non-Spirit parser + return read_graphviz_new(std::string(user_first, user_last), graph, dp, node_id); +#endif +} + +// Parse the passed stream as a GraphViz dot file. +template +bool read_graphviz(std::istream& in, MutableGraph& graph, + dynamic_properties& dp, + std::string const& node_id = "node_id") +{ + typedef std::istream_iterator is_t; + in >> std::noskipws; + return read_graphviz(is_t(in), is_t(), graph, dp, node_id); +} + +} // namespace boost + +#ifdef BOOST_GRAPH_USE_MPI +# include +#endif + +#endif // BOOST_GRAPHVIZ_HPP diff --git a/patched/graph/parallel/distribution.hpp b/patched/graph/parallel/distribution.hpp new file mode 100644 index 0000000..1275af7 --- /dev/null +++ b/patched/graph/parallel/distribution.hpp @@ -0,0 +1,615 @@ +// Copyright 2004 The Trustees of Indiana University. + +// Use, modification and distribution is subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// Authors: Douglas Gregor +// Peter Gottschling +// Andrew Lumsdaine +#ifndef BOOST_PARALLEL_DISTRIBUTION_HPP +#define BOOST_PARALLEL_DISTRIBUTION_HPP + +#ifndef BOOST_GRAPH_USE_MPI +#error "Parallel BGL files should not be included unless has been included" +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { namespace parallel { + +template +class variant_distribution +{ +public: + typedef typename ProcessGroup::process_id_type process_id_type; + typedef typename ProcessGroup::process_size_type process_size_type; + typedef SizeType size_type; + +private: + struct basic_distribution + { + virtual ~basic_distribution() {} + virtual size_type block_size(process_id_type, size_type) const = 0; + virtual process_id_type in_process(size_type) const = 0; + virtual size_type local(size_type) const = 0; + virtual size_type global(size_type) const = 0; + virtual size_type global(process_id_type, size_type) const = 0; + virtual void* address() = 0; + virtual const void* address() const = 0; + virtual type_index type() const = 0; + }; + + template + struct poly_distribution : public basic_distribution + { + explicit poly_distribution(const Distribution& distribution) + : distribution_(distribution) { } + + virtual size_type block_size(process_id_type id, size_type n) const + { return distribution_.block_size(id, n); } + + virtual process_id_type in_process(size_type i) const + { return distribution_(i); } + + virtual size_type local(size_type i) const + { return distribution_.local(i); } + + virtual size_type global(size_type n) const + { return distribution_.global(n); } + + virtual size_type global(process_id_type id, size_type n) const + { return distribution_.global(id, n); } + + virtual void* address() { return &distribution_; } + virtual const void* address() const { return &distribution_; } + virtual type_index type() const { return type_id(); } + + private: + Distribution distribution_; + }; + +public: + variant_distribution() { } + + template + variant_distribution(const Distribution& distribution) + : distribution_(new poly_distribution(distribution)) { } + + size_type block_size(process_id_type id, size_type n) const + { return distribution_->block_size(id, n); } + + process_id_type operator()(size_type i) const + { return distribution_->in_process(i); } + + size_type local(size_type i) const + { return distribution_->local(i); } + + size_type global(size_type n) const + { return distribution_->global(n); } + + size_type global(process_id_type id, size_type n) const + { return distribution_->global(id, n); } + + operator bool() const { return distribution_; } + + void clear() { distribution_.reset(); } + + template + T* as() + { + if (distribution_->type() == type_id()) + return static_cast(distribution_->address()); + else + return 0; + } + + template + const T* as() const + { + if (distribution_->type() == type_id()) + return static_cast(distribution_->address()); + else + return 0; + } + +private: + shared_ptr distribution_; +}; + +struct block +{ + template + explicit block(const LinearProcessGroup& pg, std::size_t n) + : id(process_id(pg)), p(num_processes(pg)), n(n) { } + + // If there are n elements in the distributed data structure, returns the number of elements stored locally. + template + SizeType block_size(SizeType n) const + { return (n / p) + ((std::size_t)(n % p) > id? 1 : 0); } + + // If there are n elements in the distributed data structure, returns the number of elements stored on processor ID + template + SizeType block_size(ProcessID id, SizeType n) const + { return (n / p) + ((ProcessID)(n % p) > id? 1 : 0); } + + // Returns the processor on which element with global index i is stored + template + SizeType operator()(SizeType i) const + { + SizeType cutoff_processor = n % p; + SizeType cutoff = cutoff_processor * (n / p + 1); + + if (i < cutoff) return i / (n / p + 1); + else return cutoff_processor + (i - cutoff) / (n / p); + } + + // Find the starting index for processor with the given id + template + std::size_t start(ID id) const + { + std::size_t estimate = id * (n / p + 1); + ID cutoff_processor = n % p; + if (id < cutoff_processor) return estimate; + else return estimate - (id - cutoff_processor); + } + + // Find the local index for the ith global element + template + SizeType local(SizeType i) const + { + SizeType owner = (*this)(i); + return i - start(owner); + } + + // Returns the global index of local element i + template + SizeType global(SizeType i) const + { return global(id, i); } + + // Returns the global index of the ith local element on processor id + template + SizeType global(ProcessID id, SizeType i) const + { return i + start(id); } + + private: + std::size_t id; //< The ID number of this processor + std::size_t p; //< The number of processors + std::size_t n; //< The size of the problem space +}; + +// Block distribution with arbitrary block sizes +struct uneven_block +{ + typedef std::vector size_vector; + + template + explicit uneven_block(const LinearProcessGroup& pg, const std::vector& local_sizes) + : id(process_id(pg)), p(num_processes(pg)), local_sizes(local_sizes) + { + BOOST_ASSERT(local_sizes.size() == p); + local_starts.resize(p + 1); + local_starts[0] = 0; + std::partial_sum(local_sizes.begin(), local_sizes.end(), &local_starts[1]); + n = local_starts[p]; + } + + // To do maybe: enter local size in each process and gather in constructor (much handier) + // template + // explicit uneven_block(const LinearProcessGroup& pg, std::size_t my_local_size) + + // If there are n elements in the distributed data structure, returns the number of elements stored locally. + template + SizeType block_size(SizeType) const + { return local_sizes[id]; } + + // If there are n elements in the distributed data structure, returns the number of elements stored on processor ID + template + SizeType block_size(ProcessID id, SizeType) const + { return local_sizes[id]; } + + // Returns the processor on which element with global index i is stored + template + SizeType operator()(SizeType i) const + { + BOOST_ASSERT (i >= (SizeType) 0 && i < (SizeType) n); // check for valid range + size_vector::const_iterator lb = std::lower_bound(local_starts.begin(), local_starts.end(), (std::size_t) i); + return ((SizeType)(*lb) == i ? lb : --lb) - local_starts.begin(); + } + + // Find the starting index for processor with the given id + template + std::size_t start(ID id) const + { + return local_starts[id]; + } + + // Find the local index for the ith global element + template + SizeType local(SizeType i) const + { + SizeType owner = (*this)(i); + return i - start(owner); + } + + // Returns the global index of local element i + template + SizeType global(SizeType i) const + { return global(id, i); } + + // Returns the global index of the ith local element on processor id + template + SizeType global(ProcessID id, SizeType i) const + { return i + start(id); } + + private: + std::size_t id; //< The ID number of this processor + std::size_t p; //< The number of processors + std::size_t n; //< The size of the problem space + std::vector local_sizes; //< The sizes of all blocks + std::vector local_starts; //< Lowest global index of each block +}; + + +struct oned_block_cyclic +{ + template + explicit oned_block_cyclic(const LinearProcessGroup& pg, std::size_t size) + : id(process_id(pg)), p(num_processes(pg)), size(size) { } + + template + SizeType block_size(SizeType n) const + { + return block_size(id, n); + } + + template + SizeType block_size(ProcessID id, SizeType n) const + { + SizeType all_blocks = n / size; + SizeType extra_elements = n % size; + SizeType everyone_gets = all_blocks / p; + SizeType extra_blocks = all_blocks % p; + SizeType my_blocks = everyone_gets + (p < extra_blocks? 1 : 0); + SizeType my_elements = my_blocks * size + + (p == extra_blocks? extra_elements : 0); + return my_elements; + } + + template + SizeType operator()(SizeType i) const + { + return (i / size) % p; + } + + template + SizeType local(SizeType i) const + { + return ((i / size) / p) * size + i % size; + } + + template + SizeType global(SizeType i) const + { return global(id, i); } + + template + SizeType global(ProcessID id, SizeType i) const + { + return ((i / size) * p + id) * size + i % size; + } + + private: + std::size_t id; //< The ID number of this processor + std::size_t p; //< The number of processors + std::size_t size; //< Block size +}; + +struct twod_block_cyclic +{ + template + explicit twod_block_cyclic(const LinearProcessGroup& pg, + std::size_t block_rows, std::size_t block_columns, + std::size_t data_columns_per_row) + : id(process_id(pg)), p(num_processes(pg)), + block_rows(block_rows), block_columns(block_columns), + data_columns_per_row(data_columns_per_row) + { } + + template + SizeType block_size(SizeType n) const + { + return block_size(id, n); + } + + template + SizeType block_size(ProcessID id, SizeType n) const + { + // TBD: This is really lame :) + int result = -1; + while (n > 0) { + --n; + if ((*this)(n) == id && (int)local(n) > result) result = local(n); + } + ++result; + + // std::cerr << "Block size of id " << id << " is " << result << std::endl; + return result; + } + + template + SizeType operator()(SizeType i) const + { + SizeType result = get_block_num(i) % p; + // std::cerr << "Item " << i << " goes on processor " << result << std::endl; + return result; + } + + template + SizeType local(SizeType i) const + { + // Compute the start of the block + std::size_t block_num = get_block_num(i); + // std::cerr << "Item " << i << " is in block #" << block_num << std::endl; + + std::size_t local_block_num = block_num / p; + std::size_t block_start = local_block_num * block_rows * block_columns; + + // Compute the offset into the block + std::size_t data_row = i / data_columns_per_row; + std::size_t data_col = i % data_columns_per_row; + std::size_t block_offset = (data_row % block_rows) * block_columns + + (data_col % block_columns); + + // std::cerr << "Item " << i << " maps to local index " << block_start+block_offset << std::endl; + return block_start + block_offset; + } + + template + SizeType global(SizeType i) const + { + // Compute the (global) block in which this element resides + SizeType local_block_num = i / (block_rows * block_columns); + SizeType block_offset = i % (block_rows * block_columns); + SizeType block_num = local_block_num * p + id; + + // Compute the position of the start of the block (globally) + SizeType block_start = block_num * block_rows * block_columns; + + std::cerr << "Block " << block_num << " starts at index " << block_start + << std::endl; + + // Compute the row and column of this block + SizeType block_row = block_num / (data_columns_per_row / block_columns); + SizeType block_col = block_num % (data_columns_per_row / block_columns); + + SizeType row_in_block = block_offset / block_columns; + SizeType col_in_block = block_offset % block_columns; + + std::cerr << "Local index " << i << " is in block at row " << block_row + << ", column " << block_col << ", in-block row " << row_in_block + << ", in-block col " << col_in_block << std::endl; + + SizeType result = block_row * block_rows + block_col * block_columns + + row_in_block * block_rows + col_in_block; + + + std::cerr << "global(" << i << "@" << id << ") = " << result + << " =? " << local(result) << std::endl; + BOOST_ASSERT(i == local(result)); + return result; + } + + private: + template + std::size_t get_block_num(SizeType i) const + { + std::size_t data_row = i / data_columns_per_row; + std::size_t data_col = i % data_columns_per_row; + std::size_t block_row = data_row / block_rows; + std::size_t block_col = data_col / block_columns; + std::size_t blocks_in_row = data_columns_per_row / block_columns; + std::size_t block_num = block_col * blocks_in_row + block_row; + return block_num; + } + + std::size_t id; //< The ID number of this processor + std::size_t p; //< The number of processors + std::size_t block_rows; //< The # of rows in each block + std::size_t block_columns; //< The # of columns in each block + std::size_t data_columns_per_row; //< The # of columns per row of data +}; + +class twod_random +{ + template + struct random_int + { + explicit random_int(RandomNumberGen& gen) : gen(gen) { } + + template + T operator()(T n) const + { + uniform_int distrib(0, n-1); + return distrib(gen); + } + + private: + RandomNumberGen& gen; + }; + + public: + template + explicit twod_random(const LinearProcessGroup& pg, + std::size_t block_rows, std::size_t block_columns, + std::size_t data_columns_per_row, + std::size_t n, + RandomNumberGen& gen) + : id(process_id(pg)), p(num_processes(pg)), + block_rows(block_rows), block_columns(block_columns), + data_columns_per_row(data_columns_per_row), + global_to_local(n / (block_rows * block_columns)) + { + std::copy(make_counting_iterator(std::size_t(0)), + make_counting_iterator(global_to_local.size()), + global_to_local.begin()); + + random_int rand(gen); + std::random_shuffle(global_to_local.begin(), global_to_local.end(), rand); + } + + template + SizeType block_size(SizeType n) const + { + return block_size(id, n); + } + + template + SizeType block_size(ProcessID id, SizeType n) const + { + // TBD: This is really lame :) + int result = -1; + while (n > 0) { + --n; + if ((*this)(n) == id && (int)local(n) > result) result = local(n); + } + ++result; + + // std::cerr << "Block size of id " << id << " is " << result << std::endl; + return result; + } + + template + SizeType operator()(SizeType i) const + { + SizeType result = get_block_num(i) % p; + // std::cerr << "Item " << i << " goes on processor " << result << std::endl; + return result; + } + + template + SizeType local(SizeType i) const + { + // Compute the start of the block + std::size_t block_num = get_block_num(i); + // std::cerr << "Item " << i << " is in block #" << block_num << std::endl; + + std::size_t local_block_num = block_num / p; + std::size_t block_start = local_block_num * block_rows * block_columns; + + // Compute the offset into the block + std::size_t data_row = i / data_columns_per_row; + std::size_t data_col = i % data_columns_per_row; + std::size_t block_offset = (data_row % block_rows) * block_columns + + (data_col % block_columns); + + // std::cerr << "Item " << i << " maps to local index " << block_start+block_offset << std::endl; + return block_start + block_offset; + } + + private: + template + std::size_t get_block_num(SizeType i) const + { + std::size_t data_row = i / data_columns_per_row; + std::size_t data_col = i % data_columns_per_row; + std::size_t block_row = data_row / block_rows; + std::size_t block_col = data_col / block_columns; + std::size_t blocks_in_row = data_columns_per_row / block_columns; + std::size_t block_num = block_col * blocks_in_row + block_row; + return global_to_local[block_num]; + } + + std::size_t id; //< The ID number of this processor + std::size_t p; //< The number of processors + std::size_t block_rows; //< The # of rows in each block + std::size_t block_columns; //< The # of columns in each block + std::size_t data_columns_per_row; //< The # of columns per row of data + std::vector global_to_local; +}; + +class random_distribution +{ + template + struct random_int + { + explicit random_int(RandomNumberGen& gen) : gen(gen) { } + + template + T operator()(T n) const + { + uniform_int distrib(0, n-1); + return distrib(gen); + } + + private: + RandomNumberGen& gen; + }; + + public: + template + random_distribution(const LinearProcessGroup& pg, RandomNumberGen& gen, + std::size_t n) + : base(pg, n), local_to_global(n), global_to_local(n) + { + std::copy(make_counting_iterator(std::size_t(0)), + make_counting_iterator(n), + local_to_global.begin()); + + random_int rand(gen); + std::random_shuffle(local_to_global.begin(), local_to_global.end(), rand); + + + for (std::vector::size_type i = 0; i < n; ++i) + global_to_local[local_to_global[i]] = i; + } + + template + SizeType block_size(SizeType n) const + { return base.block_size(n); } + + template + SizeType block_size(ProcessID id, SizeType n) const + { return base.block_size(id, n); } + + template + SizeType operator()(SizeType i) const + { + return base(global_to_local[i]); + } + + template + SizeType local(SizeType i) const + { + return base.local(global_to_local[i]); + } + + template + SizeType global(ProcessID p, SizeType i) const + { + return local_to_global[base.global(p, i)]; + } + + template + SizeType global(SizeType i) const + { + return local_to_global[base.global(i)]; + } + + private: + block base; + std::vector local_to_global; + std::vector global_to_local; +}; + +} } // end namespace boost::parallel + +#endif // BOOST_PARALLEL_DISTRIBUTION_HPP + diff --git a/patched/math/constants/info.hpp b/patched/math/constants/info.hpp new file mode 100644 index 0000000..579ad4d --- /dev/null +++ b/patched/math/constants/info.hpp @@ -0,0 +1,163 @@ +// Copyright John Maddock 2010. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifdef _MSC_VER +# pragma once +#endif + +#ifndef BOOST_MATH_CONSTANTS_INFO_INCLUDED +#define BOOST_MATH_CONSTANTS_INFO_INCLUDED + +#include +#include +#include +#include + +namespace boost{ namespace math{ namespace constants{ + + namespace detail{ + + template + const char* nameof(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) + { + return type_id().name_demangled(); + } + template <> + const char* nameof(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(float)) + { + return "float"; + } + template <> + const char* nameof(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(double)) + { + return "double"; + } + template <> + const char* nameof(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(long double)) + { + return "long double"; + } + + } + +template +void print_info_on_type(std::ostream& os = std::cout BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T) BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(Policy)) +{ + using detail::nameof; +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + os << + "Information on the Implementation and Handling of \n" + "Mathematical Constants for Type " << nameof() << + "\n\n" + "Checking for std::numeric_limits<" << nameof() << "> specialisation: " << + (std::numeric_limits::is_specialized ? "yes" : "no") << std::endl; + if(std::numeric_limits::is_specialized) + { + os << + "std::numeric_limits<" << nameof() << ">::digits reports that the radix is " << std::numeric_limits::radix << ".\n"; + if (std::numeric_limits::radix == 2) + { + os << + "std::numeric_limits<" << nameof() << ">::digits reports that the precision is \n" << std::numeric_limits::digits << " binary digits.\n"; + } + else if (std::numeric_limits::radix == 10) + { + os << + "std::numeric_limits<" << nameof() << ">::digits reports that the precision is \n" << std::numeric_limits::digits10 << " decimal digits.\n"; + os << + "std::numeric_limits<" << nameof() << ">::digits reports that the precision is \n" + << std::numeric_limits::digits * 1000L /301L << " binary digits.\n"; // divide by log2(10) - about 3 bits per decimal digit. + } + else + { + os << "Unknown radix = " << std::numeric_limits::radix << "\n"; + } + } + typedef typename boost::math::policies::precision::type precision_type; + if(precision_type::value) + { + if (std::numeric_limits::radix == 2) + { + os << + "boost::math::policies::precision<" << nameof() << ", " << nameof() << " reports that the compile time precision is \n" << precision_type::value << " binary digits.\n"; + } + else if (std::numeric_limits::radix == 10) + { + os << + "boost::math::policies::precision<" << nameof() << ", " << nameof() << " reports that the compile time precision is \n" << precision_type::value << " binary digits.\n"; + } + else + { + os << "Unknown radix = " << std::numeric_limits::radix << "\n"; + } + } + else + { + os << + "boost::math::policies::precision<" << nameof() << ", Policy> \n" + "reports that there is no compile type precision available.\n" + "boost::math::tools::digits<" << nameof() << ">() \n" + "reports that the current runtime precision is \n" << + boost::math::tools::digits() << " binary digits.\n"; + } + + typedef typename construction_traits::type construction_type; + + switch(construction_type::value) + { + case 0: + os << + "No compile time precision is available, the construction method \n" + "will be decided at runtime and results will not be cached \n" + "- this may lead to poor runtime performance.\n" + "Current runtime precision indicates that\n"; + if(boost::math::tools::digits() > max_string_digits) + { + os << "the constant will be recalculated on each call.\n"; + } + else + { + os << "the constant will be constructed from a string on each call.\n"; + } + break; + case 1: + os << + "The constant will be constructed from a float.\n"; + break; + case 2: + os << + "The constant will be constructed from a double.\n"; + break; + case 3: + os << + "The constant will be constructed from a long double.\n"; + break; + case 4: + os << + "The constant will be constructed from a string (and the result cached).\n"; + break; + default: + os << + "The constant will be calculated (and the result cached).\n"; + break; + } + os << std::endl; +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif +} + +template +void print_info_on_type(std::ostream& os = std::cout BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) +{ + print_info_on_type >(os); +} + +}}} // namespaces + +#endif // BOOST_MATH_CONSTANTS_INFO_INCLUDED diff --git a/patched/math/distributions/detail/hypergeometric_pdf.hpp b/patched/math/distributions/detail/hypergeometric_pdf.hpp new file mode 100644 index 0000000..bf22fc9 --- /dev/null +++ b/patched/math/distributions/detail/hypergeometric_pdf.hpp @@ -0,0 +1,488 @@ +// Copyright 2008 Gautam Sewani +// Copyright 2008 John Maddock +// +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP +#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP + +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_MATH_INSTRUMENT +#include +#endif + +namespace boost{ namespace math{ namespace detail{ + +template +void bubble_down_one(T* first, T* last, Func f) +{ + using std::swap; + T* next = first; + ++next; + while((next != last) && (!f(*first, *next))) + { + swap(*first, *next); + ++first; + ++next; + } +} + +template +struct sort_functor +{ + sort_functor(const T* exponents) : m_exponents(exponents){} + bool operator()(int i, int j) + { + return m_exponents[i] > m_exponents[j]; + } +private: + const T* m_exponents; +}; + +template +T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&) +{ + BOOST_MATH_STD_USING + + BOOST_MATH_INSTRUMENT_FPU + BOOST_MATH_INSTRUMENT_VARIABLE(x); + BOOST_MATH_INSTRUMENT_VARIABLE(r); + BOOST_MATH_INSTRUMENT_VARIABLE(n); + BOOST_MATH_INSTRUMENT_VARIABLE(N); + BOOST_MATH_INSTRUMENT_VARIABLE(type_id().name_demangled()); + + T bases[9] = { + T(n) + Lanczos::g() + 0.5f, + T(r) + Lanczos::g() + 0.5f, + T(N - n) + Lanczos::g() + 0.5f, + T(N - r) + Lanczos::g() + 0.5f, + 1 / (T(N) + Lanczos::g() + 0.5f), + 1 / (T(x) + Lanczos::g() + 0.5f), + 1 / (T(n - x) + Lanczos::g() + 0.5f), + 1 / (T(r - x) + Lanczos::g() + 0.5f), + 1 / (T(N - n - r + x) + Lanczos::g() + 0.5f) + }; + T exponents[9] = { + n + T(0.5f), + r + T(0.5f), + N - n + T(0.5f), + N - r + T(0.5f), + N + T(0.5f), + x + T(0.5f), + n - x + T(0.5f), + r - x + T(0.5f), + N - n - r + x + T(0.5f) + }; + int base_e_factors[9] = { + -1, -1, -1, -1, 1, 1, 1, 1, 1 + }; + int sorted_indexes[9] = { + 0, 1, 2, 3, 4, 5, 6, 7, 8 + }; +#ifdef BOOST_MATH_INSTRUMENT + BOOST_MATH_INSTRUMENT_FPU + for(unsigned i = 0; i < 9; ++i) + { + BOOST_MATH_INSTRUMENT_VARIABLE(i); + BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); + } +#endif + std::sort(sorted_indexes, sorted_indexes + 9, sort_functor(exponents)); +#ifdef BOOST_MATH_INSTRUMENT + BOOST_MATH_INSTRUMENT_FPU + for(unsigned i = 0; i < 9; ++i) + { + BOOST_MATH_INSTRUMENT_VARIABLE(i); + BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); + } +#endif + + do{ + exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]]; + bases[sorted_indexes[1]] *= bases[sorted_indexes[0]]; + if((bases[sorted_indexes[1]] < tools::min_value()) && (exponents[sorted_indexes[1]] != 0)) + { + return 0; + } + base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]]; + bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor(exponents)); + +#ifdef BOOST_MATH_INSTRUMENT + for(unsigned i = 0; i < 9; ++i) + { + BOOST_MATH_INSTRUMENT_VARIABLE(i); + BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); + } +#endif + }while(exponents[sorted_indexes[1]] > 1); + + // + // Combine equal powers: + // + int j = 8; + while(exponents[sorted_indexes[j]] == 0) --j; + while(j) + { + while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]])) + { + bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]]; + exponents[sorted_indexes[j]] = 0; + base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]]; + bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor(exponents)); + --j; + } + --j; + +#ifdef BOOST_MATH_INSTRUMENT + BOOST_MATH_INSTRUMENT_VARIABLE(j); + for(unsigned i = 0; i < 9; ++i) + { + BOOST_MATH_INSTRUMENT_VARIABLE(i); + BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); + } +#endif + } + +#ifdef BOOST_MATH_INSTRUMENT + BOOST_MATH_INSTRUMENT_FPU + for(unsigned i = 0; i < 9; ++i) + { + BOOST_MATH_INSTRUMENT_VARIABLE(i); + BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); + BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); + } +#endif + + T result; + BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast(base_e_factors[sorted_indexes[0]]))); + BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]); + { + BOOST_FPU_EXCEPTION_GUARD + result = pow(bases[sorted_indexes[0]] * exp(static_cast(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]); + } + BOOST_MATH_INSTRUMENT_VARIABLE(result); + for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i) + { + BOOST_FPU_EXCEPTION_GUARD + if(result < tools::min_value()) + return 0; // short circuit further evaluation + if(exponents[sorted_indexes[i]] == 1) + result *= bases[sorted_indexes[i]] * exp(static_cast(base_e_factors[sorted_indexes[i]])); + else if(exponents[sorted_indexes[i]] == 0.5f) + result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast(base_e_factors[sorted_indexes[i]]))); + else + result *= pow(bases[sorted_indexes[i]] * exp(static_cast(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]); + + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + + result *= Lanczos::lanczos_sum_expG_scaled(static_cast(n + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(r + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(N - n + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(N - r + 1)) + / + ( Lanczos::lanczos_sum_expG_scaled(static_cast(N + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(x + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(n - x + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(r - x + 1)) + * Lanczos::lanczos_sum_expG_scaled(static_cast(N - n - r + x + 1))); + + BOOST_MATH_INSTRUMENT_VARIABLE(result); + return result; +} + +template +T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol) +{ + BOOST_MATH_STD_USING + return exp( + boost::math::lgamma(T(n + 1), pol) + + boost::math::lgamma(T(r + 1), pol) + + boost::math::lgamma(T(N - n + 1), pol) + + boost::math::lgamma(T(N - r + 1), pol) + - boost::math::lgamma(T(N + 1), pol) + - boost::math::lgamma(T(x + 1), pol) + - boost::math::lgamma(T(n - x + 1), pol) + - boost::math::lgamma(T(r - x + 1), pol) + - boost::math::lgamma(T(N - n - r + x + 1), pol)); +} + +template +inline T integer_power(const T& x, int ex) +{ + if(ex < 0) + return 1 / integer_power(x, -ex); + switch(ex) + { + case 0: + return 1; + case 1: + return x; + case 2: + return x * x; + case 3: + return x * x * x; + case 4: + return boost::math::pow<4>(x); + case 5: + return boost::math::pow<5>(x); + case 6: + return boost::math::pow<6>(x); + case 7: + return boost::math::pow<7>(x); + case 8: + return boost::math::pow<8>(x); + } + BOOST_MATH_STD_USING +#ifdef __SUNPRO_CC + return pow(x, T(ex)); +#else + return pow(x, ex); +#endif +} +template +struct hypergeometric_pdf_prime_loop_result_entry +{ + T value; + const hypergeometric_pdf_prime_loop_result_entry* next; +}; + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4510 4512 4610) +#endif + +struct hypergeometric_pdf_prime_loop_data +{ + const unsigned x; + const unsigned r; + const unsigned n; + const unsigned N; + unsigned prime_index; + unsigned current_prime; +}; + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +template +T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry& result) +{ + while(data.current_prime <= data.N) + { + unsigned base = data.current_prime; + int prime_powers = 0; + while(base <= data.N) + { + prime_powers += data.n / base; + prime_powers += data.r / base; + prime_powers += (data.N - data.n) / base; + prime_powers += (data.N - data.r) / base; + prime_powers -= data.N / base; + prime_powers -= data.x / base; + prime_powers -= (data.n - data.x) / base; + prime_powers -= (data.r - data.x) / base; + prime_powers -= (data.N - data.n - data.r + data.x) / base; + base *= data.current_prime; + } + if(prime_powers) + { + T p = integer_power(data.current_prime, prime_powers); + if((p > 1) && (tools::max_value() / p < result.value)) + { + // + // The next calculation would overflow, use recursion + // to sidestep the issue: + // + hypergeometric_pdf_prime_loop_result_entry t = { p, &result }; + data.current_prime = prime(++data.prime_index); + return hypergeometric_pdf_prime_loop_imp(data, t); + } + if((p < 1) && (tools::min_value() / p > result.value)) + { + // + // The next calculation would underflow, use recursion + // to sidestep the issue: + // + hypergeometric_pdf_prime_loop_result_entry t = { p, &result }; + data.current_prime = prime(++data.prime_index); + return hypergeometric_pdf_prime_loop_imp(data, t); + } + result.value *= p; + } + data.current_prime = prime(++data.prime_index); + } + // + // When we get to here we have run out of prime factors, + // the overall result is the product of all the partial + // results we have accumulated on the stack so far, these + // are in a linked list starting with "data.head" and ending + // with "result". + // + // All that remains is to multiply them together, taking + // care not to overflow or underflow. + // + // Enumerate partial results >= 1 in variable i + // and partial results < 1 in variable j: + // + hypergeometric_pdf_prime_loop_result_entry const *i, *j; + i = &result; + while(i && i->value < 1) + i = i->next; + j = &result; + while(j && j->value >= 1) + j = j->next; + + T prod = 1; + + while(i || j) + { + while(i && ((prod <= 1) || (j == 0))) + { + prod *= i->value; + i = i->next; + while(i && i->value < 1) + i = i->next; + } + while(j && ((prod >= 1) || (i == 0))) + { + prod *= j->value; + j = j->next; + while(j && j->value >= 1) + j = j->next; + } + } + + return prod; +} + +template +inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&) +{ + hypergeometric_pdf_prime_loop_result_entry result = { 1, 0 }; + hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) }; + return hypergeometric_pdf_prime_loop_imp(data, result); +} + +template +T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&) +{ + BOOST_MATH_STD_USING + BOOST_ASSERT(N < boost::math::max_factorial::value); + T result = boost::math::unchecked_factorial(n); + T num[3] = { + boost::math::unchecked_factorial(r), + boost::math::unchecked_factorial(N - n), + boost::math::unchecked_factorial(N - r) + }; + T denom[5] = { + boost::math::unchecked_factorial(N), + boost::math::unchecked_factorial(x), + boost::math::unchecked_factorial(n - x), + boost::math::unchecked_factorial(r - x), + boost::math::unchecked_factorial(N - n - r + x) + }; + int i = 0; + int j = 0; + while((i < 3) || (j < 5)) + { + while((j < 5) && ((result >= 1) || (i >= 3))) + { + result /= denom[j]; + ++j; + } + while((i < 3) && ((result <= 1) || (j >= 5))) + { + result *= num[i]; + ++i; + } + } + return result; +} + + +template +inline typename tools::promote_args::type + hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + value_type result; + if(N <= boost::math::max_factorial::value) + { + // + // If N is small enough then we can evaluate the PDF via the factorials + // directly: table lookup of the factorials gives the best performance + // of the methods available: + // + result = detail::hypergeometric_pdf_factorial_imp(x, r, n, N, forwarding_policy()); + } + else if(N <= boost::math::prime(boost::math::max_prime - 1)) + { + // + // If N is no larger than the largest prime number in our lookup table + // (104729) then we can use prime factorisation to evaluate the PDF, + // this is slow but accurate: + // + result = detail::hypergeometric_pdf_prime_imp(x, r, n, N, forwarding_policy()); + } + else + { + // + // Catch all case - use the lanczos approximation - where available - + // to evaluate the ratio of factorials. This is reasonably fast + // (almost as quick as using logarithmic evaluation in terms of lgamma) + // but only a few digits better in accuracy than using lgamma: + // + result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy()); + } + + if(result > 1) + { + result = 1; + } + if(result < 0) + { + result = 0; + } + + return policies::checked_narrowing_cast(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)"); +} + +}}} // namespaces + +#endif + diff --git a/patched/math/policies/error_handling.hpp b/patched/math/policies/error_handling.hpp new file mode 100644 index 0000000..862071f --- /dev/null +++ b/patched/math/policies/error_handling.hpp @@ -0,0 +1,692 @@ +// Copyright John Maddock 2007. +// Copyright Paul A. Bristow 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_POLICY_ERROR_HANDLING_HPP +#define BOOST_MATH_POLICY_ERROR_HANDLING_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef BOOST_MSVC +# pragma warning(push) // Quiet warnings in boost/format.hpp +# pragma warning(disable: 4996) // _SCL_SECURE_NO_DEPRECATE +# pragma warning(disable: 4512) // assignment operator could not be generated. +// And warnings in error handling: +# pragma warning(disable: 4702) // unreachable code +// Note that this only occurs when the compiler can deduce code is unreachable, +// for example when policy macros are used to ignore errors rather than throw. +#endif +#include +#include + +namespace boost{ namespace math{ + +class evaluation_error : public std::runtime_error +{ +public: + evaluation_error(const std::string& s) : std::runtime_error(s){} +}; + +class rounding_error : public std::runtime_error +{ +public: + rounding_error(const std::string& s) : std::runtime_error(s){} +}; + +namespace policies{ +// +// Forward declarations of user error handlers, +// it's up to the user to provide the definition of these: +// +template +T user_domain_error(const char* function, const char* message, const T& val); +template +T user_pole_error(const char* function, const char* message, const T& val); +template +T user_overflow_error(const char* function, const char* message, const T& val); +template +T user_underflow_error(const char* function, const char* message, const T& val); +template +T user_denorm_error(const char* function, const char* message, const T& val); +template +T user_evaluation_error(const char* function, const char* message, const T& val); +template +T user_rounding_error(const char* function, const char* message, const T& val, const TargetType& t); +template +T user_indeterminate_result_error(const char* function, const char* message, const T& val); + +namespace detail +{ +// +// Helper function to avoid binding rvalue to non-const-reference, +// in other words a warning suppression mechansim: +// +template +inline std::string do_format(Formatter f, const Group& g) +{ + return (f % g).str(); +} + +template +void raise_error(const char* function, const char* message) +{ + if(function == 0) + function = "Unknown function operating on type %1%"; + if(message == 0) + message = "Cause unknown"; + + std::string msg("Error in function "); + msg += (boost::format(function) % type_id().name_demangled()).str(); + msg += ": "; + msg += message; + + E e(msg); + boost::throw_exception(e); +} + +template +void raise_error(const char* function, const char* message, const T& val) +{ + if(function == 0) + function = "Unknown function operating on type %1%"; + if(message == 0) + message = "Cause unknown: error caused by bad argument with value %1%"; + + std::string msg("Error in function "); + msg += (boost::format(function) % type_id().name_demangled()).str(); + msg += ": "; + msg += message; + + int prec = 2 + (boost::math::policies::digits >() * 30103UL) / 100000UL; + msg = do_format(boost::format(msg), boost::io::group(std::setprecision(prec), val)); + + E e(msg); + boost::throw_exception(e); +} + +template +inline T raise_domain_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::domain_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message, val); + // we never get here: + return std::numeric_limits::quiet_NaN(); +} + +template +inline T raise_domain_error( + const char* , + const char* , + const T& , + const ::boost::math::policies::domain_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return std::numeric_limits::quiet_NaN(); +} + +template +inline T raise_domain_error( + const char* , + const char* , + const T& , + const ::boost::math::policies::domain_error< ::boost::math::policies::errno_on_error>&) +{ + errno = EDOM; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return std::numeric_limits::quiet_NaN(); +} + +template +inline T raise_domain_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::domain_error< ::boost::math::policies::user_error>&) +{ + return user_domain_error(function, message, val); +} + +template +inline T raise_pole_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::pole_error< ::boost::math::policies::throw_on_error>&) +{ + return boost::math::policies::detail::raise_domain_error(function, message, val, ::boost::math::policies::domain_error< ::boost::math::policies::throw_on_error>()); +} + +template +inline T raise_pole_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::pole_error< ::boost::math::policies::ignore_error>&) +{ + return ::boost::math::policies::detail::raise_domain_error(function, message, val, ::boost::math::policies::domain_error< ::boost::math::policies::ignore_error>()); +} + +template +inline T raise_pole_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::pole_error< ::boost::math::policies::errno_on_error>&) +{ + return ::boost::math::policies::detail::raise_domain_error(function, message, val, ::boost::math::policies::domain_error< ::boost::math::policies::errno_on_error>()); +} + +template +inline T raise_pole_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::pole_error< ::boost::math::policies::user_error>&) +{ + return user_pole_error(function, message, val); +} + +template +inline T raise_overflow_error( + const char* function, + const char* message, + const ::boost::math::policies::overflow_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message ? message : "numeric overflow"); + // we never get here: + return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(); +} + +template +inline T raise_overflow_error( + const char* , + const char* , + const ::boost::math::policies::overflow_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(); +} + +template +inline T raise_overflow_error( + const char* , + const char* , + const ::boost::math::policies::overflow_error< ::boost::math::policies::errno_on_error>&) +{ + errno = ERANGE; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(); +} + +template +inline T raise_overflow_error( + const char* function, + const char* message, + const ::boost::math::policies::overflow_error< ::boost::math::policies::user_error>&) +{ + return user_overflow_error(function, message, std::numeric_limits::infinity()); +} + +template +inline T raise_underflow_error( + const char* function, + const char* message, + const ::boost::math::policies::underflow_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message ? message : "numeric underflow"); + // we never get here: + return 0; +} + +template +inline T raise_underflow_error( + const char* , + const char* , + const ::boost::math::policies::underflow_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return T(0); +} + +template +inline T raise_underflow_error( + const char* /* function */, + const char* /* message */, + const ::boost::math::policies::underflow_error< ::boost::math::policies::errno_on_error>&) +{ + errno = ERANGE; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return T(0); +} + +template +inline T raise_underflow_error( + const char* function, + const char* message, + const ::boost::math::policies::underflow_error< ::boost::math::policies::user_error>&) +{ + return user_underflow_error(function, message, T(0)); +} + +template +inline T raise_denorm_error( + const char* function, + const char* message, + const T& /* val */, + const ::boost::math::policies::denorm_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message ? message : "denormalised result"); + // we never get here: + return T(0); +} + +template +inline T raise_denorm_error( + const char* , + const char* , + const T& val, + const ::boost::math::policies::denorm_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return val; +} + +template +inline T raise_denorm_error( + const char* , + const char* , + const T& val, + const ::boost::math::policies::denorm_error< ::boost::math::policies::errno_on_error>&) +{ + errno = ERANGE; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return val; +} + +template +inline T raise_denorm_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::denorm_error< ::boost::math::policies::user_error>&) +{ + return user_denorm_error(function, message, val); +} + +template +inline T raise_evaluation_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::evaluation_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message, val); + // we never get here: + return T(0); +} + +template +inline T raise_evaluation_error( + const char* , + const char* , + const T& val, + const ::boost::math::policies::evaluation_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return val; +} + +template +inline T raise_evaluation_error( + const char* , + const char* , + const T& val, + const ::boost::math::policies::evaluation_error< ::boost::math::policies::errno_on_error>&) +{ + errno = EDOM; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return val; +} + +template +inline T raise_evaluation_error( + const char* function, + const char* message, + const T& val, + const ::boost::math::policies::evaluation_error< ::boost::math::policies::user_error>&) +{ + return user_evaluation_error(function, message, val); +} + +template +inline T raise_rounding_error( + const char* function, + const char* message, + const T& val, + const TargetType&, + const ::boost::math::policies::rounding_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message, val); + // we never get here: + return T(0); +} + +template +inline T raise_rounding_error( + const char* , + const char* , + const T& val, + const TargetType&, + const ::boost::math::policies::rounding_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return std::numeric_limits::is_specialized ? (val > 0 ? (std::numeric_limits::max)() : -(std::numeric_limits::max)()): val; +} + +template +inline T raise_rounding_error( + const char* , + const char* , + const T& val, + const TargetType&, + const ::boost::math::policies::rounding_error< ::boost::math::policies::errno_on_error>&) +{ + errno = ERANGE; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return std::numeric_limits::is_specialized ? (val > 0 ? (std::numeric_limits::max)() : -(std::numeric_limits::max)()): val; +} + +template +inline T raise_rounding_error( + const char* function, + const char* message, + const T& val, + const TargetType& t, + const ::boost::math::policies::rounding_error< ::boost::math::policies::user_error>&) +{ + return user_rounding_error(function, message, val, t); +} + +template +inline T raise_indeterminate_result_error( + const char* function, + const char* message, + const T& val, + const R& , + const ::boost::math::policies::indeterminate_result_error< ::boost::math::policies::throw_on_error>&) +{ + raise_error(function, message, val); + // we never get here: + return std::numeric_limits::quiet_NaN(); +} + +template +inline T raise_indeterminate_result_error( + const char* , + const char* , + const T& , + const R& result, + const ::boost::math::policies::indeterminate_result_error< ::boost::math::policies::ignore_error>&) +{ + // This may or may not do the right thing, but the user asked for the error + // to be ignored so here we go anyway: + return result; +} + +template +inline T raise_indeterminate_result_error( + const char* , + const char* , + const T& , + const R& result, + const ::boost::math::policies::indeterminate_result_error< ::boost::math::policies::errno_on_error>&) +{ + errno = EDOM; + // This may or may not do the right thing, but the user asked for the error + // to be silent so here we go anyway: + return result; +} + +template +inline T raise_indeterminate_result_error( + const char* function, + const char* message, + const T& val, + const R& , + const ::boost::math::policies::indeterminate_result_error< ::boost::math::policies::user_error>&) +{ + return user_indeterminate_result_error(function, message, val); +} + +} // namespace detail + +template +inline T raise_domain_error(const char* function, const char* message, const T& val, const Policy&) +{ + typedef typename Policy::domain_error_type policy_type; + return detail::raise_domain_error( + function, message ? message : "Domain Error evaluating function at %1%", + val, policy_type()); +} + +template +inline T raise_pole_error(const char* function, const char* message, const T& val, const Policy&) +{ + typedef typename Policy::pole_error_type policy_type; + return detail::raise_pole_error( + function, message ? message : "Evaluation of function at pole %1%", + val, policy_type()); +} + +template +inline T raise_overflow_error(const char* function, const char* message, const Policy&) +{ + typedef typename Policy::overflow_error_type policy_type; + return detail::raise_overflow_error( + function, message ? message : "Overflow Error", + policy_type()); +} + +template +inline T raise_underflow_error(const char* function, const char* message, const Policy&) +{ + typedef typename Policy::underflow_error_type policy_type; + return detail::raise_underflow_error( + function, message ? message : "Underflow Error", + policy_type()); +} + +template +inline T raise_denorm_error(const char* function, const char* message, const T& val, const Policy&) +{ + typedef typename Policy::denorm_error_type policy_type; + return detail::raise_denorm_error( + function, message ? message : "Denorm Error", + val, + policy_type()); +} + +template +inline T raise_evaluation_error(const char* function, const char* message, const T& val, const Policy&) +{ + typedef typename Policy::evaluation_error_type policy_type; + return detail::raise_evaluation_error( + function, message ? message : "Internal Evaluation Error, best value so far was %1%", + val, policy_type()); +} + +template +inline T raise_rounding_error(const char* function, const char* message, const T& val, const TargetType& t, const Policy&) +{ + typedef typename Policy::rounding_error_type policy_type; + return detail::raise_rounding_error( + function, message ? message : "Value %1% can not be represented in the target integer type.", + val, t, policy_type()); +} + +template +inline T raise_indeterminate_result_error(const char* function, const char* message, const T& val, const R& result, const Policy&) +{ + typedef typename Policy::indeterminate_result_error_type policy_type; + return detail::raise_indeterminate_result_error( + function, message ? message : "Indeterminate result with value %1%", + val, result, policy_type()); +} + +// +// checked_narrowing_cast: +// +namespace detail +{ + +template +inline bool check_overflow(T val, R* result, const char* function, const Policy& pol) +{ + BOOST_MATH_STD_USING + if(fabs(val) > tools::max_value()) + { + *result = static_cast(boost::math::policies::detail::raise_overflow_error(function, 0, pol)); + return true; + } + return false; +} +template +inline bool check_overflow(std::complex val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_overflow(val.real(), &re, function, pol) || check_overflow(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} +template +inline bool check_underflow(T val, R* result, const char* function, const Policy& pol) +{ + if((val != 0) && (static_cast(val) == 0)) + { + *result = static_cast(boost::math::policies::detail::raise_underflow_error(function, 0, pol)); + return true; + } + return false; +} +template +inline bool check_underflow(std::complex val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_underflow(val.real(), &re, function, pol) || check_underflow(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} +template +inline bool check_denorm(T val, R* result, const char* function, const Policy& pol) +{ + BOOST_MATH_STD_USING + if((fabs(val) < static_cast(tools::min_value())) && (static_cast(val) != 0)) + { + *result = static_cast(boost::math::policies::detail::raise_denorm_error(function, 0, static_cast(val), pol)); + return true; + } + return false; +} +template +inline bool check_denorm(std::complex val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_denorm(val.real(), &re, function, pol) || check_denorm(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} + +// Default instantiations with ignore_error policy. +template +inline bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error&){ return false; } +template +inline bool check_overflow(std::complex /* val */, R* /* result */, const char* /* function */, const overflow_error&){ return false; } +template +inline bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error&){ return false; } +template +inline bool check_underflow(std::complex /* val */, R* /* result */, const char* /* function */, const underflow_error&){ return false; } +template +inline bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error&){ return false; } +template +inline bool check_denorm(std::complex /* val */, R* /* result*/, const char* /* function */, const denorm_error&){ return false; } + +} // namespace detail + +template +inline R checked_narrowing_cast(T val, const char* function) +{ + typedef typename Policy::overflow_error_type overflow_type; + typedef typename Policy::underflow_error_type underflow_type; + typedef typename Policy::denorm_error_type denorm_type; + // + // Most of what follows will evaluate to a no-op: + // + R result = 0; + if(detail::check_overflow(val, &result, function, overflow_type())) + return result; + if(detail::check_underflow(val, &result, function, underflow_type())) + return result; + if(detail::check_denorm(val, &result, function, denorm_type())) + return result; + + return static_cast(val); +} + +template +inline void check_series_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol) +{ + if(max_iter >= policies::get_max_series_iterations()) + raise_evaluation_error( + function, + "Series evaluation exceeded %1% iterations, giving up now.", static_cast(static_cast(max_iter)), pol); +} + +template +inline void check_root_iterations(const char* function, boost::uintmax_t max_iter, const Policy& pol) +{ + if(max_iter >= policies::get_max_root_iterations()) + raise_evaluation_error( + function, + "Root finding evaluation exceeded %1% iterations, giving up now.", static_cast(static_cast(max_iter)), pol); +} + +} //namespace policies + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +}} // namespaces boost/math + +#endif // BOOST_MATH_POLICY_ERROR_HANDLING_HPP + diff --git a/patched/math/special_functions/erf.hpp b/patched/math/special_functions/erf.hpp new file mode 100644 index 0000000..cb2c763 --- /dev/null +++ b/patched/math/special_functions/erf.hpp @@ -0,0 +1,1151 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SPECIAL_ERF_HPP +#define BOOST_MATH_SPECIAL_ERF_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#include +#include +#include +#include +#include + +#include + +namespace boost{ namespace math{ + +namespace detail +{ + +// +// Asymptotic series for large z: +// +template +struct erf_asympt_series_t +{ + erf_asympt_series_t(T z) : xx(2 * -z * z), tk(1) + { + BOOST_MATH_STD_USING + result = -exp(-z * z) / sqrt(boost::math::constants::pi()); + result /= z; + } + + typedef T result_type; + + T operator()() + { + BOOST_MATH_STD_USING + T r = result; + result *= tk / xx; + tk += 2; + if( fabs(r) < fabs(result)) + result = 0; + return r; + } +private: + T result; + T xx; + int tk; +}; +// +// How large z has to be in order to ensure that the series converges: +// +template +inline float erf_asymptotic_limit_N(const T&) +{ + return (std::numeric_limits::max)(); +} +inline float erf_asymptotic_limit_N(const mpl::int_<24>&) +{ + return 2.8F; +} +inline float erf_asymptotic_limit_N(const mpl::int_<53>&) +{ + return 4.3F; +} +inline float erf_asymptotic_limit_N(const mpl::int_<64>&) +{ + return 4.8F; +} +inline float erf_asymptotic_limit_N(const mpl::int_<106>&) +{ + return 6.5F; +} +inline float erf_asymptotic_limit_N(const mpl::int_<113>&) +{ + return 6.8F; +} + +template +inline T erf_asymptotic_limit() +{ + typedef typename policies::precision::type precision_type; + typedef typename mpl::if_< + mpl::less_equal >, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<0>, + mpl::int_<24> + >::type, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<53>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<106>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<113>, + mpl::int_<0> + >::type + >::type + >::type + >::type + >::type tag_type; + return erf_asymptotic_limit_N(tag_type()); +} + +template +T erf_imp(T z, bool invert, const Policy& pol, const Tag& t) +{ + BOOST_MATH_STD_USING + + BOOST_MATH_INSTRUMENT_CODE("Generic erf_imp called"); + + if(z < 0) + { + if(!invert) + return -erf_imp(T(-z), invert, pol, t); + else + return 1 + erf_imp(T(-z), false, pol, t); + } + + T result; + + if(!invert && (z > detail::erf_asymptotic_limit())) + { + detail::erf_asympt_series_t s(z); + boost::uintmax_t max_iter = policies::get_max_series_iterations(); + result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter, 1); + policies::check_series_iterations("boost::math::erf<%1%>(%1%, %1%)", max_iter, pol); + } + else + { + T x = z * z; + if(x < 0.6) + { + // Compute P: + result = z * exp(-x); + result /= sqrt(boost::math::constants::pi()); + if(result != 0) + result *= 2 * detail::lower_gamma_series(T(0.5f), x, pol); + } + else if(x < 1.1f) + { + // Compute Q: + invert = !invert; + result = tgamma_small_upper_part(T(0.5f), x, pol); + result /= sqrt(boost::math::constants::pi()); + } + else + { + // Compute Q: + invert = !invert; + result = z * exp(-x); + result /= sqrt(boost::math::constants::pi()); + result *= upper_gamma_fraction(T(0.5f), x, policies::get_epsilon()); + } + } + if(invert) + result = 1 - result; + return result; +} + +template +T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<53>& t) +{ + BOOST_MATH_STD_USING + + BOOST_MATH_INSTRUMENT_CODE("53-bit precision erf_imp called"); + + if(z < 0) + { + if(!invert) + return -erf_imp(T(-z), invert, pol, t); + else if(z < -0.5) + return 2 - erf_imp(T(-z), invert, pol, t); + else + return 1 + erf_imp(T(-z), false, pol, t); + } + + T result; + + // + // Big bunch of selection statements now to pick + // which implementation to use, + // try to put most likely options first: + // + if(z < 0.5) + { + // + // We're going to calculate erf: + // + if(z < 1e-10) + { + if(z == 0) + { + result = T(0); + } + else + { + static const T c = BOOST_MATH_BIG_CONSTANT(T, 53, 0.003379167095512573896158903121545171688); + result = static_cast(z * 1.125f + z * c); + } + } + else + { + // Maximum Deviation Found: 1.561e-17 + // Expected Error Term: 1.561e-17 + // Maximum Relative Change in Control Points: 1.155e-04 + // Max Error found at double precision = 2.961182e-17 + + static const T Y = 1.044948577880859375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0834305892146531832907), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.338165134459360935041), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.0509990735146777432841), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.00772758345802133288487), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.000322780120964605683831), + }; + static const T Q[] = { + 1L, + BOOST_MATH_BIG_CONSTANT(T, 53, 0.455004033050794024546), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0875222600142252549554), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00858571925074406212772), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.000370900071787748000569), + }; + T zz = z * z; + result = z * (Y + tools::evaluate_polynomial(P, zz) / tools::evaluate_polynomial(Q, zz)); + } + } + else if(invert ? (z < 28) : (z < 5.8f)) + { + // + // We'll be calculating erfc: + // + invert = !invert; + if(z < 1.5f) + { + // Maximum Deviation Found: 3.702e-17 + // Expected Error Term: 3.702e-17 + // Maximum Relative Change in Control Points: 2.845e-04 + // Max Error found at double precision = 4.841816e-17 + static const T Y = 0.405935764312744140625f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, -0.098090592216281240205), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.178114665841120341155), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.191003695796775433986), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0888900368967884466578), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0195049001251218801359), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00180424538297014223957), + }; + static const T Q[] = { + 1L, + BOOST_MATH_BIG_CONSTANT(T, 53, 1.84759070983002217845), + BOOST_MATH_BIG_CONSTANT(T, 53, 1.42628004845511324508), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.578052804889902404909), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.12385097467900864233), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0113385233577001411017), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.337511472483094676155e-5), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 0.5)) / tools::evaluate_polynomial(Q, T(z - 0.5)); + result *= exp(-z * z) / z; + } + else if(z < 2.5f) + { + // Max Error found at double precision = 6.599585e-18 + // Maximum Deviation Found: 3.909e-18 + // Expected Error Term: 3.909e-18 + // Maximum Relative Change in Control Points: 9.886e-05 + static const T Y = 0.50672817230224609375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, -0.0243500476207698441272), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0386540375035707201728), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.04394818964209516296), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175679436311802092299), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00323962406290842133584), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.000235839115596880717416), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 1), + BOOST_MATH_BIG_CONSTANT(T, 53, 1.53991494948552447182), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.982403709157920235114), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.325732924782444448493), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0563921837420478160373), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00410369723978904575884), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 1.5)) / tools::evaluate_polynomial(Q, T(z - 1.5)); + result *= exp(-z * z) / z; + } + else if(z < 4.5f) + { + // Maximum Deviation Found: 1.512e-17 + // Expected Error Term: 1.512e-17 + // Maximum Relative Change in Control Points: 2.222e-04 + // Max Error found at double precision = 2.062515e-17 + static const T Y = 0.5405750274658203125f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00295276716530971662634), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0137384425896355332126), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00840807615555585383007), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00212825620914618649141), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.000250269961544794627958), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.113212406648847561139e-4), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 1), + BOOST_MATH_BIG_CONSTANT(T, 53, 1.04217814166938418171), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.442597659481563127003), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0958492726301061423444), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0105982906484876531489), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.000479411269521714493907), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 3.5)) / tools::evaluate_polynomial(Q, T(z - 3.5)); + result *= exp(-z * z) / z; + } + else + { + // Max Error found at double precision = 2.997958e-17 + // Maximum Deviation Found: 2.860e-17 + // Expected Error Term: 2.859e-17 + // Maximum Relative Change in Control Points: 1.357e-05 + static const T Y = 0.5579090118408203125f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 0.00628057170626964891937), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.0175389834052493308818), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.212652252872804219852), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.687717681153649930619), + BOOST_MATH_BIG_CONSTANT(T, 53, -2.5518551727311523996), + BOOST_MATH_BIG_CONSTANT(T, 53, -3.22729451764143718517), + BOOST_MATH_BIG_CONSTANT(T, 53, -2.8175401114513378771), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 1), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.79257750980575282228), + BOOST_MATH_BIG_CONSTANT(T, 53, 11.0567237927800161565), + BOOST_MATH_BIG_CONSTANT(T, 53, 15.930646027911794143), + BOOST_MATH_BIG_CONSTANT(T, 53, 22.9367376522880577224), + BOOST_MATH_BIG_CONSTANT(T, 53, 13.5064170191802889145), + BOOST_MATH_BIG_CONSTANT(T, 53, 5.48409182238641741584), + }; + result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); + result *= exp(-z * z) / z; + } + } + else + { + // + // Any value of z larger than 28 will underflow to zero: + // + result = 0; + invert = !invert; + } + + if(invert) + { + result = 1 - result; + } + + return result; +} // template T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<53>& t) + + +template +T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<64>& t) +{ + BOOST_MATH_STD_USING + + BOOST_MATH_INSTRUMENT_CODE("64-bit precision erf_imp called"); + + if(z < 0) + { + if(!invert) + return -erf_imp(T(-z), invert, pol, t); + else if(z < -0.5) + return 2 - erf_imp(T(-z), invert, pol, t); + else + return 1 + erf_imp(T(-z), false, pol, t); + } + + T result; + + // + // Big bunch of selection statements now to pick which + // implementation to use, try to put most likely options + // first: + // + if(z < 0.5) + { + // + // We're going to calculate erf: + // + if(z == 0) + { + result = 0; + } + else if(z < 1e-10) + { + static const T c = BOOST_MATH_BIG_CONSTANT(T, 64, 0.003379167095512573896158903121545171688); + result = z * 1.125 + z * c; + } + else + { + // Max Error found at long double precision = 1.623299e-20 + // Maximum Deviation Found: 4.326e-22 + // Expected Error Term: -4.326e-22 + // Maximum Relative Change in Control Points: 1.474e-04 + static const T Y = 1.044948577880859375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0834305892146531988966), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.338097283075565413695), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.0509602734406067204596), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.00904906346158537794396), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.000489468651464798669181), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.200305626366151877759e-4), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.455817300515875172439), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0916537354356241792007), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0102722652675910031202), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.000650511752687851548735), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.189532519105655496778e-4), + }; + result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z))); + } + } + else if(invert ? (z < 110) : (z < 6.4f)) + { + // + // We'll be calculating erfc: + // + invert = !invert; + if(z < 1.5) + { + // Max Error found at long double precision = 3.239590e-20 + // Maximum Deviation Found: 2.241e-20 + // Expected Error Term: -2.241e-20 + // Maximum Relative Change in Control Points: 5.110e-03 + static const T Y = 0.405935764312744140625f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -0.0980905922162812031672), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.159989089922969141329), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.222359821619935712378), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.127303921703577362312), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0384057530342762400273), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00628431160851156719325), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.000441266654514391746428), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.266689068336295642561e-7), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.03237474985469469291), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.78355454954969405222), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.867940326293760578231), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.248025606990021698392), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0396649631833002269861), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00279220237309449026796), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f)); + result *= exp(-z * z) / z; + } + else if(z < 2.5) + { + // Max Error found at long double precision = 3.686211e-21 + // Maximum Deviation Found: 1.495e-21 + // Expected Error Term: -1.494e-21 + // Maximum Relative Change in Control Points: 1.793e-04 + static const T Y = 0.50672817230224609375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -0.024350047620769840217), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0343522687935671451309), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0505420824305544949541), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0257479325917757388209), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00669349844190354356118), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00090807914416099524444), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.515917266698050027934e-4), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.71657861671930336344), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.26409634824280366218), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.512371437838969015941), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.120902623051120950935), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0158027197831887485261), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.000897871370778031611439), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f)); + result *= exp(-z * z) / z; + } + else if(z < 4.5) + { + // Maximum Deviation Found: 1.107e-20 + // Expected Error Term: -1.106e-20 + // Maximum Relative Change in Control Points: 1.709e-04 + // Max Error found at long double precision = 1.446908e-20 + static const T Y = 0.5405750274658203125f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0029527671653097284033), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0141853245895495604051), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0104959584626432293901), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00343963795976100077626), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00059065441194877637899), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.523435380636174008685e-4), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.189896043050331257262e-5), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.19352160185285642574), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.603256964363454392857), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.165411142458540585835), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0259729870946203166468), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00221657568292893699158), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.804149464190309799804e-4), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 3.5f)) / tools::evaluate_polynomial(Q, T(z - 3.5f)); + result *= exp(-z * z) / z; + } + else + { + // Max Error found at long double precision = 7.961166e-21 + // Maximum Deviation Found: 6.677e-21 + // Expected Error Term: 6.676e-21 + // Maximum Relative Change in Control Points: 2.319e-05 + static const T Y = 0.55825519561767578125f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 0.00593438793008050214106), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280666231009089713937), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.141597835204583050043), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.978088201154300548842), + BOOST_MATH_BIG_CONSTANT(T, 64, -5.47351527796012049443), + BOOST_MATH_BIG_CONSTANT(T, 64, -13.8677304660245326627), + BOOST_MATH_BIG_CONSTANT(T, 64, -27.1274948720539821722), + BOOST_MATH_BIG_CONSTANT(T, 64, -29.2545152747009461519), + BOOST_MATH_BIG_CONSTANT(T, 64, -16.8865774499799676937), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.72948911186645394541), + BOOST_MATH_BIG_CONSTANT(T, 64, 23.6750543147695749212), + BOOST_MATH_BIG_CONSTANT(T, 64, 60.0021517335693186785), + BOOST_MATH_BIG_CONSTANT(T, 64, 131.766251645149522868), + BOOST_MATH_BIG_CONSTANT(T, 64, 178.167924971283482513), + BOOST_MATH_BIG_CONSTANT(T, 64, 182.499390505915222699), + BOOST_MATH_BIG_CONSTANT(T, 64, 104.365251479578577989), + BOOST_MATH_BIG_CONSTANT(T, 64, 30.8365511891224291717), + }; + result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); + result *= exp(-z * z) / z; + } + } + else + { + // + // Any value of z larger than 110 will underflow to zero: + // + result = 0; + invert = !invert; + } + + if(invert) + { + result = 1 - result; + } + + return result; +} // template T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<64>& t) + + +template +T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t) +{ + BOOST_MATH_STD_USING + + BOOST_MATH_INSTRUMENT_CODE("113-bit precision erf_imp called"); + + if(z < 0) + { + if(!invert) + return -erf_imp(T(-z), invert, pol, t); + else if(z < -0.5) + return 2 - erf_imp(T(-z), invert, pol, t); + else + return 1 + erf_imp(T(-z), false, pol, t); + } + + T result; + + // + // Big bunch of selection statements now to pick which + // implementation to use, try to put most likely options + // first: + // + if(z < 0.5) + { + // + // We're going to calculate erf: + // + if(z == 0) + { + result = 0; + } + else if(z < 1e-20) + { + static const T c = BOOST_MATH_BIG_CONSTANT(T, 113, 0.003379167095512573896158903121545171688); + result = z * 1.125 + z * c; + } + else + { + // Max Error found at long double precision = 2.342380e-35 + // Maximum Deviation Found: 6.124e-36 + // Expected Error Term: -6.124e-36 + // Maximum Relative Change in Control Points: 3.492e-10 + static const T Y = 1.0841522216796875f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0442269454158250738961589031215451778), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.35549265736002144875335323556961233), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0582179564566667896225454670863270393), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0112694696904802304229950538453123925), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.000805730648981801146251825329609079099), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.566304966591936566229702842075966273e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.169655010425186987820201021510002265e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.344448249920445916714548295433198544e-7), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.466542092785657604666906909196052522), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.100005087012526447295176964142107611), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0128341535890117646540050072234142603), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00107150448466867929159660677016658186), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.586168368028999183607733369248338474e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.196230608502104324965623171516808796e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.313388521582925207734229967907890146e-7), + }; + result = z * (Y + tools::evaluate_polynomial(P, T(z * z)) / tools::evaluate_polynomial(Q, T(z * z))); + } + } + else if(invert ? (z < 110) : (z < 8.65f)) + { + // + // We'll be calculating erfc: + // + invert = !invert; + if(z < 1) + { + // Max Error found at long double precision = 3.246278e-35 + // Maximum Deviation Found: 1.388e-35 + // Expected Error Term: 1.387e-35 + // Maximum Relative Change in Control Points: 6.127e-05 + static const T Y = 0.371877193450927734375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0640320213544647969396032886581290455), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.200769874440155895637857443946706731), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.378447199873537170666487408805779826), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.30521399466465939450398642044975127), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.146890026406815277906781824723458196), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0464837937749539978247589252732769567), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00987895759019540115099100165904822903), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00137507575429025512038051025154301132), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0001144764551085935580772512359680516), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.436544865032836914773944382339900079e-5), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.47651182872457465043733800302427977), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.78706486002517996428836400245547955), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.87295924621659627926365005293130693), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.829375825174365625428280908787261065), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.251334771307848291593780143950311514), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0522110268876176186719436765734722473), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00718332151250963182233267040106902368), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000595279058621482041084986219276392459), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.226988669466501655990637599399326874e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.270666232259029102353426738909226413e-10), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 0.5f)) / tools::evaluate_polynomial(Q, T(z - 0.5f)); + result *= exp(-z * z) / z; + } + else if(z < 1.5) + { + // Max Error found at long double precision = 2.215785e-35 + // Maximum Deviation Found: 1.539e-35 + // Expected Error Term: 1.538e-35 + // Maximum Relative Change in Control Points: 6.104e-05 + static const T Y = 0.45658016204833984375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0289965858925328393392496555094848345), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0868181194868601184627743162571779226), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.169373435121178901746317404936356745), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.13350446515949251201104889028133486), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0617447837290183627136837688446313313), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0185618495228251406703152962489700468), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00371949406491883508764162050169531013), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000485121708792921297742105775823900772), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.376494706741453489892108068231400061e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.133166058052466262415271732172490045e-5), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.32970330146503867261275580968135126), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.46325715420422771961250513514928746), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.55307882560757679068505047390857842), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.644274289865972449441174485441409076), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.182609091063258208068606847453955649), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0354171651271241474946129665801606795), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00454060370165285246451879969534083997), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000349871943711566546821198612518656486), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.123749319840299552925421880481085392e-4), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 1.0f)) / tools::evaluate_polynomial(Q, T(z - 1.0f)); + result *= exp(-z * z) / z; + } + else if(z < 2.25) + { + // Maximum Deviation Found: 1.418e-35 + // Expected Error Term: 1.418e-35 + // Maximum Relative Change in Control Points: 1.316e-04 + // Max Error found at long double precision = 1.998462e-35 + static const T Y = 0.50250148773193359375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0201233630504573402185161184151016606), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0331864357574860196516686996302305002), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0716562720864787193337475444413405461), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0545835322082103985114927569724880658), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0236692635189696678976549720784989593), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00656970902163248872837262539337601845), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00120282643299089441390490459256235021), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000142123229065182650020762792081622986), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.991531438367015135346716277792989347e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.312857043762117596999398067153076051e-6), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.13506082409097783827103424943508554), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.06399257267556230937723190496806215), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.18678481279932541314830499880691109), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.447733186643051752513538142316799562), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.11505680005657879437196953047542148), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.020163993632192726170219663831914034), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00232708971840141388847728782209730585), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000160733201627963528519726484608224112), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.507158721790721802724402992033269266e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.18647774409821470950544212696270639e-12), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 1.5f)) / tools::evaluate_polynomial(Q, T(z - 1.5f)); + result *= exp(-z * z) / z; + } + else if (z < 3) + { + // Maximum Deviation Found: 3.575e-36 + // Expected Error Term: 3.575e-36 + // Maximum Relative Change in Control Points: 7.103e-05 + // Max Error found at long double precision = 5.794737e-36 + static const T Y = 0.52896785736083984375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.00902152521745813634562524098263360074), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0145207142776691539346923710537580927), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0301681239582193983824211995978678571), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0215548540823305814379020678660434461), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00864683476267958365678294164340749949), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00219693096885585491739823283511049902), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000364961639163319762492184502159894371), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.388174251026723752769264051548703059e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.241918026931789436000532513553594321e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.676586625472423508158937481943649258e-7), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.93669171363907292305550231764920001), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.69468476144051356810672506101377494), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.880023580986436640372794392579985511), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.299099106711315090710836273697708402), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0690593962363545715997445583603382337), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0108427016361318921960863149875360222), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00111747247208044534520499324234317695), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.686843205749767250666787987163701209e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.192093541425429248675532015101904262e-5), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 2.25f)) / tools::evaluate_polynomial(Q, T(z - 2.25f)); + result *= exp(-z * z) / z; + } + else if(z < 3.5) + { + // Maximum Deviation Found: 8.126e-37 + // Expected Error Term: -8.126e-37 + // Maximum Relative Change in Control Points: 1.363e-04 + // Max Error found at long double precision = 1.747062e-36 + static const T Y = 0.54037380218505859375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0033703486408887424921155540591370375), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0104948043110005245215286678898115811), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0148530118504000311502310457390417795), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00816693029245443090102738825536188916), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00249716579989140882491939681805594585), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0004655591010047353023978045800916647), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.531129557920045295895085236636025323e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.343526765122727069515775194111741049e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.971120407556888763695313774578711839e-7), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.59911256167540354915906501335919317), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.136006830764025173864831382946934), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.468565867990030871678574840738423023), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.122821824954470343413956476900662236), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0209670914950115943338996513330141633), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00227845718243186165620199012883547257), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000144243326443913171313947613547085553), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.407763415954267700941230249989140046e-5), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 3.0f)) / tools::evaluate_polynomial(Q, T(z - 3.0f)); + result *= exp(-z * z) / z; + } + else if(z < 5.5) + { + // Maximum Deviation Found: 5.804e-36 + // Expected Error Term: -5.803e-36 + // Maximum Relative Change in Control Points: 2.475e-05 + // Max Error found at long double precision = 1.349545e-35 + static const T Y = 0.55000019073486328125f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00118142849742309772151454518093813615), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0072201822885703318172366893469382745), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0078782276276860110721875733778481505), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00418229166204362376187593976656261146), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00134198400587769200074194304298642705), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000283210387078004063264777611497435572), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.405687064094911866569295610914844928e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.39348283801568113807887364414008292e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.248798540917787001526976889284624449e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.929502490223452372919607105387474751e-8), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.156161469668275442569286723236274457e-9), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.52955245103668419479878456656709381), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.06263944820093830054635017117417064), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.441684612681607364321013134378316463), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.121665258426166960049773715928906382), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0232134512374747691424978642874321434), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00310778180686296328582860464875562636), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000288361770756174705123674838640161693), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.177529187194133944622193191942300132e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.655068544833064069223029299070876623e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.11005507545746069573608988651927452e-7), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 4.5f)) / tools::evaluate_polynomial(Q, T(z - 4.5f)); + result *= exp(-z * z) / z; + } + else if(z < 7.5) + { + // Maximum Deviation Found: 1.007e-36 + // Expected Error Term: 1.007e-36 + // Maximum Relative Change in Control Points: 1.027e-03 + // Max Error found at long double precision = 2.646420e-36 + static const T Y = 0.5574436187744140625f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000293236907400849056269309713064107674), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00225110719535060642692275221961480162), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00190984458121502831421717207849429799), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000747757733460111743833929141001680706), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000170663175280949889583158597373928096), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.246441188958013822253071608197514058e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.229818000860544644974205957895688106e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.134886977703388748488480980637704864e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.454764611880548962757125070106650958e-8), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.673002744115866600294723141176820155e-10), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.12843690320861239631195353379313367), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.569900657061622955362493442186537259), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.169094404206844928112348730277514273), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0324887449084220415058158657252147063), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00419252877436825753042680842608219552), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00036344133176118603523976748563178578), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.204123895931375107397698245752850347e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.674128352521481412232785122943508729e-6), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.997637501418963696542159244436245077e-8), + }; + result = Y + tools::evaluate_polynomial(P, T(z - 6.5f)) / tools::evaluate_polynomial(Q, T(z - 6.5f)); + result *= exp(-z * z) / z; + } + else if(z < 11.5) + { + // Maximum Deviation Found: 8.380e-36 + // Expected Error Term: 8.380e-36 + // Maximum Relative Change in Control Points: 2.632e-06 + // Max Error found at long double precision = 9.849522e-36 + static const T Y = 0.56083202362060546875f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000282420728751494363613829834891390121), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00175387065018002823433704079355125161), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0021344978564889819420775336322920375), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00124151356560137532655039683963075661), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000423600733566948018555157026862139644), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.914030340865175237133613697319509698e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.126999927156823363353809747017945494e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.110610959842869849776179749369376402e-5), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.55075079477173482096725348704634529e-7), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.119735694018906705225870691331543806e-8), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.69889613396167354566098060039549882), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.28824647372749624464956031163282674), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.572297795434934493541628008224078717), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.164157697425571712377043857240773164), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0315311145224594430281219516531649562), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00405588922155632380812945849777127458), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000336929033691445666232029762868642417), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.164033049810404773469413526427932109e-4), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.356615210500531410114914617294694857e-6), + }; + result = Y + tools::evaluate_polynomial(P, T(z / 2 - 4.75f)) / tools::evaluate_polynomial(Q, T(z / 2 - 4.75f)); + result *= exp(-z * z) / z; + } + else + { + // Maximum Deviation Found: 1.132e-35 + // Expected Error Term: -1.132e-35 + // Maximum Relative Change in Control Points: 4.674e-04 + // Max Error found at long double precision = 1.162590e-35 + static const T Y = 0.5632686614990234375f; + static const T P[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000920922048732849448079451574171836943), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00321439044532288750501700028748922439), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.250455263029390118657884864261823431), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.906807635364090342031792404764598142), + BOOST_MATH_BIG_CONSTANT(T, 113, -8.92233572835991735876688745989985565), + BOOST_MATH_BIG_CONSTANT(T, 113, -21.7797433494422564811782116907878495), + BOOST_MATH_BIG_CONSTANT(T, 113, -91.1451915251976354349734589601171659), + BOOST_MATH_BIG_CONSTANT(T, 113, -144.1279109655993927069052125017673), + BOOST_MATH_BIG_CONSTANT(T, 113, -313.845076581796338665519022313775589), + BOOST_MATH_BIG_CONSTANT(T, 113, -273.11378811923343424081101235736475), + BOOST_MATH_BIG_CONSTANT(T, 113, -271.651566205951067025696102600443452), + BOOST_MATH_BIG_CONSTANT(T, 113, -60.0530577077238079968843307523245547), + }; + static const T Q[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.49040448075464744191022350947892036), + BOOST_MATH_BIG_CONSTANT(T, 113, 34.3563592467165971295915749548313227), + BOOST_MATH_BIG_CONSTANT(T, 113, 84.4993232033879023178285731843850461), + BOOST_MATH_BIG_CONSTANT(T, 113, 376.005865281206894120659401340373818), + BOOST_MATH_BIG_CONSTANT(T, 113, 629.95369438888946233003926191755125), + BOOST_MATH_BIG_CONSTANT(T, 113, 1568.35771983533158591604513304269098), + BOOST_MATH_BIG_CONSTANT(T, 113, 1646.02452040831961063640827116581021), + BOOST_MATH_BIG_CONSTANT(T, 113, 2299.96860633240298708910425594484895), + BOOST_MATH_BIG_CONSTANT(T, 113, 1222.73204392037452750381340219906374), + BOOST_MATH_BIG_CONSTANT(T, 113, 799.359797306084372350264298361110448), + BOOST_MATH_BIG_CONSTANT(T, 113, 72.7415265778588087243442792401576737), + }; + result = Y + tools::evaluate_polynomial(P, T(1 / z)) / tools::evaluate_polynomial(Q, T(1 / z)); + result *= exp(-z * z) / z; + } + } + else + { + // + // Any value of z larger than 110 will underflow to zero: + // + result = 0; + invert = !invert; + } + + if(invert) + { + result = 1 - result; + } + + return result; +} // template T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<113>& t) + +template +struct erf_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + static void do_init(const mpl::int_<0>&){} + static void do_init(const mpl::int_<53>&) + { + boost::math::erf(static_cast(1e-12), Policy()); + boost::math::erf(static_cast(0.25), Policy()); + boost::math::erf(static_cast(1.25), Policy()); + boost::math::erf(static_cast(2.25), Policy()); + boost::math::erf(static_cast(4.25), Policy()); + boost::math::erf(static_cast(5.25), Policy()); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::erf(static_cast(1e-12), Policy()); + boost::math::erf(static_cast(0.25), Policy()); + boost::math::erf(static_cast(1.25), Policy()); + boost::math::erf(static_cast(2.25), Policy()); + boost::math::erf(static_cast(4.25), Policy()); + boost::math::erf(static_cast(5.25), Policy()); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::erf(static_cast(1e-22), Policy()); + boost::math::erf(static_cast(0.25), Policy()); + boost::math::erf(static_cast(1.25), Policy()); + boost::math::erf(static_cast(2.125), Policy()); + boost::math::erf(static_cast(2.75), Policy()); + boost::math::erf(static_cast(3.25), Policy()); + boost::math::erf(static_cast(5.25), Policy()); + boost::math::erf(static_cast(7.25), Policy()); + boost::math::erf(static_cast(11.25), Policy()); + boost::math::erf(static_cast(12.5), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template +const typename erf_initializer::init erf_initializer::initializer; + +} // namespace detail + +template +inline typename tools::promote_args::type erf(T z, const Policy& /* pol */) +{ + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::precision::type precision_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + BOOST_MATH_INSTRUMENT_CODE("result_type = " << type_id().name_demangled()); + BOOST_MATH_INSTRUMENT_CODE("value_type = " << type_id().name_demangled()); + BOOST_MATH_INSTRUMENT_CODE("precision_type = " << type_id().name_demangled()); + + typedef typename mpl::if_< + mpl::less_equal >, + mpl::int_<0>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<53>, // double + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, // 80-bit long double + typename mpl::if_< + mpl::less_equal >, + mpl::int_<113>, // 128-bit long double + mpl::int_<0> // too many bits, use generic version. + >::type + >::type + >::type + >::type tag_type; + + BOOST_MATH_INSTRUMENT_CODE("tag_type = " << type_id().name_demangled()); + + detail::erf_initializer::force_instantiate(); // Force constants to be initialized before main + + return policies::checked_narrowing_cast(detail::erf_imp( + static_cast(z), + false, + forwarding_policy(), + tag_type()), "boost::math::erf<%1%>(%1%, %1%)"); +} + +template +inline typename tools::promote_args::type erfc(T z, const Policy& /* pol */) +{ + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::precision::type precision_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + BOOST_MATH_INSTRUMENT_CODE("result_type = " << type_id().name_demangled()); + BOOST_MATH_INSTRUMENT_CODE("value_type = " << type_id().name_demangled()); + BOOST_MATH_INSTRUMENT_CODE("precision_type = " << type_id().name_demangled()); + + typedef typename mpl::if_< + mpl::less_equal >, + mpl::int_<0>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<53>, // double + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, // 80-bit long double + typename mpl::if_< + mpl::less_equal >, + mpl::int_<113>, // 128-bit long double + mpl::int_<0> // too many bits, use generic version. + >::type + >::type + >::type + >::type tag_type; + + BOOST_MATH_INSTRUMENT_CODE("tag_type = " << type_id().name_demangled()); + + detail::erf_initializer::force_instantiate(); // Force constants to be initialized before main + + return policies::checked_narrowing_cast(detail::erf_imp( + static_cast(z), + true, + forwarding_policy(), + tag_type()), "boost::math::erfc<%1%>(%1%, %1%)"); +} + +template +inline typename tools::promote_args::type erf(T z) +{ + return boost::math::erf(z, policies::policy<>()); +} + +template +inline typename tools::promote_args::type erfc(T z) +{ + return boost::math::erfc(z, policies::policy<>()); +} + +} // namespace math +} // namespace boost + +#include + +#endif // BOOST_MATH_SPECIAL_ERF_HPP + + + + diff --git a/patched/math/special_functions/gamma.hpp b/patched/math/special_functions/gamma.hpp new file mode 100644 index 0000000..5322f18 --- /dev/null +++ b/patched/math/special_functions/gamma.hpp @@ -0,0 +1,1659 @@ + +// Copyright John Maddock 2006-7. +// Copyright Paul A. Bristow 2007. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SF_GAMMA_HPP +#define BOOST_MATH_SF_GAMMA_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4127 4701) +// // For lexical_cast, until fixed in 1.35? +// // conditional expression is constant & +// // Potentially uninitialized local variable 'name' used +#endif +#include +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifdef BOOST_MATH_INSTRUMENT +#include +#include +#include +#endif + +#ifdef BOOST_MSVC +# pragma warning(push) +# pragma warning(disable: 4702) // unreachable code (return after domain_error throw). +# pragma warning(disable: 4127) // conditional expression is constant. +# pragma warning(disable: 4100) // unreferenced formal parameter. +// Several variables made comments, +// but some difficulty as whether referenced on not may depend on macro values. +// So to be safe, 4100 warnings suppressed. +// TODO - revisit this? +#endif + +namespace boost{ namespace math{ + +namespace detail{ + +template +inline bool is_odd(T v, const boost::true_type&) +{ + int i = static_cast(v); + return i&1; +} +template +inline bool is_odd(T v, const boost::false_type&) +{ + // Oh dear can't cast T to int! + BOOST_MATH_STD_USING + T modulus = v - 2 * floor(v/2); + return static_cast(modulus != 0); +} +template +inline bool is_odd(T v) +{ + return is_odd(v, ::boost::is_convertible()); +} + +template +T sinpx(T z) +{ + // Ad hoc function calculates x * sin(pi * x), + // taking extra care near when x is near a whole number. + BOOST_MATH_STD_USING + int sign = 1; + if(z < 0) + { + z = -z; + } + else + { + sign = -sign; + } + T fl = floor(z); + T dist; + if(is_odd(fl)) + { + fl += 1; + dist = fl - z; + sign = -sign; + } + else + { + dist = z - fl; + } + BOOST_ASSERT(fl >= 0); + if(dist > 0.5) + dist = 1 - dist; + T result = sin(dist*boost::math::constants::pi()); + return sign*z*result; +} // template T sinpx(T z) +// +// tgamma(z), with Lanczos support: +// +template +T gamma_imp(T z, const Policy& pol, const Lanczos& l) +{ + BOOST_MATH_STD_USING + + T result = 1; + +#ifdef BOOST_MATH_INSTRUMENT + static bool b = false; + if(!b) + { + std::cout << "tgamma_imp called with " << type_id().name_demangled() << " " << type_id().name_demangled() << std::endl; + b = true; + } +#endif + static const char* function = "boost::math::tgamma<%1%>(%1%)"; + + if(z <= 0) + { + if(floor(z) == z) + return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); + if(z <= -20) + { + result = gamma_imp(T(-z), pol, l) * sinpx(z); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + if((fabs(result) < 1) && (tools::max_value() * fabs(result) < boost::math::constants::pi())) + return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + result = -boost::math::constants::pi() / result; + if(result == 0) + return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); + if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL) + return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", result, pol); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + return result; + } + + // shift z to > 1: + while(z < 0) + { + result /= z; + z += 1; + } + } + BOOST_MATH_INSTRUMENT_VARIABLE(result); + if((floor(z) == z) && (z < max_factorial::value)) + { + result *= unchecked_factorial(itrunc(z, pol) - 1); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + else + { + result *= Lanczos::lanczos_sum(z); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value()); + if(z * log(z) > tools::log_max_value()) + { + // we're going to overflow unless this is done with care: + T zgh = (z + static_cast(Lanczos::g()) - boost::math::constants::half()); + BOOST_MATH_INSTRUMENT_VARIABLE(zgh); + if(log(zgh) * z / 2 > tools::log_max_value()) + return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + T hp = pow(zgh, (z / 2) - T(0.25)); + BOOST_MATH_INSTRUMENT_VARIABLE(hp); + result *= hp / exp(zgh); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + if(tools::max_value() / hp < result) + return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + result *= hp; + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + else + { + T zgh = (z + static_cast(Lanczos::g()) - boost::math::constants::half()); + BOOST_MATH_INSTRUMENT_VARIABLE(zgh); + BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half())); + BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh)); + result *= pow(zgh, z - boost::math::constants::half()) / exp(zgh); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + } + return result; +} +// +// lgamma(z) with Lanczos support: +// +template +T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0) +{ +#ifdef BOOST_MATH_INSTRUMENT + static bool b = false; + if(!b) + { + std::cout << "lgamma_imp called with " << type_id().name_demangled() << " " << type_id().name_demangled() << std::endl; + b = true; + } +#endif + + BOOST_MATH_STD_USING + + static const char* function = "boost::math::lgamma<%1%>(%1%)"; + + T result = 0; + int sresult = 1; + if(z <= 0) + { + // reflection formula: + if(floor(z) == z) + return policies::raise_pole_error(function, "Evaluation of lgamma at a negative integer %1%.", z, pol); + + T t = sinpx(z); + z = -z; + if(t < 0) + { + t = -t; + } + else + { + sresult = -sresult; + } + result = log(boost::math::constants::pi()) - lgamma_imp(z, pol, l) - log(t); + } + else if(z < 15) + { + typedef typename policies::precision::type precision_type; + typedef typename mpl::if_< + mpl::and_< + mpl::less_equal >, + mpl::greater > + >, + mpl::int_<64>, + typename mpl::if_< + mpl::and_< + mpl::less_equal >, + mpl::greater > + >, + mpl::int_<113>, mpl::int_<0> >::type + >::type tag_type; + result = lgamma_small_imp(z, T(z - 1), T(z - 2), tag_type(), pol, l); + } + else if((z >= 3) && (z < 100)) + { + // taking the log of tgamma reduces the error, no danger of overflow here: + result = log(gamma_imp(z, pol, l)); + } + else + { + // regular evaluation: + T zgh = static_cast(z + Lanczos::g() - boost::math::constants::half()); + result = log(zgh) - 1; + result *= z - 0.5f; + result += log(Lanczos::lanczos_sum_expG_scaled(z)); + } + + if(sign) + *sign = sresult; + return result; +} + +// +// Incomplete gamma functions follow: +// +template +struct upper_incomplete_gamma_fract +{ +private: + T z, a; + int k; +public: + typedef std::pair result_type; + + upper_incomplete_gamma_fract(T a1, T z1) + : z(z1-a1+1), a(a1), k(0) + { + } + + result_type operator()() + { + ++k; + z += 2; + return result_type(k * (a - k), z); + } +}; + +template +inline T upper_gamma_fraction(T a, T z, T eps) +{ + // Multiply result by z^a * e^-z to get the full + // upper incomplete integral. Divide by tgamma(z) + // to normalise. + upper_incomplete_gamma_fract f(a, z); + return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)); +} + +template +struct lower_incomplete_gamma_series +{ +private: + T a, z, result; +public: + typedef T result_type; + lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){} + + T operator()() + { + T r = result; + a += 1; + result *= z/a; + return r; + } +}; + +template +inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0) +{ + // Multiply result by ((z^a) * (e^-z) / a) to get the full + // lower incomplete integral. Then divide by tgamma(a) + // to get the normalised value. + lower_incomplete_gamma_series s(a, z); + boost::uintmax_t max_iter = policies::get_max_series_iterations(); + T factor = policies::get_epsilon(); + T result = boost::math::tools::sum_series(s, factor, max_iter, init_value); + policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol); + return result; +} + +// +// Fully generic tgamma and lgamma use the incomplete partial +// sums added together: +// +template +T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l) +{ + static const char* function = "boost::math::tgamma<%1%>(%1%)"; + BOOST_MATH_STD_USING + if((z <= 0) && (floor(z) == z)) + return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); + if(z <= -20) + { + T result = gamma_imp(T(-z), pol, l) * sinpx(z); + if((fabs(result) < 1) && (tools::max_value() * fabs(result) < boost::math::constants::pi())) + return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + result = -boost::math::constants::pi() / result; + if(result == 0) + return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); + if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL) + return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", result, pol); + return result; + } + // + // The upper gamma fraction is *very* slow for z < 6, actually it's very + // slow to converge everywhere but recursing until z > 6 gets rid of the + // worst of it's behaviour. + // + T prefix = 1; + while(z < 6) + { + prefix /= z; + z += 1; + } + BOOST_MATH_INSTRUMENT_CODE(prefix); + if((floor(z) == z) && (z < max_factorial::value)) + { + prefix *= unchecked_factorial(itrunc(z, pol) - 1); + } + else + { + prefix = prefix * pow(z / boost::math::constants::e(), z); + BOOST_MATH_INSTRUMENT_CODE(prefix); + T sum = detail::lower_gamma_series(z, z, pol) / z; + BOOST_MATH_INSTRUMENT_CODE(sum); + sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon()); + BOOST_MATH_INSTRUMENT_CODE(sum); + if(fabs(tools::max_value() / prefix) < fabs(sum)) + return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + BOOST_MATH_INSTRUMENT_CODE((sum * prefix)); + return sum * prefix; + } + return prefix; +} + +template +T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign) +{ + BOOST_MATH_STD_USING + + static const char* function = "boost::math::lgamma<%1%>(%1%)"; + T result = 0; + int sresult = 1; + if(z <= 0) + { + if(floor(z) == z) + return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); + T t = detail::sinpx(z); + z = -z; + if(t < 0) + { + t = -t; + } + else + { + sresult = -sresult; + } + result = log(boost::math::constants::pi()) - lgamma_imp(z, pol, l, 0) - log(t); + } + else if((z != 1) && (z != 2)) + { + T limit = (std::max)(T(z+1), T(10)); + T prefix = z * log(limit) - limit; + T sum = detail::lower_gamma_series(z, limit, pol) / z; + sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon()); + result = log(sum) + prefix; + } + if(sign) + *sign = sresult; + return result; +} +// +// This helper calculates tgamma(dz+1)-1 without cancellation errors, +// used by the upper incomplete gamma with z < 1: +// +template +T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l) +{ + BOOST_MATH_STD_USING + + typedef typename policies::precision::type precision_type; + + typedef typename mpl::if_< + mpl::or_< + mpl::less_equal >, + mpl::greater > + >, + typename mpl::if_< + is_same, + mpl::int_<113>, + mpl::int_<0> + >::type, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, mpl::int_<113> >::type + >::type tag_type; + + T result; + if(dz < 0) + { + if(dz < -0.5) + { + // Best method is simply to subtract 1 from tgamma: + result = boost::math::tgamma(1+dz, pol) - 1; + BOOST_MATH_INSTRUMENT_CODE(result); + } + else + { + // Use expm1 on lgamma: + result = boost::math::expm1(-boost::math::log1p(dz, pol) + + lgamma_small_imp(dz+2, dz + 1, dz, tag_type(), pol, l)); + BOOST_MATH_INSTRUMENT_CODE(result); + } + } + else + { + if(dz < 2) + { + // Use expm1 on lgamma: + result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag_type(), pol, l), pol); + BOOST_MATH_INSTRUMENT_CODE(result); + } + else + { + // Best method is simply to subtract 1 from tgamma: + result = boost::math::tgamma(1+dz, pol) - 1; + BOOST_MATH_INSTRUMENT_CODE(result); + } + } + + return result; +} + +template +inline T tgammap1m1_imp(T dz, Policy const& pol, + const ::boost::math::lanczos::undefined_lanczos& l) +{ + BOOST_MATH_STD_USING // ADL of std names + // + // There should be a better solution than this, but the + // algebra isn't easy for the general case.... + // Start by subracting 1 from tgamma: + // + T result = gamma_imp(T(1 + dz), pol, l) - 1; + BOOST_MATH_INSTRUMENT_CODE(result); + // + // Test the level of cancellation error observed: we loose one bit + // for each power of 2 the result is less than 1. If we would get + // more bits from our most precise lgamma rational approximation, + // then use that instead: + // + BOOST_MATH_INSTRUMENT_CODE((dz > -0.5)); + BOOST_MATH_INSTRUMENT_CODE((dz < 2)); + BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits()) * fabs(result) < 1e34)); + if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits()) * fabs(result) < 1e34)) + { + result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113()); + BOOST_MATH_INSTRUMENT_CODE(result); + } + return result; +} + +// +// Series representation for upper fraction when z is small: +// +template +struct small_gamma2_series +{ + typedef T result_type; + + small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){} + + T operator()() + { + T r = result / (apn); + result *= x; + result /= ++n; + apn += 1; + return r; + } + +private: + T result, x, apn; + int n; +}; +// +// calculate power term prefix (z^a)(e^-z) used in the non-normalised +// incomplete gammas: +// +template +T full_igamma_prefix(T a, T z, const Policy& pol) +{ + BOOST_MATH_STD_USING + + T prefix; + T alz = a * log(z); + + if(z >= 1) + { + if((alz < tools::log_max_value()) && (-z > tools::log_min_value())) + { + prefix = pow(z, a) * exp(-z); + } + else if(a >= 1) + { + prefix = pow(z / exp(z/a), a); + } + else + { + prefix = exp(alz - z); + } + } + else + { + if(alz > tools::log_min_value()) + { + prefix = pow(z, a) * exp(-z); + } + else if(z/a < tools::log_max_value()) + { + prefix = pow(z / exp(z/a), a); + } + else + { + prefix = exp(alz - z); + } + } + // + // This error handling isn't very good: it happens after the fact + // rather than before it... + // + if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE) + policies::raise_overflow_error("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol); + + return prefix; +} +// +// Compute (z^a)(e^-z)/tgamma(a) +// most if the error occurs in this function: +// +template +T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l) +{ + BOOST_MATH_STD_USING + T agh = a + static_cast(Lanczos::g()) - T(0.5); + T prefix; + T d = ((z - a) - static_cast(Lanczos::g()) + T(0.5)) / agh; + + if(a < 1) + { + // + // We have to treat a < 1 as a special case because our Lanczos + // approximations are optimised against the factorials with a > 1, + // and for high precision types especially (128-bit reals for example) + // very small values of a can give rather eroneous results for gamma + // unless we do this: + // + // TODO: is this still required? Lanczos approx should be better now? + // + if(z <= tools::log_min_value()) + { + // Oh dear, have to use logs, should be free of cancellation errors though: + return exp(a * log(z) - z - lgamma_imp(a, pol, l)); + } + else + { + // direct calculation, no danger of overflow as gamma(a) < 1/a + // for small a. + return pow(z, a) * exp(-z) / gamma_imp(a, pol, l); + } + } + else if((fabs(d*d*a) <= 100) && (a > 150)) + { + // special case for large a and a ~ z. + prefix = a * boost::math::log1pmx(d, pol) + z * static_cast(0.5 - Lanczos::g()) / agh; + prefix = exp(prefix); + } + else + { + // + // general case. + // direct computation is most accurate, but use various fallbacks + // for different parts of the problem domain: + // + T alz = a * log(z / agh); + T amz = a - z; + if(((std::min)(alz, amz) <= tools::log_min_value()) || ((std::max)(alz, amz) >= tools::log_max_value())) + { + T amza = amz / a; + if(((std::min)(alz, amz)/2 > tools::log_min_value()) && ((std::max)(alz, amz)/2 < tools::log_max_value())) + { + // compute square root of the result and then square it: + T sq = pow(z / agh, a / 2) * exp(amz / 2); + prefix = sq * sq; + } + else if(((std::min)(alz, amz)/4 > tools::log_min_value()) && ((std::max)(alz, amz)/4 < tools::log_max_value()) && (z > a)) + { + // compute the 4th root of the result then square it twice: + T sq = pow(z / agh, a / 4) * exp(amz / 4); + prefix = sq * sq; + prefix *= prefix; + } + else if((amza > tools::log_min_value()) && (amza < tools::log_max_value())) + { + prefix = pow((z * exp(amza)) / agh, a); + } + else + { + prefix = exp(alz + amz); + } + } + else + { + prefix = pow(z / agh, a) * exp(amz); + } + } + prefix *= sqrt(agh / boost::math::constants::e()) / Lanczos::lanczos_sum_expG_scaled(a); + return prefix; +} +// +// And again, without Lanczos support: +// +template +T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&) +{ + BOOST_MATH_STD_USING + + T limit = (std::max)(T(10), a); + T sum = detail::lower_gamma_series(a, limit, pol) / a; + sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon()); + + if(a < 10) + { + // special case for small a: + T prefix = pow(z / 10, a); + prefix *= exp(10-z); + if(0 == prefix) + { + prefix = pow((z * exp((10-z)/a)) / 10, a); + } + prefix /= sum; + return prefix; + } + + T zoa = z / a; + T amz = a - z; + T alzoa = a * log(zoa); + T prefix; + if(((std::min)(alzoa, amz) <= tools::log_min_value()) || ((std::max)(alzoa, amz) >= tools::log_max_value())) + { + T amza = amz / a; + if((amza <= tools::log_min_value()) || (amza >= tools::log_max_value())) + { + prefix = exp(alzoa + amz); + } + else + { + prefix = pow(zoa * exp(amza), a); + } + } + else + { + prefix = pow(zoa, a) * exp(amz); + } + prefix /= sum; + return prefix; +} +// +// Upper gamma fraction for very small a: +// +template +inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0) +{ + BOOST_MATH_STD_USING // ADL of std functions. + // + // Compute the full upper fraction (Q) when a is very small: + // + T result; + result = boost::math::tgamma1pm1(a, pol); + if(pgam) + *pgam = (result + 1) / a; + T p = boost::math::powm1(x, a, pol); + result -= p; + result /= a; + detail::small_gamma2_series s(a, x); + boost::uintmax_t max_iter = policies::get_max_series_iterations() - 10; + p += 1; + if(pderivative) + *pderivative = p / (*pgam * exp(x)); + T init_value = invert ? *pgam : 0; + result = -p * tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, (init_value - result) / p); + policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol); + if(invert) + result = -result; + return result; +} +// +// Upper gamma fraction for integer a: +// +template +inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0) +{ + // + // Calculates normalised Q when a is an integer: + // + BOOST_MATH_STD_USING + T e = exp(-x); + T sum = e; + if(sum != 0) + { + T term = sum; + for(unsigned n = 1; n < a; ++n) + { + term /= n; + term *= x; + sum += term; + } + } + if(pderivative) + { + *pderivative = e * pow(x, a) / boost::math::unchecked_factorial(itrunc(T(a - 1), pol)); + } + return sum; +} +// +// Upper gamma fraction for half integer a: +// +template +T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol) +{ + // + // Calculates normalised Q when a is a half-integer: + // + BOOST_MATH_STD_USING + T e = boost::math::erfc(sqrt(x), pol); + if((e != 0) && (a > 1)) + { + T term = exp(-x) / sqrt(constants::pi() * x); + term *= x; + static const T half = T(1) / 2; + term /= half; + T sum = term; + for(unsigned n = 2; n < a; ++n) + { + term /= n - half; + term *= x; + sum += term; + } + e += sum; + if(p_derivative) + { + *p_derivative = 0; + } + } + else if(p_derivative) + { + // We'll be dividing by x later, so calculate derivative * x: + *p_derivative = sqrt(x) * exp(-x) / constants::root_pi(); + } + return e; +} +// +// Main incomplete gamma entry point, handles all four incomplete gamma's: +// +template +T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, + const Policy& pol, T* p_derivative) +{ + static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)"; + if(a <= 0) + policies::raise_domain_error(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); + if(x < 0) + policies::raise_domain_error(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); + + BOOST_MATH_STD_USING + + typedef typename lanczos::lanczos::type lanczos_type; + + T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used + + BOOST_ASSERT((p_derivative == 0) || (normalised == true)); + + bool is_int, is_half_int; + bool is_small_a = (a < 30) && (a <= x + 1); + if(is_small_a) + { + T fa = floor(a); + is_int = (fa == a); + is_half_int = is_int ? false : (fabs(fa - a) == 0.5f); + } + else + { + is_int = is_half_int = false; + } + + int eval_method; + + if(is_int && (x > 0.6)) + { + // calculate Q via finite sum: + invert = !invert; + eval_method = 0; + } + else if(is_half_int && (x > 0.2)) + { + // calculate Q via finite sum for half integer a: + invert = !invert; + eval_method = 1; + } + else if(x < 0.5) + { + // + // Changeover criterion chosen to give a changeover at Q ~ 0.33 + // + if(-0.4 / log(x) < a) + { + eval_method = 2; + } + else + { + eval_method = 3; + } + } + else if(x < 1.1) + { + // + // Changover here occurs when P ~ 0.75 or Q ~ 0.25: + // + if(x * 0.75f < a) + { + eval_method = 2; + } + else + { + eval_method = 3; + } + } + else + { + // + // Begin by testing whether we're in the "bad" zone + // where the result will be near 0.5 and the usual + // series and continued fractions are slow to converge: + // + bool use_temme = false; + if(normalised && std::numeric_limits::is_specialized && (a > 20)) + { + T sigma = fabs((x-a)/a); + if((a > 200) && (policies::digits() <= 113)) + { + // + // This limit is chosen so that we use Temme's expansion + // only if the result would be larger than about 10^-6. + // Below that the regular series and continued fractions + // converge OK, and if we use Temme's method we get increasing + // errors from the dominant erfc term as it's (inexact) argument + // increases in magnitude. + // + if(20 / a > sigma * sigma) + use_temme = true; + } + else if(policies::digits() <= 64) + { + // Note in this zone we can't use Temme's expansion for + // types longer than an 80-bit real: + // it would require too many terms in the polynomials. + if(sigma < 0.4) + use_temme = true; + } + } + if(use_temme) + { + eval_method = 5; + } + else + { + // + // Regular case where the result will not be too close to 0.5. + // + // Changeover here occurs at P ~ Q ~ 0.5 + // Note that series computation of P is about x2 faster than continued fraction + // calculation of Q, so try and use the CF only when really necessary, especially + // for small x. + // + if(x - (1 / (3 * x)) < a) + { + eval_method = 2; + } + else + { + eval_method = 4; + invert = !invert; + } + } + } + + switch(eval_method) + { + case 0: + { + result = finite_gamma_q(a, x, pol, p_derivative); + if(normalised == false) + result *= boost::math::tgamma(a, pol); + break; + } + case 1: + { + result = finite_half_gamma_q(a, x, p_derivative, pol); + if(normalised == false) + result *= boost::math::tgamma(a, pol); + if(p_derivative && (*p_derivative == 0)) + *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); + break; + } + case 2: + { + // Compute P: + result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); + if(p_derivative) + *p_derivative = result; + if(result != 0) + { + T init_value = 0; + if(invert) + { + init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result; + } + result *= detail::lower_gamma_series(a, x, pol, init_value) / a; + if(invert) + { + invert = false; + result = -result; + } + } + break; + } + case 3: + { + // Compute Q: + invert = !invert; + T g; + result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative); + invert = false; + if(normalised) + result /= g; + break; + } + case 4: + { + // Compute Q: + result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); + if(p_derivative) + *p_derivative = result; + if(result != 0) + result *= upper_gamma_fraction(a, x, policies::get_epsilon()); + break; + } + case 5: + { + // + // Use compile time dispatch to the appropriate + // Temme asymptotic expansion. This may be dead code + // if T does not have numeric limits support, or has + // too many digits for the most precise version of + // these expansions, in that case we'll be calling + // an empty function. + // + typedef typename policies::precision::type precision_type; + + typedef typename mpl::if_< + mpl::or_ >, + mpl::greater > >, + mpl::int_<0>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<53>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, + mpl::int_<113> + >::type + >::type + >::type tag_type; + + result = igamma_temme_large(a, x, pol, static_cast(0)); + if(x >= a) + invert = !invert; + if(p_derivative) + *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); + break; + } + } + + if(normalised && (result > 1)) + result = 1; + if(invert) + { + T gam = normalised ? 1 : boost::math::tgamma(a, pol); + result = gam - result; + } + if(p_derivative) + { + // + // Need to convert prefix term to derivative: + // + if((x < 1) && (tools::max_value() * x < *p_derivative)) + { + // overflow, just return an arbitrarily large value: + *p_derivative = tools::max_value() / 2; + } + + *p_derivative /= x; + } + + return result; +} + +// +// Ratios of two gamma functions: +// +template +T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&) +{ + BOOST_MATH_STD_USING + T zgh = z + Lanczos::g() - constants::half(); + T result; + if(fabs(delta) < 10) + { + result = exp((constants::half() - z) * boost::math::log1p(delta / zgh, pol)); + } + else + { + result = pow(zgh / (zgh + delta), z - constants::half()); + } + result *= pow(constants::e() / (zgh + delta), delta); + result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta)); + return result; +} +// +// And again without Lanczos support this time: +// +template +T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&) +{ + BOOST_MATH_STD_USING + // + // The upper gamma fraction is *very* slow for z < 6, actually it's very + // slow to converge everywhere but recursing until z > 6 gets rid of the + // worst of it's behaviour. + // + T prefix = 1; + T zd = z + delta; + while((zd < 6) && (z < 6)) + { + prefix /= z; + prefix *= zd; + z += 1; + zd += 1; + } + if(delta < 10) + { + prefix *= exp(-z * boost::math::log1p(delta / z, pol)); + } + else + { + prefix *= pow(z / zd, z); + } + prefix *= pow(constants::e() / zd, delta); + T sum = detail::lower_gamma_series(z, z, pol) / z; + sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon()); + T sumd = detail::lower_gamma_series(zd, zd, pol) / zd; + sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon()); + sum /= sumd; + if(fabs(tools::max_value() / prefix) < fabs(sum)) + return policies::raise_overflow_error("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol); + return sum * prefix; +} + +template +T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol) +{ + BOOST_MATH_STD_USING + + if(z <= 0) + policies::raise_domain_error("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", z, pol); + if(z+delta <= 0) + policies::raise_domain_error("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta, pol); + + if(floor(delta) == delta) + { + if(floor(z) == z) + { + // + // Both z and delta are integers, see if we can just use table lookup + // of the factorials to get the result: + // + if((z <= max_factorial::value) && (z + delta <= max_factorial::value)) + { + return unchecked_factorial((unsigned)itrunc(z, pol) - 1) / unchecked_factorial((unsigned)itrunc(T(z + delta), pol) - 1); + } + } + if(fabs(delta) < 20) + { + // + // delta is a small integer, we can use a finite product: + // + if(delta == 0) + return 1; + if(delta < 0) + { + z -= 1; + T result = z; + while(0 != (delta += 1)) + { + z -= 1; + result *= z; + } + return result; + } + else + { + T result = 1 / z; + while(0 != (delta -= 1)) + { + z += 1; + result /= z; + } + return result; + } + } + } + typedef typename lanczos::lanczos::type lanczos_type; + return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type()); +} + +template +T gamma_p_derivative_imp(T a, T x, const Policy& pol) +{ + // + // Usual error checks first: + // + if(a <= 0) + policies::raise_domain_error("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); + if(x < 0) + policies::raise_domain_error("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); + // + // Now special cases: + // + if(x == 0) + { + return (a > 1) ? 0 : + (a == 1) ? 1 : policies::raise_overflow_error("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol); + } + // + // Normal case: + // + typedef typename lanczos::lanczos::type lanczos_type; + T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type()); + if((x < 1) && (tools::max_value() * x < f1)) + { + // overflow: + return policies::raise_overflow_error("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol); + } + + f1 /= x; + + return f1; +} + +template +inline typename tools::promote_args::type + tgamma(T z, const Policy& /* pol */, const mpl::true_) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast(detail::gamma_imp(static_cast(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)"); +} + +template +struct igamma_initializer +{ + struct init + { + init() + { + typedef typename policies::precision::type precision_type; + + typedef typename mpl::if_< + mpl::or_ >, + mpl::greater > >, + mpl::int_<0>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<53>, + typename mpl::if_< + mpl::less_equal >, + mpl::int_<64>, + mpl::int_<113> + >::type + >::type + >::type tag_type; + + do_init(tag_type()); + } + template + static void do_init(const mpl::int_&) + { + boost::math::gamma_p(static_cast(400), static_cast(400), Policy()); + } + static void do_init(const mpl::int_<53>&){} + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template +const typename igamma_initializer::init igamma_initializer::initializer; + +template +struct lgamma_initializer +{ + struct init + { + init() + { + typedef typename policies::precision::type precision_type; + typedef typename mpl::if_< + mpl::and_< + mpl::less_equal >, + mpl::greater > + >, + mpl::int_<64>, + typename mpl::if_< + mpl::and_< + mpl::less_equal >, + mpl::greater > + >, + mpl::int_<113>, mpl::int_<0> >::type + >::type tag_type; + do_init(tag_type()); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::lgamma(static_cast(2.5), Policy()); + boost::math::lgamma(static_cast(1.25), Policy()); + boost::math::lgamma(static_cast(1.75), Policy()); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::lgamma(static_cast(2.5), Policy()); + boost::math::lgamma(static_cast(1.25), Policy()); + boost::math::lgamma(static_cast(1.5), Policy()); + boost::math::lgamma(static_cast(1.75), Policy()); + } + static void do_init(const mpl::int_<0>&) + { + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template +const typename lgamma_initializer::init lgamma_initializer::initializer; + +template +inline typename tools::promote_args::type + tgamma(T1 a, T2 z, const Policy&, const mpl::false_) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + igamma_initializer::force_instantiate(); + + return policies::checked_narrowing_cast( + detail::gamma_incomplete_imp(static_cast(a), + static_cast(z), false, true, + forwarding_policy(), static_cast(0)), "boost::math::tgamma<%1%>(%1%, %1%)"); +} + +template +inline typename tools::promote_args::type + tgamma(T1 a, T2 z, const mpl::false_ tag) +{ + return tgamma(a, z, policies::policy<>(), tag); +} + + +} // namespace detail + +template +inline typename tools::promote_args::type + tgamma(T z) +{ + return tgamma(z, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + lgamma(T z, int* sign, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + detail::lgamma_initializer::force_instantiate(); + + return policies::checked_narrowing_cast(detail::lgamma_imp(static_cast(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)"); +} + +template +inline typename tools::promote_args::type + lgamma(T z, int* sign) +{ + return lgamma(z, sign, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + lgamma(T x, const Policy& pol) +{ + return ::boost::math::lgamma(x, 0, pol); +} + +template +inline typename tools::promote_args::type + lgamma(T x) +{ + return ::boost::math::lgamma(x, 0, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + tgamma1pm1(T z, const Policy& /* pol */) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)"); +} + +template +inline typename tools::promote_args::type + tgamma1pm1(T z) +{ + return tgamma1pm1(z, policies::policy<>()); +} + +// +// Full upper incomplete gamma: +// +template +inline typename tools::promote_args::type + tgamma(T1 a, T2 z) +{ + // + // Type T2 could be a policy object, or a value, select the + // right overload based on T2: + // + typedef typename policies::is_policy::type maybe_policy; + return detail::tgamma(a, z, maybe_policy()); +} +template +inline typename tools::promote_args::type + tgamma(T1 a, T2 z, const Policy& pol) +{ + return detail::tgamma(a, z, pol, mpl::false_()); +} +// +// Full lower incomplete gamma: +// +template +inline typename tools::promote_args::type + tgamma_lower(T1 a, T2 z, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + detail::igamma_initializer::force_instantiate(); + + return policies::checked_narrowing_cast( + detail::gamma_incomplete_imp(static_cast(a), + static_cast(z), false, false, + forwarding_policy(), static_cast(0)), "tgamma_lower<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + tgamma_lower(T1 a, T2 z) +{ + return tgamma_lower(a, z, policies::policy<>()); +} +// +// Regularised upper incomplete gamma: +// +template +inline typename tools::promote_args::type + gamma_q(T1 a, T2 z, const Policy& /* pol */) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + detail::igamma_initializer::force_instantiate(); + + return policies::checked_narrowing_cast( + detail::gamma_incomplete_imp(static_cast(a), + static_cast(z), true, true, + forwarding_policy(), static_cast(0)), "gamma_q<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + gamma_q(T1 a, T2 z) +{ + return gamma_q(a, z, policies::policy<>()); +} +// +// Regularised lower incomplete gamma: +// +template +inline typename tools::promote_args::type + gamma_p(T1 a, T2 z, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + detail::igamma_initializer::force_instantiate(); + + return policies::checked_narrowing_cast( + detail::gamma_incomplete_imp(static_cast(a), + static_cast(z), true, false, + forwarding_policy(), static_cast(0)), "gamma_p<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + gamma_p(T1 a, T2 z) +{ + return gamma_p(a, z, policies::policy<>()); +} + +// ratios of gamma functions: +template +inline typename tools::promote_args::type + tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast(detail::tgamma_delta_ratio_imp(static_cast(z), static_cast(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + tgamma_delta_ratio(T1 z, T2 delta) +{ + return tgamma_delta_ratio(z, delta, policies::policy<>()); +} +template +inline typename tools::promote_args::type + tgamma_ratio(T1 a, T2 b, const Policy&) +{ + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast(detail::tgamma_delta_ratio_imp(static_cast(a), static_cast(static_cast(b) - static_cast(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + tgamma_ratio(T1 a, T2 b) +{ + return tgamma_ratio(a, b, policies::policy<>()); +} + +template +inline typename tools::promote_args::type + gamma_p_derivative(T1 a, T2 x, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast(detail::gamma_p_derivative_imp(static_cast(a), static_cast(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)"); +} +template +inline typename tools::promote_args::type + gamma_p_derivative(T1 a, T2 x) +{ + return gamma_p_derivative(a, x, policies::policy<>()); +} + +} // namespace math +} // namespace boost + +#ifdef BOOST_MSVC +# pragma warning(pop) +#endif + +#include +#include +#include + +#endif // BOOST_MATH_SF_GAMMA_HPP + + + + diff --git a/patched/msm/back/tools.hpp b/patched/msm/back/tools.hpp new file mode 100644 index 0000000..641e625 --- /dev/null +++ b/patched/msm/back/tools.hpp @@ -0,0 +1,69 @@ +// Copyright 2008 Christophe Henry +// henry UNDERSCORE christophe AT hotmail DOT com +// This is an extended version of the state machine available in the boost::mpl library +// Distributed under the same license as the original. +// Copyright for the original version: +// Copyright 2005 David Abrahams and Aleksey Gurtovoy. Distributed +// under the Boost Software License, Version 1.0. (See accompanying +// file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MSM_BACK_TOOLS_H +#define BOOST_MSM_BACK_TOOLS_H + + +#include +#include +#include +#include +#include + +namespace boost { namespace msm { namespace back +{ + +// fills the array passed in with the state names in the correct order +// the array must be big enough. To know the needed size, use mpl::size +// on fsm::generate_state_set +template +struct fill_state_names +{ + fill_state_names(char const** names):m_names(names){} + template + void operator()(boost::msm::wrap const&) + { + m_names[get_state_id::value]= type_id().name_demangled(); + } +private: + char const** m_names; +}; + +// fills the typeid-generated name of the given state in the string passed as argument +template +struct get_state_name +{ + get_state_name(std::string& name_to_fill, int state_id):m_name(name_to_fill),m_state_id(state_id){} + template + void operator()(boost::msm::wrap const&) + { + if (get_state_id::value == m_state_id) + { + m_name = type_id().name_demangled(); + } + } +private: + std::string& m_name; + int m_state_id; +}; + +// displays the typeid of the given Type +struct display_type +{ + template + void operator()(boost::msm::wrap const&) + { + std::cout << type_id().name_demangled() << std::endl; + } +}; + +} } }//boost::msm::back +#endif //BOOST_MSM_BACK_TOOLS_H diff --git a/patched/property_map/dynamic_property_map.hpp b/patched/property_map/dynamic_property_map.hpp new file mode 100644 index 0000000..28c07d7 --- /dev/null +++ b/patched/property_map/dynamic_property_map.hpp @@ -0,0 +1,345 @@ +#ifndef DYNAMIC_PROPERTY_MAP_RG09302004_HPP +#define DYNAMIC_PROPERTY_MAP_RG09302004_HPP + +// Copyright 2004-5 The Trustees of Indiana University. + +// Use, modification and distribution is subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +// dynamic_property_map.hpp - +// Support for runtime-polymorphic property maps. This header is factored +// out of Doug Gregor's routines for reading GraphML files for use in reading +// GraphViz graph files. + +// Authors: Doug Gregor +// Ronald Garcia +// + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + +namespace detail { + + // read_value - + // A wrapper around lexical_cast, which does not behave as + // desired for std::string types. + template + inline Value read_value(const std::string& value) + { return boost::lexical_cast(value); } + + template<> + inline std::string read_value(const std::string& value) + { return value; } + +} + + +// dynamic_property_map - +// This interface supports polymorphic manipulation of property maps. +class dynamic_property_map +{ +public: + virtual ~dynamic_property_map() { } + + virtual boost::any get(const any& key) = 0; + virtual std::string get_string(const any& key) = 0; + virtual void put(const any& key, const any& value) = 0; + virtual type_index key() const = 0; + virtual type_index value() const = 0; +}; + + +////////////////////////////////////////////////////////////////////// +// Property map exceptions +////////////////////////////////////////////////////////////////////// + +struct dynamic_property_exception : public std::exception { + virtual ~dynamic_property_exception() throw() {} + virtual const char* what() const throw() = 0; +}; + +struct property_not_found : public dynamic_property_exception { + std::string property; + mutable std::string statement; + property_not_found(const std::string& property) : property(property) {} + virtual ~property_not_found() throw() {} + + const char* what() const throw() { + if(statement.empty()) + statement = + std::string("Property not found: ") + property + "."; + + return statement.c_str(); + } +}; + +struct dynamic_get_failure : public dynamic_property_exception { + std::string property; + mutable std::string statement; + dynamic_get_failure(const std::string& property) : property(property) {} + virtual ~dynamic_get_failure() throw() {} + + const char* what() const throw() { + if(statement.empty()) + statement = + std::string( + "dynamic property get cannot retrieve value for property: ") + + property + "."; + + return statement.c_str(); + } +}; + +struct dynamic_const_put_error : public dynamic_property_exception { + virtual ~dynamic_const_put_error() throw() {} + + const char* what() const throw() { + return "Attempt to put a value into a const property map: "; + } +}; + + +namespace detail { + +// Trying to work around VC++ problem that seems to relate to having too many +// functions named "get" +template +typename boost::property_traits::reference +get_wrapper_xxx(const PMap& pmap, const Key& key) { + using boost::get; + return get(pmap, key); +} + +// +// dynamic_property_map_adaptor - +// property-map adaptor to support runtime polymorphism. +template +class dynamic_property_map_adaptor : public dynamic_property_map +{ + typedef typename property_traits::key_type key_type; + typedef typename property_traits::value_type value_type; + typedef typename property_traits::category category; + + // do_put - overloaded dispatches from the put() member function. + // Attempts to "put" to a property map that does not model + // WritablePropertyMap result in a runtime exception. + + // in_value must either hold an object of value_type or a string that + // can be converted to value_type via iostreams. + void do_put(const any& in_key, const any& in_value, mpl::bool_) + { + using boost::put; + + key_type key = any_cast(in_key); + if (in_value.type() == type_id()) { + put(property_map_, key, any_cast(in_value)); + } else { + // if in_value is an empty string, put a default constructed value_type. + std::string v = any_cast(in_value); + if (v.empty()) { + put(property_map_, key, value_type()); + } else { + put(property_map_, key, detail::read_value(v)); + } + } + } + + void do_put(const any&, const any&, mpl::bool_) + { + BOOST_THROW_EXCEPTION(dynamic_const_put_error()); + } + +public: + explicit dynamic_property_map_adaptor(const PropertyMap& property_map_) + : property_map_(property_map_) { } + + virtual boost::any get(const any& key) + { + return get_wrapper_xxx(property_map_, any_cast(key)); + } + + virtual std::string get_string(const any& key) + { + std::ostringstream out; + out << get_wrapper_xxx(property_map_, any_cast(key)); + return out.str(); + } + + virtual void put(const any& in_key, const any& in_value) + { + do_put(in_key, in_value, + mpl::bool_<(is_convertible::value)>()); + } + + virtual type_index key() const { return type_id(); } + virtual type_index value() const { return type_id(); } + + PropertyMap& base() { return property_map_; } + const PropertyMap& base() const { return property_map_; } + +private: + PropertyMap property_map_; +}; + +} // namespace detail + +// +// dynamic_properties - +// container for dynamic property maps +// +struct dynamic_properties +{ + typedef std::multimap > + property_maps_type; + typedef boost::function3, + const std::string&, + const boost::any&, + const boost::any&> generate_fn_type; +public: + + typedef property_maps_type::iterator iterator; + typedef property_maps_type::const_iterator const_iterator; + + dynamic_properties() : generate_fn() { } + dynamic_properties(const generate_fn_type& g) : generate_fn(g) {} + + ~dynamic_properties() {} + + template + dynamic_properties& + property(const std::string& name, PropertyMap property_map_) + { + boost::shared_ptr pm( + boost::make_shared >(property_map_)); + property_maps.insert(property_maps_type::value_type(name, pm)); + + return *this; + } + + iterator begin() { return property_maps.begin(); } + const_iterator begin() const { return property_maps.begin(); } + iterator end() { return property_maps.end(); } + const_iterator end() const { return property_maps.end(); } + + iterator lower_bound(const std::string& name) + { return property_maps.lower_bound(name); } + + const_iterator lower_bound(const std::string& name) const + { return property_maps.lower_bound(name); } + + void + insert(const std::string& name, boost::shared_ptr pm) + { + property_maps.insert(property_maps_type::value_type(name, pm)); + } + + template + boost::shared_ptr + generate(const std::string& name, const Key& key, const Value& value) + { + if(!generate_fn) { + BOOST_THROW_EXCEPTION(property_not_found(name)); + } else { + return generate_fn(name,key,value); + } + } + +private: + property_maps_type property_maps; + generate_fn_type generate_fn; +}; + +template +bool +put(const std::string& name, dynamic_properties& dp, const Key& key, + const Value& value) +{ + for (dynamic_properties::iterator i = dp.lower_bound(name); + i != dp.end() && i->first == name; ++i) { + if (i->second->key() == type_id()) { + i->second->put(key, value); + return true; + } + } + + boost::shared_ptr new_map = dp.generate(name, key, value); + if (new_map.get()) { + new_map->put(key, value); + dp.insert(name, new_map); + return true; + } else { + return false; + } +} + +#ifndef BOOST_NO_EXPLICIT_FUNCTION_TEMPLATE_ARGUMENTS +template +Value +get(const std::string& name, const dynamic_properties& dp, const Key& key) +{ + for (dynamic_properties::const_iterator i = dp.lower_bound(name); + i != dp.end() && i->first == name; ++i) { + if (i->second->key() == type_id()) + return any_cast(i->second->get(key)); + } + + BOOST_THROW_EXCEPTION(dynamic_get_failure(name)); +} +#endif + +template +Value +get(const std::string& name, const dynamic_properties& dp, const Key& key, type) +{ + for (dynamic_properties::const_iterator i = dp.lower_bound(name); + i != dp.end() && i->first == name; ++i) { + if (i->second->key() == type_id()) + return any_cast(i->second->get(key)); + } + + BOOST_THROW_EXCEPTION(dynamic_get_failure(name)); +} + +template +std::string +get(const std::string& name, const dynamic_properties& dp, const Key& key) +{ + for (dynamic_properties::const_iterator i = dp.lower_bound(name); + i != dp.end() && i->first == name; ++i) { + if (i->second->key() == type_id()) + return i->second->get_string(key); + } + + BOOST_THROW_EXCEPTION(dynamic_get_failure(name)); +} + +// The easy way to ignore properties. +inline +boost::shared_ptr +ignore_other_properties(const std::string&, + const boost::any&, + const boost::any&) { + return boost::shared_ptr(); +} + +} // namespace boost + +#endif // DYNAMIC_PROPERTY_MAP_RG09302004_HPP diff --git a/patched/property_tree/detail/ptree_implementation.hpp b/patched/property_tree/detail/ptree_implementation.hpp new file mode 100644 index 0000000..bf40c72 --- /dev/null +++ b/patched/property_tree/detail/ptree_implementation.hpp @@ -0,0 +1,915 @@ +// ---------------------------------------------------------------------------- +// Copyright (C) 2002-2006 Marcin Kalicinski +// Copyright (C) 2009 Sebastian Redl +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// +// For more information, see www.boost.org +// ---------------------------------------------------------------------------- +#ifndef BOOST_PROPERTY_TREE_DETAIL_PTREE_IMPLEMENTATION_HPP_INCLUDED +#define BOOST_PROPERTY_TREE_DETAIL_PTREE_IMPLEMENTATION_HPP_INCLUDED + +#include +#include +#include +#include +#include +#include + +#if (defined(BOOST_MSVC) && \ + (_MSC_FULL_VER >= 160000000 && _MSC_FULL_VER < 170000000)) || \ + (defined(BOOST_INTEL_WIN) && \ + defined(BOOST_DINKUMWARE_STDLIB)) +#define BOOST_PROPERTY_TREE_PAIR_BUG +#endif + +namespace boost { namespace property_tree +{ + template + struct basic_ptree::subs + { + struct by_name {}; + // The actual child container. +#if defined(BOOST_PROPERTY_TREE_PAIR_BUG) + // MSVC 10 has moved std::pair's members to a base + // class. Unfortunately this does break the interface. + BOOST_STATIC_CONSTANT(unsigned, + first_offset = offsetof(value_type, first)); +#endif + typedef multi_index_container, + multi_index::ordered_non_unique, +#if defined(BOOST_PROPERTY_TREE_PAIR_BUG) + multi_index::member_offset, +#else + multi_index::member, +#endif + key_compare + > + > + > base_container; + + // The by-name lookup index. + typedef typename base_container::template index::type + by_name_index; + + // Access functions for getting to the children of a tree. + static base_container& ch(self_type *s) { + return *static_cast(s->m_children); + } + static const base_container& ch(const self_type *s) { + return *static_cast(s->m_children); + } + static by_name_index& assoc(self_type *s) { + return ch(s).BOOST_NESTED_TEMPLATE get(); + } + static const by_name_index& assoc(const self_type *s) { + return ch(s).BOOST_NESTED_TEMPLATE get(); + } + }; + template + class basic_ptree::iterator : public boost::iterator_adaptor< + iterator, typename subs::base_container::iterator, value_type> + { + friend class boost::iterator_core_access; + typedef boost::iterator_adaptor< + iterator, typename subs::base_container::iterator, value_type> + baset; + public: + typedef typename baset::reference reference; + iterator() {} + explicit iterator(typename iterator::base_type b) + : iterator::iterator_adaptor_(b) + {} + reference dereference() const + { + // multi_index doesn't allow modification of its values, because + // indexes could sort by anything, and modification screws that up. + // However, we only sort by the key, and it's protected against + // modification in the value_type, so this const_cast is safe. + return const_cast(*this->base_reference()); + } + }; + template + class basic_ptree::const_iterator : public boost::iterator_adaptor< + const_iterator, typename subs::base_container::const_iterator> + { + public: + const_iterator() {} + explicit const_iterator(typename const_iterator::base_type b) + : const_iterator::iterator_adaptor_(b) + {} + const_iterator(iterator b) + : const_iterator::iterator_adaptor_(b.base()) + {} + }; + template + class basic_ptree::reverse_iterator + : public boost::reverse_iterator + { + public: + reverse_iterator() {} + explicit reverse_iterator(iterator b) + : boost::reverse_iterator(b) + {} + }; + template + class basic_ptree::const_reverse_iterator + : public boost::reverse_iterator + { + public: + const_reverse_iterator() {} + explicit const_reverse_iterator(const_iterator b) + : boost::reverse_iterator(b) + {} + const_reverse_iterator( + typename basic_ptree::reverse_iterator b) + : boost::reverse_iterator(b) + {} + }; + template + class basic_ptree::assoc_iterator + : public boost::iterator_adaptor + { + friend class boost::iterator_core_access; + typedef boost::iterator_adaptor + baset; + public: + typedef typename baset::reference reference; + assoc_iterator() {} + explicit assoc_iterator(typename assoc_iterator::base_type b) + : assoc_iterator::iterator_adaptor_(b) + {} + reference dereference() const + { + return const_cast(*this->base_reference()); + } + }; + template + class basic_ptree::const_assoc_iterator + : public boost::iterator_adaptor + { + public: + const_assoc_iterator() {} + explicit const_assoc_iterator( + typename const_assoc_iterator::base_type b) + : const_assoc_iterator::iterator_adaptor_(b) + {} + const_assoc_iterator(assoc_iterator b) + : const_assoc_iterator::iterator_adaptor_(b.base()) + {} + }; + + + // Big five + + // Perhaps the children collection could be created on-demand only, to + // reduce heap traffic. But that's a lot more work to implement. + + template inline + basic_ptree::basic_ptree() + : m_children(new typename subs::base_container) + { + } + + template inline + basic_ptree::basic_ptree(const data_type &d) + : m_data(d), m_children(new typename subs::base_container) + { + } + + template inline + basic_ptree::basic_ptree(const basic_ptree &rhs) + : m_data(rhs.m_data), + m_children(new typename subs::base_container(subs::ch(&rhs))) + { + } + + template + basic_ptree & + basic_ptree::operator =(const basic_ptree &rhs) + { + self_type(rhs).swap(*this); + return *this; + } + + template + basic_ptree::~basic_ptree() + { + delete &subs::ch(this); + } + + template inline + void basic_ptree::swap(basic_ptree &rhs) + { + boost::swap(m_data, rhs.m_data); + // Void pointers, no ADL necessary + std::swap(m_children, rhs.m_children); + } + + // Container view + + template inline + typename basic_ptree::size_type + basic_ptree::size() const + { + return subs::ch(this).size(); + } + + template inline + typename basic_ptree::size_type + basic_ptree::max_size() const + { + return subs::ch(this).max_size(); + } + + template inline + bool basic_ptree::empty() const + { + return subs::ch(this).empty(); + } + + template inline + typename basic_ptree::iterator + basic_ptree::begin() + { + return iterator(subs::ch(this).begin()); + } + + template inline + typename basic_ptree::const_iterator + basic_ptree::begin() const + { + return const_iterator(subs::ch(this).begin()); + } + + template inline + typename basic_ptree::iterator + basic_ptree::end() + { + return iterator(subs::ch(this).end()); + } + + template inline + typename basic_ptree::const_iterator + basic_ptree::end() const + { + return const_iterator(subs::ch(this).end()); + } + + template inline + typename basic_ptree::reverse_iterator + basic_ptree::rbegin() + { + return reverse_iterator(this->end()); + } + + template inline + typename basic_ptree::const_reverse_iterator + basic_ptree::rbegin() const + { + return const_reverse_iterator(this->end()); + } + + template inline + typename basic_ptree::reverse_iterator + basic_ptree::rend() + { + return reverse_iterator(this->begin()); + } + + template inline + typename basic_ptree::const_reverse_iterator + basic_ptree::rend() const + { + return const_reverse_iterator(this->begin()); + } + + template inline + typename basic_ptree::value_type & + basic_ptree::front() + { + return const_cast(subs::ch(this).front()); + } + + template inline + const typename basic_ptree::value_type & + basic_ptree::front() const + { + return subs::ch(this).front(); + } + + template inline + typename basic_ptree::value_type & + basic_ptree::back() + { + return const_cast(subs::ch(this).back()); + } + + template inline + const typename basic_ptree::value_type & + basic_ptree::back() const + { + return subs::ch(this).back(); + } + + template inline + typename basic_ptree::iterator + basic_ptree::insert(iterator where, const value_type &value) + { + return iterator(subs::ch(this).insert(where.base(), value).first); + } + + template + template inline + void basic_ptree::insert(iterator where, It first, It last) + { + subs::ch(this).insert(where.base(), first, last); + } + + template inline + typename basic_ptree::iterator + basic_ptree::erase(iterator where) + { + return iterator(subs::ch(this).erase(where.base())); + } + + template inline + typename basic_ptree::iterator + basic_ptree::erase(iterator first, iterator last) + { + return iterator(subs::ch(this).erase(first.base(), last.base())); + } + + template inline + typename basic_ptree::iterator + basic_ptree::push_front(const value_type &value) + { + return iterator(subs::ch(this).push_front(value).first); + } + + template inline + typename basic_ptree::iterator + basic_ptree::push_back(const value_type &value) + { + return iterator(subs::ch(this).push_back(value).first); + } + + template inline + void basic_ptree::pop_front() + { + subs::ch(this).pop_front(); + } + + template inline + void basic_ptree::pop_back() + { + subs::ch(this).pop_back(); + } + + template inline + void basic_ptree::reverse() + { + subs::ch(this).reverse(); + } + + namespace impl + { + struct by_first + { + template + bool operator ()(const P& lhs, const P& rhs) const { + return lhs.first < rhs.first; + } + }; + } + + template inline + void basic_ptree::sort() + { + sort(impl::by_first()); + } + + template + template inline + void basic_ptree::sort(Compare comp) + { + subs::ch(this).sort(comp); + } + + // Equality + + template inline + bool basic_ptree::operator ==( + const basic_ptree &rhs) const + { + // The size test is cheap, so add it as an optimization + return size() == rhs.size() && data() == rhs.data() && + subs::ch(this) == subs::ch(&rhs); + } + + template inline + bool basic_ptree::operator !=( + const basic_ptree &rhs) const + { + return !(*this == rhs); + } + + // Associative view + + template inline + typename basic_ptree::assoc_iterator + basic_ptree::ordered_begin() + { + return assoc_iterator(subs::assoc(this).begin()); + } + + template inline + typename basic_ptree::const_assoc_iterator + basic_ptree::ordered_begin() const + { + return const_assoc_iterator(subs::assoc(this).begin()); + } + + template inline + typename basic_ptree::assoc_iterator + basic_ptree::not_found() + { + return assoc_iterator(subs::assoc(this).end()); + } + + template inline + typename basic_ptree::const_assoc_iterator + basic_ptree::not_found() const + { + return const_assoc_iterator(subs::assoc(this).end()); + } + + template inline + typename basic_ptree::assoc_iterator + basic_ptree::find(const key_type &key) + { + return assoc_iterator(subs::assoc(this).find(key)); + } + + template inline + typename basic_ptree::const_assoc_iterator + basic_ptree::find(const key_type &key) const + { + return const_assoc_iterator(subs::assoc(this).find(key)); + } + + template inline + std::pair< + typename basic_ptree::assoc_iterator, + typename basic_ptree::assoc_iterator + > basic_ptree::equal_range(const key_type &key) + { + std::pair r( + subs::assoc(this).equal_range(key)); + return std::pair( + assoc_iterator(r.first), assoc_iterator(r.second)); + } + + template inline + std::pair< + typename basic_ptree::const_assoc_iterator, + typename basic_ptree::const_assoc_iterator + > basic_ptree::equal_range(const key_type &key) const + { + std::pair r( + subs::assoc(this).equal_range(key)); + return std::pair( + const_assoc_iterator(r.first), const_assoc_iterator(r.second)); + } + + template inline + typename basic_ptree::size_type + basic_ptree::count(const key_type &key) const + { + return subs::assoc(this).count(key); + } + + template inline + typename basic_ptree::size_type + basic_ptree::erase(const key_type &key) + { + return subs::assoc(this).erase(key); + } + + template inline + typename basic_ptree::iterator + basic_ptree::to_iterator(assoc_iterator ai) + { + return iterator(subs::ch(this). + BOOST_NESTED_TEMPLATE project<0>(ai.base())); + } + + template inline + typename basic_ptree::const_iterator + basic_ptree::to_iterator(const_assoc_iterator ai) const + { + return const_iterator(subs::ch(this). + BOOST_NESTED_TEMPLATE project<0>(ai.base())); + } + + // Property tree view + + template inline + typename basic_ptree::data_type & + basic_ptree::data() + { + return m_data; + } + + template inline + const typename basic_ptree::data_type & + basic_ptree::data() const + { + return m_data; + } + + template inline + void basic_ptree::clear() + { + m_data = data_type(); + subs::ch(this).clear(); + } + + template + basic_ptree & + basic_ptree::get_child(const path_type &path) + { + path_type p(path); + self_type *n = walk_path(p); + if (!n) { + BOOST_PROPERTY_TREE_THROW(ptree_bad_path("No such node", path)); + } + return *n; + } + + template inline + const basic_ptree & + basic_ptree::get_child(const path_type &path) const + { + return const_cast(this)->get_child(path); + } + + template inline + basic_ptree & + basic_ptree::get_child(const path_type &path, + self_type &default_value) + { + path_type p(path); + self_type *n = walk_path(p); + return n ? *n : default_value; + } + + template inline + const basic_ptree & + basic_ptree::get_child(const path_type &path, + const self_type &default_value) const + { + return const_cast(this)->get_child(path, + const_cast(default_value)); + } + + + template + optional &> + basic_ptree::get_child_optional(const path_type &path) + { + path_type p(path); + self_type *n = walk_path(p); + if (!n) { + return optional(); + } + return *n; + } + + template + optional &> + basic_ptree::get_child_optional(const path_type &path) const + { + path_type p(path); + self_type *n = walk_path(p); + if (!n) { + return optional(); + } + return *n; + } + + template + basic_ptree & + basic_ptree::put_child(const path_type &path, + const self_type &value) + { + path_type p(path); + self_type &parent = force_path(p); + // Got the parent. Now get the correct child. + key_type fragment = p.reduce(); + assoc_iterator el = parent.find(fragment); + // If the new child exists, replace it. + if(el != parent.not_found()) { + return el->second = value; + } else { + return parent.push_back(value_type(fragment, value))->second; + } + } + + template + basic_ptree & + basic_ptree::add_child(const path_type &path, + const self_type &value) + { + path_type p(path); + self_type &parent = force_path(p); + // Got the parent. + key_type fragment = p.reduce(); + return parent.push_back(value_type(fragment, value))->second; + } + + template + template + typename boost::enable_if, Type>::type + basic_ptree::get_value(Translator tr) const + { + if(boost::optional o = get_value_optional(tr)) { + return *o; + } + BOOST_PROPERTY_TREE_THROW(ptree_bad_data( + std::string("conversion of data to type \"") + + type_id().name_demangled() + "\" failed", data())); + } + + template + template inline + Type basic_ptree::get_value() const + { + return get_value( + typename translator_between::type()); + } + + template + template inline + Type basic_ptree::get_value(const Type &default_value, + Translator tr) const + { + return get_value_optional(tr).get_value_or(default_value); + } + + template + template + typename boost::enable_if< + detail::is_character, + std::basic_string + >::type + basic_ptree::get_value(const Ch *default_value, Translator tr)const + { + return get_value, Translator>(default_value, tr); + } + + template + template inline + typename boost::disable_if, Type>::type + basic_ptree::get_value(const Type &default_value) const + { + return get_value(default_value, + typename translator_between::type()); + } + + template + template + typename boost::enable_if< + detail::is_character, + std::basic_string + >::type + basic_ptree::get_value(const Ch *default_value) const + { + return get_value< std::basic_string >(default_value); + } + + template + template inline + optional basic_ptree::get_value_optional( + Translator tr) const + { + return tr.get_value(data()); + } + + template + template inline + optional basic_ptree::get_value_optional() const + { + return get_value_optional( + typename translator_between::type()); + } + + template + template inline + typename boost::enable_if, Type>::type + basic_ptree::get(const path_type &path, + Translator tr) const + { + return get_child(path).BOOST_NESTED_TEMPLATE get_value(tr); + } + + template + template inline + Type basic_ptree::get(const path_type &path) const + { + return get_child(path).BOOST_NESTED_TEMPLATE get_value(); + } + + template + template inline + Type basic_ptree::get(const path_type &path, + const Type &default_value, + Translator tr) const + { + return get_optional(path, tr).get_value_or(default_value); + } + + template + template + typename boost::enable_if< + detail::is_character, + std::basic_string + >::type + basic_ptree::get( + const path_type &path, const Ch *default_value, Translator tr) const + { + return get, Translator>(path, default_value, tr); + } + + template + template inline + typename boost::disable_if, Type>::type + basic_ptree::get(const path_type &path, + const Type &default_value) const + { + return get_optional(path).get_value_or(default_value); + } + + template + template + typename boost::enable_if< + detail::is_character, + std::basic_string + >::type + basic_ptree::get( + const path_type &path, const Ch *default_value) const + { + return get< std::basic_string >(path, default_value); + } + + template + template + optional basic_ptree::get_optional(const path_type &path, + Translator tr) const + { + if (optional child = get_child_optional(path)) + return child.get(). + BOOST_NESTED_TEMPLATE get_value_optional(tr); + else + return optional(); + } + + template + template + optional basic_ptree::get_optional( + const path_type &path) const + { + if (optional child = get_child_optional(path)) + return child.get().BOOST_NESTED_TEMPLATE get_value_optional(); + else + return optional(); + } + + template + template + void basic_ptree::put_value(const Type &value, Translator tr) + { + if(optional o = tr.put_value(value)) { + data() = *o; + } else { + BOOST_PROPERTY_TREE_THROW(ptree_bad_data( + std::string("conversion of type \"") + type_id().name_demangled() + + "\" to data failed", boost::any())); + } + } + + template + template inline + void basic_ptree::put_value(const Type &value) + { + put_value(value, typename translator_between::type()); + } + + template + template + basic_ptree & basic_ptree::put( + const path_type &path, const Type &value, Translator tr) + { + if(optional child = get_child_optional(path)) { + child.get().put_value(value, tr); + return *child; + } else { + self_type &child2 = put_child(path, self_type()); + child2.put_value(value, tr); + return child2; + } + } + + template + template inline + basic_ptree & basic_ptree::put( + const path_type &path, const Type &value) + { + return put(path, value, + typename translator_between::type()); + } + + template + template inline + basic_ptree & basic_ptree::add( + const path_type &path, const Type &value, Translator tr) + { + self_type &child = add_child(path, self_type()); + child.put_value(value, tr); + return child; + } + + template + template inline + basic_ptree & basic_ptree::add( + const path_type &path, const Type &value) + { + return add(path, value, + typename translator_between::type()); + } + + + template + basic_ptree * + basic_ptree::walk_path(path_type &p) const + { + if(p.empty()) { + // I'm the child we're looking for. + return const_cast(this); + } + // Recurse down the tree to find the path. + key_type fragment = p.reduce(); + const_assoc_iterator el = find(fragment); + if(el == not_found()) { + // No such child. + return 0; + } + // Not done yet, recurse. + return el->second.walk_path(p); + } + + template + basic_ptree & basic_ptree::force_path(path_type &p) + { + BOOST_ASSERT(!p.empty() && "Empty path not allowed for put_child."); + if(p.single()) { + // I'm the parent we're looking for. + return *this; + } + key_type fragment = p.reduce(); + assoc_iterator el = find(fragment); + // If we've found an existing child, go down that path. Else + // create a new one. + self_type& child = el == not_found() ? + push_back(value_type(fragment, self_type()))->second : el->second; + return child.force_path(p); + } + + // Free functions + + template + inline void swap(basic_ptree &pt1, basic_ptree &pt2) + { + pt1.swap(pt2); + } + +} } + +#if defined(BOOST_PROPERTY_TREE_PAIR_BUG) +#undef BOOST_PROPERTY_TREE_PAIR_BUG +#endif + +#endif