diff --git a/include/boost/math/algebraic_traits.hpp b/include/boost/math/algebraic_traits.hpp
new file mode 100644
index 000000000..9b3578317
--- /dev/null
+++ b/include/boost/math/algebraic_traits.hpp
@@ -0,0 +1,180 @@
+//  (C) Copyright Jeremy William Murphy 2016.
+
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_ALGEBRAIC_TRAITS
+#define BOOST_MATH_ALGEBRAIC_TRAITS
+
+#include <boost/config.hpp>
+#include <boost/core/enable_if.hpp>
+#include <boost/mpl/bool.hpp>
+#include <boost/mpl/logical.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/type_traits/is_floating_point.hpp>
+
+#include <limits>
+#include <functional>
+#include <string>
+
+/*                                      PURPOSE
+ * 
+ * A facility to provide convenient access to algebraic traits of types.
+ * 
+ * For example:
+ * 
+ *      algebraic_traits<std::multiplies, int>::zero_element() -> 0
+ *      algebraic_traits<std::multiplies, int>::identity_element() -> 1
+ *      algebraic_traits<std::plus, std::string>::identity_element() -> ""
+ * 
+ *      algebraic_traits<std::multiplies, int>::is_associative() -> true
+ *      algebraic_traits<std::multiplies, double>::is_associative() -> false
+ *      algebraic_traits<std::plus, std::string>::is_associative() -> true
+ *      algebraic_traits<std::multiplies, double>::is_commutative() -> true
+ *      algebraic_traits<std::plus, std::string>::is_commutative() -> false
+ * 
+ * Or alternatively:
+ * 
+ *      std::multiplies<int> op1;
+ *      std::multiplies<double> op2;
+ *      std::plus<std::string> op3;
+ *      algebraic_traits::zero_element(op1) -> 0
+ *      algebraic_traits::zero_element(op2) -> 0.0
+ *      algebraic_traits::identity_element(op3) -> ""
+ * 
+ * 
+ * Properties to consider:
+ * 
+ * Medial
+ * Distributive
+ * Commutative
+ * Idempotent
+ * Associative
+ * Unital
+ * Cancellative
+ */
+
+namespace boost {
+namespace math {
+namespace detail {
+    
+struct algebraic_traits_base
+{
+    BOOST_STATIC_CONSTANT(bool, is_specialized = false);
+};
+
+} // detail
+
+
+template <template <typename> class BinaryOperation, typename T, typename Enable = void>
+struct algebraic_traits : public detail::algebraic_traits_base
+{
+    // To provide default implementations or not?
+};
+
+
+// associativity
+
+template <template <typename> class BinaryOperation, typename T, typename Enable = void>
+struct algebraic_traits_associativity
+{};
+
+
+template <template <typename> class BinaryOperation, typename T>
+struct algebraic_traits_associativity<BinaryOperation, T, 
+typename enable_if< typename mpl::or_<is_integral<T>, mpl::bool_<std::numeric_limits<T>::is_exact> >::type >::type >
+{
+    BOOST_STATIC_CONSTEXPR
+    bool is_associative() { return true; }
+};
+
+
+template <template <typename> class BinaryOperation, typename T>
+struct algebraic_traits_associativity<BinaryOperation, T, 
+typename enable_if< typename mpl::or_<is_floating_point<T>, mpl::bool_<!std::numeric_limits<T>::is_exact> >::type >::type >
+{
+    BOOST_STATIC_CONSTEXPR
+    bool is_associative() { return false; }
+};
+
+
+// commutativity
+
+
+// associativity
+
+template <template <typename> class BinaryOperation, typename T, typename Enable = void>
+struct algebraic_traits_commutativity
+{};
+
+
+template <typename T>
+struct algebraic_traits_commutativity<std::plus, T>
+{
+    BOOST_STATIC_CONSTEXPR
+    bool is_commutative() { return true; }
+};
+
+
+template <>
+struct algebraic_traits_commutativity<std::plus, std::string>
+{
+    BOOST_STATIC_CONSTEXPR
+    bool is_commutative() { return false; }
+};
+
+
+//  multiplication
+
+template <typename T>
+struct algebraic_traits< std::multiplies, T >
+{
+    BOOST_STATIC_CONSTANT(bool, is_specialized = true);
+    
+    BOOST_STATIC_CONSTEXPR
+    T zero_element() { return T(0); }
+    
+    BOOST_STATIC_CONSTEXPR
+    T identity_element() { return T(1); }
+};
+
+
+
+// What makes a type associative in programming: exactness?
+
+
+//  plus
+
+// The inefficient way to multiply inherit?
+template <typename T>
+struct algebraic_traits<std::plus, T> : public algebraic_traits_commutativity<std::plus, T>, public algebraic_traits_associativity<std::plus, T>
+{
+    BOOST_STATIC_CONSTANT(bool, is_specialized = true);
+
+    BOOST_STATIC_CONSTEXPR
+    T identity_element() { return T(0); }
+};
+
+
+template <>
+struct algebraic_traits<std::plus, std::string>
+{
+    BOOST_STATIC_CONSTANT(bool, is_specialized = true);
+
+    // TODO: Can be constexpr? String literal in C++14?
+    BOOST_STATIC_CONSTEXPR
+    std::string identity_element() { return std::string(); }
+    
+    BOOST_STATIC_CONSTEXPR
+    bool is_commutative() { return false; }
+
+    BOOST_STATIC_CONSTEXPR
+    bool is_associative() { return true; }
+};
+
+
+} // math
+} // boost
+
+#endif
diff --git a/include/boost/math/algorithm.hpp b/include/boost/math/algorithm.hpp
new file mode 100644
index 000000000..ee76dc16c
--- /dev/null
+++ b/include/boost/math/algorithm.hpp
@@ -0,0 +1,62 @@
+//  (C) Copyright Jeremy William Murphy 2016.
+
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_ALGORITHM
+#define BOOST_MATH_ALGORITHM
+
+#include <boost/math/algebraic_traits.hpp>
+#include <boost/concept/assert.hpp>
+#include <boost/concept_check.hpp>
+#include <boost/bind.hpp>
+
+#include <functional>
+#include <numeric>
+#include <utility>
+
+namespace boost {
+namespace math {
+    namespace detail
+    {
+        struct not_equal_to_zero
+        {
+            typedef bool result_type;
+            
+            // Slightly awkward.
+            template <template <typename> class Op, typename T>
+            bool operator()(Op<T> const &, T const &x) const
+            {
+                return x != algebraic_traits<Op, T>::zero_element();
+            }
+        };
+    }
+    
+    
+    template<typename I, typename T, typename BinaryOperation, typename Predicate>
+    inline std::pair<T, I>
+    accumulate_if(I first, I last, T init, BinaryOperation op, Predicate p)
+    {
+        BOOST_CONCEPT_ASSERT((InputIterator<I>));
+        
+        for (; first != last && p(init); ++first )
+            init = op ( init, *first );
+        return std::make_pair(init, first);
+    }
+    
+
+    template<typename I, typename T, typename BinaryOperation>
+    inline std::pair<T, I>
+    accumulate_if(I first, I last, T init, BinaryOperation op)
+    {
+        using boost::placeholders::_1;
+        
+        detail::not_equal_to_zero const pred;
+        return accumulate_if(first, last, init, op, bind(pred, op, _1));
+    }
+} // math
+} // boost
+
+#endif
+
diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp
index 1dda577b4..0680e7e21 100644
--- a/include/boost/math/common_factor_rt.hpp
+++ b/include/boost/math/common_factor_rt.hpp
@@ -11,11 +11,13 @@
 #include <boost/core/enable_if.hpp>
 #include <boost/mpl/and.hpp>
 #include <boost/type_traits.hpp>
-
+#include <boost/math/algorithm.hpp>
 #include <boost/config.hpp>  // for BOOST_NESTED_TEMPLATE, etc.
 #include <boost/limits.hpp>  // for std::numeric_limits
 #include <climits>           // for CHAR_MIN
+#include <iterator>
 #include <boost/detail/workaround.hpp>
+#include <boost/math/algebraic_traits.hpp>
 #include <iterator>
 #include <algorithm>
 #include <limits>
@@ -379,12 +381,38 @@ inline Integer gcd(Integer const &a, Integer const &b)
     return detail::optimal_gcd_select(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
 }
 
+
+// Proof-of-concept only.
+template <typename T>
+struct gcd_evaluator
+{
+    T operator()(T const &x, T const &y) const
+    {
+        return gcd(x, y);
+    }
+};
+
+
 template <typename Integer>
 inline Integer lcm(Integer const &a, Integer const &b)
 {
    return detail::lcm_imp(static_cast<Integer>(gcd_traits<Integer>::abs(a)), static_cast<Integer>(gcd_traits<Integer>::abs(b)));
 }
 
+template <typename T>
+struct algebraic_traits<gcd_evaluator, T>
+{
+    BOOST_STATIC_CONSTEXPR
+    T zero_element() { return T(1); }
+    
+    BOOST_STATIC_CONSTEXPR
+    bool is_commutative() { return true; }
+    
+    BOOST_STATIC_CONSTEXPR
+    bool is_associative() { return true; }
+};
+
+
 /**
  * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
  * Chapter 4.5.2, Algorithm C: Greatest common divisor of n integers.
@@ -404,14 +432,8 @@ gcd_range(I first, I last)
 {
     BOOST_ASSERT(first != last);
     typedef typename std::iterator_traits<I>::value_type T;
-    
-    T d = *first++;
-    while (d != T(1) && first != last)
-    {
-        d = gcd(d, *first);
-        first++;
-    }
-    return std::make_pair(d, first);
+    I const init(first++);
+    return accumulate_if(first, last, *init, gcd_evaluator<T>());
 }
 
 }  // namespace math