diff --git a/include/boost/math/big_number.hpp b/include/boost/math/big_number.hpp
index f9a7c98b..434ddb79 100644
--- a/include/boost/math/big_number.hpp
+++ b/include/boost/math/big_number.hpp
@@ -17,156 +17,10 @@
 #include <boost/type_traits/is_signed.hpp>
 #include <boost/type_traits/is_unsigned.hpp>
 #include <boost/type_traits/is_floating_point.hpp>
+#include <boost/math/big_number/default_ops.hpp>
 
 namespace boost{ namespace math{
 
-template <class Backend>
-class big_number;
-
-namespace detail{
-
-// Forward-declare an expression wrapper
-template<typename Expr>
-struct big_number_exp;
-//
-// Declare our grammars:
-//
-struct big_number_grammar
-  : proto::or_<
-        proto::terminal< proto::_ >
-      , proto::plus< big_number_grammar, big_number_grammar >
-      , proto::multiplies< big_number_grammar, big_number_grammar >
-      , proto::minus< big_number_grammar, big_number_grammar >
-      , proto::divides< big_number_grammar, big_number_grammar >
-      , proto::unary_plus< big_number_grammar >
-      , proto::negate< big_number_grammar >
-      , proto::modulus<big_number_grammar, big_number_grammar>
-    >
-{};
-
-// Define a calculator domain. Expression within
-// the calculator domain will be wrapped in the
-// calculator<> expression wrapper.
-struct big_number_domain
-  : proto::domain< proto::generator<big_number_exp>, big_number_grammar>
-{};
-
-template<typename Expr>
-struct big_number_exp
-  : proto::extends<Expr, big_number_exp<Expr>, big_number_domain>
-{
-    typedef
-        proto::extends<Expr, big_number_exp<Expr>, big_number_domain> base_type;
-
-    big_number_exp(Expr const &expr = Expr())
-      : base_type(expr)
-    {}
-    template <class Other>
-    big_number_exp(const Other& o)
-       : base_type(o)
-    {}
-};
-
-struct CalcDepth
-  : proto::or_<
-        proto::when< proto::terminal<proto::_>,
-            mpl::int_<0>()
-        >
-      , proto::when< proto::unary_expr<proto::_, CalcDepth>,
-            CalcDepth(proto::_child)
-        >
-      , proto::when< proto::binary_expr<proto::_, CalcDepth, CalcDepth>,
-            mpl::plus<mpl::max<CalcDepth(proto::_left),
-                     CalcDepth(proto::_right)>, mpl::int_<1> >()
-        >
-    >
-{};
-
-template <int b>
-struct has_enough_bits
-{
-   template <class T>
-   struct type : public mpl::bool_<std::numeric_limits<T>::digits >= b>{};
-};
-
-template <class Val, class Backend, class Tag>
-struct canonical_imp
-{
-   typedef Val type;
-};
-template <class Val, class Backend>
-struct canonical_imp<Val, Backend, mpl::int_<0> >
-{
-   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
-   typedef typename mpl::find_if<
-      typename Backend::signed_types,
-      pred_type
-   >::type iter_type;
-   typedef typename mpl::deref<iter_type>::type type;
-};
-template <class Val, class Backend>
-struct canonical_imp<Val, Backend, mpl::int_<1> >
-{
-   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
-   typedef typename mpl::find_if<
-      typename Backend::unsigned_types,
-      pred_type
-   >::type iter_type;
-   typedef typename mpl::deref<iter_type>::type type;
-};
-template <class Val, class Backend>
-struct canonical_imp<Val, Backend, mpl::int_<2> >
-{
-   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
-   typedef typename mpl::find_if<
-      typename Backend::real_types,
-      pred_type
-   >::type iter_type;
-   typedef typename mpl::deref<iter_type>::type type;
-};
-template <class Val, class Backend>
-struct canonical_imp<Val, Backend, mpl::int_<3> >
-{
-   typedef const char* type;
-};
-
-template <class Val, class Backend>
-struct canonical
-{
-   typedef typename mpl::if_<
-      is_signed<Val>,
-      mpl::int_<0>,
-      typename mpl::if_<
-         is_unsigned<Val>,
-         mpl::int_<1>,
-         typename mpl::if_<
-            is_floating_point<Val>,
-            mpl::int_<2>,
-            typename mpl::if_<
-               mpl::or_<
-                  is_convertible<Val, const char*>,
-                  is_same<Val, std::string>
-               >,
-               mpl::int_<3>,
-               mpl::int_<4>
-            >::type
-         >::type
-      >::type
-   >::type tag_type;
-
-   typedef typename canonical_imp<Val, Backend, tag_type>::type type;
-};
-
-} // namespace detail
-
-//
-// Traits class, lets us know whether a backend is an integer type, otherwise assumed to be a real number type:
-//
-template <class Num>
-struct is_extended_integer : public mpl::false_ {};
-template <class Backend>
-struct is_extended_integer<big_number<Backend> > : public is_extended_integer<Backend>{};
-
 template <class Backend>
 class big_number : public detail::big_number_exp<typename proto::terminal<big_number<Backend>*>::type >
 {
@@ -207,7 +61,7 @@ public:
    template <class Exp>
    big_number& operator=(const detail::big_number_exp<Exp>& e)
    {
-      do_assign(e, typename proto::tag_of<Exp>::type());
+      do_assign(e, typename detail::assign_and_eval<Exp>::type());
       return *this;
    }
 
@@ -230,7 +84,7 @@ public:
    {
       proto::value(*this) = this;
       BOOST_ASSERT(proto::value(*this) == this);
-      do_assign(e, typename proto::tag_of<Exp>::type());
+      do_assign(e, typename detail::assign_and_eval<Exp>::type());
    }
 
 #ifndef BOOST_NO_RVALUE_REFERENCES
@@ -267,7 +121,8 @@ public:
    typename enable_if<boost::is_arithmetic<V>, big_number<Backend>& >::type 
       operator+=(const V& v)
    {
-      do_add_value(canonical_value(v), mpl::false_());
+      using big_num_default_ops::add;
+      add(m_backend, canonical_value(v));
       return *this;
    }
 
@@ -291,7 +146,8 @@ public:
    typename enable_if<boost::is_arithmetic<V>, big_number<Backend>& >::type 
       operator-=(const V& v)
    {
-      do_subtract_value(canonical_value(v), mpl::false_());
+      using big_num_default_ops::subtract;
+      subtract(m_backend, canonical_value(v));
       return *this;
    }
 
@@ -316,7 +172,8 @@ public:
    typename enable_if<boost::is_arithmetic<V>, big_number<Backend>& >::type 
       operator*=(const V& v)
    {
-      do_multiplies_value(canonical_value(v), mpl::false_());
+      using big_num_default_ops::multiply;
+      multiply(m_backend, canonical_value(v));
       return *this;
    }
 
@@ -342,7 +199,8 @@ public:
       operator%=(const V& v)
    {
       BOOST_STATIC_ASSERT_MSG(is_extended_integer<Backend>::value, "The modulus operation is only valid for integer types");
-      do_modulus_value(canonical_value(v), mpl::false_());
+      using big_num_default_ops::modulus;
+      modulus(m_backend, canonical_value(v));
       return *this;
    }
 
@@ -366,16 +224,17 @@ public:
    typename enable_if<boost::is_arithmetic<V>, big_number<Backend>& >::type 
       operator/=(const V& v)
    {
-      do_divide_value(canonical_value(v), mpl::false_());
+      using big_num_default_ops::divide;
+      divide(m_backend, canonical_value(v));
       return *this;
    }
 
    //
    // String conversion functions:
    //
-   std::string str(unsigned digits = 0)const
+   std::string str(unsigned digits = 0, bool scientific = true)const
    {
-      return m_backend.str(digits);
+      return m_backend.str(digits, scientific);
    }
    //
    // Default precision:
@@ -417,18 +276,103 @@ public:
       return m_backend;
    }
 private:
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::add_immediates&)
+   {
+      using big_num_default_ops::add;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      add(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::add_and_negate_immediates&)
+   {
+      using big_num_default_ops::add;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      add(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+      m_backend.negate();
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::subtract_immediates&)
+   {
+      using big_num_default_ops::subtract;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      subtract(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::subtract_and_negate_immediates&)
+   {
+      using big_num_default_ops::subtract;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      subtract(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+      m_backend.negate();
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::multiply_immediates&)
+   {
+      using big_num_default_ops::multiply;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      multiply(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::multiply_and_negate_immediates&)
+   {
+      using big_num_default_ops::multiply;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      multiply(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+      m_backend.negate();
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::divide_immediates&)
+   {
+      using big_num_default_ops::divide;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      divide(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::divide_and_negate_immediates&)
+   {
+      using big_num_default_ops::divide;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      divide(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+      m_backend.negate();
+   }
+
+   template <class Exp>
+   void do_assign(const Exp& e, const detail::modulus_immediates&)
+   {
+      using big_num_default_ops::modulus;
+      typedef typename proto::tag_of<typename proto::result_of::left<Exp>::type>::type left_tag;
+      typedef typename proto::tag_of<typename proto::result_of::right<Exp>::type>::type right_tag;
+      modulus(m_backend, canonical_value(underlying_value(proto::left(e), left_tag())), canonical_value(underlying_value(proto::right(e), right_tag())));
+   }
+
    template <class Exp>
    void do_assign(const Exp& e, const proto::tag::unary_plus&)
    {
       typedef typename proto::result_of::left<Exp>::type left_type;
-      do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+      do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
    }
 
    template <class Exp>
    void do_assign(const Exp& e, const proto::tag::negate&)
    {
       typedef typename proto::result_of::left<Exp>::type left_type;
-      do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+      do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
       m_backend.negate();
    }
 
@@ -461,12 +405,12 @@ private:
       }
       else if(left_depth >= right_depth)
       {
-         do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+         do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
          do_add(proto::right(e), typename proto::tag_of<right_type>::type());
       }
       else
       {
-         do_assign(proto::right(e), typename proto::tag_of<right_type>::type());
+         do_assign(proto::right(e), typename detail::assign_and_eval<right_type>::type());
          do_add(proto::left(e), typename proto::tag_of<left_type>::type());
       }
    }
@@ -500,12 +444,12 @@ private:
       }
       else if(left_depth >= right_depth)
       {
-         do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+         do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
          do_subtract(proto::right(e), typename proto::tag_of<right_type>::type());
       }
       else
       {
-         do_assign(proto::right(e), typename proto::tag_of<right_type>::type());
+         do_assign(proto::right(e), typename detail::assign_and_eval<right_type>::type());
          do_subtract(proto::left(e), typename proto::tag_of<left_type>::type());
          m_backend.negate();
       }
@@ -539,12 +483,12 @@ private:
       }
       else if(left_depth >= right_depth)
       {
-         do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+         do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
          do_multiplies(proto::right(e), typename proto::tag_of<right_type>::type());
       }
       else
       {
-         do_assign(proto::right(e), typename proto::tag_of<right_type>::type());
+         do_assign(proto::right(e), typename detail::assign_and_eval<right_type>::type());
          do_multiplies(proto::left(e), typename proto::tag_of<left_type>::type());
       }
    }
@@ -572,7 +516,7 @@ private:
       }
       else
       {
-         do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+         do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
          do_divide(proto::right(e), typename proto::tag_of<right_type>::type());
       }
    }
@@ -605,7 +549,7 @@ private:
       }
       else
       {
-         do_assign(proto::left(e), typename proto::tag_of<left_type>::type());
+         do_assign(proto::left(e), typename detail::assign_and_eval<left_type>::type());
          do_modulus(proto::right(e), typename proto::tag_of<right_type>::type());
       }
    }
@@ -617,28 +561,55 @@ private:
          m_backend = canonical_value(proto::value(e));
       }
    }
+   template <class Exp>
+   void do_assign(const Exp& e, const proto::tag::function&)
+   {
+      typedef typename proto::arity_of<Exp>::type tag_type;
+      do_assign_function(e, tag_type());
+   }
+   template <class Exp>
+   void do_assign_function(const Exp& e, const mpl::long_<1>&)
+   {
+      proto::value(proto::left(e))(&m_backend);
+   }
+   template <class Exp>
+   void do_assign_function(const Exp& e, const mpl::long_<2>&)
+   {
+      typedef typename proto::result_of::right<Exp>::type right_type;
+      typedef typename proto::tag_of<right_type>::type tag_type;
+      do_assign_function_1(proto::value(proto::left(e)), proto::right(e), tag_type());
+   }
+   template <class F, class Exp>
+   void do_assign_function_1(const F& f, const Exp& val, const proto::tag::terminal&)
+   {
+      f(&m_backend, canonical_value(proto::value(val)));
+   }
+   template <class F, class Exp, class Tag>
+   void do_assign_function_1(const F& f, const Exp& val, const Tag&)
+   {
+      big_number t(val);
+      f(&m_backend, t.backend());
+   }
+   template <class Exp>
+   void do_assign_function(const Exp& e, const mpl::long_<3>&)
+   {
+      typedef typename proto::result_of::right<Exp>::type right_type;
+      typedef typename proto::tag_of<right_type>::type tag_type;
+      typedef typename proto::result_of::child_c<Exp, 2>::type end_type;
+      typedef typename proto::tag_of<end_type>::type end_tag;
+      do_assign_function_2(proto::value(proto::left(e)), proto::right(e), proto::child_c<2>(e), tag_type(), end_tag());
+   }
+   template <class F, class Exp1, class Exp2>
+   void do_assign_function_2(const F& f, const Exp1& val1, const Exp2& val2, const proto::tag::terminal&, const proto::tag::terminal&)
+   {
+      f(&m_backend, canonical_value(proto::value(val1)), canonical_value(proto::value(val2)));
+   }
 
    template <class Exp>
    void do_add(const Exp& e, const proto::tag::terminal&)
    {
-      typedef typename proto::result_of::value<Exp>::type t1;
-      typedef typename remove_reference<t1>::type t2;
-      typedef typename remove_cv<t2>::type t3;
-      typedef typename detail::canonical<t3, Backend>::type t4;
-      typedef typename is_convertible<t4, const char*>::type tag;
-      do_add_value(canonical_value(proto::value(e)), tag());
-   }
-
-   template <class V>
-   void do_add_value(const V& v, const mpl::false_&)
-   {
-      m_backend += v;
-   }
-   template <class V>
-   void do_add_value(const V& v, const mpl::true_&)
-   {
-      self_type temp(v);
-      m_backend += temp.m_backend;
+      using big_num_default_ops::add;
+      add(m_backend, canonical_value(proto::value(e)));
    }
 
    template <class Exp>
@@ -683,24 +654,8 @@ private:
    template <class Exp>
    void do_subtract(const Exp& e, const proto::tag::terminal&)
    {
-      typedef typename proto::result_of::value<Exp>::type t1;
-      typedef typename remove_reference<t1>::type t2;
-      typedef typename remove_cv<t2>::type t3;
-      typedef typename detail::canonical<t3, Backend>::type t4;
-      typedef typename is_convertible<t4, const char*>::type tag;
-      do_subtract_value(canonical_value(proto::value(e)), tag());
-   }
-
-   template <class V>
-   void do_subtract_value(const V& v, const mpl::false_&)
-   {
-      m_backend -= v;
-   }
-   template <class V>
-   void do_subtract_value(const V& v, const mpl::true_&)
-   {
-      self_type temp(v);
-      m_backend -= temp.m_backend;
+      using big_num_default_ops::subtract;
+      subtract(m_backend, canonical_value(proto::value(e)));
    }
 
    template <class Exp>
@@ -745,25 +700,8 @@ private:
    template <class Exp>
    void do_multiplies(const Exp& e, const proto::tag::terminal&)
    {
-      typedef typename proto::result_of::value<Exp>::type t1;
-      typedef typename remove_reference<t1>::type t2;
-      typedef typename remove_cv<t2>::type t3;
-      typedef typename detail::canonical<t3, Backend>::type t4;
-      typedef typename is_convertible<t4, const char*>::type tag;
-      do_multiplies_value(canonical_value(proto::value(e)), tag());
-   }
-
-   template <class Val>
-   void do_multiplies_value(const Val& v, const mpl::false_&)
-   {
-      m_backend *= v;
-   }
-
-   template <class Val>
-   void do_multiplies_value(const Val& e, const mpl::true_&)
-   {
-      self_type temp(e);
-      m_backend *= temp.m_backend;
+      using big_num_default_ops::multiply;
+      multiply(m_backend, canonical_value(proto::value(e)));
    }
 
    template <class Exp>
@@ -802,32 +740,16 @@ private:
    template <class Exp, class unknown>
    void do_multiplies(const Exp& e, const unknown&)
    {
+      using big_num_default_ops::multiply;
       self_type temp(e);
-      m_backend *= temp.m_backend;
+      multiply(m_backend, temp.m_backend);
    }
 
    template <class Exp>
    void do_divide(const Exp& e, const proto::tag::terminal&)
    {
-      typedef typename proto::result_of::value<Exp>::type t1;
-      typedef typename remove_reference<t1>::type t2;
-      typedef typename remove_cv<t2>::type t3;
-      typedef typename detail::canonical<t3, Backend>::type t4;
-      typedef typename is_convertible<t4, const char*>::type tag;
-      do_divide_value(canonical_value(proto::value(e)), tag());
-   }
-
-   template <class Val>
-   void do_divide_value(const Val& v, const mpl::false_&)
-   {
-      m_backend /= v;
-   }
-
-   template <class Val>
-   void do_divide_value(const Val& e, const mpl::true_&)
-   {
-      self_type temp(e);
-      m_backend /= temp.m_backend;
+      using big_num_default_ops::divide;
+      divide(m_backend, canonical_value(proto::value(e)));
    }
 
    template <class Exp>
@@ -866,39 +788,24 @@ private:
    template <class Exp, class unknown>
    void do_divide(const Exp& e, const unknown&)
    {
+      using big_num_default_ops::multiply;
       self_type temp(e);
-      m_backend /= temp.m_backend;
+      divide(m_backend, temp.m_backend);
    }
 
    template <class Exp>
    void do_modulus(const Exp& e, const proto::tag::terminal&)
    {
-      typedef typename proto::result_of::value<Exp>::type t1;
-      typedef typename remove_reference<t1>::type t2;
-      typedef typename remove_cv<t2>::type t3;
-      typedef typename detail::canonical<t3, Backend>::type t4;
-      typedef typename is_convertible<t4, const char*>::type tag;
-      do_modulus_value(canonical_value(proto::value(e)), tag());
-   }
-
-   template <class Val>
-   void do_modulus_value(const Val& v, const mpl::false_&)
-   {
-      m_backend %= v;
-   }
-
-   template <class Val>
-   void do_modulus_value(const Val& e, const mpl::true_&)
-   {
-      self_type temp(e);
-      m_backend %= temp.m_backend;
+      using big_num_default_ops::modulus;
+      modulus(m_backend, canonical_value(proto::value(e)));
    }
 
    template <class Exp, class Unknown>
    void do_modulus(const Exp& e, const Unknown&)
    {
+      using big_num_default_ops::modulus;
       self_type temp(e);
-      do_modulus_value(canonical_value(proto::value(temp)), mpl::false_());
+      modulus(m_backend, canonical_value(proto::value(temp)));
    }
 
    // Tests if the expression contains a reference to *this:
@@ -953,6 +860,19 @@ private:
    { 
       return v == this; 
    }
+   template <class Exp>
+   static typename detail::underlying_result<Exp>::type underlying_value(const big_number_exp<Exp>& e, const proto::tag::terminal&)
+   {
+      return proto::value(e);
+   }
+   template <class Exp, class tag>
+   static typename detail::underlying_result<Exp>::type 
+      underlying_value(const big_number_exp<Exp>& e, const tag&)
+   {
+      typedef typename proto::result_of::left<Exp>::type left_type;
+      typedef typename proto::tag_of<left_type>::type tag_type;
+      return underlying_value(proto::left(e), tag_type());
+   }
 
    static const Backend& canonical_value(const self_type& v){  return v.m_backend;  }
    static const Backend& canonical_value(const self_type* v){  return v->m_backend;  }
@@ -967,72 +887,6 @@ private:
 namespace detail
 {
 
-template <class Exp1, class Exp2>
-struct combine_expression_type
-{
-   typedef void type;
-};
-
-template <class Backend>
-struct combine_expression_type<boost::math::big_number<Backend>, boost::math::big_number<Backend> >
-{
-   typedef boost::math::big_number<Backend> type;
-};
-
-template <class Backend, class Exp>
-struct combine_expression_type<boost::math::big_number<Backend>, Exp>
-{
-   typedef boost::math::big_number<Backend> type;
-};
-
-template <class Backend, class Exp>
-struct combine_expression_type<Exp, boost::math::big_number<Backend> >
-{
-   typedef boost::math::big_number<Backend> type;
-};
-
-template <class T>
-struct is_big_number : public mpl::false_{};
-template <class T>
-struct is_big_number<boost::math::big_number<T> > : public mpl::true_{};
-template <class T>
-struct is_big_number_exp : public mpl::false_{};
-template <class T>
-struct is_big_number_exp<boost::math::detail::big_number_exp<T> > : public mpl::true_{};
-
-
-template <class Exp, int arity>
-struct expression_type_imp;
-
-template <class Exp>
-struct expression_type_imp<Exp, 0>
-{
-   typedef typename remove_pointer<typename proto::result_of::value<Exp>::type>::type type;
-};
-
-template <class Exp>
-struct expression_type_imp<Exp, 1>
-{
-   typedef typename proto::result_of::left<Exp>::type nested_type;
-   typedef typename expression_type_imp<nested_type, proto::arity_of<nested_type>::value>::type type;
-};
-
-template <class Exp>
-struct expression_type_imp<Exp, 2>
-{
-   typedef typename proto::result_of::left<Exp>::type left_type;
-   typedef typename proto::result_of::right<Exp>::type right_type;
-   typedef typename expression_type_imp<left_type, proto::arity_of<left_type>::value>::type left_result;
-   typedef typename expression_type_imp<right_type, proto::arity_of<right_type>::value>::type right_result;
-   typedef typename combine_expression_type<left_result, right_result>::type type;
-};
-
-template <class Exp>
-struct expression_type
-{
-   typedef typename expression_type_imp<Exp, proto::arity_of<Exp>::value>::type type;
-};
-
 template <class Backend>
 inline int big_number_compare(const big_number<Backend>& a, const big_number<Backend>& b)
 {
@@ -1164,13 +1018,20 @@ inline typename boost::enable_if<detail::is_valid_comparison<Exp1, Exp2>, bool>:
 }
 
 template <class Backend>
-std::ostream& operator << (std::ostream& os, const big_number<Backend>& r)
+inline std::ostream& operator << (std::ostream& os, const big_number<Backend>& r)
 {
-   return os << r.str(static_cast<unsigned>(os.precision()));
+   return os << r.str(static_cast<unsigned>(os.precision(), os.flags() & os.scientific));
+}
+template <class Exp>
+inline std::ostream& operator << (std::ostream& os, const detail::big_number_exp<Exp>& r)
+{
+   typedef typename detail::expression_type<Exp>::type value_type;
+   value_type temp(r);
+   return os << temp;
 }
 
 template <class Backend>
-std::istream& operator >> (std::istream& is, big_number<Backend>& r)
+inline std::istream& operator >> (std::istream& is, big_number<Backend>& r)
 {
    std::string s;
    is >> s;
@@ -1178,47 +1039,6 @@ std::istream& operator >> (std::istream& is, big_number<Backend>& r)
    return is;
 }
 
-//
-// Non-member functions accepting an expression-template as argument:
-//
-#undef sqrt
-template <class Exp>
-typename boost::math::detail::expression_type<Exp>::type sqrt(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return sqrt(result_type(val));
-}
-template <class Exp>
-typename detail::expression_type<Exp>::type abs(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return abs(result_type(val));
-}
-template <class Exp>
-typename detail::expression_type<Exp>::type fabs(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return fabs(result_type(val));
-}
-template <class Exp>
-typename detail::expression_type<Exp>::type ceil(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return ceil(result_type(val));
-}
-template <class Exp>
-typename detail::expression_type<Exp>::type floor(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return floor(result_type(val));
-}
-template <class Exp>
-typename detail::expression_type<Exp>::type trunc(const detail::big_number_exp<Exp>& val)
-{
-   typedef typename detail::expression_type<Exp>::type result_type;
-   return trunc(result_type(val));
-}
-
 }} // namespaces
 
 #endif
diff --git a/include/boost/math/big_number/big_number_base.hpp b/include/boost/math/big_number/big_number_base.hpp
new file mode 100644
index 00000000..b3caed30
--- /dev/null
+++ b/include/boost/math/big_number/big_number_base.hpp
@@ -0,0 +1,463 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_BIG_NUM_BASE_HPP
+#define BOOST_MATH_BIG_NUM_BASE_HPP
+
+namespace boost{ namespace math{
+
+template <class Backend>
+class big_number;
+
+namespace detail{
+
+// Forward-declare an expression wrapper
+template<typename Expr>
+struct big_number_exp;
+//
+// Declare our grammars:
+//
+struct big_number_grammar
+  : proto::or_<
+        proto::terminal< proto::_ >
+      , proto::plus< big_number_grammar, big_number_grammar >
+      , proto::multiplies< big_number_grammar, big_number_grammar >
+      , proto::minus< big_number_grammar, big_number_grammar >
+      , proto::divides< big_number_grammar, big_number_grammar >
+      , proto::unary_plus< big_number_grammar >
+      , proto::negate< big_number_grammar >
+      , proto::modulus<big_number_grammar, big_number_grammar>
+    >
+{};
+
+// Define a calculator domain. Expression within
+// the calculator domain will be wrapped in the
+// calculator<> expression wrapper.
+struct big_number_domain
+  : proto::domain< proto::generator<big_number_exp>, big_number_grammar>
+{};
+
+template<typename Expr>
+struct big_number_exp
+  : proto::extends<Expr, big_number_exp<Expr>, big_number_domain>
+{
+    typedef
+        proto::extends<Expr, big_number_exp<Expr>, big_number_domain> base_type;
+
+    big_number_exp(Expr const &expr = Expr())
+      : base_type(expr)
+    {}
+    template <class Other>
+    big_number_exp(const Other& o)
+       : base_type(o)
+    {}
+};
+
+struct CalcDepth
+  : proto::or_<
+        proto::when< proto::terminal<proto::_>,
+            mpl::int_<0>()
+        >
+      , proto::when< proto::unary_expr<proto::_, CalcDepth>,
+            CalcDepth(proto::_child)
+        >
+      , proto::when< proto::binary_expr<proto::_, CalcDepth, CalcDepth>,
+            mpl::plus<mpl::max<CalcDepth(proto::_left),
+                     CalcDepth(proto::_right)>, mpl::int_<1> >()
+        >
+    >
+{};
+
+template <int b>
+struct has_enough_bits
+{
+   template <class T>
+   struct type : public mpl::bool_<std::numeric_limits<T>::digits >= b>{};
+};
+
+template <class Val, class Backend, class Tag>
+struct canonical_imp
+{
+   typedef Val type;
+};
+template <class Val, class Backend>
+struct canonical_imp<Val, Backend, mpl::int_<0> >
+{
+   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
+   typedef typename mpl::find_if<
+      typename Backend::signed_types,
+      pred_type
+   >::type iter_type;
+   typedef typename mpl::deref<iter_type>::type type;
+};
+template <class Val, class Backend>
+struct canonical_imp<Val, Backend, mpl::int_<1> >
+{
+   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
+   typedef typename mpl::find_if<
+      typename Backend::unsigned_types,
+      pred_type
+   >::type iter_type;
+   typedef typename mpl::deref<iter_type>::type type;
+};
+template <class Val, class Backend>
+struct canonical_imp<Val, Backend, mpl::int_<2> >
+{
+   typedef typename has_enough_bits<std::numeric_limits<Val>::digits>::template type<mpl::_> pred_type;
+   typedef typename mpl::find_if<
+      typename Backend::real_types,
+      pred_type
+   >::type iter_type;
+   typedef typename mpl::deref<iter_type>::type type;
+};
+template <class Val, class Backend>
+struct canonical_imp<Val, Backend, mpl::int_<3> >
+{
+   typedef const char* type;
+};
+
+template <class Val, class Backend>
+struct canonical
+{
+   typedef typename mpl::if_<
+      is_signed<Val>,
+      mpl::int_<0>,
+      typename mpl::if_<
+         is_unsigned<Val>,
+         mpl::int_<1>,
+         typename mpl::if_<
+            is_floating_point<Val>,
+            mpl::int_<2>,
+            typename mpl::if_<
+               mpl::or_<
+                  is_convertible<Val, const char*>,
+                  is_same<Val, std::string>
+               >,
+               mpl::int_<3>,
+               mpl::int_<4>
+            >::type
+         >::type
+      >::type
+   >::type tag_type;
+
+   typedef typename canonical_imp<Val, Backend, tag_type>::type type;
+};
+
+template <class Exp, class tag>
+struct assign_and_eval_imp
+{
+   typedef tag type;
+};
+
+struct add_immediates{};
+struct add_and_negate_immediates{};
+struct subtract_immediates{};
+struct subtract_and_negate_immediates{};
+struct multiply_immediates{};
+struct multiply_and_negate_immediates{};
+struct divide_immediates{};
+struct divide_and_negate_immediates{};
+struct modulus_immediates{};
+
+struct immediate{};
+struct negative_immediate{};
+
+template <class Exp, class tag>
+struct immediate_type
+{
+   typedef tag type;
+};
+template <class Exp>
+struct immediate_type<Exp, proto::tag::terminal>
+{
+   typedef immediate type;
+};
+template <class Exp>
+struct immediate_type<Exp, proto::tag::unary_plus>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename mpl::if_<
+      is_same<left_tag, proto::tag::terminal>,
+      immediate,
+      left_tag
+      >::type type;
+};
+template <class Exp>
+struct immediate_type<Exp, proto::tag::negate>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename immediate_type<left_type, left_tag>::type tag;
+   typedef typename mpl::if_<
+      is_same<immediate, tag>,
+      negative_immediate,
+      tag
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval_imp<Exp, proto::tag::plus>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename proto::tag_of<right_type>::type right_tag;
+   typedef typename immediate_type<left_type, left_tag>::type left_imm;
+   typedef typename immediate_type<right_type, right_tag>::type right_imm;
+   typedef typename mpl::if_<
+      mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, immediate> >,
+      add_immediates,
+      typename mpl::if_<
+         mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, negative_immediate> >,
+         subtract_immediates,
+         typename mpl::if_<
+            mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, immediate> >,
+            subtract_and_negate_immediates,
+            typename mpl::if_<
+               mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, negative_immediate> >,
+               add_and_negate_immediates,
+               proto::tag::plus
+            >::type
+         >::type
+      >::type
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval_imp<Exp, proto::tag::minus>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename proto::tag_of<right_type>::type right_tag;
+   typedef typename immediate_type<left_type, left_tag>::type left_imm;
+   typedef typename immediate_type<right_type, right_tag>::type right_imm;
+   typedef typename mpl::if_<
+      mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, immediate> >,
+      subtract_immediates,
+      typename mpl::if_<
+         mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, negative_immediate> >,
+         add_immediates,
+         typename mpl::if_<
+            mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, immediate> >,
+            add_and_negate_immediates,
+            typename mpl::if_<
+               mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, negative_immediate> >,
+               subtract_and_negate_immediates,
+               proto::tag::minus
+            >::type
+         >::type
+      >::type
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval_imp<Exp, proto::tag::multiplies>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename proto::tag_of<right_type>::type right_tag;
+   typedef typename immediate_type<left_type, left_tag>::type left_imm;
+   typedef typename immediate_type<right_type, right_tag>::type right_imm;
+   typedef typename mpl::if_<
+      mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, immediate> >,
+      multiply_immediates,
+      typename mpl::if_<
+         mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, negative_immediate> >,
+         multiply_and_negate_immediates,
+         typename mpl::if_<
+            mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, immediate> >,
+            multiply_and_negate_immediates,
+            typename mpl::if_<
+               mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, negative_immediate> >,
+               multiply_immediates,
+               proto::tag::multiplies
+            >::type
+         >::type
+      >::type
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval_imp<Exp, proto::tag::divides>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename proto::tag_of<right_type>::type right_tag;
+   typedef typename immediate_type<left_type, left_tag>::type left_imm;
+   typedef typename immediate_type<right_type, right_tag>::type right_imm;
+   typedef typename mpl::if_<
+      mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, immediate> >,
+      divide_immediates,
+      typename mpl::if_<
+         mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, negative_immediate> >,
+         divide_and_negate_immediates,
+         typename mpl::if_<
+            mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, immediate> >,
+            divide_and_negate_immediates,
+            typename mpl::if_<
+               mpl::and_<is_same<left_imm, negative_immediate>, is_same<right_imm, negative_immediate> >,
+               divide_immediates,
+               proto::tag::divides
+            >::type
+         >::type
+      >::type
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval_imp<Exp, proto::tag::modulus>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::tag_of<left_type>::type left_tag;
+   typedef typename proto::tag_of<right_type>::type right_tag;
+   typedef typename immediate_type<left_type, left_tag>::type left_imm;
+   typedef typename immediate_type<right_type, right_tag>::type right_imm;
+   typedef typename mpl::if_<
+      mpl::and_<is_same<left_imm, immediate>, is_same<right_imm, immediate> >,
+      modulus_immediates,
+      proto::tag::modulus
+   >::type type;
+};
+
+template <class Exp>
+struct assign_and_eval
+{
+   typedef typename proto::tag_of<Exp>::type tag_type;
+   typedef typename assign_and_eval_imp<Exp, tag_type>::type type;
+};
+
+template <class Exp, class tag>
+struct underlying_result_imp
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::tag_of<left_type>::type tag_type;
+   typedef typename underlying_result_imp<left_type, tag_type>::type type;
+};
+
+template <class Exp>
+struct underlying_result_imp<Exp, proto::tag::terminal>
+{
+   typedef typename proto::result_of::value<Exp const &>::type type;
+};
+
+template <class Exp>
+struct underlying_result
+{
+   typedef typename proto::tag_of<Exp>::type tag_type;
+   typedef typename underlying_result_imp<Exp, tag_type>::type type;
+};
+
+template <class Exp1, class Exp2>
+struct combine_expression_type
+{
+   typedef void type;
+};
+
+template <class Backend>
+struct combine_expression_type<boost::math::big_number<Backend>, boost::math::big_number<Backend> >
+{
+   typedef boost::math::big_number<Backend> type;
+};
+
+template <class Backend, class Exp>
+struct combine_expression_type<boost::math::big_number<Backend>, Exp>
+{
+   typedef boost::math::big_number<Backend> type;
+};
+
+template <class Backend, class Exp>
+struct combine_expression_type<Exp, boost::math::big_number<Backend> >
+{
+   typedef boost::math::big_number<Backend> type;
+};
+
+template <class T>
+struct is_big_number : public mpl::false_{};
+template <class T>
+struct is_big_number<boost::math::big_number<T> > : public mpl::true_{};
+template <class T>
+struct is_big_number_exp : public mpl::false_{};
+template <class T>
+struct is_big_number_exp<boost::math::detail::big_number_exp<T> > : public mpl::true_{};
+
+
+template <class Exp, int arity>
+struct expression_type_imp;
+
+template <class Exp>
+struct expression_type_imp<Exp, 0>
+{
+   typedef typename remove_pointer<typename proto::result_of::value<Exp>::type>::type type;
+};
+
+template <class Exp>
+struct expression_type_imp<Exp, 1>
+{
+   typedef typename proto::result_of::left<Exp>::type nested_type;
+   typedef typename expression_type_imp<nested_type, proto::arity_of<nested_type>::value>::type type;
+};
+
+template <class Exp>
+struct expression_type_imp<Exp, 2>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename expression_type_imp<left_type, proto::arity_of<left_type>::value>::type left_result;
+   typedef typename expression_type_imp<right_type, proto::arity_of<right_type>::value>::type right_result;
+   typedef typename combine_expression_type<left_result, right_result>::type type;
+};
+
+template <class Exp>
+struct expression_type_imp<Exp, 3>
+{
+   typedef typename proto::result_of::left<Exp>::type left_type;
+   typedef typename proto::result_of::right<Exp>::type right_type;
+   typedef typename proto::result_of::child_c<Exp, 2>::type end_type;
+   typedef typename expression_type_imp<left_type, proto::arity_of<left_type>::value>::type left_result;
+   typedef typename expression_type_imp<right_type, proto::arity_of<right_type>::value>::type right_result;
+   typedef typename expression_type_imp<end_type, proto::arity_of<end_type>::value>::type end_result;
+   typedef typename combine_expression_type<left_result, typename combine_expression_type<right_result, end_type>::type>::type type;
+};
+
+template <class Exp>
+struct expression_type
+{
+   typedef typename expression_type_imp<Exp, proto::arity_of<Exp>::value>::type type;
+};
+
+template <class Exp>
+struct backend_type
+{
+   typedef typename expression_type<Exp>::type num_type;
+   typedef typename backend_type<num_type>::type type;
+};
+
+template <class Backend>
+struct backend_type<boost::math::big_number<Backend> >
+{
+   typedef Backend type;
+};
+
+} // namespace detail
+
+//
+// Traits class, lets us know whether a backend is an integer type, otherwise assumed to be a real number type:
+//
+template <class Num>
+struct is_extended_integer : public mpl::false_ {};
+template <class Backend>
+struct is_extended_integer<big_number<Backend> > : public is_extended_integer<Backend>{};
+
+
+}} // namespaces
+
+#endif // BOOST_MATH_BIG_NUM_BASE_HPP
+
+
diff --git a/include/boost/math/big_number/default_ops.hpp b/include/boost/math/big_number/default_ops.hpp
new file mode 100644
index 00000000..569a32f6
--- /dev/null
+++ b/include/boost/math/big_number/default_ops.hpp
@@ -0,0 +1,302 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2011 John Maddock. Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_BIG_NUM_DEF_OPS
+#define BOOST_MATH_BIG_NUM_DEF_OPS
+
+#include <boost/math/big_number/big_number_base.hpp>
+
+namespace boost{ namespace math{ namespace big_num_default_ops{
+
+//
+// Default versions of mixed arithmetic, these just construct a temporary
+// from the arithmetic value and then do the arithmetic on that:
+//
+template <class T, class V>
+inline typename enable_if<mpl::or_<is_arithmetic<V>, is_convertible<V, const char*>, is_same<V, std::string> > >::type 
+   add(T& result, V const& v)
+{
+   T t;
+   t = v;
+   add(result, t);
+}
+template <class T, class V>
+inline typename enable_if<mpl::or_<is_arithmetic<V>, is_convertible<V, const char*>, is_same<V, std::string> > >::type
+   subtract(T& result, V const& v)
+{
+   T t;
+   t = v;
+   subtract(result, t);
+}
+template <class T, class V>
+inline typename enable_if<mpl::or_<is_arithmetic<V>, is_convertible<V, const char*>, is_same<V, std::string> > >::type
+   multiply(T& result, V const& v)
+{
+   T t;
+   t = v;
+   multiply(result, t);
+}
+template <class T, class V>
+inline typename enable_if<mpl::or_<is_arithmetic<V>, is_convertible<V, const char*>, is_same<V, std::string> > >::type
+   divide(T& result, V const& v)
+{
+   T t;
+   t = v;
+   divide(result, t);
+}
+template <class T, class V>
+inline typename enable_if<mpl::or_<is_arithmetic<V>, is_convertible<V, const char*>, is_same<V, std::string> > >::type
+   modulus(T& result, V const& v)
+{
+   T t;
+   t = v;
+   mudulus(result, t);
+}
+
+template <class T>
+inline bool is_same_object(const T& u, const T&v)
+{  return &u == &v;  }
+template <class T, class U>
+inline bool is_same_object(const T& u, const U&v)
+{  return false;  }
+
+//
+// Default versions of 3-arg arithmetic functions, these just forward to the 2 arg versions:
+//
+template <class T, class U, class V>
+inline void add(T& t, const U& u, const V& v)
+{
+   if(is_same_object(t, v))
+   {
+      add(t, u);
+   }
+   else if(is_same_object(t, u))
+   {
+      add(t, v);
+   }
+   else
+   {
+      t = u;
+      add(t, v);
+   }
+}
+template <class T, class U, class V>
+inline void subtract(T& t, const U& u, const V& v)
+{
+   if(is_same_object(t, u))
+      subtract(t, v);
+   else if(is_same_object(t, v))
+   {
+      subtract(t, u);
+      t.negate();
+   }
+   else
+   {
+      t = u;
+      subtract(t, v);
+   }
+}
+template <class T, class U, class V>
+inline void multiply(T& t, const U& u, const V& v)
+{
+   if(is_same_object(t, u))
+      multiply(t, v);
+   else if(is_same_object(t, v))
+      multiply(t, u);
+   else
+   {
+      t = u;
+      multiply(t, v);
+   }
+}
+template <class T, class U, class V>
+inline void divide(T& t, const U& u, const V& v)
+{
+   if(is_same_object(t, u))
+      divide(t, v);
+   else if(is_same_object(t, v))
+   {
+      T temp = t;
+      divide(temp, u, v);
+      temp.swap(t);
+   }
+   else
+   {
+      t = u;
+      divide(t, v);
+   }
+}
+template <class T, class U, class V>
+inline void modulus(T& t, const U& u, const V& v)
+{
+   if(is_same_object(t, u))
+      modulus(t, v);
+   else if(is_same_object(t, v))
+   {
+      T temp;
+      modulus(temp, u, v);
+      temp.swap(t);
+   }
+   else
+   {
+      t = u;
+      modulus(t, v);
+   }
+}
+
+//
+// Functions:
+//
+template <class T>
+void abs(T* result, const T& arg)
+{
+   typedef typename T::signed_types type_list;
+   typedef typename mpl::front<type_list>::type front;
+   *result = arg;
+   if(arg.compare(front(0)) < 0)
+      result->negate();
+}
+template <class T>
+void fabs(T* result, const T& arg)
+{
+   typedef typename T::signed_types type_list;
+   typedef typename mpl::front<type_list>::type front;
+   *result = arg;
+   if(arg.compare(front(0)) < 0)
+      result->negate();
+}
+
+//
+// These have to implemented by the backend, declared here so that our macro generated code compiles OK.
+//
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type floor();
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type ceil();
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type trunc();
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type sqrt();
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type ldexp();
+template <class T>
+typename enable_if_c<sizeof(T) == 0>::type frexp();
+
+}
+
+
+template <class Backend>
+class big_number;
+
+namespace detail{
+
+template<typename Expr>
+struct big_number_exp;
+
+}
+
+#define UNARY_OP_FUNCTOR(func)\
+namespace detail{\
+template <class Backend>\
+struct BOOST_JOIN(func, _funct)\
+{\
+   void operator()(Backend* result, const Backend& arg)const\
+   {\
+      using big_num_default_ops::func;\
+      func(result, arg);\
+   }\
+};\
+\
+}\
+\
+template <class Exp>\
+typename proto::result_of::make_expr<\
+    proto::tag::function\
+  , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<Exp>::type>\
+  , detail::big_number_exp<Exp>\
+>::type const \
+func(const detail::big_number_exp<Exp>& arg)\
+{\
+    return proto::make_expr<proto::tag::function>(\
+        detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<Exp>::type>() \
+      , arg   \
+    );\
+}\
+template <class Backend>\
+typename proto::result_of::make_expr<\
+    proto::tag::function\
+  , detail::BOOST_JOIN(func, _funct)<Backend>\
+  , detail::big_number_exp<typename proto::terminal<big_number<Backend>*>::type>\
+>::type const \
+func(const big_number<Backend>& arg)\
+{\
+    return proto::make_expr<proto::tag::function>(\
+        detail::BOOST_JOIN(func, _funct)<Backend>() \
+      , static_cast<const detail::big_number_exp<typename proto::terminal<big_number<Backend>*>::type>&>(arg)   \
+    );\
+}
+
+#define BINARY_OP_FUNCTOR(func)\
+namespace detail{\
+template <class Backend>\
+struct BOOST_JOIN(func, _funct)\
+{\
+   template <class A2>\
+   void operator()(Backend* result, const Backend& arg, const A2& a)const\
+   {\
+      using big_num_default_ops::func;\
+      func(result, arg, a);\
+   }\
+};\
+\
+}\
+\
+template <class Exp, class A2>\
+typename proto::result_of::make_expr<\
+    proto::tag::function\
+  , detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<Exp>::type>\
+  , detail::big_number_exp<Exp>\
+  , typename proto::result_of::as_child<const A2&>::type\
+>::type const \
+func(const detail::big_number_exp<Exp>& arg, const A2& a)\
+{\
+    return proto::make_expr<proto::tag::function>(\
+        detail::BOOST_JOIN(func, _funct)<typename detail::backend_type<Exp>::type>() \
+      , arg, proto::as_child(a)   \
+    );\
+}\
+template <class Backend, class A2>\
+typename proto::result_of::make_expr<\
+    proto::tag::function\
+  , detail::BOOST_JOIN(func, _funct)<Backend>\
+  , detail::big_number_exp<typename proto::terminal<big_number<Backend>*>::type>\
+  , typename proto::result_of::as_child<const A2&>::type\
+>::type const \
+func(const big_number<Backend>& arg, const A2& a)\
+{\
+    return proto::make_expr<proto::tag::function>(\
+        detail::BOOST_JOIN(func, _funct)<Backend>() \
+      , static_cast<const detail::big_number_exp<typename proto::terminal<big_number<Backend>*>::type>&>(arg),\
+      proto::as_child(a)\
+    );\
+}
+
+UNARY_OP_FUNCTOR(abs)
+UNARY_OP_FUNCTOR(fabs)
+UNARY_OP_FUNCTOR(sqrt)
+UNARY_OP_FUNCTOR(floor)
+UNARY_OP_FUNCTOR(ceil)
+UNARY_OP_FUNCTOR(trunc)
+
+BINARY_OP_FUNCTOR(ldexp)
+BINARY_OP_FUNCTOR(frexp)
+
+#undef BINARY_OP_FUNCTOR
+#undef UNARY_OP_FUNCTOR
+
+}} // namespaces
+
+#endif
+
diff --git a/include/boost/math/big_number/gmp.hpp b/include/boost/math/big_number/gmp.hpp
index 95fddb04..20dc9203 100644
--- a/include/boost/math/big_number/gmp.hpp
+++ b/include/boost/math/big_number/gmp.hpp
@@ -145,130 +145,45 @@ struct gmp_real_imp
       mpf_set_str(m_data, s, 10);
       return *this;
    }
-   gmp_real_imp& operator += (const gmp_real<digits10>& o)
-   {
-      mpf_add(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_real_imp& operator += (V v)
-   {
-      gmp_real<digits10> d;
-      d = v;
-      return *this += d;
-   }
-   gmp_real_imp& operator -= (const gmp_real<digits10>& o)
-   {
-      mpf_sub(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_real_imp& operator -= (V v)
-   {
-      gmp_real<digits10> d;
-      d = v;
-      return *this -= d;
-   }
-   gmp_real_imp& operator *= (const gmp_real<digits10>& o)
-   {
-      mpf_mul(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_real_imp& operator *= (V v)
-   {
-      gmp_real<digits10> d;
-      d = v;
-      return *this *= d;
-   }
-   gmp_real_imp& operator /= (const gmp_real<digits10>& o)
-   {
-      mpf_div(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_real_imp& operator /= (V v)
-   {
-      gmp_real<digits10> d;
-      d = v;
-      return *this /= d;
-   }
-   gmp_real_imp& operator += (unsigned long i)
-   {
-      mpf_add_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_real_imp& operator -= (unsigned long i)
-   {
-      mpf_sub_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_real_imp& operator *= (unsigned long i)
-   {
-      mpf_mul_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_real_imp& operator /= (unsigned long i)
-   {
-      mpf_div_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_real_imp& operator += (long i)
-   {
-      if(i > 0)
-         mpf_add_ui(m_data, m_data, i);
-      else
-         mpf_sub_ui(m_data, m_data, std::abs(i));
-      return *this;
-   }
-   gmp_real_imp& operator -= (long i)
-   {
-      if(i > 0)
-         mpf_sub_ui(m_data, m_data, i);
-      else
-         mpf_add_ui(m_data, m_data, std::abs(i));
-      return *this;
-   }
-   gmp_real_imp& operator *= (long i)
-   {
-      mpf_mul_ui(m_data, m_data, std::abs(i));
-      if(i < 0)
-         mpf_neg(m_data, m_data);
-      return *this;
-   }
-   gmp_real_imp& operator /= (long i)
-   {
-      mpf_div_ui(m_data, m_data, std::abs(i));
-      if(i < 0)
-         mpf_neg(m_data, m_data);
-      return *this;
-   }
    void swap(gmp_real_imp& o)
    {
       mpf_swap(m_data, o.m_data);
    }
-   std::string str(unsigned digits)const
+   std::string str(unsigned digits, bool scientific)const
    {
+      std::string result;
       mp_exp_t e;
       void *(*alloc_func_ptr) (size_t);
       void *(*realloc_func_ptr) (void *, size_t, size_t);
       void (*free_func_ptr) (void *, size_t);
       const char* ps = mpf_get_str (0, &e, 10, digits, m_data);
-      std::string s("0.");
-      if(ps[0] == '-')
+      std::ptrdiff_t sl = std::strlen(ps);
+      if(sl == 0)
+         return "0";
+      if(*ps == '-')
+         --sl; // number of digits excluding sign.
+      if(!scientific
+         && (sl <= std::numeric_limits<boost::uintmax_t>::digits10 + 1)
+         && (e >= sl)
+         && (sl <= std::numeric_limits<boost::uintmax_t>::digits10 + 1))
       {
-         s.insert(0, ps, 1);
-         s += ps + 1;
+         result = ps;
+         result.append(e-sl, '0');
       }
       else
       {
-         s += ps;
+         result = ps;
+         if(ps[0] == '-')
+            result.insert(2, 1, '.');
+         else
+            result.insert(1, 1, '.');
+         --e;
+         if(e)
+            result += "e" + lexical_cast<std::string>(e);
       }
-      s += "e";
-      s += boost::lexical_cast<std::string>(e);
       mp_get_memory_functions(&alloc_func_ptr, &realloc_func_ptr, &free_func_ptr);
       (*free_func_ptr)((void*)ps, std::strlen(ps) + 1);
-      return s;
+      return result;
    }
    ~gmp_real_imp()
    {
@@ -400,74 +315,280 @@ private:
    }
 };
 
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& result, const gmp_real<digits10>& o)
+{
+   mpf_add(result.data(), result.data(), o.data());
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& result, const gmp_real<digits10>& o)
+{
+   mpf_sub(result.data(), result.data(), o.data());
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& result, const gmp_real<digits10>& o)
+{
+   mpf_mul(result.data(), result.data(), o.data());
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& result, const gmp_real<digits10>& o)
+{
+   mpf_div(result.data(), result.data(), o.data());
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& result, unsigned long i)
+{
+   mpf_add_ui(result.data(), result.data(), i);
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& result, unsigned long i)
+{
+   mpf_sub_ui(result.data(), result.data(), i);
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& result, unsigned long i)
+{
+   mpf_mul_ui(result.data(), result.data(), i);
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& result, unsigned long i)
+{
+   mpf_div_ui(result.data(), result.data(), i);
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& result, long i)
+{
+   if(i > 0)
+      mpf_add_ui(result.data(), result.data(), i);
+   else
+      mpf_sub_ui(result.data(), result.data(), std::abs(i));
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& result, long i)
+{
+   if(i > 0)
+      mpf_sub_ui(result.data(), result.data(), i);
+   else
+      mpf_add_ui(result.data(), result.data(), std::abs(i));
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& result, long i)
+{
+   mpf_mul_ui(result.data(), result.data(), std::abs(i));
+   if(i < 0)
+      mpf_neg(result.data(), result.data());
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& result, long i)
+{
+   mpf_div_ui(result.data(), result.data(), std::abs(i));
+   if(i < 0)
+      mpf_neg(result.data(), result.data());
+}
+//
+// Specialised 3 arg versions of the basic operators:
+//
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& a, const gmp_real<digits10>& x, const gmp_real<digits10>& y)
+{
+   mpf_add(a.data(), x.data(), y.data());
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& a, const gmp_real<digits10>& x, unsigned long y)
+{
+   mpf_add_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& a, const gmp_real<digits10>& x, long y)
+{
+   if(y < 0)
+      mpf_sub_ui(a.data(), x.data(), -y);
+   else
+      mpf_add_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& a, unsigned long x, const gmp_real<digits10>& y)
+{
+   mpf_add_ui(a.data(), y.data(), x);
+}
+template <unsigned digits10>
+inline void add(gmp_real<digits10>& a, long x, const gmp_real<digits10>& y)
+{
+   if(x < 0)
+   {
+      mpf_ui_sub(a.data(), -x, y.data());
+      mpf_neg(a.data(), a.data());
+   }
+   else
+      mpf_add_ui(a.data(), y.data(), x);
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& a, const gmp_real<digits10>& x, const gmp_real<digits10>& y)
+{
+   mpf_sub(a.data(), x.data(), y.data());
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& a, const gmp_real<digits10>& x, unsigned long y)
+{
+   mpf_sub_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& a, const gmp_real<digits10>& x, long y)
+{
+   if(y < 0)
+      mpf_add_ui(a.data(), x.data(), -y);
+   else
+      mpf_sub_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& a, unsigned long x, const gmp_real<digits10>& y)
+{
+   mpf_ui_sub(a.data(), x, y.data());
+}
+template <unsigned digits10>
+inline void subtract(gmp_real<digits10>& a, long x, const gmp_real<digits10>& y)
+{
+   if(x < 0)
+   {
+      mpf_add_ui(a.data(), y.data(), -x);
+      mpf_neg(a.data(), a.data());
+   }
+   else
+      mpf_ui_sub(a.data(), x, y.data());
+}
+
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& a, const gmp_real<digits10>& x, const gmp_real<digits10>& y)
+{
+   mpf_mul(a.data(), x.data(), y.data());
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& a, const gmp_real<digits10>& x, unsigned long y)
+{
+   mpf_mul_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& a, const gmp_real<digits10>& x, long y)
+{
+   if(y < 0)
+   {
+      mpf_mul_ui(a.data(), x.data(), -y);
+      a.negate();
+   }
+   else
+      mpf_mul_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& a, unsigned long x, const gmp_real<digits10>& y)
+{
+   mpf_mul_ui(a.data(), y.data(), x);
+}
+template <unsigned digits10>
+inline void multiply(gmp_real<digits10>& a, long x, const gmp_real<digits10>& y)
+{
+   if(x < 0)
+   {
+      mpf_mul_ui(a.data(), y.data(), -x);
+      mpf_neg(a.data(), a.data());
+   }
+   else
+      mpf_mul_ui(a.data(), y.data(), x);
+}
+
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& a, const gmp_real<digits10>& x, const gmp_real<digits10>& y)
+{
+   mpf_div(a.data(), x.data(), y.data());
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& a, const gmp_real<digits10>& x, unsigned long y)
+{
+   mpf_div_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& a, const gmp_real<digits10>& x, long y)
+{
+   if(y < 0)
+   {
+      mpf_div_ui(a.data(), x.data(), -y);
+      a.negate();
+   }
+   else
+      mpf_div_ui(a.data(), x.data(), y);
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& a, unsigned long x, const gmp_real<digits10>& y)
+{
+   mpf_ui_div(a.data(), x, y.data());
+}
+template <unsigned digits10>
+inline void divide(gmp_real<digits10>& a, long x, const gmp_real<digits10>& y)
+{
+   if(x < 0)
+   {
+      mpf_ui_div(a.data(), -x, y.data());
+      mpf_neg(a.data(), a.data());
+   }
+   else
+      mpf_ui_div(a.data(), x, y.data());
+}
+
 //
 // Native non-member operations:
 //
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > sqrt(const big_number<gmp_real<Digits10> >& val)
+inline void sqrt(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
 {
-   big_number<gmp_real<Digits10> > result;
-   mpf_sqrt(result.backend().data(), val.backend().data());
-   return result;
+   mpf_sqrt(result->data(), val.data());
+}
+
+template <unsigned Digits10>
+inline void abs(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
+{
+   mpf_abs(result->data(), val.data());
+}
+
+template <unsigned Digits10>
+inline void fabs(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
+{
+   mpf_abs(result->data(), val.data());
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > abs(const big_number<gmp_real<Digits10> >& val)
+inline void ceil(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
 {
-   big_number<gmp_real<Digits10> > result;
-   mpf_abs(result.backend().data(), val.backend().data());
-   return result;
+   mpf_ceil(result->data(), val.data());
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > fabs(const big_number<gmp_real<Digits10> >& val)
+inline void floor(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
 {
-   big_number<gmp_real<Digits10> > result;
-   mpf_abs(result.backend().data(), val.backend().data());
-   return result;
+   mpf_floor(result->data(), val.data());
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > ceil(const big_number<gmp_real<Digits10> >& val)
+inline void trunc(gmp_real<Digits10>* result, const gmp_real<Digits10>& val)
 {
-   big_number<gmp_real<Digits10> > result;
-   mpf_ceil(result.backend().data(), val.backend().data());
-   return result;
+   mpf_trunc(result->data(), val.data());
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > floor(const big_number<gmp_real<Digits10> >& val)
+inline void ldexp(gmp_real<Digits10>* result, const gmp_real<Digits10>& val, long e)
 {
-   big_number<gmp_real<Digits10> > result;
-   mpf_floor(result.backend().data(), val.backend().data());
-   return result;
-}
-template <unsigned Digits10>
-big_number<gmp_real<Digits10> > trunc(const big_number<gmp_real<Digits10> >& val)
-{
-   big_number<gmp_real<Digits10> > result;
-   mpf_trunc(result.backend().data(), val.backend().data());
-   return result;
-}
-template <unsigned Digits10>
-big_number<gmp_real<Digits10> > ldexp(const big_number<gmp_real<Digits10> >& val, long e)
-{
-   big_number<gmp_real<Digits10> > result;
    if(e > 0)
-      mpf_mul_2exp(result.backend().data(), val.backend().data(), e);
+      mpf_mul_2exp(result->data(), val.data(), e);
    else if(e < 0)
-      mpf_div_2exp(result.backend().data(), val.backend().data(), -e);
-   return result;
+      mpf_div_2exp(result->data(), val.data(), -e);
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > frexp(const big_number<gmp_real<Digits10> >& val, int* e)
+inline void frexp(gmp_real<Digits10>* result, const gmp_real<Digits10>& val, int* e)
 {
    long v;
-   mpf_get_d_2exp(&v, val.backend().data());
+   mpf_get_d_2exp(&v, val.data());
    *e = v;
-   return ldexp(val, -v);
+   return ldexp(result, val, -v);
 }
 template <unsigned Digits10>
-big_number<gmp_real<Digits10> > frexp(const big_number<gmp_real<Digits10> >& val, long* e)
+inline void frexp(gmp_real<Digits10>* result, const gmp_real<Digits10>& val, long* e)
 {
-   mpf_get_d_2exp(e, val.backend().data());
-   return ldexp(val, -*e);
+   mpf_get_d_2exp(e, val.data());
+   return ldexp(result, val, -*e);
 }
 
 struct gmp_int
@@ -582,163 +703,11 @@ struct gmp_int
       mpz_set_str(m_data, s, 10);
       return *this;
    }
-   gmp_int& operator += (const gmp_int& o)
-   {
-      mpz_add(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_int& operator += (V v)
-   {
-      gmp_int d;
-      d = v;
-      return *this += d;
-   }
-   gmp_int& operator -= (const gmp_int& o)
-   {
-      mpz_sub(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_int& operator -= (V v)
-   {
-      gmp_int d;
-      d = v;
-      return *this -= d;
-   }
-   gmp_int& operator *= (const gmp_int& o)
-   {
-      mpz_mul(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_int& operator *= (V v)
-   {
-      gmp_int d;
-      d = v;
-      return *this *= d;
-   }
-   gmp_int& operator /= (const gmp_int& o)
-   {
-      mpz_div(m_data, m_data, o.m_data);
-      return *this;
-   }
-   template <class V>
-   gmp_int& operator /= (V v)
-   {
-      gmp_int d;
-      d = v;
-      return *this /= d;
-   }
-   gmp_int& operator %= (const gmp_int& o)
-   {
-      bool neg = mpz_sgn(m_data) < 0;
-      bool neg2 = mpz_sgn(o.m_data) < 0;
-      mpz_mod(m_data, m_data, o.m_data);
-      if(neg)
-      {
-         if(!neg2)
-            negate();
-         mpz_add(m_data, m_data, o.m_data);
-         if(!neg2)
-            negate();
-      }
-      return *this;
-   }
-   template <class V>
-   gmp_int& operator %= (V v)
-   {
-      gmp_int d;
-      d = v;
-      return *this %= d;
-   }
-   gmp_int& operator += (unsigned long i)
-   {
-      mpz_add_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_int& operator -= (unsigned long i)
-   {
-      mpz_sub_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_int& operator *= (unsigned long i)
-   {
-      mpz_mul_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_int& operator %= (unsigned long i)
-   {
-      bool neg = mpz_sgn(m_data) < 0;
-      mpz_mod_ui(m_data, m_data, i);
-      if(neg)
-      {
-         negate();
-         mpz_add_ui(m_data, m_data, i);
-         negate();
-      }
-      return *this;
-   }
-   gmp_int& operator /= (unsigned long i)
-   {
-      mpz_div_ui(m_data, m_data, i);
-      return *this;
-   }
-   gmp_int& operator += (long i)
-   {
-      if(i > 0)
-         mpz_add_ui(m_data, m_data, i);
-      else
-         mpz_sub_ui(m_data, m_data, std::abs(i));
-      return *this;
-   }
-   gmp_int& operator -= (long i)
-   {
-      if(i > 0)
-         mpz_sub_ui(m_data, m_data, i);
-      else
-         mpz_add_ui(m_data, m_data, std::abs(i));
-      return *this;
-   }
-   gmp_int& operator *= (long i)
-   {
-      mpz_mul_ui(m_data, m_data, std::abs(i));
-      if(i < 0)
-         mpz_neg(m_data, m_data);
-      return *this;
-   }
-   gmp_int& operator %= (long i)
-   {
-      bool neg = mpz_sgn(m_data) < 0;
-      bool neg2 = i < 0;
-      mpz_mod_ui(m_data, m_data, std::abs(i));
-      if(neg)
-      {
-         if(!neg2)
-         {
-            negate();
-            mpz_add_ui(m_data, m_data, std::abs(i));
-            negate();
-         }
-         else
-         {
-            mpz_sub_ui(m_data, m_data, std::abs(i));
-         }
-      }
-      return *this;
-   }
-   gmp_int& operator /= (long i)
-   {
-      mpz_div_ui(m_data, m_data, std::abs(i));
-      if(i < 0)
-         mpz_neg(m_data, m_data);
-      return *this;
-   }
    void swap(gmp_int& o)
    {
       mpz_swap(m_data, o.m_data);
    }
-   std::string str(unsigned)const
+   std::string str(unsigned /*digits*/, bool /*scientific*/)const
    {
       void *(*alloc_func_ptr) (size_t);
       void *(*realloc_func_ptr) (void *, size_t, size_t);
@@ -776,10 +745,218 @@ struct gmp_int
       d = v;
       return compare(d);
    }
+   mpz_t& data() { return m_data; }
+   const mpz_t& data()const { return m_data; }
 protected:
    mpz_t m_data;
 };
 
+inline void add(gmp_int& t, const gmp_int& o)
+{
+   mpz_add(t.data(), t.data(), o.data());
+}
+inline void subtract(gmp_int& t, const gmp_int& o)
+{
+   mpz_sub(t.data(), t.data(), o.data());
+}
+inline void multiply(gmp_int& t, const gmp_int& o)
+{
+   mpz_mul(t.data(), t.data(), o.data());
+}
+inline void divide(gmp_int& t, const gmp_int& o)
+{
+   mpz_div(t.data(), t.data(), o.data());
+}
+inline void modulus(gmp_int& t, const gmp_int& o)
+{
+   bool neg = mpz_sgn(t.data()) < 0;
+   bool neg2 = mpz_sgn(o.data()) < 0;
+   mpz_mod(t.data(), t.data(), o.data());
+   if(neg)
+   {
+      if(!neg2)
+         t.negate();
+      mpz_add(t.data(), t.data(), o.data());
+      if(!neg2)
+         t.negate();
+   }
+}
+inline void add(gmp_int& t, unsigned long i)
+{
+   mpz_add_ui(t.data(), t.data(), i);
+}
+inline void subtract(gmp_int& t, unsigned long i)
+{
+   mpz_sub_ui(t.data(), t.data(), i);
+}
+inline void multiply(gmp_int& t, unsigned long i)
+{
+   mpz_mul_ui(t.data(), t.data(), i);
+}
+inline void modulus(gmp_int& t, unsigned long i)
+{
+   bool neg = mpz_sgn(t.data()) < 0;
+   mpz_mod_ui(t.data(), t.data(), i);
+   if(neg)
+   {
+      t.negate();
+      mpz_add_ui(t.data(), t.data(), i);
+      t.negate();
+   }
+}
+inline void divide(gmp_int& t, unsigned long i)
+{
+   mpz_div_ui(t.data(), t.data(), i);
+}
+inline void add(gmp_int& t, long i)
+{
+   if(i > 0)
+      mpz_add_ui(t.data(), t.data(), i);
+   else
+      mpz_sub_ui(t.data(), t.data(), -i);
+}
+inline void subtract(gmp_int& t, long i)
+{
+   if(i > 0)
+      mpz_sub_ui(t.data(), t.data(), i);
+   else
+      mpz_add_ui(t.data(), t.data(), -i);
+}
+inline void multiply(gmp_int& t, long i)
+{
+   mpz_mul_ui(t.data(), t.data(), std::abs(i));
+   if(i < 0)
+      mpz_neg(t.data(), t.data());
+}
+inline void modulus(gmp_int& t, long i)
+{
+   bool neg = mpz_sgn(t.data()) < 0;
+   bool neg2 = i < 0;
+   mpz_mod_ui(t.data(), t.data(), std::abs(i));
+   if(neg)
+   {
+      if(!neg2)
+      {
+         t.negate();
+         mpz_add_ui(t.data(), t.data(), std::abs(i));
+         t.negate();
+      }
+      else
+      {
+         mpz_sub_ui(t.data(), t.data(), std::abs(i));
+      }
+   }
+}
+inline void divide(gmp_int& t, long i)
+{
+   mpz_div_ui(t.data(), t.data(), std::abs(i));
+   if(i < 0)
+      mpz_neg(t.data(), t.data());
+}
+
+inline void add(gmp_int& t, const gmp_int& p, const gmp_int& o)
+{
+   mpz_add(t.data(), p.data(), o.data());
+}
+inline void subtract(gmp_int& t, const gmp_int& p, const gmp_int& o)
+{
+   mpz_sub(t.data(), p.data(), o.data());
+}
+inline void multiply(gmp_int& t, const gmp_int& p, const gmp_int& o)
+{
+   mpz_mul(t.data(), p.data(), o.data());
+}
+inline void divide(gmp_int& t, const gmp_int& p, const gmp_int& o)
+{
+   mpz_div(t.data(), p.data(), o.data());
+}
+inline void modulus(gmp_int& t, const gmp_int& p, const gmp_int& o)
+{
+   bool neg = mpz_sgn(p.data()) < 0;
+   bool neg2 = mpz_sgn(o.data()) < 0;
+   mpz_mod(t.data(), p.data(), o.data());
+   if(neg)
+   {
+      if(!neg2)
+         t.negate();
+      mpz_add(t.data(), t.data(), o.data());
+      if(!neg2)
+         t.negate();
+   }
+}
+inline void add(gmp_int& t, const gmp_int& p, unsigned long i)
+{
+   mpz_add_ui(t.data(), p.data(), i);
+}
+inline void subtract(gmp_int& t, const gmp_int& p, unsigned long i)
+{
+   mpz_sub_ui(t.data(), p.data(), i);
+}
+inline void multiply(gmp_int& t, const gmp_int& p, unsigned long i)
+{
+   mpz_mul_ui(t.data(), p.data(), i);
+}
+inline void modulus(gmp_int& t, const gmp_int& p, unsigned long i)
+{
+   bool neg = mpz_sgn(p.data()) < 0;
+   mpz_mod_ui(t.data(), p.data(), i);
+   if(neg)
+   {
+      t.negate();
+      mpz_add_ui(t.data(), t.data(), i);
+      t.negate();
+   }
+}
+inline void divide(gmp_int& t, const gmp_int& p, unsigned long i)
+{
+   mpz_div_ui(t.data(), p.data(), i);
+}
+inline void add(gmp_int& t, const gmp_int& p, long i)
+{
+   if(i > 0)
+      mpz_add_ui(t.data(), p.data(), i);
+   else
+      mpz_sub_ui(t.data(), p.data(), -i);
+}
+inline void subtract(gmp_int& t, const gmp_int& p, long i)
+{
+   if(i > 0)
+      mpz_sub_ui(t.data(), p.data(), i);
+   else
+      mpz_add_ui(t.data(), p.data(), -i);
+}
+inline void multiply(gmp_int& t, const gmp_int& p, long i)
+{
+   mpz_mul_ui(t.data(), p.data(), std::abs(i));
+   if(i < 0)
+      mpz_neg(t.data(), t.data());
+}
+inline void modulus(gmp_int& t, const gmp_int& p, long i)
+{
+   bool neg = mpz_sgn(p.data()) < 0;
+   bool neg2 = i < 0;
+   mpz_mod_ui(t.data(), p.data(), std::abs(i));
+   if(neg)
+   {
+      if(!neg2)
+      {
+         t.negate();
+         mpz_add_ui(t.data(), t.data(), std::abs(i));
+         t.negate();
+      }
+      else
+      {
+         mpz_sub_ui(t.data(), t.data(), std::abs(i));
+      }
+   }
+}
+inline void divide(gmp_int& t, const gmp_int& p, long i)
+{
+   mpz_div_ui(t.data(), p.data(), std::abs(i));
+   if(i < 0)
+      mpz_neg(t.data(), t.data());
+}
+   
 template<>
 struct is_extended_integer<gmp_int> : public mpl::true_ {};
 
diff --git a/include/boost/math/concepts/big_number_architypes.hpp b/include/boost/math/concepts/big_number_architypes.hpp
index 31b4d992..6067b6e6 100644
--- a/include/boost/math/concepts/big_number_architypes.hpp
+++ b/include/boost/math/concepts/big_number_architypes.hpp
@@ -18,9 +18,9 @@ namespace concepts{
 
 struct big_number_backend_real_architype
 {
-   typedef mpl::list<int, long long>                 signed_types;
-   typedef mpl::list<unsigned, unsigned long long>   unsigned_types;
-   typedef mpl::list<double, long double>            real_types;
+   typedef mpl::list<long long>                 signed_types;
+   typedef mpl::list<unsigned long long>        unsigned_types;
+   typedef mpl::list<long double>               real_types;
 
    big_number_backend_real_architype()
    {
@@ -49,24 +49,6 @@ struct big_number_backend_real_architype
       std::cout << "Int Assignment (" << i << ")" << std::endl;
       return *this;
    }
-   big_number_backend_real_architype& operator = (unsigned i)
-   {
-      m_value = i;
-      std::cout << "UInt Assignment (" << i << ")" << std::endl;
-      return *this;
-   }
-   big_number_backend_real_architype& operator = (int i)
-   {
-      m_value = i;
-      std::cout << "Int Assignment (" << i << ")" << std::endl;
-      return *this;
-   }
-   big_number_backend_real_architype& operator = (double d)
-   {
-      m_value = d;
-      std::cout << "double Assignment (" << d << ")" << std::endl;
-      return *this;
-   }
    big_number_backend_real_architype& operator = (long double d)
    {
       m_value = d;
@@ -79,182 +61,29 @@ struct big_number_backend_real_architype
       std::cout << "const char* Assignment (" << s << ")" << std::endl;
       return *this;
    }
-   big_number_backend_real_architype& operator += (const big_number_backend_real_architype& o)
-   {
-      std::cout << "Addition (" << m_value << " += " << o.m_value << ")" << std::endl;
-      m_value += o.m_value;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (const big_number_backend_real_architype& o)
-   {
-      std::cout << "Subtraction (" << m_value << " -= " << o.m_value << ")" << std::endl;
-      m_value -= o.m_value;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (const big_number_backend_real_architype& o)
-   {
-      std::cout << "Multiplication (" << m_value << " *= " << o.m_value << ")" << std::endl;
-      m_value *= o.m_value;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (const big_number_backend_real_architype& o)
-   {
-      std::cout << "Division (" << m_value << " /= " << o.m_value << ")" << std::endl;
-      m_value /= o.m_value;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (boost::uintmax_t i)
-   {
-      std::cout << "UIntmax_t Addition (" << m_value << " += " << i << ")" << std::endl;
-      m_value += i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (boost::uintmax_t i)
-   {
-      std::cout << "UIntmax_t Subtraction (" << m_value << " -= " << i << ")" << std::endl;
-      m_value -= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (boost::uintmax_t i)
-   {
-      std::cout << "UIntmax_t Multiplication (" << m_value << " *= " << i << ")" << std::endl;
-      m_value *= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (boost::uintmax_t i)
-   {
-      std::cout << "UIntmax_t Division (" << m_value << " /= " << i << ")" << std::endl;
-      m_value /= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (boost::intmax_t i)
-   {
-      std::cout << "Intmax_t Addition (" << m_value << " += " << i << ")" << std::endl;
-      m_value += i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (boost::intmax_t i)
-   {
-      std::cout << "Intmax_t Subtraction (" << m_value << " -= " << i << ")" << std::endl;
-      m_value -= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (boost::intmax_t i)
-   {
-      std::cout << "Intmax_t Multiplication (" << m_value << " *= " << i << ")" << std::endl;
-      m_value *= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (boost::intmax_t i)
-   {
-      std::cout << "Intmax_t Division (" << m_value << " /= " << i << ")" << std::endl;
-      m_value /= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (unsigned i)
-   {
-      std::cout << "UInt Addition (" << m_value << " += " << i << ")" << std::endl;
-      m_value += i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (unsigned i)
-   {
-      std::cout << "UInt Subtraction (" << m_value << " -= " << i << ")" << std::endl;
-      m_value -= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (unsigned i)
-   {
-      std::cout << "UInt Multiplication (" << m_value << " *= " << i << ")" << std::endl;
-      m_value *= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (unsigned i)
-   {
-      std::cout << "UInt Division (" << m_value << " /= " << i << ")" << std::endl;
-      m_value /= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (int i)
-   {
-      std::cout << "Int Addition (" << m_value << " += " << i << ")" << std::endl;
-      m_value += i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (int i)
-   {
-      std::cout << "Int Subtraction (" << m_value << " -= " << i << ")" << std::endl;
-      m_value -= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (int i)
-   {
-      std::cout << "Int Multiplication (" << m_value << " *= " << i << ")" << std::endl;
-      m_value *= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (int i)
-   {
-      std::cout << "Int Division (" << m_value << " /= " << i << ")" << std::endl;
-      m_value /= i;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (double d)
-   {
-      std::cout << "double Addition (" << m_value << " += " << d << ")" << std::endl;
-      m_value += d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (double d)
-   {
-      std::cout << "double Subtraction (" << m_value << " -= " << d << ")" << std::endl;
-      m_value -= d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (double d)
-   {
-      std::cout << "double Multiplication (" << m_value << " *= " << d << ")" << std::endl;
-      m_value *= d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (double d)
-   {
-      std::cout << "double Division (" << m_value << " /= " << d << ")" << std::endl;
-      m_value /= d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator += (long double d)
-   {
-      std::cout << "long double Addition (" << m_value << " += " << d << ")" << std::endl;
-      m_value += d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator -= (long double d)
-   {
-      std::cout << "long double Subtraction (" << m_value << " -= " << d << ")" << std::endl;
-      m_value -= d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator *= (long double d)
-   {
-      std::cout << "long double Multiplication (" << m_value << " *= " << d << ")" << std::endl;
-      m_value *= d;
-      return *this;
-   }
-   big_number_backend_real_architype& operator /= (long double d)
-   {
-      std::cout << "long double Division (" << m_value << " /= " << d << ")" << std::endl;
-      m_value /= d;
-      return *this;
-   }
    void swap(big_number_backend_real_architype& o)
    {
       std::cout << "Swapping (" << m_value << " with " << o.m_value << ")" << std::endl;
       std::swap(m_value, o.m_value);
    }
-   std::string str(unsigned)const
+   std::string str(unsigned digits, bool scientific)const
    {
-      std::string s(boost::lexical_cast<std::string>(m_value));
+      std::stringstream ss;
+      if(scientific)
+         ss.setf(ss.scientific);
+      if(digits)
+         ss.precision(digits);
+      else
+         ss.precision(std::numeric_limits<long double>::digits10 + 2);
+      boost::intmax_t i = m_value;
+      boost::uintmax_t u = m_value;
+      if(!scientific && m_value == i)
+         ss << i;
+      else if(!scientific && m_value == u)
+         ss << u;
+      else
+         ss << m_value;
+      std::string s = ss.str();
       std::cout << "Converting to string (" << s << ")" << std::endl;
       return s;
    }
@@ -268,16 +97,6 @@ struct big_number_backend_real_architype
       std::cout << "Comparison" << std::endl;
       return m_value > o.m_value ? 1 : (m_value < o.m_value ? -1 : 0);
    }
-   int compare(int i)const
-   {
-      std::cout << "Comparison with int" << std::endl;
-      return m_value > i ? 1 : (m_value < i ? -1 : 0);
-   }
-   int compare(unsigned i)const
-   {
-      std::cout << "Comparison with unsigned" << std::endl;
-      return m_value > i ? 1 : (m_value < i ? -1 : 0);
-   }
    int compare(boost::intmax_t i)const
    {
       std::cout << "Comparison with int" << std::endl;
@@ -288,20 +107,35 @@ struct big_number_backend_real_architype
       std::cout << "Comparison with unsigned" << std::endl;
       return m_value > i ? 1 : (m_value < i ? -1 : 0);
    }
-   int compare(double d)const
-   {
-      std::cout << "Comparison with double" << std::endl;
-      return m_value > d ? 1 : (m_value < d ? -1 : 0);
-   }
    int compare(long double d)const
    {
       std::cout << "Comparison with long double" << std::endl;
       return m_value > d ? 1 : (m_value < d ? -1 : 0);
    }
-private:
    long double m_value;
 };
 
+inline void add(big_number_backend_real_architype& result, const big_number_backend_real_architype& o)
+{
+   std::cout << "Addition (" << result.m_value << " += " << o.m_value << ")" << std::endl;
+   result.m_value += o.m_value;
+}
+inline void subtract(big_number_backend_real_architype& result, const big_number_backend_real_architype& o)
+{
+   std::cout << "Subtraction (" << result.m_value << " -= " << o.m_value << ")" << std::endl;
+   result.m_value -= o.m_value;
+}
+inline void multiply(big_number_backend_real_architype& result, const big_number_backend_real_architype& o)
+{
+   std::cout << "Multiplication (" << result.m_value << " *= " << o.m_value << ")" << std::endl;
+   result.m_value *= o.m_value;
+}
+inline void divide(big_number_backend_real_architype& result, const big_number_backend_real_architype& o)
+{
+   std::cout << "Division (" << result.m_value << " /= " << o.m_value << ")" << std::endl;
+   result.m_value /= o.m_value;
+}
+
 typedef boost::math::big_number<big_number_backend_real_architype> big_number_real_architype;
 
 }}} // namespaces
diff --git a/math/test/linpack-benchmark.cpp b/math/test/linpack-benchmark.cpp
index bf5300fb..89ba6099 100644
--- a/math/test/linpack-benchmark.cpp
+++ b/math/test/linpack-benchmark.cpp
@@ -18,6 +18,7 @@ http://www.netlib.org/f2c/libf2c.zip
 #include <iostream>
 #include <iomanip>
 #include <cmath>
+
 #ifdef TEST_BIG_NUMBER
 #include <boost/math/big_number/gmp.hpp>
 typedef boost::math::mpf_real_100 real_type;
diff --git a/math/test/test_arithmetic.cpp b/math/test/test_arithmetic.cpp
index 99436325..9aeefdb2 100644
--- a/math/test/test_arithmetic.cpp
+++ b/math/test/test_arithmetic.cpp
@@ -50,6 +50,16 @@ void test_integer_ops(const boost::mpl::true_&)
    BOOST_TEST(a % -7 == 20 % -7);
    BOOST_TEST(-a % 7u == -20 % 7);
 
+   b = -b;
+   BOOST_TEST(a % b == 20 % -7);
+   a = -a;
+   BOOST_TEST(a % b == -20 % -7);
+   BOOST_TEST(a % -7 == -20 % -7);
+   b = 7;
+   BOOST_TEST(a % b == -20 % 7);
+   BOOST_TEST(a % 7 == -20 % 7);
+   BOOST_TEST(a % 7u == -20 % 7);
+
    a = 20;
    a %= b;
    BOOST_TEST(a == 20 % 7);
@@ -83,6 +93,7 @@ void test_real_ops(const boost::mpl::false_&){}
 template <class Real>
 void test_real_ops(const boost::mpl::true_&)
 {
+#if defined(TEST_MPF) || defined(TEST_MPF_50)
    std::cout << "Root2 = " << sqrt(Real(2)) << std::endl;
    BOOST_TEST(abs(Real(2)) == 2);
    BOOST_TEST(abs(Real(-2)) == 2);
@@ -92,12 +103,8 @@ void test_real_ops(const boost::mpl::true_&)
    BOOST_TEST(ceil(Real(5) / 2) == 3);
    BOOST_TEST(floor(Real(-5) / 2) == -3);
    BOOST_TEST(ceil(Real(-5) / 2) == -2);
-#ifndef TEST_E_FLOAT
    BOOST_TEST(trunc(Real(5) / 2) == 2);
    BOOST_TEST(trunc(Real(-5) / 2) == -2);
-#endif
-
-#ifndef TEST_E_FLOAT
    //
    // ldexp and frexp, these pretty much have to implemented by each backend:
    //
@@ -108,6 +115,7 @@ void test_real_ops(const boost::mpl::true_&)
    Real r = frexp(v, &exp);
    BOOST_TEST(r == 0.5);
    BOOST_TEST(exp == 10);
+   BOOST_TEST(v == 512);
    v = 1 / v;
    r = frexp(v, &exp);
    BOOST_TEST(r == 0.5);
@@ -115,9 +123,24 @@ void test_real_ops(const boost::mpl::true_&)
 #endif
 }
 
+template <class T>
+struct lexical_cast_target_type
+{
+   typedef typename boost::mpl::if_<
+      boost::is_signed<T>,
+      boost::intmax_t,
+      typename boost::mpl::if_<
+         boost::is_unsigned<T>,
+         boost::uintmax_t,
+         T
+      >::type
+   >::type type;
+};
+
 template <class Real, class Num>
 void test_negative_mixed(boost::mpl::true_ const&)
 {
+   typedef typename lexical_cast_target_type<Num>::type target_type;
    std::cout << "Testing mixed arithmetic with type: " << typeid(Real).name()  << " and " << typeid(Num).name() << std::endl;
    Num n1 = -static_cast<Num>(1uLL << (std::numeric_limits<Num>::digits - 1));
    Num n2 = -1;
@@ -133,10 +156,10 @@ void test_negative_mixed(boost::mpl::true_ const&)
    BOOST_TEST(n2 == Real(n2));
    BOOST_TEST(n3 == Real(n3));
    BOOST_TEST(n4 == Real(n4));
-   BOOST_TEST(n1 == boost::lexical_cast<long double>(Real(n1)));
-   BOOST_TEST(n2 == boost::lexical_cast<long double>(Real(n2)));
-   BOOST_TEST(n3 == boost::lexical_cast<long double>(Real(n3)));
-   BOOST_TEST(n4 == boost::lexical_cast<long double>(Real(n4)));
+   BOOST_TEST(n1 == boost::lexical_cast<target_type>(Real(n1).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n2 == boost::lexical_cast<target_type>(Real(n2).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n3 == boost::lexical_cast<target_type>(Real(n3).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n4 == boost::lexical_cast<target_type>(Real(n4).str(0, boost::is_floating_point<Num>::value)));
    // Assignment:
    Real r(0);
    BOOST_TEST(r != n1);
@@ -172,6 +195,21 @@ void test_negative_mixed(boost::mpl::true_ const&)
    BOOST_TEST(Real(r / n5) == n1 / n5);
    r /= n5;
    BOOST_TEST(r == n1 / n5);
+   //
+   // Extra cases for full coverage:
+   //
+   r = Real(n4) + n5;
+   BOOST_TEST(r == n4 + n5);
+   r = n4 + Real(n5);
+   BOOST_TEST(r == n4 + n5);
+   r = Real(n4) - n5;
+   BOOST_TEST(r == n4 - n5);
+   r = n4 - Real(n5);
+   BOOST_TEST(r == n4 - n5);
+   r = n4 * Real(n5);
+   BOOST_TEST(r == n4 * n5);
+   r = (4 * n4) / Real(4);
+   BOOST_TEST(r == n4);
 }
 
 template <class Real, class Num>
@@ -182,6 +220,7 @@ void test_negative_mixed(boost::mpl::false_ const&)
 template <class Real, class Num>
 void test_mixed()
 {
+   typedef typename lexical_cast_target_type<Num>::type target_type;
    std::cout << "Testing mixed arithmetic with type: " << typeid(Real).name()  << " and " << typeid(Num).name() << std::endl;
    Num n1 = static_cast<Num>(1uLL << (std::numeric_limits<Num>::digits - 1));
    Num n2 = 1;
@@ -197,10 +236,10 @@ void test_mixed()
    BOOST_TEST(n2 == Real(n2));
    BOOST_TEST(n3 == Real(n3));
    BOOST_TEST(n4 == Real(n4));
-   BOOST_TEST(n1 == boost::lexical_cast<long double>(Real(n1)));
-   BOOST_TEST(n2 == boost::lexical_cast<long double>(Real(n2)));
-   BOOST_TEST(n3 == boost::lexical_cast<long double>(Real(n3)));
-   BOOST_TEST(n4 == boost::lexical_cast<long double>(Real(n4)));
+   BOOST_TEST(n1 == boost::lexical_cast<target_type>(Real(n1).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n2 == boost::lexical_cast<target_type>(Real(n2).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n3 == boost::lexical_cast<target_type>(Real(n3).str(0, boost::is_floating_point<Num>::value)));
+   BOOST_TEST(n4 == boost::lexical_cast<target_type>(Real(n4).str(0, boost::is_floating_point<Num>::value)));
    // Assignment:
    Real r(0);
    BOOST_TEST(r != n1);
@@ -237,6 +276,18 @@ void test_mixed()
    r /= n5;
    BOOST_TEST(r == n1 / n5);
 
+   //
+   // special cases for full coverage:
+   //
+   r = n5 + Real(n4);
+   BOOST_TEST(r == n4 + n5);
+   r = n4 - Real(n5);
+   BOOST_TEST(r == n4 - n5);
+   r = n4 * Real(n5);
+   BOOST_TEST(r == n4 * n5);
+   r = (4 * n4) / Real(4);
+   BOOST_TEST(r == n4);
+
    test_negative_mixed<Real, Num>(boost::mpl::bool_<std::numeric_limits<Num>::is_signed>());
 }
 
@@ -439,6 +490,28 @@ void test()
    ac = a + std::string("8");
    BOOST_TEST(ac == 16);
 #endif
+   //
+   // simple tests with immediate values, these calls can be optimised in many backends:
+   //
+   ac = a + b;
+   BOOST_TEST(ac == 72);
+   ac = a + +b;
+   BOOST_TEST(ac == 72);
+   ac = +a + b;
+   BOOST_TEST(ac == 72);
+   ac = +a + +b;
+   BOOST_TEST(ac == 72);
+   ac = a + -b;
+   BOOST_TEST(ac == 8 - 64);
+   ac = -a + b;
+   BOOST_TEST(ac == -8+64);
+   ac = -a + -b;
+   BOOST_TEST(ac == -72);
+   ac = a + - + -b; // lots of unary operators!!
+   BOOST_TEST(ac == 72);
+   ac = a;
+   ac = b / ac;
+   BOOST_TEST(ac == b / a);
    //
    // Comparisons:
    //