1 files changed, 150 insertions, 0 deletions
diff --git a/src/boost/libs/math/example/continued_fractions.cpp b/src/boost/libs/math/example/continued_fractions.cpp
new file mode 100644
index 00000000..23f479fc
--- /dev/null
+++ b/src/boost/libs/math/example/continued_fractions.cpp
@@ -0,0 +1,150 @@
+//  (C) Copyright John Maddock 2018.
+//  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)
+
+#include <boost/math/tools/fraction.hpp>
+#include <iostream>
+#include <complex>
+#include <boost/multiprecision/cpp_complex.hpp>
+
+//[golden_ratio_1
+template <class T>
+struct golden_ratio_fraction
+{
+   typedef T result_type;
+
+   result_type operator()()
+   {
+      return 1;
+   }
+};
+//]
+
+//[cf_tan_fraction
+template <class T>
+struct tan_fraction
+{
+private:
+   T a, b;
+public:
+   tan_fraction(T v)
+      : a(-v * v), b(-1)
+   {}
+
+   typedef std::pair<T, T> result_type;
+
+   std::pair<T, T> operator()()
+   {
+      b += 2;
+      return std::make_pair(a, b);
+   }
+};
+//]
+//[cf_tan
+template <class T>
+T tan(T a)
+{
+   tan_fraction<T> fract(a);
+   return a / continued_fraction_b(fract, std::numeric_limits<T>::epsilon());
+}
+//]
+//[cf_expint_fraction
+template <class T>
+struct expint_fraction
+{
+   typedef std::pair<T, T> result_type;
+   expint_fraction(unsigned n_, T z_) : b(z_ + T(n_)), i(-1), n(n_) {}
+   std::pair<T, T> operator()()
+   {
+      std::pair<T, T> result = std::make_pair(-static_cast<T>((i + 1) * (n + i)), b);
+      b += 2;
+      ++i;
+      return result;
+   }
+private:
+   T b;
+   int i;
+   unsigned n;
+};
+//]
+//[cf_expint
+template <class T>
+inline std::complex<T> expint_as_fraction(unsigned n, std::complex<T> const& z)
+{
+   boost::uintmax_t max_iter = 1000;
+   expint_fraction<std::complex<T> > f(n, z);
+   std::complex<T> result = boost::math::tools::continued_fraction_b(
+      f,
+      std::complex<T>(std::numeric_limits<T>::epsilon()),
+      max_iter);
+   result = exp(-z) / result;
+   return result;
+}
+//]
+//[cf_upper_gamma_fraction
+template <class T>
+struct upper_incomplete_gamma_fract
+{
+private:
+   typedef typename T::value_type scalar_type;
+   T z, a;
+   int k;
+public:
+   typedef std::pair<T, T> result_type;
+
+   upper_incomplete_gamma_fract(T a1, T z1)
+      : z(z1 - a1 + scalar_type(1)), a(a1), k(0)
+   {
+   }
+
+   result_type operator()()
+   {
+      ++k;
+      z += scalar_type(2);
+      return result_type(scalar_type(k) * (a - scalar_type(k)), z);
+   }
+};
+//]
+//[cf_gamma_Q
+template <class T>
+inline std::complex<T> gamma_Q_as_fraction(const std::complex<T>& a, const std::complex<T>& z)
+{
+   upper_incomplete_gamma_fract<std::complex<T> > f(a, z);
+   std::complex<T> eps(std::numeric_limits<T>::epsilon());
+   return pow(z, a) / (exp(z) *(z - a + T(1) + boost::math::tools::continued_fraction_a(f, eps)));
+}
+//]
+inline boost::multiprecision::cpp_complex_50 gamma_Q_as_fraction(const boost::multiprecision::cpp_complex_50& a, const boost::multiprecision::cpp_complex_50& z)
+{
+   upper_incomplete_gamma_fract<boost::multiprecision::cpp_complex_50> f(a, z);
+   boost::multiprecision::cpp_complex_50 eps(std::numeric_limits<boost::multiprecision::cpp_complex_50::value_type>::epsilon());
+   return pow(z, a) / (exp(z) * (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)));
+}
+
+
+int main()
+{
+   using namespace boost::math::tools;
+
+   //[cf_gr
+   golden_ratio_fraction<double> func;
+   double gr = continued_fraction_a(
+      func,
+      std::numeric_limits<double>::epsilon());
+   std::cout << "The golden ratio is: " << gr << std::endl;
+   //]
+
+   std::cout << tan(0.5) << std::endl;
+
+   std::complex<double> arg(3, 2);
+   std::cout << expint_as_fraction(5, arg) << std::endl;
+
+   std::complex<double> a(3, 3), z(3, 2);
+   std::cout << gamma_Q_as_fraction(a, z) << std::endl;
+
+   boost::multiprecision::cpp_complex_50 am(3, 3), zm(3, 2);
+   std::cout << gamma_Q_as_fraction(am, zm) << std::endl;
+
+   return 0;
+}