1 files changed, 237 insertions, 0 deletions
diff --git a/src/boost/libs/math/example/polynomial_arithmetic.cpp b/src/boost/libs/math/example/polynomial_arithmetic.cpp
new file mode 100644
index 00000000..78efd06d
--- /dev/null
+++ b/src/boost/libs/math/example/polynomial_arithmetic.cpp
@@ -0,0 +1,237 @@
+// 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)
+
+// Copyright Jeremy W. Murphy 2015.
+
+// This file is written to be included from a Quickbook .qbk document.
+// It can be compiled by the C++ compiler, and run. Any output can
+// also be added here as comment or included or pasted in elsewhere.
+// Caution: this file contains Quickbook markup as well as code
+// and comments: don't change any of the special comment markups!
+
+//[polynomial_arithmetic_0
+/*`First include the essential polynomial header (and others) to make the example:
+*/
+#include <boost/math/tools/polynomial.hpp>
+//] [polynomial_arithmetic_0
+
+#include <boost/array.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/assert.hpp>
+
+#include <iostream>
+#include <stdexcept>
+#include <cmath>
+#include <string>
+#include <utility>
+
+//[polynomial_arithmetic_1
+/*`and some using statements are convenient:
+*/
+
+using std::string;
+using std::exception;
+using std::cout;
+using std::abs;
+using std::pair;
+
+using namespace boost::math;
+using namespace boost::math::tools; // for polynomial
+using boost::lexical_cast;
+
+//] [/polynomial_arithmetic_1]
+
+template <typename T>
+string sign_str(T const &x)
+{
+  return x < 0 ? "-" : "+";
+}
+
+template <typename T>
+string inner_coefficient(T const &x)
+{
+  string result(" " + sign_str(x) + " ");
+  if (abs(x) != T(1))
+      result += lexical_cast<string>(abs(x));
+  return result;
+}
+
+/*! Output in formula format.
+For example: from a polynomial in Boost container storage  [ 10, -6, -4, 3 ]
+show as human-friendly formula notation: 3x^3 - 4x^2 - 6x + 10.
+*/
+template <typename T>
+string formula_format(polynomial<T> const &a)
+{
+  string result;
+  if (a.size() == 0)
+      result += lexical_cast<string>(T(0));
+  else
+  {
+    // First one is a special case as it may need unary negate.
+    unsigned i = a.size() - 1;
+    if (a[i] < 0)
+        result += "-";
+    if (abs(a[i]) != T(1))
+        result += lexical_cast<string>(abs(a[i]));
+
+    if (i > 0)
+    {
+      result += "x";
+      if (i > 1)
+      {
+          result += "^" + lexical_cast<string>(i);
+          i--;
+          for (; i != 1; i--)
+              if (a[i])
+                result += inner_coefficient(a[i]) + "x^" + lexical_cast<string>(i);
+
+          if (a[i])
+            result += inner_coefficient(a[i]) + "x";
+      }
+      i--;
+
+      if (a[i])
+        result += " " + sign_str(a[i]) + " " + lexical_cast<string>(abs(a[i]));
+    }
+  }
+  return result;
+} // string formula_format(polynomial<T> const &a)
+
+
+int main()
+{
+  cout << "Example: Polynomial arithmetic.\n\n";
+
+  try
+  {
+//[polynomial_arithmetic_2
+/*`Store the coefficients in a convenient way to access them,
+then create some polynomials using construction from an iterator range,
+and finally output in a 'pretty' formula format.
+
+[tip Although we might conventionally write a polynomial from left to right
+in descending order of degree, Boost.Math stores in [*ascending order of degree].]
+
+  Read/write for humans:    3x^3 - 4x^2 - 6x + 10
+  Boost polynomial storage: [ 10, -6, -4, 3 ]
+*/
+  boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
+  polynomial<double> const a(d3a.begin(), d3a.end());
+
+  // With C++11 and later, you can also use initializer_list construction.
+  polynomial<double> const b{{-2.0, 1.0}};
+
+  // formula_format() converts from Boost storage to human notation.
+  cout << "a = " << formula_format(a)
+  << "\nb = " << formula_format(b) << "\n\n";
+
+//] [/polynomial_arithmetic_2]
+
+//[polynomial_arithmetic_3
+  // Now we can do arithmetic with the usual infix operators: + - * / and %.
+  polynomial<double> s = a + b;
+  cout << "a + b = " << formula_format(s) << "\n";
+  polynomial<double> d = a - b;
+  cout << "a - b = " << formula_format(d) << "\n";
+  polynomial<double> p = a * b;
+  cout << "a * b = " << formula_format(p) << "\n";
+  polynomial<double> q = a / b;
+  cout << "a / b = " << formula_format(q) << "\n";
+  polynomial<double> r = a % b;
+  cout << "a % b = " << formula_format(r) << "\n";
+//] [/polynomial_arithmetic_3]
+
+//[polynomial_arithmetic_4
+/*`
+Division is a special case where you can calculate two for the price of one.
+
+Actually, quotient and remainder are always calculated together due to the nature
+of the algorithm: the infix operators return one result and throw the other
+away.
+
+If you are doing a lot of division and want both the quotient and remainder, then
+you don't want to do twice the work necessary.
+
+In that case you can call the underlying function, [^quotient_remainder],
+to get both results together as a pair.
+*/
+  pair< polynomial<double>, polynomial<double> > result;
+  result = quotient_remainder(a, b);
+// Reassure ourselves that the result is the same.
+  BOOST_ASSERT(result.first == q);
+  BOOST_ASSERT(result.second == r);
+//] [/polynomial_arithmetic_4]
+//[polynomial_arithmetic_5
+  /* 
+We can use the right and left shift operators to add and remove a factor of x.
+This has the same semantics as left and right shift for integers where it is a 
+factor of 2. x is the smallest prime factor of a polynomial as is 2 for integers.
+*/
+    cout << "Right and left shift operators.\n";
+    cout << "\n" << formula_format(p) << "\n";
+    cout << "... right shift by 1 ...\n";
+    p >>= 1;
+    cout << formula_format(p) << "\n";
+    cout << "... left shift by 2 ...\n";
+    p <<= 2;
+    cout << formula_format(p) << "\n";    
+  
+/*
+We can also give a meaning to odd and even for a polynomial that is consistent
+with these operations: a polynomial is odd if it has a non-zero constant value, 
+even otherwise. That is:
+    x^2 + 1     odd
+    x^2         even    
+   */
+    cout << std::boolalpha;
+    cout << "\nPrint whether a polynomial is odd.\n";
+    cout << formula_format(s) << "   odd? " << odd(s) << "\n";
+    // We cheekily use the internal details to subtract the constant, making it even.
+    s -= s.data().front();
+    cout << formula_format(s) << "   odd? " << odd(s) << "\n";
+    // And of course you can check if it is even:
+    cout << formula_format(s) << "   even? " << even(s) << "\n";
+    
+    
+    //] [/polynomial_arithmetic_5]
+    //[polynomial_arithmetic_6]
+    /* For performance and convenience, we can test whether a polynomial is zero 
+     * by implicitly converting to bool with the same semantics as int.    */
+    polynomial<double> zero; // Default construction is 0.
+    cout << "zero: " << (zero ? "not zero" : "zero") << "\n";
+    cout << "r: " << (r ? "not zero" : "zero") << "\n";
+    /* We can also set a polynomial to zero without needing a another zero 
+     * polynomial to assign to it. */
+    r.set_zero();
+    cout << "r: " << (r ? "not zero" : "zero") << "\n";    
+    //] [/polynomial_arithmetic_6]
+}
+catch (exception const &e)
+{
+  cout << "\nMessage from thrown exception was:\n   " << e.what() << "\n";
+}
+} // int main()
+
+/*
+//[polynomial_output_1
+
+a = 3x^3 - 4x^2 - 6x + 10
+b = x - 2
+
+//] [/polynomial_output_1]
+
+
+//[polynomial_output_2
+
+a + b = 3x^3 - 4x^2 - 5x + 8
+a - b = 3x^3 - 4x^2 - 7x + 12
+a * b = 3x^4 - 10x^3 + 2x^2 + 22x - 20
+a / b = 3x^2 + 2x - 2
+a % b = 6
+
+//] [/polynomial_output_2]
+
+*/