///////////////////////////////////////////////////////////////////////////////

// Copyright Christopher Kormanyos 2002 - 2013.

// Copyright 2011 -2013 John Maddock. 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)

//

// This work is based on an earlier work:

// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",

// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469

//

// Note that there are no "noexcept" specifications on the functions in this file: there are too many

// calls to lexical_cast (and similar) to easily analyse the code for correctness. So until compilers

// can detect noexcept misuse at compile time, the only realistic option is to simply not use it here.

//

#ifndef BOOST_MP_CPP_DEC_FLOAT_BACKEND_HPP

#define BOOST_MP_CPP_DEC_FLOAT_BACKEND_HPP

#include <boost/config.hpp>

#include <boost/cstdint.hpp>

#include <limits>

#ifndef BOOST_NO_CXX11_HDR_ARRAY

#include <array>

#else

#include <boost/array.hpp>

#endif

#include <boost/cstdint.hpp>

#include <boost/functional/hash_fwd.hpp>

#include <boost/multiprecision/number.hpp>

#include <boost/multiprecision/detail/big_lanczos.hpp>

#include <boost/multiprecision/detail/dynamic_array.hpp>

//

// Headers required for Boost.Math integration:

//

#include <boost/math/policies/policy.hpp>

//

// Some includes we need from Boost.Math, since we rely on that library to provide these functions:

//

#include <boost/math/special_functions/asinh.hpp>

#include <boost/math/special_functions/acosh.hpp>

#include <boost/math/special_functions/atanh.hpp>

#include <boost/math/special_functions/cbrt.hpp>

#include <boost/math/special_functions/expm1.hpp>

#include <boost/math/special_functions/gamma.hpp>

#ifdef BOOST_MSVC

#pragma warning(push)

#pragma warning(disable:6326) // comparison of two constants

#endif

namespace boost{

namespace multiprecision{

namespace backends{

template <unsigned Digits10, class ExponentType = boost::int32_t, class Allocator = void>

class cpp_dec_float;

} // namespace

template <unsigned Digits10, class ExponentType, class Allocator>

struct number_category<backends::cpp_dec_float<Digits10, ExponentType, Allocator> > : public mpl::int_<number_kind_floating_point>{};

namespace backends{

template <unsigned Digits10, class ExponentType, class Allocator>

class cpp_dec_float

{

private:

static const boost::int32_t cpp_dec_float_digits10_setting = Digits10;

// We need at least 16-bits in the exponent type to do anything sensible:

BOOST_STATIC_ASSERT_MSG(boost::is_signed<ExponentType>::value, "ExponentType must be a signed built in integer type.");

BOOST_STATIC_ASSERT_MSG(sizeof(ExponentType) > 1, "ExponentType is too small.");

public:

typedef mpl::list<boost::long_long_type> signed_types;

typedef mpl::list<boost::ulong_long_type> unsigned_types;

typedef mpl::list<long double> float_types;

typedef ExponentType exponent_type;

static const boost::int32_t cpp_dec_float_radix = 10L;

static const boost::int32_t cpp_dec_float_digits10_limit_lo = 9L;

static const boost::int32_t cpp_dec_float_digits10_limit_hi = boost::integer_traits<boost::int32_t>::const_max - 100;

static const boost::int32_t cpp_dec_float_digits10 = ((cpp_dec_float_digits10_setting < cpp_dec_float_digits10_limit_lo) ? cpp_dec_float_digits10_limit_lo : ((cpp_dec_float_digits10_setting > cpp_dec_float_digits10_limit_hi) ? cpp_dec_float_digits10_limit_hi : cpp_dec_float_digits10_setting));

static const ExponentType cpp_dec_float_max_exp10 = (static_cast<ExponentType>(1) << (std::numeric_limits<ExponentType>::digits - 5));

static const ExponentType cpp_dec_float_min_exp10 = -cpp_dec_float_max_exp10;

static const ExponentType cpp_dec_float_max_exp = cpp_dec_float_max_exp10;

static const ExponentType cpp_dec_float_min_exp = cpp_dec_float_min_exp10;

BOOST_STATIC_ASSERT((cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp10 == -cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp10));

private:

static const boost::int32_t cpp_dec_float_elem_digits10 = 8L;

static const boost::int32_t cpp_dec_float_elem_mask = 100000000L;

BOOST_STATIC_ASSERT(0 == cpp_dec_float_max_exp10 % cpp_dec_float_elem_digits10);

// There are three guard limbs.

// 1) The first limb has 'play' from 1...8 decimal digits.

// 2) The last limb also has 'play' from 1...8 decimal digits.

// 3) One limb can get lost when justifying after multiply,

// as only half of the triangle is multiplied and a carry

// from below is missing.

static const boost::int32_t cpp_dec_float_elem_number_request = static_cast<boost::int32_t>((cpp_dec_float_digits10 / cpp_dec_float_elem_digits10) + (((cpp_dec_float_digits10 % cpp_dec_float_elem_digits10) != 0) ? 1 : 0));

// The number of elements needed (with a minimum of two) plus three added guard limbs.

static const boost::int32_t cpp_dec_float_elem_number = static_cast<boost::int32_t>(((cpp_dec_float_elem_number_request < 2L) ? 2L : cpp_dec_float_elem_number_request) + 3L);

public:

static const boost::int32_t cpp_dec_float_total_digits10 = static_cast<boost::int32_t>(cpp_dec_float_elem_number * cpp_dec_float_elem_digits10);

private:

typedef enum enum_fpclass_type

{

cpp_dec_float_finite,

cpp_dec_float_inf,

cpp_dec_float_NaN

}

fpclass_type;

#ifndef BOOST_NO_CXX11_HDR_ARRAY

typedef typename mpl::if_<is_void<Allocator>,

std::array<boost::uint32_t, cpp_dec_float_elem_number>,

detail::dynamic_array<boost::uint32_t, cpp_dec_float_elem_number, Allocator>

>::type array_type;

#else

typedef typename mpl::if_<is_void<Allocator>,

boost::array<boost::uint32_t, cpp_dec_float_elem_number>,

detail::dynamic_array<boost::uint32_t, cpp_dec_float_elem_number, Allocator>

>::type array_type;

#endif

array_type data;

ExponentType exp;

bool neg;

fpclass_type fpclass;

boost::int32_t prec_elem;

//

// Special values constructor:

//

cpp_dec_float(fpclass_type c) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (c),

prec_elem(cpp_dec_float_elem_number) { }

//

// Static data initializer:

//

struct initializer

{

initializer()

{

cpp_dec_float<Digits10, ExponentType, Allocator>::nan();

cpp_dec_float<Digits10, ExponentType, Allocator>::inf();

(cpp_dec_float<Digits10, ExponentType, Allocator>::min)();

(cpp_dec_float<Digits10, ExponentType, Allocator>::max)();

cpp_dec_float<Digits10, ExponentType, Allocator>::zero();

cpp_dec_float<Digits10, ExponentType, Allocator>::one();

cpp_dec_float<Digits10, ExponentType, Allocator>::two();

cpp_dec_float<Digits10, ExponentType, Allocator>::half();

cpp_dec_float<Digits10, ExponentType, Allocator>::double_min();

cpp_dec_float<Digits10, ExponentType, Allocator>::double_max();

cpp_dec_float<Digits10, ExponentType, Allocator>::long_double_max();

cpp_dec_float<Digits10, ExponentType, Allocator>::long_double_min();

cpp_dec_float<Digits10, ExponentType, Allocator>::long_long_max();

cpp_dec_float<Digits10, ExponentType, Allocator>::long_long_min();

cpp_dec_float<Digits10, ExponentType, Allocator>::ulong_long_max();

cpp_dec_float<Digits10, ExponentType, Allocator>::eps();

cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(0);

}

void do_nothing(){}

};

static initializer init;

public:

// Constructors

cpp_dec_float() BOOST_MP_NOEXCEPT_IF(noexcept(array_type())) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (cpp_dec_float_finite),

prec_elem(cpp_dec_float_elem_number) { }

cpp_dec_float(const char* s) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (cpp_dec_float_finite),

prec_elem(cpp_dec_float_elem_number)

{

*this = s;

}

template<class I>

cpp_dec_float(I i, typename enable_if<is_unsigned<I> >::type* = 0) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (cpp_dec_float_finite),

prec_elem(cpp_dec_float_elem_number)

{

from_unsigned_long_long(i);

}

template <class I>

cpp_dec_float(I i, typename enable_if<is_signed<I> >::type* = 0) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (cpp_dec_float_finite),

prec_elem(cpp_dec_float_elem_number)

{

if(i < 0)

{

from_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(i));

negate();

}

else

from_unsigned_long_long(i);

}

cpp_dec_float(const cpp_dec_float& f) BOOST_MP_NOEXCEPT_IF(noexcept(array_type(std::declval<const array_type&>()))) :

data (f.data),

exp (f.exp),

neg (f.neg),

fpclass (f.fpclass),

prec_elem(f.prec_elem) { }

template <unsigned D, class ET, class A>

cpp_dec_float(const cpp_dec_float<D, ET, A>& f, typename enable_if_c<D <= Digits10>::type* = 0) :

data(),

exp (f.exp),

neg (f.neg),

fpclass (static_cast<fpclass_type>(static_cast<int>(f.fpclass))),

prec_elem(cpp_dec_float_elem_number)

{

std::copy(f.data.begin(), f.data.begin() + f.prec_elem, data.begin());

}

template <unsigned D, class ET, class A>

explicit cpp_dec_float(const cpp_dec_float<D, ET, A>& f, typename disable_if_c<D <= Digits10>::type* = 0) :

data(),

exp (f.exp),

neg (f.neg),

fpclass (static_cast<fpclass_type>(static_cast<int>(f.fpclass))),

prec_elem(cpp_dec_float_elem_number)

{

// TODO: this doesn't round!

std::copy(f.data.begin(), f.data.begin() + prec_elem, data.begin());

}

template <class F>

cpp_dec_float(const F val, typename enable_if_c<is_floating_point<F>::value

#ifdef BOOST_HAS_FLOAT128

&& !boost::is_same<F, __float128>::value

#endif

>::type* = 0) :

data(),

exp (static_cast<ExponentType>(0)),

neg (false),

fpclass (cpp_dec_float_finite),

prec_elem(cpp_dec_float_elem_number)

{

*this = val;

}

cpp_dec_float(const double mantissa, const ExponentType exponent);

std::size_t hash()const

{

std::size_t result = 0;

for(int i = 0; i < prec_elem; ++i)

boost::hash_combine(result, data[i]);

boost::hash_combine(result, exp);

boost::hash_combine(result, neg);

boost::hash_combine(result, fpclass);

return result;

}

// Specific special values.

static const cpp_dec_float& nan()

{

static const cpp_dec_float val(cpp_dec_float_NaN);

init.do_nothing();

return val;

}

static const cpp_dec_float& inf()

{

static const cpp_dec_float val(cpp_dec_float_inf);

init.do_nothing();

return val;

}

static const cpp_dec_float& (max)()

{

init.do_nothing();

static cpp_dec_float val_max = std::string("1.0e" + boost::lexical_cast<std::string>(cpp_dec_float_max_exp10)).c_str();

return val_max;

}

static const cpp_dec_float& (min)()

{

init.do_nothing();

static cpp_dec_float val_min = std::string("1.0e" + boost::lexical_cast<std::string>(cpp_dec_float_min_exp10)).c_str();

return val_min;

}

static const cpp_dec_float& zero()

{

init.do_nothing();

static cpp_dec_float val(static_cast<boost::ulong_long_type>(0u));

return val;

}

static const cpp_dec_float& one()

{

init.do_nothing();

static cpp_dec_float val(static_cast<boost::ulong_long_type>(1u));

return val;

}

static const cpp_dec_float& two()

{

init.do_nothing();

static cpp_dec_float val(static_cast<boost::ulong_long_type>(2u));

return val;

}

static const cpp_dec_float& half()

{

init.do_nothing();

static cpp_dec_float val(0.5L);

return val;

}

static const cpp_dec_float& double_min()

{

init.do_nothing();

static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::min)()));

return val;

}

static const cpp_dec_float& double_max()

{

init.do_nothing();

static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::max)()));

return val;

}

static const cpp_dec_float& long_double_min()

{

init.do_nothing();

#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS

static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::min)()));

#else

static cpp_dec_float val((std::numeric_limits<long double>::min)());

#endif

return val;

}

static const cpp_dec_float& long_double_max()

{

init.do_nothing();

#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS

static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::max)()));

#else

static cpp_dec_float val((std::numeric_limits<long double>::max)());

#endif

return val;

}

static const cpp_dec_float& long_long_max()

{

init.do_nothing();

static cpp_dec_float val((std::numeric_limits<boost::long_long_type>::max)());

return val;

}

static const cpp_dec_float& long_long_min()

{

init.do_nothing();

static cpp_dec_float val((std::numeric_limits<boost::long_long_type>::min)());

return val;

}

static const cpp_dec_float& ulong_long_max()

{

init.do_nothing();

static cpp_dec_float val((std::numeric_limits<boost::ulong_long_type>::max)());

return val;

}

static const cpp_dec_float& eps()

{

init.do_nothing();

static cpp_dec_float val(1.0, 1 - static_cast<int>(cpp_dec_float_digits10));

return val;

}

// Basic operations.

cpp_dec_float& operator=(const cpp_dec_float& v) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<array_type&>() = std::declval<const array_type&>()))

{

data = v.data;

exp = v.exp;

neg = v.neg;

fpclass = v.fpclass;

prec_elem = v.prec_elem;

return *this;

}

template <unsigned D>

cpp_dec_float& operator=(const cpp_dec_float<D>& f)

{

exp = f.exp;

neg = f.neg;

fpclass = static_cast<enum_fpclass_type>(static_cast<int>(f.fpclass));

unsigned elems = (std::min)(f.prec_elem, cpp_dec_float_elem_number);

std::copy(f.data.begin(), f.data.begin() + elems, data.begin());

std::fill(data.begin() + elems, data.end(), 0);

prec_elem = cpp_dec_float_elem_number;

return *this;

}

cpp_dec_float& operator=(boost::long_long_type v)

{

if(v < 0)

{

from_unsigned_long_long(1u - boost::ulong_long_type(v + 1)); // Avoid undefined behaviour in negation of minimum value for long long

negate();

}

else

from_unsigned_long_long(v);

return *this;

}

cpp_dec_float& operator=(boost::ulong_long_type v)

{

from_unsigned_long_long(v);

return *this;

}

cpp_dec_float& operator=(long double v);

cpp_dec_float& operator=(const char* v)

{

rd_string(v);

return *this;

}

cpp_dec_float& operator+=(const cpp_dec_float& v);

cpp_dec_float& operator-=(const cpp_dec_float& v);

cpp_dec_float& operator*=(const cpp_dec_float& v);

cpp_dec_float& operator/=(const cpp_dec_float& v);

cpp_dec_float& add_unsigned_long_long(const boost::ulong_long_type n)

{

cpp_dec_float t;

t.from_unsigned_long_long(n);

return *this += t;

}

cpp_dec_float& sub_unsigned_long_long(const boost::ulong_long_type n)

{

cpp_dec_float t;

t.from_unsigned_long_long(n);

return *this -= t;

}

cpp_dec_float& mul_unsigned_long_long(const boost::ulong_long_type n);

cpp_dec_float& div_unsigned_long_long(const boost::ulong_long_type n);

// Elementary primitives.

cpp_dec_float& calculate_inv ();

cpp_dec_float& calculate_sqrt();

void negate()

{

if(!iszero())

neg = !neg;

}

// Comparison functions

bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_NaN); }

bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_inf); }

bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_finite); }

bool iszero () const

{

return ((fpclass == cpp_dec_float_finite) && (data[0u] == 0u));

}

bool isone () const;

bool isint () const;

bool isneg () const { return neg; }

// Operators pre-increment and pre-decrement

cpp_dec_float& operator++()

{

return *this += one();

}

cpp_dec_float& operator--()

{

return *this -= one();

}

std::string str(boost::intmax_t digits, std::ios_base::fmtflags f)const;

int compare(const cpp_dec_float& v)const;

template <class V>

int compare(const V& v)const

{

cpp_dec_float<Digits10, ExponentType, Allocator> t;

t = v;

return compare(t);

}

void swap(cpp_dec_float& v)

{

data.swap(v.data);

std::swap(exp, v.exp);

std::swap(neg, v.neg);

std::swap(fpclass, v.fpclass);

std::swap(prec_elem, v.prec_elem);

}

double extract_double() const;

long double extract_long_double() const;

boost::long_long_type extract_signed_long_long() const;

boost::ulong_long_type extract_unsigned_long_long() const;

void extract_parts(double& mantissa, ExponentType& exponent) const;

cpp_dec_float extract_integer_part() const;

void precision(const boost::int32_t prec_digits)

{

if(prec_digits >= cpp_dec_float_total_digits10)

{

prec_elem = cpp_dec_float_elem_number;

}

else

{

const boost::int32_t elems = static_cast<boost::int32_t>( static_cast<boost::int32_t>( (prec_digits + (cpp_dec_float_elem_digits10 / 2)) / cpp_dec_float_elem_digits10)

+ static_cast<boost::int32_t>(((prec_digits % cpp_dec_float_elem_digits10) != 0) ? 1 : 0));

prec_elem = (std::min)(cpp_dec_float_elem_number, (std::max)(elems, static_cast<boost::int32_t>(2)));

}

}

static cpp_dec_float pow2(boost::long_long_type i);

ExponentType order()const

{

const bool bo_order_is_zero = ((!(isfinite)()) || (data[0] == static_cast<boost::uint32_t>(0u)));

//

// Binary search to find the order of the leading term:

//

ExponentType prefix = 0;

if(data[0] >= 100000UL)

{

if(data[0] >= 10000000UL)

{

if(data[0] >= 100000000UL)

{

if(data[0] >= 1000000000UL)

prefix = 9;

else

prefix = 8;

}

else

prefix = 7;

}

else

{

if(data[0] >= 1000000UL)

prefix = 6;

else

prefix = 5;

}

}

else

{

if(data[0] >= 1000UL)

{

if(data[0] >= 10000UL)

prefix = 4;

else

prefix = 3;

}

else

{

if(data[0] >= 100)

prefix = 2;

else if(data[0] >= 10)

prefix = 1;

}

}

return (bo_order_is_zero ? static_cast<ExponentType>(0) : static_cast<ExponentType>(exp + prefix));

}

template<class Archive>

void serialize(Archive & ar, const unsigned int /*version*/)

{

for(unsigned i = 0; i < data.size(); ++i)

ar & boost::serialization::make_nvp("digit", data[i]);

ar & boost::serialization::make_nvp("exponent", exp);

ar & boost::serialization::make_nvp("sign", neg);

ar & boost::serialization::make_nvp("class-type", fpclass);

ar & boost::serialization::make_nvp("precision", prec_elem);

}

private:

static bool data_elem_is_non_zero_predicate(const boost::uint32_t& d) { return (d != static_cast<boost::uint32_t>(0u)); }

static bool data_elem_is_non_nine_predicate(const boost::uint32_t& d) { return (d != static_cast<boost::uint32_t>(cpp_dec_float::cpp_dec_float_elem_mask - 1)); }

static bool char_is_nonzero_predicate(const char& c) { return (c != static_cast<char>('0')); }

void from_unsigned_long_long(const boost::ulong_long_type u);

int cmp_data(const array_type& vd) const;

static boost::uint32_t mul_loop_uv(boost::uint32_t* const u, const boost::uint32_t* const v, const boost::int32_t p);

static boost::uint32_t mul_loop_n (boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p);

static boost::uint32_t div_loop_n (boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p);

bool rd_string(const char* const s);

template <unsigned D, class ET, class A>

friend class cpp_dec_float;

};

template <unsigned Digits10, class ExponentType, class Allocator>

typename cpp_dec_float<Digits10, ExponentType, Allocator>::initializer cpp_dec_float<Digits10, ExponentType, Allocator>::init;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_radix;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_setting;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_limit_lo;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_limit_hi;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10;

template <unsigned Digits10, class ExponentType, class Allocator>

const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp;

template <unsigned Digits10, class ExponentType, class Allocator>

const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp;

template <unsigned Digits10, class ExponentType, class Allocator>

const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp10;

template <unsigned Digits10, class ExponentType, class Allocator>

const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp10;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_digits10;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_number_request;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_number;

template <unsigned Digits10, class ExponentType, class Allocator>

const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_mask;

template <unsigned Digits10, class ExponentType, class Allocator>

cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator+=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)

{

if((isnan)())

{

return *this;

}

if((isinf)())

{

if((v.isinf)() && (isneg() != v.isneg()))

{

*this = nan();

}

return *this;

}

if(iszero())

{

return operator=(v);

}

if((v.isnan)() || (v.isinf)())

{

*this = v;

return *this;

}

// Get the offset for the add/sub operation.

static const ExponentType max_delta_exp = static_cast<ExponentType>((cpp_dec_float_elem_number - 1) * cpp_dec_float_elem_digits10);

const ExponentType ofs_exp = static_cast<ExponentType>(exp - v.exp);

// Check if the operation is out of range, requiring special handling.

if(v.iszero() || (ofs_exp > max_delta_exp))

{

// Result is *this unchanged since v is negligible compared to *this.

return *this;

}

else if(ofs_exp < -max_delta_exp)

{

// Result is *this = v since *this is negligible compared to v.

return operator=(v);

}

// Do the add/sub operation.

typename array_type::iterator p_u = data.begin();

typename array_type::const_iterator p_v = v.data.begin();

bool b_copy = false;

const boost::int32_t ofs = static_cast<boost::int32_t>(static_cast<boost::int32_t>(ofs_exp) / cpp_dec_float_elem_digits10);

array_type n_data;

if(neg == v.neg)

{

// Add v to *this, where the data array of either *this or v

// might have to be treated with a positive, negative or zero offset.

// The result is stored in *this. The data are added one element

// at a time, each element with carry.

if(ofs >= static_cast<boost::int32_t>(0))

{

std::copy(v.data.begin(), v.data.end() - static_cast<size_t>(ofs), n_data.begin() + static_cast<size_t>(ofs));

std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(ofs), static_cast<boost::uint32_t>(0u));

p_v = n_data.begin();

}

else

{

std::copy(data.begin(), data.end() - static_cast<size_t>(-ofs), n_data.begin() + static_cast<size_t>(-ofs));

std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(-ofs), static_cast<boost::uint32_t>(0u));

p_u = n_data.begin();

b_copy = true;

}

// Addition algorithm

boost::uint32_t carry = static_cast<boost::uint32_t>(0u);

for(boost::int32_t j = static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)); j >= static_cast<boost::int32_t>(0); j--)

{

boost::uint32_t t = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(p_u[j] + p_v[j]) + carry);

carry = t / static_cast<boost::uint32_t>(cpp_dec_float_elem_mask);

p_u[j] = static_cast<boost::uint32_t>(t - static_cast<boost::uint32_t>(carry * static_cast<boost::uint32_t>(cpp_dec_float_elem_mask)));

}

if(b_copy)

{

data = n_data;

exp = v.exp;

}

// There needs to be a carry into the element -1 of the array data

if(carry != static_cast<boost::uint32_t>(0u))

{

std::copy_backward(data.begin(), data.end() - static_cast<std::size_t>(1u), data.end());

data[0] = carry;

exp += static_cast<ExponentType>(cpp_dec_float_elem_digits10);

}

}

else

{

// Subtract v from *this, where the data array of either *this or v

// might have to be treated with a positive, negative or zero offset.

if((ofs > static_cast<boost::int32_t>(0))

|| ( (ofs == static_cast<boost::int32_t>(0))

&& (cmp_data(v.data) > static_cast<boost::int32_t>(0)))

)

{

// In this case, |u| > |v| and ofs is positive.

// Copy the data of v, shifted down to a lower value

// into the data array m_n. Set the operand pointer p_v

// to point to the copied, shifted data m_n.

std::copy(v.data.begin(), v.data.end() - static_cast<size_t>(ofs), n_data.begin() + static_cast<size_t>(ofs));

std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(ofs), static_cast<boost::uint32_t>(0u));

p_v = n_data.begin();

}

else

{

if(ofs != static_cast<boost::int32_t>(0))

{

// In this case, |u| < |v| and ofs is negative.

// Shift the data of u down to a lower value.

std::copy_backward(data.begin(), data.end() - static_cast<size_t>(-ofs), data.end());

std::fill(data.begin(), data.begin() + static_cast<size_t>(-ofs), static_cast<boost::uint32_t>(0u));

}

// Copy the data of v into the data array n_data.

// Set the u-pointer p_u to point to m_n and the

// operand pointer p_v to point to the shifted

// data m_data.

n_data = v.data;

p_u = n_data.begin();

p_v = data.begin();

b_copy = true;

}

boost::int32_t j;

// Subtraction algorithm

boost::int32_t borrow = static_cast<boost::int32_t>(0);

for(j = static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)); j >= static_cast<boost::int32_t>(0); j--)

{

boost::int32_t t = static_cast<boost::int32_t>(static_cast<boost::int32_t>( static_cast<boost::int32_t>(p_u[j])

- static_cast<boost::int32_t>(p_v[j])) - borrow);

// Underflow? Borrow?

if(t < static_cast<boost::int32_t>(0))

{

// Yes, underflow and borrow

t += static_cast<boost::int32_t>(cpp_dec_float_elem_mask);

borrow = static_cast<boost::int32_t>(1);

}

else

{

borrow = static_cast<boost::int32_t>(0);

}

p_u[j] = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(t) % static_cast<boost::uint32_t>(cpp_dec_float_elem_mask));

}

if(b_copy)

{

data = n_data;

exp = v.exp;

neg = v.neg;

}

// Is it necessary to justify the data?

const typename array_type::const_iterator first_nonzero_elem = std::find_if(data.begin(), data.end(), data_elem_is_non_zero_predicate);

if(first_nonzero_elem != data.begin())

{

if(first_nonzero_elem == data.end())

{

// This result of the subtraction is exactly zero.

// Reset the sign and the exponent.

neg = false;

exp = static_cast<ExponentType>(0);

}

else

{

// Justify the data

const std::size_t sj = static_cast<std::size_t>(std::distance<typename array_type::const_iterator>(data.begin(), first_nonzero_elem));

std::copy(data.begin() + static_cast<std::size_t>(sj), data.end(), data.begin());

std::fill(data.end() - sj, data.end(), static_cast<boost::uint32_t>(0u));

exp -= static_cast<ExponentType>(sj * static_cast<std::size_t>(cpp_dec_float_elem_digits10));

}

}

}

// Handle underflow.

if(iszero())

return (*this = zero());

// Check for potential overflow.

const bool b_result_might_overflow = (exp >= static_cast<ExponentType>(cpp_dec_float_max_exp10));

// Handle overflow.

if(b_result_might_overflow)

{

const bool b_result_is_neg = neg;

neg = false;

if(compare((cpp_dec_float::max)()) > 0)

*this = inf();

neg = b_result_is_neg;

}

return *this;

}

template <unsigned Digits10, class ExponentType, class Allocator>

cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator-=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)

{

// Use *this - v = -(-*this + v).

negate();

*this += v;

negate();

return *this;

}

template <unsigned Digits10, class ExponentType, class Allocator>

cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator*=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)

{

// Evaluate the sign of the result.

const bool b_result_is_neg = (neg != v.neg);

// Artificially set the sign of the result to be positive.

neg = false;

// Handle special cases like zero, inf and NaN.

const bool b_u_is_inf = (isinf)();

const bool b_v_is_inf = (v.isinf)();

const bool b_u_is_zero = iszero();

const bool b_v_is_zero = v.iszero();

if( ((isnan)() || (v.isnan)())

|| (b_u_is_inf && b_v_is_zero)

|| (b_v_is_inf && b_u_is_zero)

)

{

*this = nan();

return *this;

}

if(b_u_is_inf || b_v_is_inf)

{

*this = inf();

if(b_result_is_neg)

negate();

return *this;

}

if(b_u_is_zero || b_v_is_zero)

{

return *this = zero();

}

// Check for potential overflow or underflow.

const bool b_result_might_overflow = ((exp + v.exp) >= static_cast<ExponentType>(cpp_dec_float_max_exp10));

const bool b_result_might_underflow = ((exp + v.exp) <= static_cast<ExponentType>(cpp_dec_float_min_exp10));

// Set the exponent of the result.

exp += v.exp;

const boost::int32_t prec_mul = (std::min)(prec_elem, v.prec_elem);

const boost::uint32_t carry = mul_loop_uv(data.data(), v.data.data(), prec_mul);

// Handle a potential carry.

if(carry != static_cast<boost::uint32_t>(0u))

{

exp += cpp_dec_float_elem_digits10;

// Shift the result of the multiplication one element to the right...

std::copy_backward(data.begin(),

data.begin() + static_cast<std::size_t>(prec_elem - static_cast<boost::int32_t>(1)),

data.begin() + static_cast<std::size_t>(prec_elem));

// ... And insert the carry.

data.front() = carry;

}

// Handle overflow.

if(b_result_might_overflow && (compare((cpp_dec_float::max)()) > 0))

{

*this = inf();

}

// Handle underflow.

if(b_result_might_underflow && (compare((cpp_dec_float::min)()) < 0))

{

*this = zero();

return *this;

}

// Set the sign of the result.

neg = b_result_is_neg;

return *this;

}

template <unsigned Digits10, class ExponentType, class Allocator>

cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator/=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)

{

if(iszero())

{

if((v.isnan)())

{

return *this = v;

}

else if(v.iszero())

{

return *this = nan();

}

}

const bool u_and_v_are_finite_and_identical = ( (isfinite)()

&& (fpclass == v.fpclass)

&& (exp == v.exp)

&& (cmp_data(v.data) == static_cast<boost::int32_t>(0)));

if(u_and_v_are_finite_and_identical)

{

if(neg != v.neg)

{

*this = one();

negate();

}

else