[Boost-commit] svn:boost r75121 - in sandbox/multiprecision/libs/multiprecision: build example

Subject: [Boost-commit] svn:boost r75121 - in sandbox/multiprecision/libs/multiprecision: build example
From: e_float_at_[hidden]
Date: 2011-10-26 13:25:29

Author: christopher_kormanyos
Date: 2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
New Revision: 75121
URL: http://svn.boost.org/trac/boost/changeset/75121

Log:
- Added an example of generic mathematical programming using orthogonal polynomials.
Added:
   sandbox/multiprecision/libs/multiprecision/example/example_010_orthogonal_polynomials.cpp   (contents, props changed)
Text files modified: 
   sandbox/multiprecision/libs/multiprecision/build/MakefileFiles.gmk              |     4 +++-                                    
   sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj         |     1 +                                       
   sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj.filters |     3 +++                                     
   sandbox/multiprecision/libs/multiprecision/example/examples.h                   |    28 ++++++++++++++++++++++++++++            
   4 files changed, 35 insertions(+), 1 deletions(-)

Modified: sandbox/multiprecision/libs/multiprecision/build/MakefileFiles.gmk
==============================================================================
--- sandbox/multiprecision/libs/multiprecision/build/MakefileFiles.gmk	(original)
+++ sandbox/multiprecision/libs/multiprecision/build/MakefileFiles.gmk	2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
@@ -50,7 +50,9 @@
 FILES_EXAMPLES  = ../example/example_001_basic_usage_real                       \
                   ../example/example_002_basic_usage_imag                       \
                   ../example/example_005_recursive_trapezoid_integral           \
-                  ../example/example_008_gauss_laguerre
+                  ../example/example_008_gauss_laguerre                         \
+                  ../example/example_010_orthogonal_polynomials
+
 
 #
 # Test files

Modified: sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj
==============================================================================
--- sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj	(original)
+++ sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj	2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
@@ -771,6 +771,7 @@
     <ClCompile Include="..\example\example_002_basic_usage_imag.cpp" />
     <ClCompile Include="..\example\example_005_recursive_trapezoid_integral.cpp" />
     <ClCompile Include="..\example\example_008_gauss_laguerre.cpp" />
+    <ClCompile Include="..\example\example_010_orthogonal_polynomials.cpp" />
     <ClCompile Include="..\src\backends\float\efx\mp_float_efx.cpp">
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='debug-gmp|Win32'">true</ExcludedFromBuild>
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='debug-mpfr|Win32'">true</ExcludedFromBuild>

Modified: sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj.filters
==============================================================================
--- sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj.filters	(original)
+++ sandbox/multiprecision/libs/multiprecision/build/multiprecision.vcxproj.filters	2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
@@ -388,6 +388,9 @@
     <ClCompile Include="..\src\utility\util_timer.cpp">
       <Filter>libs\multiprecision\src\utility</Filter>
     </ClCompile>
+    <ClCompile Include="..\example\example_010_orthogonal_polynomials.cpp">
+      <Filter>libs\multiprecision\example</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\src\functions\tables\tables.h">

Added: sandbox/multiprecision/libs/multiprecision/example/example_010_orthogonal_polynomials.cpp
==============================================================================
--- (empty file)
+++ sandbox/multiprecision/libs/multiprecision/example/example_010_orthogonal_polynomials.cpp	2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
@@ -0,0 +1,133 @@
+
+//          Copyright Christopher Kormanyos 2002 - 2011.
+// 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)
+
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#include <boost/multiprecision/mp_float.hpp>
+#include <boost/multiprecision/mp_float_functions.hpp>
+
+#include "examples.h"
+
+namespace
+{
+  template<typename T, const int N>
+  struct integer_value_template
+  {
+    static const T& value(void)
+    {
+      static const T t(N);
+      return t;
+    }
+  };
+
+  inline std::complex<double> operator+(const boost::int32_t n, const std::complex<double>& x) { return std::complex<double>(x.real() + n, x.imag()); }
+  inline std::complex<double> operator-(const boost::int32_t n, const std::complex<double>& x) { return std::complex<double>(n - x.real(), x.imag()); }
+  inline std::complex<double> operator*(const boost::int32_t n, const std::complex<double>& x) { return std::complex<double>(x.real() * n, x.imag() * n); }
+  inline std::complex<double> operator/(const boost::int32_t n, const std::complex<double>& x) { return std::complex<double>(x.real() * n, -x.imag() * n) / ((x.real() * x.real()) + (x.imag() * x.imag())); }
+
+  inline std::complex<double> operator+(const std::complex<double>& x, const boost::int32_t n) { return std::complex<double>(x.real() + n, x.imag()); }
+  inline std::complex<double> operator-(const std::complex<double>& x, const boost::int32_t n) { return std::complex<double>(x.real() - n, x.imag()); }
+  inline std::complex<double> operator*(const std::complex<double>& x, const boost::int32_t n) { return std::complex<double>(x.real() * n, x.imag() * n); }
+  inline std::complex<double> operator/(const std::complex<double>& x, const boost::int32_t n) { return std::complex<double>(x.real() / n, x.imag() / n); }
+
+  typedef enum enum_orthogonal_polynomial_type
+  {
+    orthogonal_polynomial_type_chebyshev_t = 1,
+    orthogonal_polynomial_type_chebyshev_u = 2,
+    orthogonal_polynomial_type_laguerre_l  = 3,
+    orthogonal_polynomial_type_hermite_h   = 4
+  }
+  orthogonal_polynomial_type;
+
+  template<typename T>
+  static inline void orthogonal_polynomial_template(const T& x,
+                                                    const boost::uint32_t n,
+                                                    const orthogonal_polynomial_type polynomial_type,
+                                                    std::vector<T>& results)
+  {
+    // Compute the value of an orthogonal polynomial of one of the following types:
+    // Chebyshev 1st, Chebyshev 2nd, Laguerre, or Hermite
+
+    results.clear();
+    results.reserve(static_cast<std::size_t>(n + 1u));
+
+    T y0 = integer_value_template<T, 1>::value();
+
+    results.push_back(y0);
+
+    if(n == static_cast<boost::uint32_t>(0u)) { return; }
+
+    T y1;
+
+    if(polynomial_type == orthogonal_polynomial_type_chebyshev_t)
+    {
+      y1 = x;
+    }
+    else if(polynomial_type == orthogonal_polynomial_type_laguerre_l)
+    {
+      y1 = static_cast<boost::int32_t>(1) - x;
+    }
+    else
+    {
+      y1 = x * static_cast<boost::int32_t>(2);
+    }
+
+    results.push_back(y1);
+
+    if(n == static_cast<boost::uint32_t>(1u)) { return; }
+
+    T a = integer_value_template<T, 2>::value();
+    T b = integer_value_template<T, 0>::value();
+    T c = integer_value_template<T, 1>::value();
+
+    T yk;
+
+    // Calculate higher orders using the recurrence relation.
+    // The direction of stability is upward recurrence.
+    for(boost::int32_t k = static_cast<boost::int32_t>(2); k <= static_cast<boost::int32_t>(n); k++)
+    {
+      if(polynomial_type == orthogonal_polynomial_type_laguerre_l)
+      {
+        a = integer_value_template<T, -1>::value() / k;
+        b = integer_value_template<T, 2>::value() + a;
+        c = integer_value_template<T, 1>::value() + a;
+      }
+      else if(polynomial_type == orthogonal_polynomial_type_hermite_h)
+      {
+        c = integer_value_template<T, 2>::value() * static_cast<boost::int32_t>(k - static_cast<boost::int32_t>(1));
+      }
+
+      yk = (((a * x) + b) * y1) - (c * y0);
+
+      y0 = y1;
+      y1 = yk;
+
+      results.push_back(yk);
+    }
+  }
+}
+
+void examples::nr_010::chebyshev_t(const boost::uint32_t n, const double& x, std::vector<double>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::chebyshev_u(const boost::uint32_t n, const double& x, std::vector<double>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::hermite    (const boost::uint32_t n, const double& x, std::vector<double>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_laguerre_l,  results); }
+void examples::nr_010::laguerre   (const boost::uint32_t n, const double& x, std::vector<double>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_hermite_h,   results); }
+
+void examples::nr_010::chebyshev_t(const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::chebyshev_u(const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::hermite    (const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_laguerre_l,  results); }
+void examples::nr_010::laguerre   (const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_hermite_h,   results); }
+
+void examples::nr_010::chebyshev_t(const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::chebyshev_u(const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::hermite    (const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_laguerre_l,  results); }
+void examples::nr_010::laguerre   (const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_hermite_h,   results); }
+
+void examples::nr_010::chebyshev_t(const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::chebyshev_u(const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_chebyshev_t, results); }
+void examples::nr_010::hermite    (const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_laguerre_l,  results); }
+void examples::nr_010::laguerre   (const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results) { ::orthogonal_polynomial_template(x, n, ::orthogonal_polynomial_type_hermite_h,   results); }

Modified: sandbox/multiprecision/libs/multiprecision/example/examples.h
==============================================================================
--- sandbox/multiprecision/libs/multiprecision/example/examples.h	(original)
+++ sandbox/multiprecision/libs/multiprecision/example/examples.h	2011-10-26 13:25:27 EDT (Wed, 26 Oct 2011)
@@ -11,6 +11,11 @@
 #ifndef _EXAMPLES_2010_01_02_H_
   #define _EXAMPLES_2010_01_02_H_
 
+  #include <vector>
+  #include <complex>
+
+  #include <boost/cstdint.hpp>
+  #include <boost/multiprecision/mp_float.hpp>
   #include <boost/multiprecision/mp_float_functions.hpp>
 
   namespace examples
@@ -35,6 +40,29 @@
     {
       boost::multiprecision::mp_float gauss_laguerre_airy_a(const boost::multiprecision::mp_float& x);
     }
+
+    namespace nr_010
+    {
+      void chebyshev_t(const boost::uint32_t n, const double& x, std::vector<double>& results);
+      void chebyshev_u(const boost::uint32_t n, const double& x, std::vector<double>& results);
+      void hermite    (const boost::uint32_t n, const double& x, std::vector<double>& results);
+      void laguerre   (const boost::uint32_t n, const double& x, std::vector<double>& results);
+
+      void chebyshev_t(const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results);
+      void chebyshev_u(const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results);
+      void hermite    (const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results);
+      void laguerre   (const boost::uint32_t n, const std::complex<double>& x, std::vector<std::complex<double> >& results);
+
+      void chebyshev_t(const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results);
+      void chebyshev_u(const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results);
+      void hermite    (const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results);
+      void laguerre   (const boost::uint32_t n, const boost::multiprecision::mp_float& x, std::vector<boost::multiprecision::mp_float>& results);
+
+      void chebyshev_t(const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results);
+      void chebyshev_u(const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results);
+      void hermite    (const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results);
+      void laguerre   (const boost::uint32_t n, const boost::multiprecision::mp_complex& x, std::vector<boost::multiprecision::mp_complex>& results);
+    }
   }
 
 #endif // _EXAMPLES_2010_01_02_H_