code:
///////////////////////////////////////////////////////////////
//  Copyright 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_

#ifndef BOOST_MP_FLOAT128_HPP
#define BOOST_MP_FLOAT128_HPP

#include <boost/config.hpp>
#include <boost/scoped_array.hpp>
#include <boost/functional/hash.hpp>
#include <boost/multiprecision/number.hpp>

#if defined(BOOST_INTEL) && !defined(BOOST_MP_USE_FLOAT128) && !defined(BOOST_MP_USE_QUAD)
#  if defined(BOOST_INTEL_CXX_VERSION) && (BOOST_INTEL_CXX_VERSION >= 1310) && defined(__GNUC__)
#    if (__GNUC__ > 4) || ((__GNUC__ == 4) && (__GNUC_MINOR__ >= 6))
#      define BOOST_MP_USE_FLOAT128
#    endif
#  endif

#  ifndef BOOST_MP_USE_FLOAT128
#    define BOOST_MP_USE_QUAD
#  endif
#endif

#if defined(__GNUC__) && !defined(BOOST_MP_USE_FLOAT128) && !defined(BOOST_MP_USE_QUAD)
#  define BOOST_MP_USE_FLOAT128
#endif

#if !defined(BOOST_MP_USE_FLOAT128) && !defined(BOOST_MP_USE_QUAD)
#  error "Sorry compiler is neither GCC, not Intel, don't know how to configure this header."
#endif
#if defined(BOOST_MP_USE_FLOAT128) && defined(BOOST_MP_USE_QUAD)
#  error "Oh dear, both BOOST_MP_USE_FLOAT128 and BOOST_MP_USE_QUAD are defined, which one should I be using?"
#endif

#if defined(BOOST_MP_USE_FLOAT128)

extern "C" {
#include <quadmath.h>
}

typedef __float128 float128_type;

#elif defined(BOOST_MP_USE_QUAD)

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

typedef _Quad float128_type;

extern "C" {
_Quad __ldexpq(_Quad, int);
_Quad __frexpq(_Quad, int*);
_Quad __fabsq(_Quad);
_Quad __floorq(_Quad);
_Quad __ceilq(_Quad);
_Quad __sqrtq(_Quad);
_Quad __truncq(_Quad);
_Quad __expq(_Quad);
_Quad __powq(_Quad, _Quad);
_Quad __logq(_Quad);
_Quad __log10q(_Quad);
_Quad __sinq(_Quad);
_Quad __cosq(_Quad);
_Quad __tanq(_Quad);
_Quad __asinq(_Quad);
_Quad __acosq(_Quad);
_Quad __atanq(_Quad);
_Quad __sinhq(_Quad);
_Quad __coshq(_Quad);
_Quad __tanhq(_Quad);
_Quad __fmodq(_Quad, _Quad);
_Quad __atan2q(_Quad, _Quad);

#define ldexpq __ldexpq
#define frexpq __frexpq
#define fabsq __fabsq
#define floorq __floorq
#define ceilq __ceilq
#define sqrtq __sqrtq
#define truncq __truncq
#define expq __expq
#define powq __powq
#define logq __logq
#define log10q __log10q
#define sinq __sinq
#define cosq __cosq
#define tanq __tanq
#define asinq __asinq
#define acosq __acosq
#define atanq __atanq
#define sinhq __sinhq
#define coshq __coshq
#define tanhq __tanhq
#define fmodq __fmodq
#define atan2q __atan2q
}

inline _Quad isnanq(_Quad v)
{
   return v != v;
}
inline _Quad isinfq(_Quad v)
{
   return __fabsq(v) > 1.18973149535723176508575932662800702e4932Q;
}

#endif

namespace boost{
namespace multiprecision{
namespace backends{

struct float128_backend;

}

using backends::float128_backend;

template<>
struct number_category<backends::float128_backend> : public mpl::int_<number_kind_floating_point> {};
template<>
struct number_category<float128_type> : public mpl::int_<number_kind_floating_point> {};

typedef number<float128_backend, et_off> float128;

namespace backends{

struct float128_backend
{
   typedef mpl::list<signed char, short, int, long, boost::long_long_type>   signed_types;
   typedef mpl::list<unsigned char, unsigned short,
      unsigned int, unsigned long, boost::ulong_long_type>           unsigned_types;
   typedef mpl::list<float, double, long double>                 float_types;
   typedef int                                                   exponent_type;

private:
   float128_type m_value;
public:
   BOOST_CONSTEXPR float128_backend() BOOST_NOEXCEPT : m_value(0) {}
   BOOST_CONSTEXPR float128_backend(const float128_backend& o) BOOST_NOEXCEPT : m_value(o.m_value) {}
   float128_backend& operator = (const float128_backend& o) BOOST_NOEXCEPT
   {
      m_value = o.m_value;
      return *this;
   }
   template <class T>
   BOOST_CONSTEXPR float128_backend(const T& i, const typename enable_if_c<is_convertible<T, float128_type>::value>::type* = 0) BOOST_NOEXCEPT_IF(noexcept(std::declval<float128_type&>() = std::declval<const T&>()))
      : m_value(i) {}
   template <class T>
   typename enable_if_c<is_arithmetic<T>::value || is_convertible<T, float128_type>::value, float128_backend&>::type operator = (const T& i) BOOST_NOEXCEPT_IF(noexcept(std::declval<float128_type&>() = std::declval<const T&>()))
   {
      m_value = i;
      return *this;
   }
   float128_backend(long double const& f)
   {
      if(boost::math::isinf(f))
         m_value = (f < 0) ? -1.0Q / 0.0Q : 1.0Q / 0.0Q;
      else
         m_value = f;
   }
   float128_backend& operator=(long double const& f)
   {
      if(boost::math::isinf(f))
         m_value = (f < 0) ? -1.0Q / 0.0Q : 1.0Q / 0.0Q;
      else
         m_value = f;
      return *this;
   }
   float128_backend& operator = (const char* s)
   {
#ifndef BOOST_MP_USE_QUAD
      char* p_end;
      m_value = strtoflt128(s, &p_end);
      if(p_end - s != (std::ptrdiff_t)std::strlen(s))
      {
         BOOST_THROW_EXCEPTION(std::runtime_error("Unable to interpret input string as a floating point value"));
      }
#else
      boost::multiprecision::detail::convert_from_string(*this, s);
#endif
      return *this;
   }
   void swap(float128_backend& o) BOOST_NOEXCEPT
   {
      std::swap(m_value, o.value());
   }
   std::string str(std::streamsize digits, std::ios_base::fmtflags f)const
   {
#ifndef BOOST_MP_USE_QUAD
      char buf[100];
      boost::scoped_array<char> buf2;
      std::string format = "%";
      if(f & std::ios_base::showpos)
         format += "+";
      if(f & std::ios_base::showpoint)
         format += "#";
      format += ".*";
      if(digits == 0)
         digits = 36;
      format += "Q";
      if(f & std::ios_base::scientific)
         format += "e";
      else if(f & std::ios_base::fixed)
         format += "f";
      else
         format += "g";

      int v = quadmath_snprintf (buf, 100, format.c_str(), digits, m_value);

      if((v < 0) || (v >= 99))
      {
         int v_max = v;
         buf2.reset(new char[v+3]);
         v = quadmath_snprintf (&buf2[0], v_max + 3, format.c_str(), digits, m_value);
         if(v >= v_max + 3)
         {
            BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of float128_type failed."));
         }
         return &buf2[0];
      }
      return buf;
#else
      return boost::multiprecision::detail::convert_to_string(*this, digits ? digits : 37, f);
#endif
   }
   void negate() BOOST_NOEXCEPT
   {
      m_value = -m_value;
   }
   int compare(const float128_backend& o)const
   {
      return m_value == o.m_value ? 0 : m_value < o.m_value ? -1 : 1;
   }
   template <class T>
   int compare(const T& i)const
   {
      return m_value == i ? 0 : m_value < i ? -1 : 1;
   }
   float128_type& value()
   {
      return m_value;
   }
   const float128_type& value()const
   {
      return m_value;
   }
};

inline void eval_add(float128_backend& result, const float128_backend& a)
{
   result.value() += a.value();
}
template <class A>
inline void eval_add(float128_backend& result, const A& a)
{
   result.value() += a;
}
inline void eval_subtract(float128_backend& result, const float128_backend& a)
{
   result.value() -= a.value();
}
template <class A>
inline void eval_subtract(float128_backend& result, const A& a)
{
   result.value() -= a;
}
inline void eval_multiply(float128_backend& result, const float128_backend& a)
{
   result.value() *= a.value();
}
template <class A>
inline void eval_multiply(float128_backend& result, const A& a)
{
   result.value() *= a;
}
inline void eval_divide(float128_backend& result, const float128_backend& a)
{
   result.value() /= a.value();
}
template <class A>
inline void eval_divide(float128_backend& result, const A& a)
{
   result.value() /= a;
}

inline void eval_add(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = a.value() + b.value();
}
template <class A>
inline void eval_add(float128_backend& result, const float128_backend& a, const A& b)
{
   result.value() = a.value() + b;
}
inline void eval_subtract(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = a.value() - b.value();
}
template <class A>
inline void eval_subtract(float128_backend& result, const float128_backend& a, const A& b)
{
   result.value() = a.value() - b;
}
template <class A>
inline void eval_subtract(float128_backend& result, const A& a, const float128_backend& b)
{
   result.value() = a - b.value();
}
inline void eval_multiply(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = a.value() * b.value();
}
template <class A>
inline void eval_multiply(float128_backend& result, const float128_backend& a, const A& b)
{
   result.value() = a.value() * b;
}
inline void eval_divide(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = a.value() / b.value();
}

template <class R>
inline void eval_convert_to(R* result, const float128_backend& val)
{
   *result = static_cast<R>(val.value());
}

inline void eval_frexp(float128_backend& result, const float128_backend& arg, int* exp)
{
   result.value() = frexpq(arg.value(), exp);
}

inline void eval_ldexp(float128_backend& result, const float128_backend& arg, int exp)
{
   result.value() = ldexpq(arg.value(), exp);
}

inline void eval_floor(float128_backend& result, const float128_backend& arg)
{
   result.value() = floorq(arg.value());
}
inline void eval_ceil(float128_backend& result, const float128_backend& arg)
{
   result.value() = ceilq(arg.value());
}
inline void eval_sqrt(float128_backend& result, const float128_backend& arg)
{
   result.value() = sqrtq(arg.value());
}
inline int eval_fpclassify(const float128_backend& arg)
{
   if(isnanq(arg.value()))
      return FP_NAN;
   else if(isinfq(arg.value()))
      return FP_INFINITE;
   else if(arg.value() == 0)
      return FP_ZERO;

   float128_backend t(arg);
   if(t.value() < 0)
      t.negate();
   if(t.value() < 3.36210314311209350626267781732175260e-4932Q)
      return FP_SUBNORMAL;
   return FP_NORMAL;
}

inline void eval_increment(float128_backend& arg)
{
   ++arg.value();
}
inline void eval_decrement(float128_backend& arg)
{
   --arg.value();
}

/*********************************************************************
*
* abs/fabs:
*
*********************************************************************/

inline void eval_abs(float128_backend& result, const float128_backend& arg)
{
   result.value() = fabsq(arg.value());
}
inline void eval_fabs(float128_backend& result, const float128_backend& arg)
{
   result.value() = fabsq(arg.value());
}

/*********************************************************************
*
* Floating point functions:
*
*********************************************************************/

inline void eval_trunc(float128_backend& result, const float128_backend& arg)
{
   if(isnanq(arg.value()) || isinfq(arg.value()))
   {
      result = boost::math::policies::raise_rounding_error(
            "boost::multiprecision::trunc<%1%>(%1%)", 0,
            number<float128_backend, et_off>(arg),
            number<float128_backend, et_off>(arg),
            boost::math::policies::policy<>()).backend();
      return;
   }
   result.value() = truncq(arg.value());
}
/*
//
// This doesn't actually work... rely on our own default version instead.
//
inline void eval_round(float128_backend& result, const float128_backend& arg)
{
   if(isnanq(arg.value()) || isinf(arg.value()))
   {
      result = boost::math::policies::raise_rounding_error(
            "boost::multiprecision::trunc<%1%>(%1%)", 0,
            number<float128_backend, et_off>(arg),
            number<float128_backend, et_off>(arg),
            boost::math::policies::policy<>()).backend();
      return;
   }
   result.value() = roundq(arg.value());
}
*/

inline void eval_exp(float128_backend& result, const float128_backend& arg)
{
   result.value() = expq(arg.value());
}
inline void eval_log(float128_backend& result, const float128_backend& arg)
{
   result.value() = logq(arg.value());
}
inline void eval_log10(float128_backend& result, const float128_backend& arg)
{
   result.value() = log10q(arg.value());
}
inline void eval_sin(float128_backend& result, const float128_backend& arg)
{
   result.value() = sinq(arg.value());
}
inline void eval_cos(float128_backend& result, const float128_backend& arg)
{
   result.value() = cosq(arg.value());
}
inline void eval_tan(float128_backend& result, const float128_backend& arg)
{
   result.value() = tanq(arg.value());
}
inline void eval_asin(float128_backend& result, const float128_backend& arg)
{
   result.value() = asinq(arg.value());
}
inline void eval_acos(float128_backend& result, const float128_backend& arg)
{
   result.value() = acosq(arg.value());
}
inline void eval_atan(float128_backend& result, const float128_backend& arg)
{
   result.value() = atanq(arg.value());
}
inline void eval_sinh(float128_backend& result, const float128_backend& arg)
{
   result.value() = sinhq(arg.value());
}
inline void eval_cosh(float128_backend& result, const float128_backend& arg)
{
   result.value() = coshq(arg.value());
}
inline void eval_tanh(float128_backend& result, const float128_backend& arg)
{
   result.value() = tanhq(arg.value());
}
inline void eval_fmod(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = fmodq(a.value(), b.value());
}
inline void eval_pow(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = powq(a.value(), b.value());
}
inline void eval_atan2(float128_backend& result, const float128_backend& a, const float128_backend& b)
{
   result.value() = atan2q(a.value(), b.value());
}
inline void eval_multiply_add(float128_backend& result, const float128_backend& a, const float128_backend& b, const float128_backend& c)
{
   result.value() = fmaq(a.value(), b.value(), c.value());
}

inline std::size_t hash_value(const float128_backend& val)
{
   return  boost::hash_value(static_cast<double>(val.value()));
}

} // namespace backends

   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> asinh BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return asinhq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> acosh BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return acoshq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> atanh BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return atanhq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> cbrt BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return cbrtq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> erf BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return erfq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> erfc BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return erfcq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> expm1 BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return expm1q(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> lgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return lgammaq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> tgamma BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return tgammaq(arg.backend().value());
   }
   template<boost::multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<float128_backend, ExpressionTemplates> log1p BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<float128_backend, ExpressionTemplates>& arg)
   {
      return log1pq(arg.backend().value());
   }

   template <multiprecision::expression_template_option ExpressionTemplates>
   inline int signbit BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& arg)
   {
      return ::signbitq(arg.backend().value());
   }

   template <multiprecision::expression_template_option ExpressionTemplates>
   inline boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> copysign BOOST_PREVENT_MACRO_SUBSTITUTION(const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& a, const boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates>& b)
   {
      return ::copysignq(a.backend().value(), b.backend().value());
   }

   inline void eval_remainder(float128_backend& result, const float128_backend& a, const float128_backend& b)
   {
      result.value() = remainderq(a.value(), b.value());
   }
   inline void eval_remainder(float128_backend& result, const float128_backend& a, const float128_backend& b, int* pi)
   {
      result.value() = remquoq(a.value(), b.value(), pi);
   }

} // namespace multiprecision

namespace math {

   using boost::multiprecision::signbit;
   using boost::multiprecision::copysign;

} // namespace math

} // namespace boost

namespace boost{
namespace archive{

class binary_oarchive;
class binary_iarchive;

}

namespace serialization{ namespace float128_detail{

template <class Archive>
void do_serialize(Archive& ar, boost::multiprecision::backends::float128_backend& val, const mpl::false_&, const mpl::false_&)
{
   // saving
   // non-binary
   std::string s(val.str(0, std::ios_base::scientific));
   ar & s;
}
template <class Archive>
void do_serialize(Archive& ar, boost::multiprecision::backends::float128_backend& val, const mpl::true_&, const mpl::false_&)
{
   // loading
   // non-binary
   std::string s;
   ar & s;
   val = s.c_str();
}

template <class Archive>
void do_serialize(Archive& ar, boost::multiprecision::backends::float128_backend& val, const mpl::false_&, const mpl::true_&)
{
   // saving
   // binary
   ar.save_binary(&val, sizeof(val));
}
template <class Archive>
void do_serialize(Archive& ar, boost::multiprecision::backends::float128_backend& val, const mpl::true_&, const mpl::true_&)
{
   // loading
   // binary
   ar.load_binary(&val, sizeof(val));
}

} // detail

template <class Archive>
void serialize(Archive& ar, boost::multiprecision::backends::float128_backend& val, unsigned int /*version*/)
{
   typedef typename Archive::is_loading load_tag;
   typedef typename mpl::bool_<boost::is_same<Archive, boost::archive::binary_oarchive>::value || boost::is_same<Archive, boost::archive::binary_iarchive>::value> binary_tag;

   float128_detail::do_serialize(ar, val, load_tag(), binary_tag());
}

} // namepsace archive

} // namespace boost

namespace std{

template <boost::multiprecision::expression_template_option ExpressionTemplates>
class numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >
{
   typedef boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> number_type;
public:
   BOOST_STATIC_CONSTEXPR bool is_specialized = true;
   static number_type (min)() BOOST_NOEXCEPT { return 3.36210314311209350626267781732175260e-4932Q; }
   static number_type (max)() BOOST_NOEXCEPT { return 1.18973149535723176508575932662800702e4932Q; }
   static number_type lowest() BOOST_NOEXCEPT { return -(max)(); }
   BOOST_STATIC_CONSTEXPR int digits = 113;
   BOOST_STATIC_CONSTEXPR int digits10 = 33;
   BOOST_STATIC_CONSTEXPR int max_digits10 = 36;
   BOOST_STATIC_CONSTEXPR bool is_signed = true;
   BOOST_STATIC_CONSTEXPR bool is_integer = false;
   BOOST_STATIC_CONSTEXPR bool is_exact = false;
   BOOST_STATIC_CONSTEXPR int radix = 2;
   static number_type epsilon() { return 1.92592994438723585305597794258492732e-34Q; }
   static number_type round_error() { return 0.5; }
   BOOST_STATIC_CONSTEXPR int min_exponent = -16381;
   BOOST_STATIC_CONSTEXPR int min_exponent10 = min_exponent * 301L / 1000L;
   BOOST_STATIC_CONSTEXPR int max_exponent = 16384;
   BOOST_STATIC_CONSTEXPR int max_exponent10 = max_exponent * 301L / 1000L;
   BOOST_STATIC_CONSTEXPR bool has_infinity = true;
   BOOST_STATIC_CONSTEXPR bool has_quiet_NaN = true;
   BOOST_STATIC_CONSTEXPR bool has_signaling_NaN = false;
   BOOST_STATIC_CONSTEXPR float_denorm_style has_denorm = denorm_present;
   BOOST_STATIC_CONSTEXPR bool has_denorm_loss = true;
   static number_type infinity() { return 1.0q / 0.0q; }
   static number_type quiet_NaN() { return number_type("nan"); }
   static number_type signaling_NaN() { return 0; }
   static number_type denorm_min() { return 6.475175119438025110924438958227646552e-4966Q; }
   BOOST_STATIC_CONSTEXPR bool is_iec559 = true;
   BOOST_STATIC_CONSTEXPR bool is_bounded = false;
   BOOST_STATIC_CONSTEXPR bool is_modulo = false;
   BOOST_STATIC_CONSTEXPR bool traps = false;
   BOOST_STATIC_CONSTEXPR bool tinyness_before = false;
   BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest;
};

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_specialized;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::digits;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::digits10;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::max_digits10;

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_signed;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_integer;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_exact;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::radix;

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::min_exponent;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::max_exponent;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::min_exponent10;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::max_exponent10;

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::has_infinity;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::has_quiet_NaN;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::has_signaling_NaN;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::has_denorm_loss;

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_iec559;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_bounded;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::is_modulo;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::traps;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::tinyness_before;

template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::round_style;
template <boost::multiprecision::expression_template_option ExpressionTemplates>
BOOST_CONSTEXPR_OR_CONST float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::backends::float128_backend, ExpressionTemplates> >::has_denorm;

} // namespace std

#endif


code:
// Boost.Units - A C++ library for zero-overhead dimensional analysis and
// unit/quantity manipulation and conversion
//
// Copyright (C) 2003-2008 Matthias Christian Schabel
// Copyright (C) 2008 Steven Watanabe
//
// 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_UNITS_GET_DIMENSION_HPP
#define BOOST_UNITS_GET_DIMENSION_HPP

///
/// \file
/// \brief Get the dimension of a unit, absolute unit and quantity.
/// \details
///

#include <boost/units/units_fwd.hpp>

namespace boost {

namespace units {

template<class T>
struct get_dimension {};

/// Get the dimension of a unit.
template<class Dim,class System>
struct get_dimension< unit<Dim,System> >
{
    typedef Dim type;
};

/// Get the dimension of an absolute unit.
template<class Unit>
struct get_dimension< absolute<Unit> >
{
    typedef typename get_dimension<Unit>::type  type;
};

/// Get the dimension of a quantity.
template<class Unit,class Y>
struct get_dimension< quantity<Unit,Y> >
{
    typedef typename get_dimension<Unit>::type  type;
};

} // namespace units

} // namespace boost

#endif // BOOST_UNITS_GET_DIMENSION_HPP

license: mit


code:
/*
qra64.c
Encoding/decoding functions for the QRA64 mode

(c) 2016 - Nico Palermo, IV3NWV

-------------------------------------------------------------------------------

   qracodes is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
   qracodes is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with qracodes source distribution.
   If not, see <http://www.gnu.org/licenses/>.

-----------------------------------------------------------------------------

QRA code used in this sowftware release:

QRA13_64_64_IRR_E: K=13 N=64 Q=64 irregular QRA code (defined in
qra13_64_64_irr_e.h /.c)

Codes with K=13 are designed to include a CRC as the 13th information symbol
and improve the code UER (Undetected Error Rate).
The CRC symbol is not sent along the channel (the codes are punctured) and the
resulting code is a (12,63) code
*/
//----------------------------------------------------------------------------

#include <stdlib.h>
#include <string.h>

#include "qra64.h"
#include "../qracodes/qracodes.h"
#include "../qracodes/qra13_64_64_irr_e.h"
#include "../qracodes/pdmath.h"
#include "../qracodes/normrnd.h"

// Code parameters of the QRA64 mode
#define QRA64_CODE  qra_13_64_64_irr_e
#define QRA64_NMSG  218        // Must much value indicated in QRA64_CODE.NMSG

#define QRA64_KC (QRA64_K+1)   // Information symbols (crc included)
#define QRA64_NC (QRA64_N+1)   // Codeword length (as defined in the code)
#define QRA64_NITER 100	       // max number of iterations per decode

// static functions declarations ----------------------------------------------
static int  calc_crc6(const int *x, int sz);
static void ix_mask(float *dst, const float *src, const int *mask,
		    const int *x);
static int qra64_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix);
static int qra64_do_decode(int *x, const float *pix, const int *ap_mask,
			    const int *ap_x);
static float qra64_fastfading_estim_noise_std(
				const float *rxen,
				const float esnometric,
				const int submode);

static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxen,
				const float *hptr,
				const int    hlen,
				const float sigma,
				const float EsNoMetric,
				const int submode);

static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxen,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode);

// a-priori information masks for fields in JT65-like msgs --------------------

// when defined limits the AP masks to reduce the false decode rate
#define LIMIT_AP_MASKS

#ifdef LIMIT_AP_MASKS
#define MASK_CQQRZ      0xFFFFFFC
#define MASK_CALL1      0xFFFFFFC
#define MASK_CALL2      0xFFFFFFC
#define MASK_GRIDFULL	0x3FFC
#define MASK_GRIDFULL12	0x3FFC
#define MASK_GRIDBIT	0x8000
#else
#define MASK_CQQRZ      0xFFFFFFC
#define MASK_CALL1      0xFFFFFFF
#define MASK_CALL2      0xFFFFFFF
#define MASK_GRIDFULL	0xFFFF
#define MASK_GRIDFULL12	0x3FFC
#define MASK_GRIDBIT	0x8000	  // b[15] is 1 for free text, 0 otherwise
#endif

// ----------------------------------------------------------------------------

qra64codec *qra64_init(int flags)
{

  // Eb/No value for which we optimize the decoder metric
  const float EbNodBMetric = 2.8f;
  const float EbNoMetric   = (float)pow(10,EbNodBMetric/10);
  const float R = 1.0f*(QRA64_KC)/(QRA64_NC);

  qra64codec *pcodec = (qra64codec*)malloc(sizeof(qra64codec));

  if (!pcodec)
    return 0;	// can't allocate memory

  pcodec->decEsNoMetric   = 1.0f*QRA64_m*R*EbNoMetric;
  pcodec->apflags			= flags;

  memset(pcodec->apmsg_set,0,APTYPE_SIZE*sizeof(int));

  if (flags==QRA_NOAP)
    return pcodec;

  // for QRA_USERAP and QRA_AUTOAP modes we always enable [CQ/QRZ ? ?] mgs look-up.
  // encode CQ/QRZ AP messages
  // NOTE: Here we handle only CQ and QRZ msgs.
  // 'CQ nnn', 'CQ DX' and 'DE' msgs will be handled by the decoder
  // as messages with no a-priori knowledge
  qra64_apset(pcodec, CALL_CQ, 0, GRID_BLANK, APTYPE_CQQRZ);

  // initialize masks for decoding with a-priori information
  encodemsg_jt65(pcodec->apmask_cqqrz,     MASK_CQQRZ, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cqqrz_ooo, MASK_CQQRZ, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1,     MASK_CALL1, 0, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_ooo, MASK_CALL1, 0, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call2,     0, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call2_ooo, 0, MASK_CALL2, MASK_GRIDFULL);
  encodemsg_jt65(pcodec->apmask_call1_call2,     MASK_CALL1,MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_call1_call2_grid,MASK_CALL1,MASK_CALL2, MASK_GRIDFULL12);
  encodemsg_jt65(pcodec->apmask_cq_call2,     MASK_CQQRZ, MASK_CALL2, MASK_GRIDBIT);
  encodemsg_jt65(pcodec->apmask_cq_call2_ooo, MASK_CQQRZ, MASK_CALL2, MASK_GRIDFULL12);

  return pcodec;
}

void qra64_close(qra64codec *pcodec)
{
	free(pcodec);
}

int qra64_apset(qra64codec *pcodec, const int mycall, const int hiscall, const int grid, const int aptype)
{
// Set decoder a-priori knowledge accordingly to the type of the message to look up for
// arguments:
//		pcodec    = pointer to a qra64codec data structure as returned by qra64_init
//		mycall    = mycall to look for
//		hiscall   = hiscall to look for
//		grid      = grid to look for
//		aptype    = define and masks the type of AP to be set accordingly to the following:
//			APTYPE_CQQRZ     set [cq/qrz ?       ?/blank]
//			APTYPE_MYCALL    set [mycall ?       ?/blank]
//			APTYPE_HISCALL   set [?      hiscall ?/blank]
//			APTYPE_BOTHCALLS set [mycall hiscall ?]
//			APTYPE_FULL		 set [mycall hiscall grid]
//			APTYPE_CQHISCALL set [cq/qrz hiscall ?/blank] and [cq/qrz hiscall grid]
// returns:
//		0   on success
//      -1  when qra64_init was called with the QRA_NOAP flag
//		-2  invalid apytpe

	if (pcodec->apflags==QRA_NOAP)
		return -1;

	switch (aptype) {
		case APTYPE_CQQRZ:
			encodemsg_jt65(pcodec->apmsg_cqqrz,  CALL_CQ, 0, GRID_BLANK);
			break;
		case APTYPE_MYCALL:
			encodemsg_jt65(pcodec->apmsg_call1,  mycall,  0, GRID_BLANK);
			break;
		case APTYPE_HISCALL:
			encodemsg_jt65(pcodec->apmsg_call2,  0, hiscall, GRID_BLANK);
			break;
		case APTYPE_BOTHCALLS:
			encodemsg_jt65(pcodec->apmsg_call1_call2,  mycall, hiscall, GRID_BLANK);
			break;
		case APTYPE_FULL:
			encodemsg_jt65(pcodec->apmsg_call1_call2_grid,  mycall, hiscall, grid);
			break;
		case APTYPE_CQHISCALL:
			encodemsg_jt65(pcodec->apmsg_cq_call2,      CALL_CQ, hiscall, GRID_BLANK);
			encodemsg_jt65(pcodec->apmsg_cq_call2_grid, CALL_CQ, hiscall, grid);
			break;
		default:
			return -2;	// invalid ap type
		}

	  pcodec->apmsg_set[aptype]=1;	// signal the decoder to look-up for the specified type

	  return 0;
}
void qra64_apdisable(qra64codec *pcodec, const int aptype)
{
	if (pcodec->apflags==QRA_NOAP)
		return;

	if (aptype<APTYPE_CQQRZ || aptype>=APTYPE_SIZE)
		return;

	pcodec->apmsg_set[aptype] = 0;	//  signal the decoder not to look-up to the specified type
}

void qra64_encode(qra64codec *pcodec, int *y, const int *x)
{
  int encx[QRA64_KC];	// encoder input buffer
  int ency[QRA64_NC];	// encoder output buffer

  int hiscall,mycall,grid;

  memcpy(encx,x,QRA64_K*sizeof(int));		// Copy input to encoder buffer
  encx[QRA64_K]=calc_crc6(encx,QRA64_K);	// Compute and add crc symbol
  qra_encode(&QRA64_CODE, ency, encx);	 // encode msg+crc using given QRA code

  // copy codeword to output puncturing the crc symbol
  memcpy(y,ency,QRA64_K*sizeof(int));		// copy information symbols
  memcpy(y+QRA64_K,ency+QRA64_KC,QRA64_C*sizeof(int)); // copy parity symbols

  if (pcodec->apflags!=QRA_AUTOAP)
    return;

  // Here we handle the QRA_AUTOAP mode --------------------------------------------

  // When a [hiscall mycall ?] msg is detected we instruct the decoder
  // to look for [mycall hiscall ?] msgs
  // otherwise when a [cq mycall ?] msg is sent we reset the APTYPE_BOTHCALLS

  // look if the msg sent is a std type message (bit15 of grid field = 0)
  if ((x[9]&0x80)==1)
    return;	// no, it's a text message, nothing to do

  // It's a [hiscall mycall grid] message

  // We assume that mycall is our call (but we don't check it)
  // hiscall the station we are calling or a general call (CQ/QRZ/etc..)
  decodemsg_jt65(&hiscall,&mycall,&grid,x);

  if ((hiscall>=CALL_CQ && hiscall<=CALL_CQ999) || hiscall==CALL_CQDX ||
      hiscall==CALL_DE) {
	// tell the decoder to look for msgs directed to us
	qra64_apset(pcodec,mycall,0,0,APTYPE_MYCALL);
    // We are making a general call and don't know who might reply
    // Reset APTYPE_BOTHCALLS so decoder won't look for [mycall hiscall ?] msgs
    qra64_apdisable(pcodec,APTYPE_BOTHCALLS);
  } else {
    // We are replying to someone named hiscall
    // Set APTYPE_BOTHCALLS so decoder will try for [mycall hiscall ?] msgs
    qra64_apset(pcodec,mycall, hiscall, GRID_BLANK, APTYPE_BOTHCALLS);
  }

}

#define EBNO_MIN -10.0f		// minimum Eb/No value returned by the decoder (in dB)
// AWGN metric decoder
int qra64_decode(qra64codec *pcodec, float *ebno, int *x, const float *rxen)
{
  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise variance
  float msge;						// estimated message energy
  float ebnoval;					// estimated Eb/No
  int rc;

  if (QRA64_NMSG!=QRA64_CODE.NMSG)      // sanity check
    return -16;				// QRA64_NMSG define is wrong

  // compute symbols intrinsic probabilities from received energy observations
  noisestd = qra_mfskbesselmetric(ix, rxen, QRA64_m, QRA64_N,pcodec->decEsNoMetric);

  // de-puncture observations adding a uniform distribution for the crc symbol

  // move check symbols distributions one symbol towards the end
  dsttmp = PD_ROWADDR(ix,QRA64_M, QRA64_NC-1);	//Point to last symbol prob dist
  srctmp = dsttmp-QRA64_M;              // source is the previous pd
  for (k=0;k<QRA64_C;k++) {
    pd_init(dsttmp,srctmp,QRA64_M);
    dsttmp -=QRA64_M;
    srctmp -=QRA64_M;
  }
  // Initialize crc prob to a uniform distribution
  pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

  // Try to decode using all AP cases if required
  rc = qra64_decode_attempts(pcodec, xdec, ix);

  if (rc<0)
	  return rc;	// no success

  // successfull decode --------------------------------

  // copy decoded message (without crc) to output buffer
  memcpy(x,xdec,QRA64_K*sizeof(int));

  if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
	  return rc;

  // reencode message and estimate Eb/No
  qra_encode(&QRA64_CODE, ydec, xdec);
  // puncture crc
  memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));
  // compute total power of decoded message
  msge = 0;
  for (k=0;k<QRA64_N;k++) {
	  msge +=rxen[ydec[k]];	// add energy of current symbol
	  rxen+=QRA64_M;			// ptr to next symbol
	  }

  // NOTE:
  // To make a more accurate Eb/No estimation we should compute the noise variance
  // on all the rxen values but the transmitted symbols.
  // Noisestd is compute by qra_mfskbesselmetric assuming that
  // the signal power is much less than the total noise power in the QRA64_M tones
  // but this is true only if the Eb/No is low.
  // Here, in order to improve accuracy, we linearize the estimated Eb/No value empirically
  // (it gets compressed when it is very high as in this case the noise variance
  // is overestimated)

  // this would be the exact value if the noisestd were not overestimated at high Eb/No
  ebnoval = (0.5f/(QRA64_K*QRA64_m))*msge/(noisestd*noisestd)-1.0f;

  // Empirical linearization (to remove the noise variance overestimation)
  // the resulting SNR is accurate up to +20 dB (51 dB Eb/No)
  if (ebnoval>57.004f)
	  ebnoval=57.004f;
  ebnoval = ebnoval*57.03f/(57.03f-ebnoval);

  // compute value in dB
  if (ebnoval<=0) {
    ebnoval = EBNO_MIN; // assume a minimum, positive value
  }
  else {
    ebnoval = 10.0f*(float)log10(ebnoval);
    if (ebnoval<EBNO_MIN) {
      ebnoval = EBNO_MIN;
    }
  }

  *ebno = ebnoval;

  return rc;
}

//
// Fast-fading / Rayleigh channel metric decoder ----------------------------------------------
//
// Tables of fading energies coefficients for QRA64 (Ts=6912/12000)
// As the fading is assumed to be symmetric around the nominal frequency
// only the leftmost and the central coefficient are stored in the tables.
// (files have been generated with the Matlab code efgengaussenergy.m and efgenlorentzenergy.m)
#include "fadengauss.c"
#include "fadenlorentz.c"

int qra64_decode_fastfading(
				qra64codec *pcodec,		// ptr to the codec structure
				float *ebno,			// ptr to where the estimated Eb/No value will be saved
				int *x,					// ptr to decoded message
				const float *rxen,		// ptr to received symbol energies array
				const int submode,		// submode idx (0=QRA64A ... 4=QRA64E)
				const float B90,	    // spread bandwidth (90% fractional energy)
				const int fadingModel)  // 0=Gaussian 1=Lorentzian fade model

// Decode a QRA64 msg using a fast-fading metric
//
// rxen: The array of the received bin energies
//       Bins must be spaced by integer multiples of the symbol rate (1/Ts Hz)
//       The array must be an array of total length U = L x N where:
//			L: is the number of frequency bins per message symbol (see after)
//          N: is the number of symbols in a QRA64 msg (63)

//       The number of bins/symbol L depends on the selected submode accordingly to
//		 the following rule:
//			L = (64+64*2^submode+64) = 64*(2+2^submode)
//		 Tone 0 is always supposed to be at offset 64 in the array.
//		 The m-th tone nominal frequency is located at offset 64 + m*2^submode (m=0..63)
//
//		 Submode A: (2^submode = 1)
//          L = 64*3 = 196 bins/symbol
//          Total length of the energies array: U = 192*63 = 12096 floats
//
//		 Submode B: (2^submode = 2)
//          L = 64*4 = 256 bins/symbol
//          Total length of the energies array: U = 256*63 = 16128 floats
//
//		 Submode C: (2^submode = 4)
//          L = 64*6 = 384 bins/symbol
//          Total length of the energies array: U = 384*63 = 24192 floats
//
//		 Submode D: (2^submode = 8)
//          L = 64*10 = 640 bins/symbol
//          Total length of the energies array: U = 640*63 = 40320 floats
//
//		 Submode E: (2^submode = 16)
//          L = 64*18 = 1152 bins/symbol
//          Total length of the energies array: U = 1152*63 = 72576 floats
//
//		Note: The rxen array is modified and reused for internal calculations.
//
//
//	B90: spread fading bandwidth in Hz (90% fractional average energy)
//
//			B90 should be in the range 1 Hz ... 238 Hz
//			The value passed to the call is rounded to the closest value among the
//			64 available values:
//				B = 1.09^k Hz, with k=0,1,...,63
//
//			I.e. B90=27 Hz will be approximated in this way:
//				k = rnd(log(27)/log(1.09)) = 38
//              B90 = 1.09^k = 1.09^38 = 26.4 Hz
//
//          For any input value the maximum rounding error is not larger than +/- 5%
//

{

  int k;
  float *srctmp, *dsttmp;
  float ix[QRA64_NC*QRA64_M];		// (depunctured) intrisic information
  int   xdec[QRA64_KC];				// decoded message (with crc)
  int   ydec[QRA64_NC];				// re-encoded message (for snr calculations)
  float noisestd;					// estimated noise std
  float esno,ebnoval;				// estimated Eb/No
  float tempf;
  float EsNoMetric;

  int rc;
  int hidx, hlen;
  const float *hptr;

	if (QRA64_NMSG!=QRA64_CODE.NMSG)
		return -16;					// QRA64_NMSG define is wrong

	if (submode<0 || submode>4)
		return -17;				// invalid submode

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate amplitude weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian energy weighting taps
		hlen = glen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian energy weighting taps
		hlen = glen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// compute (euristically) the optimal decoder metric accordingly the given spread amount
	// We assume that the decoder threshold is:
	//		Es/No(dB) = Es/No(AWGN)(dB) + 8*log(B90)/log(240)(dB)
	// that's to say, at the maximum Doppler spread bandwidth (240 Hz) there's a ~8 dB Es/No degradation
	// over the AWGN case
	tempf = 8.0f*(float)log((float)B90)/(float)log(240.0f);
	EsNoMetric = pcodec->decEsNoMetric*(float)pow(10.0f,tempf/10.0f);

	// Step 1 -----------------------------------------------------------------------------------
	// Evaluate the noise stdev from the received energies at nominal tone frequencies
	noisestd = qra64_fastfading_estim_noise_std(rxen, EsNoMetric, submode);

	// Step 2 -----------------------------------------------------------------------------------
	// Compute message symbols probability distributions
	qra64_fastfading_intrinsics(ix, rxen, hptr, hlen, noisestd, EsNoMetric, submode);

	// Step 3 ---------------------------------------------------------------------------
	// De-puncture observations adding a uniform distribution for the crc symbol
	// Move check symbols distributions one symbol towards the end
	dsttmp = PD_ROWADDR(ix,QRA64_M, QRA64_NC-1);	//Point to last symbol prob dist
	srctmp = dsttmp-QRA64_M;              // source is the previous pd
	for (k=0;k<QRA64_C;k++) {
		pd_init(dsttmp,srctmp,QRA64_M);
		dsttmp -=QRA64_M;
		srctmp -=QRA64_M;
		}
	// Initialize crc prob to a uniform distribution
	pd_init(dsttmp,pd_uniform(QRA64_m),QRA64_M);

	// Step 4 ---------------------------------------------------------------------------
	// Attempt to decode
	rc = qra64_decode_attempts(pcodec, xdec, ix);
	if (rc<0)
	  return rc;	// no success

	// copy decoded message (without crc) to output buffer
	memcpy(x,xdec,QRA64_K*sizeof(int));

	// Step 5 ----------------------------------------------------------------------------
	// Estimate the message Eb/No

	if (ebno==0)	// null pointer indicates we are not interested in the Eb/No estimate
		return rc;

	// reencode message to estimate Eb/No
	qra_encode(&QRA64_CODE, ydec, xdec);
	// puncture crc
	memmove(ydec+QRA64_K,ydec+QRA64_KC,QRA64_C*sizeof(int));

	// compute Es/N0 of decoded message
	esno = qra64_fastfading_msg_esno(ydec,rxen,noisestd, EsNoMetric, hlen,submode);

	// as the weigthing function include about 90% of the energy
	// we could compute the unbiased esno with:
	// esno = esno/0.9;

	// Es/N0 --> Eb/N0 conversion
	ebnoval = 1.0f/(1.0f*QRA64_K/QRA64_N*QRA64_m)*esno;

	// compute value in dB
	if (ebnoval<=0) {
	  ebnoval = EBNO_MIN; // assume a minimum, positive value
    }
	else {
	  ebnoval = 10.0f*(float)log10(ebnoval);
	  if (ebnoval<EBNO_MIN) {
		  ebnoval = EBNO_MIN;
      }
    }

	*ebno = ebnoval;

  return rc;

}

/*
int qra64_fastfading_channel(float **rxen, const int *xmsg, const int submode, const float EbN0dB, const float B90, const int fadingModel)
{
	// Simulate transmission over a fading channel and non coherent detection

	// Set rxen to point to an array of bin energies formatted as required
	// by the (fast-fading) decoding routine

	// returns 0 on success or negative values on error conditions

	static float *channel_out = NULL;
	static int    channel_submode = -1;

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int n,j,hidx, hlen;
	const float *hptr;
	float *cursym,*curtone;

	float iq[2];
	float *curi, *curq;
	float sigmasig[65];	// signal standard deviation taps

	float N0, EsN0, Es, sigmanoise;

	if (rxen==NULL)
		return -1;		// rxen must be a non-null ptr

	// allocate output buffer if not yet done or if submode changed
	if (channel_out==NULL || submode!=channel_submode) {

		// unallocate previous buffer
		if (channel_out)
			free(channel_out);

		// allocate new buffer
		// we allocate twice the mem so that we can store/compute complex amplitudes
		channel_out = (float*)malloc(bpm*sizeof(float)*2);
		if (channel_out==NULL)
			return -2;	// error allocating memory

		channel_submode = submode;
		}

	if (B90<1.0f || B90>238.0f)
		return -18;				// B90 out of range

	// compute index to most appropriate energy weighting function coefficients
    hidx = (int)(log((float)B90)/log(1.09f) - 0.499f);

	if (hidx<0 || hidx > 64)
		return -19;				// index of weighting function out of range

	if (fadingModel==0) {	 // gaussian fading model
		// point to gaussian weighting taps
		hlen = glen_tab_gauss[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_gauss[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else if (fadingModel==1) {
		// point to lorentzian weighting taps
		hlen = glen_tab_lorentz[hidx];	 // hlen = (L+1)/2 (where L=(odd) number of taps of w fun)
		hptr = gptr_tab_lorentz[hidx];     // pointer to the first (L+1)/2 coefficients of w fun
		}
	else
		return -20;			// invalid fading model index

	// Compute the unfaded tone amplitudes from the Eb/No value passed to the call
	N0 = 1.0f;	// assume unitary noise PSD
	sigmanoise = (float)sqrt(N0/2);
	EsN0 = (float)pow(10.0f,EbN0dB/10.0f)*QRA64_m*QRA64_K/QRA64_N; // Es/No = m*R*Eb/No
	Es   = EsN0*N0;

	// compute signal bin sigmas
	for (n=0;n<hlen;n++)
		sigmasig[n] = (float)sqrt(hptr[n]*Es/2.0f);

	// Generate gaussian noise iq components
	normrnd_s(channel_out, bpm*2, 0 , sigmanoise);

	// Add symbols, bin by bin energies
	for (n=0;n<QRA64_N;n++) {

		cursym  = channel_out+n*bps + QRA64_M; // point to n-th symbol
		curtone = cursym+xmsg[n]*bpt;	 // point to encoded tone
		curi    = curtone-hlen+1;		 // point to real part of first bin
		curq    = curtone-hlen+1+bpm;	 // point to imag part of first bin

		// generate Rayleigh faded bins with given average energy and add to noise
		for (j=0;j<hlen;j++) {
			normrnd_s(iq, 2, 0 , sigmasig[j]);
			*curi++ += iq[0];
			*curq++ += iq[1];
			}
		for (j=hlen-2;j>=0;j--) {
			normrnd_s(iq, 2, 0 , sigmasig[j]);
			*curi++ += iq[0];
			*curq++ += iq[1];
			}

		}

	// compute total bin energies (S+N) and store in first half of buffer
	curi = channel_out;
	curq = channel_out+bpm;
	for (n=0;n<bpm;n++)
		channel_out[n] = curi[n]*curi[n] + curq[n]*curq[n];

	// set rxen to point to the channel output energies
	*rxen = channel_out;

	return 0;
}
*/

// Static functions definitions ----------------------------------------------

// fast-fading static functions --------------------------------------------------------------

static float qra64_fastfading_estim_noise_std(const float *rxen, const float esnometric, const int submode)
{
	// estimate the noise standard deviation from nominal frequency symbol bins
	// transform energies to amplitudes

	// rxen = message symbols energies (overwritten with symbols amplitudes)
	// esnometric = Es/No at nominal frequency bin for which we compute the decoder metric
	// submode = submode used (0=A...4=E)

	int bpt = (1<<submode);			// bins per tone
	int bps = QRA64_M*(2+bpt);		// total number of bins per symbols
	int bpm = bps*QRA64_N;			// total number of bins in a message
	int k;
	float sigmaest;

	// estimate noise std
	sigmaest = 0;
	for (k=0;k<bpm;k++)
		sigmaest += rxen[k];

	sigmaest = sigmaest/bpm;
	sigmaest = (float)sqrt(sigmaest/(1.0f+esnometric/bps)/2.0f);

	// Note: sigma is overestimated by the (unknown) factor sqrt((1+esno(true)/bps)/(1+esnometric/bps))

	return sigmaest;
}

static void qra64_fastfading_intrinsics(
				float *pix,
				const float *rxen,
				const float *hptr,
				const int    hlen,
				const float sigmaest,
				const float EsNoMetric,
				const int submode)
{

	// For each symbol in a message:
	// a) Compute tones loglikelihoods as a sum of products between of the expected
	// energy fading coefficient and received energies.
	// b) Compute intrinsic symbols probability distributions from symbols loglikelihoods

	int n,k,j, bps, bpt;
	const float *cursym, *curbin;
	float *curix;
	float u, maxloglh, loglh, sumix,hh;
	float w[65];
	int hhsz  = hlen-1;
	int hlast = 2*hhsz;
	float npwrest = 2.0f*sigmaest*sigmaest;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	u = EsNoMetric;
	// compute weights from energy tables
	for (j=0;j<hlen;j++) {
		hh = hptr[j]*u;
		w[j] = hh/(1+hh)/npwrest;
		}

	for (n=0;n<QRA64_N;n++) {			// for each symbol in the message
		cursym = rxen+n*bps + QRA64_M;	// point to current symbol nominal bin
		maxloglh = 0;
		curix  = pix+n*QRA64_M;
		for (k=0;k<QRA64_M;k++) {   // for each tone in the current symbol
			curbin = cursym + k*bpt -hlen+1;
			// compute tone loglikelihood (symmetric fir with given weights)
			loglh = 0.f;
			for (j=0;j<hhsz;j++)
				loglh += w[j]*(curbin[j] + curbin[hlast-j]);
			loglh += w[hhsz]*curbin[hhsz];

			if (loglh>maxloglh)		// keep track of the max loglikelihood
				maxloglh = loglh;
			curix[k]=loglh;
			}

		// scale to likelihoods
		sumix = 0.f;
		for (k=0;k<QRA64_M;k++) {
			u = (float)exp(curix[k]-maxloglh);
			curix[k]=u;
			sumix +=u;
			}
		// scale to probabilities
		sumix = 1.0f/sumix;
		for (k=0;k<QRA64_M;k++)
			curix[k] = curix[k]*sumix;
		}
}

static float qra64_fastfading_msg_esno(
			const int *ydec,
			const float *rxen,
			const float sigma,
			const float EsNoMetric,
			const int hlen,
			const int submode)
{
	// Estimate msg Es/N0

	int n,j, bps, bpt;
	const float *cursym, *curtone, *curbin;
	float u, msgsn,esno;
	int tothlen = 2*hlen-1;

	bpt = 1<<submode;				// bins per tone
	bps = QRA64_M*(2+bpt);			// bins per symbol

	msgsn = 0;
	for (n=0;n<QRA64_N;n++) {
		cursym  = rxen+n*bps + QRA64_M; // point to n-th symbol amplitudes
		curtone = cursym+ydec[n]*bpt;	 // point to decoded tone amplitudes
		curbin  = curtone-hlen+1;		 // point to first bin amplitude

		// sum bin energies
		for (j=0;j<tothlen;j++)
			msgsn += curbin[j];

		}

	msgsn =  msgsn/(QRA64_N*tothlen);	// avg msg energy per bin (noise included)

	// as sigma is overestimated (sigmatrue = sigma*sqrt((1+EsNoMetric/bps)/(1+EsNo/bps))
	// we have: msgsn = (1+x/hlen)/(1+x/bps)*2*sigma^2*(1+EsnoMetric/bps), where x = Es/N0(true)
	//
	// we can then write:
	// u = msgsn/2.0f/(sigma*sigma)/(1.0f+EsNoMetric/bps);
	// (1+x/hlen)/(1+x/bps) = u

	u = msgsn/(2.0f*sigma*sigma)/(1.0f+EsNoMetric/bps);

	// check u>1
	if (u<1)
		return 0.f;

	// check u<bps/tot hlen
	if (u>(bps/tothlen))
		return 10000.f;

	// solve for Es/No
	esno = (u-1.0f)/(1.0f/tothlen-u/bps);

	return esno;

}

#ifdef LIMIT_AP_MASKS

static int call1_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_CALL1   =   0xFFFFFFC
	int u = papmsg[4]^pdec[4];
	return (u&0x3C)==0;
}
static int call2_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_CALL2   =   0xFFFFFFC
	int u = papmsg[9]^pdec[9];
	return  (u&0x30)==0;
}
static int grid_match(const int *papmsg, const int *pdec)
{
	// assumes MASK_GRIDFULL =	0x3FFC
	int u = papmsg[11]^pdec[11];
	int rc = (u&0x03)==0;

	u = papmsg[9]^pdec[9];

	return (u&0x0C)==0 && rc;
}

#else
#define call1_match(a,b) (1)
#define call2_match(a,b) (1)
#define grid_match(a,b) (1)
#endif

// Attempt to decode given intrisic information
static int qra64_decode_attempts(qra64codec *pcodec, int *xdec, const float *ix)
{
  int rc;

  // Attempt to decode without a-priori info --------------------------------
  rc = qra64_do_decode(xdec, ix, NULL, NULL);
  if (rc>=0)
	  return 0; // successfull decode with AP0
  else
	  if (pcodec->apflags==QRA_NOAP)
		  // nothing more to do
		  return rc; // rc<0 = unsuccessful decode

  // Here we handle decoding with AP knowledge

  // Attempt to decode CQ calls
  rc = qra64_do_decode(xdec,ix,pcodec->apmask_cqqrz, pcodec->apmsg_cqqrz);
  if (rc>=0)
	  return 1;    // decoded [cq/qrz ? ?]

  rc = qra64_do_decode(xdec, ix, pcodec->apmask_cqqrz_ooo,
		       pcodec->apmsg_cqqrz);
  if (rc>=0)
	  // check that ooo really matches
	  if (grid_match(pcodec->apmsg_cqqrz,xdec))
		  return 2;    // decoded [cq/qrz ? ooo]

  // attempt to decode calls directed to us
  if (pcodec->apmsg_set[APTYPE_MYCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1,
		       pcodec->apmsg_call1);
	if (rc>=0)
		// check that mycall really matches
		if (call1_match(pcodec->apmsg_call1,xdec))
			return 3;    // decoded [mycall ? ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_ooo,
		       pcodec->apmsg_call1);
	if (rc>=0)
		// check that mycall and ooo really match
		if (call1_match(pcodec->apmsg_call1,xdec) &&
			grid_match(pcodec->apmsg_call1,xdec))
			return 4;    // decoded [mycall ? ooo]
	}

  // attempt to decode [mycall hiscall ?] msgs
  if (pcodec->apmsg_set[APTYPE_BOTHCALLS]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2,
		       pcodec->apmsg_call1_call2);
	if (rc>=0)
		// check that mycall and hiscall really match
		if (call1_match(pcodec->apmsg_call1_call2,xdec) &&
			call2_match(pcodec->apmsg_call1_call2,xdec))
			return 5;    // decoded [mycall srccall ?]
	}

  // attempt to decode [? hiscall ?/b] msgs
  if (pcodec->apmsg_set[APTYPE_HISCALL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2,
		       pcodec->apmsg_call2);
	if (rc>=0)
		// check that hiscall really match
		if (call2_match(pcodec->apmsg_call2,xdec))
			return 6;    // decoded [? hiscall ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call2_ooo,
		       pcodec->apmsg_call2);
	if (rc>=0)
		// check that hiscall and ooo match
		if (call2_match(pcodec->apmsg_call2,xdec) &&
			grid_match(pcodec->apmsg_call2,xdec))
			return 7;    // decoded [? hiscall ooo]
	}

  // attempt to decode [cq/qrz hiscall ?/b/grid] msgs
  if (pcodec->apmsg_set[APTYPE_CQHISCALL]) {

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2,
				pcodec->apmsg_cq_call2);
	if (rc>=0)
		// check that hiscall matches
		if (call2_match(pcodec->apmsg_call2,xdec))
			return 9;	// decoded [cq/qrz hiscall ?]

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2_grid,xdec, QRA64_K*sizeof(int))==0)
			return 11;		// decoded [cq/qrz hiscall grid]
		}

	rc = qra64_do_decode(xdec, ix, pcodec->apmask_cq_call2_ooo,
		       pcodec->apmsg_cq_call2);
	if (rc>=0) {
		// Full AP mask need special handling
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_cq_call2,xdec, QRA64_K*sizeof(int))==0)
			return 10;    // decoded [cq/qrz hiscall ]
		}
	}

  // attempt to decode [mycall hiscall grid]
  if (pcodec->apmsg_set[APTYPE_FULL]) {
	rc = qra64_do_decode(xdec, ix, pcodec->apmask_call1_call2_grid,
		       pcodec->apmsg_call1_call2_grid);
	if (rc>=0) {
		// Full AP mask need special handling
		// All the three msg fields were given.
		// To minimize false decodes we check the decoded message
		// with what passed in the ap_set call
		if (memcmp(pcodec->apmsg_call1_call2_grid,xdec, QRA64_K*sizeof(int))==0)
			return 8;	   // decoded [mycall hiscall grid]
		}
	}

  // all decoding attempts failed
  return -1;
}

// Decode with given a-priori information
static int qra64_do_decode(int *xdec, const float *pix, const int *ap_mask,
			   const int *ap_x)
{
  int rc;
  const float *ixsrc;
  float ix_masked[QRA64_NC*QRA64_M];  // Masked intrinsic information
  float ex[QRA64_NC*QRA64_M];	      // Extrinsic information from the decoder

  float v2cmsg[QRA64_NMSG*QRA64_M];   // buffers for the decoder messages
  float c2vmsg[QRA64_NMSG*QRA64_M];

  if (ap_mask==NULL) {   // no a-priori information
    ixsrc = pix;	 // intrinsic source is what passed as argument
  } else {
    // a-priori information provided
    // mask channel observations with a-priori
    ix_mask(ix_masked,pix,ap_mask,ap_x);
    ixsrc = ix_masked;	// intrinsic source is the masked version
  }

  // run the decoding algorithm
  rc = qra_extrinsic(&QRA64_CODE,ex,ixsrc,QRA64_NITER,v2cmsg,c2vmsg);
  if (rc<0)
    return -1;	// no convergence in given iterations

  // decode
  qra_mapdecode(&QRA64_CODE,xdec,ex,ixsrc);

  // verify crc
  if (calc_crc6(xdec,QRA64_K)!=xdec[QRA64_K]) // crc doesn't match (detected error)
    return -2;	// decoding was succesfull but crc doesn't match

  return 0;
}

// crc functions --------------------------------------------------------------
// crc-6 generator polynomial
// g(x) = x^6 + a5*x^5 + ... + a1*x + a0

// g(x) = x^6 + x + 1
#define CRC6_GEN_POL 0x30  // MSB=a0 LSB=a5

// g(x) = x^6 + x^2 + x + 1 (See:  https://users.ece.cmu.edu/~koopman/crc/)
// #define CRC6_GEN_POL 0x38  // MSB=a0 LSB=a5. Simulation results are similar

static int calc_crc6(const int *x, int sz)
{
  // todo: compute it faster using a look up table
  int k,j,t,sr = 0;
  for (k=0;k<sz;k++) {
    t = x[k];
    for (j=0;j<6;j++) {
      if ((t^sr)&0x01)
	sr = (sr>>1) ^ CRC6_GEN_POL;
      else
	sr = (sr>>1);
      t>>=1;
    }
  }
  return sr;
}

static void ix_mask(float *dst, const float *src, const int *mask,
		    const int *x)
{
  // mask intrinsic information (channel observations) with a priori knowledge

  int k,kk, smask;
  float *row;

  memcpy(dst,src,(QRA64_NC*QRA64_M)*sizeof(float));

  for (k=0;k<QRA64_K;k++) {	// we can mask only information symbols distrib
    smask = mask[k];
    row = PD_ROWADDR(dst,QRA64_M,k);
    if (smask) {
      for (kk=0;kk<QRA64_M;kk++)
	if (((kk^x[k])&smask)!=0)
	  *(row+kk) = 0.f;

      pd_norm(row,QRA64_m);
    }
  }
}

// encode/decode msgs as done in JT65
void encodemsg_jt65(int *y, const int call1, const int call2, const int grid)
{
  y[0]= (call1>>22)&0x3F;
  y[1]= (call1>>16)&0x3F;
  y[2]= (call1>>10)&0x3F;
  y[3]= (call1>>4)&0x3F;
  y[4]= (call1<<2)&0x3F;

  y[4] |= (call2>>26)&0x3F;
  y[5]= (call2>>20)&0x3F;
  y[6]= (call2>>14)&0x3F;
  y[7]= (call2>>8)&0x3F;
  y[8]= (call2>>2)&0x3F;
  y[9]= (call2<<4)&0x3F;

  y[9] |= (grid>>12)&0x3F;
  y[10]= (grid>>6)&0x3F;
  y[11]= (grid)&0x3F;

}
void decodemsg_jt65(int *call1, int *call2, int *grid, const int *x)
{
  int nc1, nc2, ng;

  nc1 = x[4]>>2;
  nc1 |= x[3]<<4;
  nc1 |= x[2]<<10;
  nc1 |= x[1]<<16;
  nc1 |= x[0]<<22;

  nc2 = x[9]>>4;
  nc2 |= x[8]<<2;
  nc2 |= x[7]<<8;
  nc2 |= x[6]<<14;
  nc2 |= x[5]<<20;
  nc2 |= (x[4]&0x03)<<26;

  ng   = x[11];
  ng  |= x[10]<<6;
  ng  |= (x[9]&0x0F)<<12;

  *call1 = nc1;
  *call2 = nc2;
  *grid  = ng;
}


code:
// Copyright Peter Dimov 2001
// Copyright Aleksey Gurtovoy 2001-2004
//
// 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)
//

// Preprocessed version of "boost/mpl/bind.hpp" header
// -- DO NOT modify by hand!

namespace boost { namespace mpl {

namespace aux {
template< bool >
struct resolve_arg_impl
{
    template<
          typename T, typename U1, typename U2, typename U3
        , typename U4, typename U5
        >
    struct result_
    {
        typedef T type;
    };
};

template<>
struct resolve_arg_impl<true>
{
    template<
          typename T, typename U1, typename U2, typename U3
        , typename U4, typename U5
        >
    struct result_
    {
        typedef typename apply_wrap5<
              T
            , U1, U2, U3, U4, U5
            >::type type;
    };
};

template< typename T > struct is_bind_template;

template<
      typename T, typename U1, typename U2, typename U3, typename U4
    , typename U5
    >
struct resolve_bind_arg
    : resolve_arg_impl< is_bind_template<T>::value >
            ::template result_< T,U1,U2,U3,U4,U5 >
{
};

template< typename T >
struct replace_unnamed_arg_impl
{
    template< typename Arg > struct result_
    {
        typedef Arg next;
        typedef T type;
    };
};

template<>
struct replace_unnamed_arg_impl< arg< -1 > >
{
    template< typename Arg > struct result_
    {
        typedef typename next<Arg>::type next;
        typedef Arg type;
    };
};

template< typename T, typename Arg >
struct replace_unnamed_arg
    : replace_unnamed_arg_impl<T>::template result_<Arg>
{
};

template< int arity_ > struct bind_chooser;

aux::no_tag is_bind_helper(...);
template< typename T > aux::no_tag is_bind_helper(protect<T>*);

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    , typename T5
    >
aux::yes_tag is_bind_helper(bind< F,T1,T2,T3,T4,T5 >*);

template< int N >
aux::yes_tag is_bind_helper(arg<N>*);

template< bool is_ref_  = true >
struct is_bind_template_impl
{
    template< typename T > struct result_
    {
        BOOST_STATIC_CONSTANT(bool, value  = false);
    };
};

template<>
struct is_bind_template_impl<false>
{
    template< typename T > struct result_
    {
        BOOST_STATIC_CONSTANT(bool, value =
              sizeof(aux::is_bind_helper(static_cast<T*>(0)))
                == sizeof(aux::yes_tag)
            );
    };
};

template< typename T > struct is_bind_template
    : is_bind_template_impl< ::boost::detail::is_reference_impl<T>::value >
        ::template result_<T>
{
};

} // namespace aux

template<
      typename F
    >
struct bind0
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
     public:
        typedef typename apply_wrap0<
              f_
            >::type type;

    };
};

namespace aux {

template<
      typename F
    >
aux::yes_tag
is_bind_helper(bind0<F>*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(1, bind0)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(1, bind0)

namespace aux {

template<>
struct bind_chooser<0>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind0<F> type;
    };
};

} // namespace aux

template<
      typename F, typename T1
    >
struct bind1
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
        typedef aux::replace_unnamed_arg< T1,n1 > r1;
        typedef typename r1::type a1;
        typedef typename r1::next n2;
        typedef aux::resolve_bind_arg< a1,U1,U2,U3,U4,U5 > t1;
        ///
     public:
        typedef typename apply_wrap1<
              f_
            , typename t1::type
            >::type type;

    };
};

namespace aux {

template<
      typename F, typename T1
    >
aux::yes_tag
is_bind_helper(bind1< F,T1 >*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(2, bind1)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(2, bind1)

namespace aux {

template<>
struct bind_chooser<1>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind1< F,T1 > type;
    };
};

} // namespace aux

template<
      typename F, typename T1, typename T2
    >
struct bind2
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
        typedef aux::replace_unnamed_arg< T1,n1 > r1;
        typedef typename r1::type a1;
        typedef typename r1::next n2;
        typedef aux::resolve_bind_arg< a1,U1,U2,U3,U4,U5 > t1;
        ///
        typedef aux::replace_unnamed_arg< T2,n2 > r2;
        typedef typename r2::type a2;
        typedef typename r2::next n3;
        typedef aux::resolve_bind_arg< a2,U1,U2,U3,U4,U5 > t2;
        ///
     public:
        typedef typename apply_wrap2<
              f_
            , typename t1::type, typename t2::type
            >::type type;

    };
};

namespace aux {

template<
      typename F, typename T1, typename T2
    >
aux::yes_tag
is_bind_helper(bind2< F,T1,T2 >*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(3, bind2)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(3, bind2)

namespace aux {

template<>
struct bind_chooser<2>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind2< F,T1,T2 > type;
    };
};

} // namespace aux

template<
      typename F, typename T1, typename T2, typename T3
    >
struct bind3
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
        typedef aux::replace_unnamed_arg< T1,n1 > r1;
        typedef typename r1::type a1;
        typedef typename r1::next n2;
        typedef aux::resolve_bind_arg< a1,U1,U2,U3,U4,U5 > t1;
        ///
        typedef aux::replace_unnamed_arg< T2,n2 > r2;
        typedef typename r2::type a2;
        typedef typename r2::next n3;
        typedef aux::resolve_bind_arg< a2,U1,U2,U3,U4,U5 > t2;
        ///
        typedef aux::replace_unnamed_arg< T3,n3 > r3;
        typedef typename r3::type a3;
        typedef typename r3::next n4;
        typedef aux::resolve_bind_arg< a3,U1,U2,U3,U4,U5 > t3;
        ///
     public:
        typedef typename apply_wrap3<
              f_
            , typename t1::type, typename t2::type, typename t3::type
            >::type type;

    };
};

namespace aux {

template<
      typename F, typename T1, typename T2, typename T3
    >
aux::yes_tag
is_bind_helper(bind3< F,T1,T2,T3 >*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(4, bind3)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(4, bind3)

namespace aux {

template<>
struct bind_chooser<3>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind3< F,T1,T2,T3 > type;
    };
};

} // namespace aux

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    >
struct bind4
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
        typedef aux::replace_unnamed_arg< T1,n1 > r1;
        typedef typename r1::type a1;
        typedef typename r1::next n2;
        typedef aux::resolve_bind_arg< a1,U1,U2,U3,U4,U5 > t1;
        ///
        typedef aux::replace_unnamed_arg< T2,n2 > r2;
        typedef typename r2::type a2;
        typedef typename r2::next n3;
        typedef aux::resolve_bind_arg< a2,U1,U2,U3,U4,U5 > t2;
        ///
        typedef aux::replace_unnamed_arg< T3,n3 > r3;
        typedef typename r3::type a3;
        typedef typename r3::next n4;
        typedef aux::resolve_bind_arg< a3,U1,U2,U3,U4,U5 > t3;
        ///
        typedef aux::replace_unnamed_arg< T4,n4 > r4;
        typedef typename r4::type a4;
        typedef typename r4::next n5;
        typedef aux::resolve_bind_arg< a4,U1,U2,U3,U4,U5 > t4;
        ///
     public:
        typedef typename apply_wrap4<
              f_
            , typename t1::type, typename t2::type, typename t3::type
            , typename t4::type
            >::type type;

    };
};

namespace aux {

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    >
aux::yes_tag
is_bind_helper(bind4< F,T1,T2,T3,T4 >*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(5, bind4)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(5, bind4)

namespace aux {

template<>
struct bind_chooser<4>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind4< F,T1,T2,T3,T4 > type;
    };
};

} // namespace aux

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    , typename T5
    >
struct bind5
{
    template<
          typename U1 = na, typename U2 = na, typename U3 = na
        , typename U4 = na, typename U5 = na
        >
    struct apply
    {
     private:
        typedef aux::replace_unnamed_arg< F, mpl::arg<1> > r0;
        typedef typename r0::type a0;
        typedef typename r0::next n1;
        typedef typename aux::resolve_bind_arg< a0,U1,U2,U3,U4,U5 >::type f_;
        ///
        typedef aux::replace_unnamed_arg< T1,n1 > r1;
        typedef typename r1::type a1;
        typedef typename r1::next n2;
        typedef aux::resolve_bind_arg< a1,U1,U2,U3,U4,U5 > t1;
        ///
        typedef aux::replace_unnamed_arg< T2,n2 > r2;
        typedef typename r2::type a2;
        typedef typename r2::next n3;
        typedef aux::resolve_bind_arg< a2,U1,U2,U3,U4,U5 > t2;
        ///
        typedef aux::replace_unnamed_arg< T3,n3 > r3;
        typedef typename r3::type a3;
        typedef typename r3::next n4;
        typedef aux::resolve_bind_arg< a3,U1,U2,U3,U4,U5 > t3;
        ///
        typedef aux::replace_unnamed_arg< T4,n4 > r4;
        typedef typename r4::type a4;
        typedef typename r4::next n5;
        typedef aux::resolve_bind_arg< a4,U1,U2,U3,U4,U5 > t4;
        ///
        typedef aux::replace_unnamed_arg< T5,n5 > r5;
        typedef typename r5::type a5;
        typedef typename r5::next n6;
        typedef aux::resolve_bind_arg< a5,U1,U2,U3,U4,U5 > t5;
        ///
     public:
        typedef typename apply_wrap5<
              f_
            , typename t1::type, typename t2::type, typename t3::type
            , typename t4::type, typename t5::type
            >::type type;

    };
};

namespace aux {

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    , typename T5
    >
aux::yes_tag
is_bind_helper(bind5< F,T1,T2,T3,T4,T5 >*);

} // namespace aux

BOOST_MPL_AUX_ARITY_SPEC(6, bind5)
BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(6, bind5)

namespace aux {

template<>
struct bind_chooser<5>
{
    template<
          typename F, typename T1, typename T2, typename T3, typename T4
        , typename T5
        >
    struct result_
    {
        typedef bind5< F,T1,T2,T3,T4,T5 > type;
    };
};

} // namespace aux

namespace aux {

template< typename T >
struct is_bind_arg
{
    BOOST_STATIC_CONSTANT(bool, value  = true);
};

template<>
struct is_bind_arg<na>
{
    BOOST_STATIC_CONSTANT(bool, value  = false);
};

template<
      typename T1, typename T2, typename T3, typename T4, typename T5
    >
struct bind_count_args
{
    BOOST_STATIC_CONSTANT(int, value =
          is_bind_arg<T1>::value + is_bind_arg<T2>::value
        + is_bind_arg<T3>::value + is_bind_arg<T4>::value
        + is_bind_arg<T5>::value
        );

};

}

template<
      typename F, typename T1, typename T2, typename T3, typename T4
    , typename T5
    >
struct bind
    : aux::bind_chooser<
          aux::bind_count_args< T1,T2,T3,T4,T5 >::value
        >::template result_< F,T1,T2,T3,T4,T5 >::type
{
};

BOOST_MPL_AUX_ARITY_SPEC(
      6
    , bind
    )

BOOST_MPL_AUX_TEMPLATE_ARITY_SPEC(
      6
    , bind
    )
}}


code:
// -*- Mode: C++ -*-
/*
 * Class to handle the formatted string as returned from the fortran decoder
 *
 * VK3ACF August 2013
 */

#ifndef DECODEDTEXT_H
#define DECODEDTEXT_H

#include <QString>

/*
012345678901234567890123456789012345678901
^    ^    ^   ^    ^  ^
2343 -11  0.8 1259 #  CQ VP2X/GM4WJS GL33
2343 -11  0.8 1259 #  CQ 999 VP2V/GM4WJS
2343 -11  0.8 1259 #  YV6BFE F6GUU R-08
2343 -19  0.3  718 #  VE6WQ SQ2NIJ -14
2343  -7  0.3  815 #  KK4DSD W7VP -16
2343 -13  0.1 3627 @  CT1FBK IK5YZT R+02

0605  Tx      1259 #  CQ VK3ACF QF22
*/

class DecodedText
{
public:
  explicit DecodedText (QString const& message, bool, QString const& my_grid);

  QString string() const { return string_; };
  QStringList messageWords () const;
  int indexOf(QString s) const { return string_.indexOf(s); };
  int indexOf(QString s, int i) const { return string_.indexOf(s,i); };
  QString mid(int f, int t) const { return string_.mid(f,t); };
  QString left(int i) const { return string_.left(i); };

  void clear() { string_.clear(); };

  QString CQersCall() const;

  bool isJT65() const;
  bool isJT9() const;
  bool isTX() const;
  bool isStandardMessage () const {return is_standard_;}
  bool isLowConfidence () const;
  int frequencyOffset() const;  // hertz offset from the tuned dial or rx frequency, aka audio frequency
  int snr() const;
  float dt() const;

  // find and extract any report. Returns true if this is a standard message
  bool report(QString const& myBaseCall, QString const& dxBaseCall, /*mod*/QString& report) const;

  // get the first message text word, usually the call
  QString call() const;

  // get the second word, most likely the de call and the third word, most likely grid
  void deCallAndGrid(/*out*/QString& call, QString& grid) const;

  unsigned timeInSeconds() const;

  // returns a string of the SNR field with a leading + or - followed by two digits
  QString report() const;

private:
  // These define the columns in the decoded text where fields are to be found.
  // We rely on these columns being the same in the fortran code (lib/decoder.f90) that formats the decoded text
  enum Columns {column_time    = 0,
      column_snr     = 5,
      column_dt      = 9,
      column_freq    = 14,
      column_mode    = 19,
      column_qsoText = 22 };

  QString string_;
  int padding_;
  bool contest_mode_;
  QString message_;
  bool is_standard_;
};

#endif // DECODEDTEXT_H


code:
// Copyright Aleksey Gurtovoy 2000-2004
//
// 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)
//

// *Preprocessed* version of the main "fold_impl.hpp" header
// -- DO NOT modify by hand!

namespace boost { namespace mpl { namespace aux {

/// forward declaration

template<
      int N
    , typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl;

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< 0,First,Last,State,ForwardOp >
{
    typedef First iter0;
    typedef State state0;
    typedef state0 state;
    typedef iter0 iterator;
};

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< 1,First,Last,State,ForwardOp >
{
    typedef First iter0;
    typedef State state0;
    typedef typename apply2< ForwardOp, state0, typename deref<iter0>::type >::type state1;
    typedef typename mpl::next<iter0>::type iter1;

    typedef state1 state;
    typedef iter1 iterator;
};

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< 2,First,Last,State,ForwardOp >
{
    typedef First iter0;
    typedef State state0;
    typedef typename apply2< ForwardOp, state0, typename deref<iter0>::type >::type state1;
    typedef typename mpl::next<iter0>::type iter1;
    typedef typename apply2< ForwardOp, state1, typename deref<iter1>::type >::type state2;
    typedef typename mpl::next<iter1>::type iter2;

    typedef state2 state;
    typedef iter2 iterator;
};

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< 3,First,Last,State,ForwardOp >
{
    typedef First iter0;
    typedef State state0;
    typedef typename apply2< ForwardOp, state0, typename deref<iter0>::type >::type state1;
    typedef typename mpl::next<iter0>::type iter1;
    typedef typename apply2< ForwardOp, state1, typename deref<iter1>::type >::type state2;
    typedef typename mpl::next<iter1>::type iter2;
    typedef typename apply2< ForwardOp, state2, typename deref<iter2>::type >::type state3;
    typedef typename mpl::next<iter2>::type iter3;

    typedef state3 state;
    typedef iter3 iterator;
};

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< 4,First,Last,State,ForwardOp >
{
    typedef First iter0;
    typedef State state0;
    typedef typename apply2< ForwardOp, state0, typename deref<iter0>::type >::type state1;
    typedef typename mpl::next<iter0>::type iter1;
    typedef typename apply2< ForwardOp, state1, typename deref<iter1>::type >::type state2;
    typedef typename mpl::next<iter1>::type iter2;
    typedef typename apply2< ForwardOp, state2, typename deref<iter2>::type >::type state3;
    typedef typename mpl::next<iter2>::type iter3;
    typedef typename apply2< ForwardOp, state3, typename deref<iter3>::type >::type state4;
    typedef typename mpl::next<iter3>::type iter4;

    typedef state4 state;
    typedef iter4 iterator;
};

template<
      int N
    , typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl
{
    typedef fold_impl<
          4
        , First
        , Last
        , State
        , ForwardOp
        > chunk_;

    typedef fold_impl<
          ( (N - 4) < 0 ? 0 : N - 4 )
        , typename chunk_::iterator
        , Last
        , typename chunk_::state
        , ForwardOp
        > res_;

    typedef typename res_::state state;
    typedef typename res_::iterator iterator;
};

template<
      typename First
    , typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< -1,First,Last,State,ForwardOp >
    : fold_impl<
          -1
        , typename mpl::next<First>::type
        , Last
        , typename apply2<ForwardOp,State, typename deref<First>::type>::type
        , ForwardOp
        >
{
};

template<
      typename Last
    , typename State
    , typename ForwardOp
    >
struct fold_impl< -1,Last,Last,State,ForwardOp >
{
    typedef State state;
    typedef Last iterator;
};

}}}


code:
//  (C) Copyright John Maddock 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)
//
//  This file is machine generated, do not edit by hand

// Unrolled polynomial evaluation using second order Horners rule
#ifndef BOOST_MATH_TOOLS_POLY_EVAL_3_HPP
#define BOOST_MATH_TOOLS_POLY_EVAL_3_HPP

namespace boost{ namespace math{ namespace tools{ namespace detail{

template <class T, class V>
inline V evaluate_polynomial_c_imp(const T*, const V&, const mpl::int_<0>*) BOOST_MATH_NOEXCEPT(V)
{
   return static_cast<V>(0);
}

template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*) BOOST_MATH_NOEXCEPT(V)
{
   return static_cast<V>(a[0]);
}

template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*) BOOST_MATH_NOEXCEPT(V)
{
   return static_cast<V>(a[1] * x + a[0]);
}

template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*) BOOST_MATH_NOEXCEPT(V)
{
   return static_cast<V>((a[2] * x + a[1]) * x + a[0]);
}

template <class T, class V>
inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*) BOOST_MATH_NOEXCEPT(V)
{
   return static_cast<V>(((a[3] * x + a[2]) * x + a[1]) * x + a[0]);
}

}}}} // namespaces

#endif // include guard


code:
/*
 [auto_generated]
 boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp

 [begin_description]
 Implementation of the Runge Kutta Cash Karp 5(4) method. It uses the generic error stepper.
 [end_description]

 Copyright 2011-2013 Mario Mulansky
 Copyright 2011-2013 Karsten Ahnert

 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_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_HPP_INCLUDED

#include <boost/fusion/container/vector.hpp>
#include <boost/fusion/container/generation/make_vector.hpp>

#include <boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp>
#include <boost/numeric/odeint/algebra/range_algebra.hpp>
#include <boost/numeric/odeint/algebra/default_operations.hpp>
#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>

#include <boost/numeric/odeint/util/state_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>

#include <boost/array.hpp>

namespace boost {
namespace numeric {
namespace odeint {

#ifndef DOXYGEN_SKIP
template< class Value = double >
struct rk54_ck_coefficients_a1 : boost::array< Value , 1 >
{
    rk54_ck_coefficients_a1( void )
    {
        (*this)[0] = static_cast< Value >( 1 )/static_cast< Value >( 5 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_a2 : boost::array< Value , 2 >
{
    rk54_ck_coefficients_a2( void )
    {
        (*this)[0] = static_cast<Value>( 3 )/static_cast<Value>( 40 );
        (*this)[1] = static_cast<Value>( 9 )/static_cast<Value>( 40 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_a3 : boost::array< Value , 3 >
{
    rk54_ck_coefficients_a3( void )
    {
        (*this)[0] = static_cast<Value>( 3 )/static_cast<Value>( 10 );
        (*this)[1] = static_cast<Value>( -9 )/static_cast<Value>( 10 );
        (*this)[2] = static_cast<Value>( 6 )/static_cast<Value>( 5 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_a4 : boost::array< Value , 4 >
{
    rk54_ck_coefficients_a4( void )
    {
        (*this)[0] = static_cast<Value>( -11 )/static_cast<Value>( 54 );
        (*this)[1] = static_cast<Value>( 5 )/static_cast<Value>( 2 );
        (*this)[2] = static_cast<Value>( -70 )/static_cast<Value>( 27 );
        (*this)[3] = static_cast<Value>( 35 )/static_cast<Value>( 27 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_a5 : boost::array< Value , 5 >
{
    rk54_ck_coefficients_a5( void )
    {
        (*this)[0] = static_cast<Value>( 1631 )/static_cast<Value>( 55296 );
        (*this)[1] = static_cast<Value>( 175 )/static_cast<Value>( 512 );
        (*this)[2] = static_cast<Value>( 575 )/static_cast<Value>( 13824 );
        (*this)[3] = static_cast<Value>( 44275 )/static_cast<Value>( 110592 );
        (*this)[4] = static_cast<Value>( 253 )/static_cast<Value>( 4096 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_b : boost::array< Value , 6 >
{
    rk54_ck_coefficients_b( void )
    {
        (*this)[0] = static_cast<Value>( 37 )/static_cast<Value>( 378 );
        (*this)[1] = static_cast<Value>( 0 );
        (*this)[2] = static_cast<Value>( 250 )/static_cast<Value>( 621 );
        (*this)[3] = static_cast<Value>( 125 )/static_cast<Value>( 594 );
        (*this)[4] = static_cast<Value>( 0 );
        (*this)[5] = static_cast<Value>( 512 )/static_cast<Value>( 1771 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_db : boost::array< Value , 6 >
{
    rk54_ck_coefficients_db( void )
    {
        (*this)[0] = static_cast<Value>( 37 )/static_cast<Value>( 378 ) - static_cast<Value>( 2825 )/static_cast<Value>( 27648 );
        (*this)[1] = static_cast<Value>( 0 );
        (*this)[2] = static_cast<Value>( 250 )/static_cast<Value>( 621 ) - static_cast<Value>( 18575 )/static_cast<Value>( 48384 );
        (*this)[3] = static_cast<Value>( 125 )/static_cast<Value>( 594 ) - static_cast<Value>( 13525 )/static_cast<Value>( 55296 );
        (*this)[4] = static_cast<Value>( -277 )/static_cast<Value>( 14336 );
        (*this)[5] = static_cast<Value>( 512 )/static_cast<Value>( 1771 ) - static_cast<Value>( 1 )/static_cast<Value>( 4 );
    }
};

template< class Value = double >
struct rk54_ck_coefficients_c : boost::array< Value , 6 >
{
    rk54_ck_coefficients_c( void )
    {
        (*this)[0] = static_cast<Value>(0);
        (*this)[1] = static_cast<Value>( 1 )/static_cast<Value>( 5 );
        (*this)[2] = static_cast<Value>( 3 )/static_cast<Value>( 10 );
        (*this)[3] = static_cast<Value>( 3 )/static_cast<Value>( 5 );
        (*this)[4] = static_cast<Value>( 1 );
        (*this)[5] = static_cast<Value>( 7 )/static_cast<Value>( 8 );
    }
};
#endif

template<
    class State ,
    class Value = double ,
    class Deriv = State ,
    class Time = Value ,
    class Algebra = typename algebra_dispatcher< State >::algebra_type ,
    class Operations = typename operations_dispatcher< State >::operations_type ,
    class Resizer = initially_resizer
    >
#ifndef DOXYGEN_SKIP
class runge_kutta_cash_karp54 : public explicit_error_generic_rk< 6 , 5 , 5 , 4 ,
        State , Value , Deriv , Time , Algebra , Operations , Resizer >
#else
class runge_kutta_cash_karp54 : public explicit_error_generic_rk
#endif
{

public:
#ifndef DOXYGEN_SKIP
    typedef explicit_error_generic_rk< 6 , 5 , 5 , 4 , State , Value , Deriv , Time ,
                               Algebra , Operations , Resizer > stepper_base_type;
#endif
    typedef typename stepper_base_type::state_type state_type;
    typedef typename stepper_base_type::value_type value_type;
    typedef typename stepper_base_type::deriv_type deriv_type;
    typedef typename stepper_base_type::time_type time_type;
    typedef typename stepper_base_type::algebra_type algebra_type;
    typedef typename stepper_base_type::operations_type operations_type;
    typedef typename stepper_base_type::resizer_type resizer_typ;

    #ifndef DOXYGEN_SKIP
    typedef typename stepper_base_type::stepper_type stepper_type;
    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
    #endif

    runge_kutta_cash_karp54( const algebra_type &algebra = algebra_type() ) : stepper_base_type(
        boost::fusion::make_vector( rk54_ck_coefficients_a1<Value>() ,
                                 rk54_ck_coefficients_a2<Value>() ,
                                 rk54_ck_coefficients_a3<Value>() ,
                                 rk54_ck_coefficients_a4<Value>() ,
                                 rk54_ck_coefficients_a5<Value>() ) ,
            rk54_ck_coefficients_b<Value>() , rk54_ck_coefficients_db<Value>() , rk54_ck_coefficients_c<Value>() ,
            algebra )
    { }
};

/********** DOXYGEN **********/

/**
 * \class runge_kutta_cash_karp54
 * \brief The Runge-Kutta Cash-Karp method.
 *
 * The Runge-Kutta Cash-Karp method is one of the standard methods for
 * solving ordinary differential equations, see
 * en.wikipedia.org/wiki/Cash-Karp_methods.
 * The method is explicit and fulfills the Error Stepper concept. Step size control
 * is provided but continuous output is not available for this method.
 *
 * This class derives from explicit_error_stepper_base and inherits its interface via CRTP (current recurring template pattern).
 * Furthermore, it derivs from explicit_error_generic_rk which is a generic Runge-Kutta algorithm with error estimation.
 * For more details see explicit_error_stepper_base and explicit_error_generic_rk.
 *
 * \tparam State The state type.
 * \tparam Value The value type.
 * \tparam Deriv The type representing the time derivative of the state.
 * \tparam Time The time representing the independent variable - the time.
 * \tparam Algebra The algebra type.
 * \tparam Operations The operations type.
 * \tparam Resizer The resizer policy type.
 */

    /**
     * \fn runge_kutta_cash_karp54::runge_kutta_cash_karp54( const algebra_type &algebra )
     * \brief Constructs the runge_kutta_cash_karp54 class. This constructor can be used as a default
     * constructor if the algebra has a default constructor.
     * \param algebra A copy of algebra is made and stored inside explicit_stepper_base.
     */
}
}
}

#endif // BOOST_NUMERIC_ODEINT_STEPPER_RUNGE_KUTTA_CASH_KARP54_HPP_INCLUDED


code:
//  Copyright 2009-2010 Vicente J. Botet Escriba

//  Distributed under the Boost Software License, Version 1.0.
//  See http://www.boost.org/LICENSE_1_0.txt

#ifndef BOOST_CHRONO_DETAIL_SYSTEM_HPP
#define BOOST_CHRONO_DETAIL_SYSTEM_HPP

#if !defined BOOST_CHRONO_DONT_PROVIDE_HYBRID_ERROR_HANDLING

#include <boost/version.hpp>
#include <boost/system/error_code.hpp>

#if ((BOOST_VERSION / 100000) < 2) && ((BOOST_VERSION / 100 % 1000) < 44)
#define BOOST_CHRONO_SYSTEM_CATEGORY boost::system::system_category
#else
#define BOOST_CHRONO_SYSTEM_CATEGORY boost::system::system_category()
#endif

#ifdef BOOST_SYSTEM_NO_DEPRECATED
#define BOOST_CHRONO_THROWS boost::throws()
#define BOOST_CHRONO_IS_THROWS(EC) (&EC==&boost::throws())
#else
#define BOOST_CHRONO_THROWS boost::system::throws
#define BOOST_CHRONO_IS_THROWS(EC) (&EC==&boost::system::throws)
#endif

#endif
#endif


code:
#if !defined(BOOST_PROTO_DONT_USE_PREPROCESSED_FILES)

    #include <boost/proto/transform/detail/preprocessed/default_function_impl.hpp>

#elif !defined(BOOST_PP_IS_ITERATING)

    #define BOOST_PROTO_DEF_FUN_INVOKE_ARG(Z, M, DATA)                                              \
        BOOST_PROTO_DEFAULT_EVAL(Z, BOOST_PP_ADD(M, 2), DATA)

    #if defined(__WAVE__) && defined(BOOST_PROTO_CREATE_PREPROCESSED_FILES)
        #pragma wave option(preserve: 2, line: 0, output: "preprocessed/default_function_impl.hpp")
    #endif

    ///////////////////////////////////////////////////////////////////////////////
    /// \file default_function_impl.hpp
    /// Contains definition of the default_function_impl, the implementation of the
    /// _default transform for function-like nodes.
    //
    //  Copyright 2008 Eric Niebler. 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)

    #if defined(__WAVE__) && defined(BOOST_PROTO_CREATE_PREPROCESSED_FILES)
        #pragma wave option(preserve: 1)
    #endif

    #define BOOST_PP_ITERATION_PARAMS_1                                                             \
        (3, (3, BOOST_PROTO_MAX_ARITY, <boost/proto/transform/detail/default_function_impl.hpp>))
    #include BOOST_PP_ITERATE()

    #if defined(__WAVE__) && defined(BOOST_PROTO_CREATE_PREPROCESSED_FILES)
        #pragma wave option(output: null)
    #endif

    #undef BOOST_PROTO_DEF_FUN_INVOKE_ARG

#else

    #define N BOOST_PP_ITERATION()

    template<typename Grammar, typename Expr, typename State, typename Data>
    struct default_function_impl<Grammar, Expr, State, Data, N>
      : transform_impl<Expr, State, Data>
    {
        BOOST_PP_REPEAT(N, BOOST_PROTO_DEFAULT_EVAL_TYPE, Expr)

        typedef
            typename proto::detail::result_of_fixup<r0>::type
        function_type;

        typedef
            typename BOOST_PROTO_RESULT_OF<
                function_type(BOOST_PP_ENUM_SHIFTED_PARAMS(N, r))
            >::type
        result_type;

        result_type operator ()(
            typename default_function_impl::expr_param e
          , typename default_function_impl::state_param s
          , typename default_function_impl::data_param d
        ) const
        {
            return this->invoke(e, s, d, is_member_function_pointer<function_type>());
        }

    private:
        result_type invoke(
            typename default_function_impl::expr_param e
          , typename default_function_impl::state_param s
          , typename default_function_impl::data_param d
          , mpl::false_
        ) const
        {
            return BOOST_PROTO_DEFAULT_EVAL(~, 0, e)(
                BOOST_PP_ENUM_SHIFTED(N, BOOST_PROTO_DEFAULT_EVAL, e)
            );
        }

        result_type invoke(
            typename default_function_impl::expr_param e
          , typename default_function_impl::state_param s
          , typename default_function_impl::data_param d
          , mpl::true_
        ) const
        {
            BOOST_PROTO_USE_GET_POINTER();
            typedef typename detail::class_member_traits<function_type>::class_type class_type;
            return (
                BOOST_PROTO_GET_POINTER(class_type, (BOOST_PROTO_DEFAULT_EVAL(~, 1, e))) ->*
                BOOST_PROTO_DEFAULT_EVAL(~, 0, e)
            )(BOOST_PP_ENUM(BOOST_PP_SUB(N, 2), BOOST_PROTO_DEF_FUN_INVOKE_ARG, e));
        }
    };

    #undef N

#endif