Subject: [Boost-commit] svn:boost r73068 - sandbox/e_float/libs/e_float/example
From: e_float_at_[hidden]
Date: 2011-07-13 16:28:02
Author: christopher_kormanyos
Date: 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
New Revision: 73068
URL: http://svn.boost.org/trac/boost/changeset/73068
Log:
- Re-added example source files.
Added:
sandbox/e_float/libs/e_float/example/ (props changed)
sandbox/e_float/libs/e_float/example/example_001_basic_usage_real.cpp (contents, props changed)
sandbox/e_float/libs/e_float/example/example_002_basic_usage_imag.cpp (contents, props changed)
sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp (contents, props changed)
sandbox/e_float/libs/e_float/example/example_008_guass_laguerre.cpp (contents, props changed)
sandbox/e_float/libs/e_float/example/examples.h (contents, props changed)
Added: sandbox/e_float/libs/e_float/example/example_001_basic_usage_real.cpp
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/example/example_001_basic_usage_real.cpp 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
@@ -0,0 +1,50 @@
+
+// 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)} © ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#include <vector>
+#include <iostream>
+#include <algorithm>
+#include <iterator>
+
+#include <e_float/e_float.hpp>
+#include <e_float/e_float_functions.hpp>
+#include "examples.h"
+#include "../src/utility/util_timer.h"
+
+void examples::nr_001::basic_usage_real(void)
+{
+ // Print 21 values of the function gamma[(222/10) + k]
+ // for 0 <= k < 21 to the standard output using e_float. Also compute the computation
+ // time for the calculation using e_float in [ms].
+ // A comparable Mathematica code is:
+ // Timing[Table[N[Gamma[(222/10) + k], 100],{k, 0, 20, 1}]].
+
+ const e_float v(222 / ef::ten());
+
+ std::vector<e_float> values(static_cast<std::size_t>(21u));
+
+ const Util::timer tm;
+ for(INT32 k = static_cast<INT32>(0); k < static_cast<INT32>(values.size()); k++)
+ {
+ values[static_cast<std::size_t>(k)] = ef::gamma(v + k);
+ }
+ const double elapsed = tm.elapsed();
+
+ std::cout << "Elapsed time: " << elapsed << "\n";
+
+ const std::streamsize original_prec = std::cout.precision(std::numeric_limits<e_float>::digits10);
+ const std::ios::fmtflags original_flag = std::cout.setf(std::ios::showpos | std::ios::scientific);
+
+ std::copy(values.begin(), values.end(), std::ostream_iterator<e_float>(std::cout, "\n"));
+
+ std::cout.precision(original_prec);
+ std::cout.unsetf(std::ios::showpos | std::ios::scientific);
+ std::cout.setf(original_flag);
+}
Added: sandbox/e_float/libs/e_float/example/example_002_basic_usage_imag.cpp
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/example/example_002_basic_usage_imag.cpp 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
@@ -0,0 +1,53 @@
+
+// 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)} © ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#include <vector>
+#include <fstream>
+#include <algorithm>
+#include <iterator>
+
+#include <e_float/e_float.hpp>
+#include <e_float/e_float_functions.hpp>
+
+#include "examples.h"
+#include "../src/utility/util_timer.h"
+
+void examples::nr_002::basic_usage_imag(void)
+{
+ // Print 21 values of the function riemann_zeta[(1/2) + (((1234/10) + k) I)]
+ // for 0 <= k < 21 to the standard output using e_float. Also compute the
+ // computation time for the calculation using e_float in [ms].
+ // A comparable Mathematica code is:
+ // Timing[Table[N[Zeta[(1/2) + (((1234/10) + k) I)], 100],{k, 0, 20, 1}]].
+
+ static const e_float y(1234 / ef::ten());
+
+ std::vector<ef_complex> values(static_cast<std::size_t>(21u));
+
+ const Util::timer tm;
+ for(INT32 k = static_cast<INT32>(0); k < static_cast<INT32>(values.size()); k++)
+ {
+ const ef_complex z(ef::half(), y + k);
+
+ values[static_cast<std::size_t>(k)] = efz::riemann_zeta(z);
+ }
+ const double elapsed = tm.elapsed();
+
+ std::cout << "Elapsed time: " << elapsed << "\n";
+
+ const std::streamsize original_prec = std::cout.precision(std::numeric_limits<e_float>::digits10);
+ const std::ios::fmtflags original_flag = std::cout.setf(std::ios::showpos | std::ios::scientific);
+
+ std::copy(values.begin(), values.end(), std::ostream_iterator<ef_complex>(std::cout, "\n"));
+
+ std::cout.precision(original_prec);
+ std::cout.unsetf(std::ios::showpos | std::ios::scientific);
+ std::cout.setf(original_flag);
+}
Added: sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/example/example_005_recursive_trapezoid_integral.cpp 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
@@ -0,0 +1,65 @@
+
+// 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)} © ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#include <e_float/e_float.hpp>
+#include <e_float/e_float_functions.hpp>
+
+#include "examples.h"
+#include "../src/utility/util_lexical_cast.h"
+#include "../src/utility/util_trapezoid.h"
+
+namespace examples
+{
+ namespace nr_005
+ {
+ class RecursiveTrapezoidJ0 : public Util::RecursiveTrapezoidRule<e_float>
+ {
+ private:
+
+ static const e_float& my_tol(void)
+ {
+ static const e_float val("1E-" + Util::lexical_cast(std::numeric_limits<e_float>::digits10 / 2));
+ return val;
+ }
+
+ private:
+
+ const e_float my_z;
+
+ public:
+
+ RecursiveTrapezoidJ0(const e_float& z) : Util::RecursiveTrapezoidRule<e_float>(ef::zero(), ef::pi(), my_tol()),
+ my_z(z) { }
+
+ virtual ~RecursiveTrapezoidJ0() { }
+
+ private:
+
+ virtual e_float my_function(const e_float& x) const
+ {
+ return ef::cos(my_z * ef::sin(x));
+ }
+ };
+ }
+}
+
+e_float examples::nr_005::recursive_trapezoid_j0(const e_float& x)
+{
+ const RecursiveTrapezoidJ0 rtj0(x);
+
+ return rtj0.operation() / ef::pi();
+}
+
+e_float examples::nr_005::recursive_trapezoid_j0_test(void)
+{
+ static const e_float x = 12 + ef::euler_gamma();
+
+ return recursive_trapezoid_j0(x);
+}
Added: sandbox/e_float/libs/e_float/example/example_008_guass_laguerre.cpp
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/example/example_008_guass_laguerre.cpp 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
@@ -0,0 +1,313 @@
+
+// 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)} © ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#include <vector>
+
+#include <e_float/e_float.hpp>
+#include <e_float/e_float_functions.hpp>
+
+#include "examples.h"
+#include "../src/utility/util_digit_scale.h"
+#include "../src/utility/util_find_root_newton_raphson.h"
+#include "../src/utility/util_lexical_cast.h"
+#include "../src/utility/util_noncopyable.h"
+#include "../src/utility/util_point.h"
+
+namespace examples
+{
+ namespace nr_008
+ {
+ class GaussQuadAbscissasAndWeights : private Util::noncopyable
+ {
+ protected:
+
+ const INT32 N;
+ std::vector<e_float> xi;
+ std::vector<e_float> wi;
+
+ protected:
+
+ GaussQuadAbscissasAndWeights(const INT32 n) : N (n),
+ xi(0u),
+ wi(0u) { }
+
+ public:
+
+ virtual ~GaussQuadAbscissasAndWeights() { }
+
+ public:
+
+ const std::vector<e_float>& abscissas(void) const { return xi; }
+ const std::vector<e_float>& weights (void) const { return wi; }
+ };
+
+ class GuassLaguerreAbscissasAndWeights : public GaussQuadAbscissasAndWeights
+ {
+ private:
+
+ const e_float alpha;
+
+ public:
+
+ GuassLaguerreAbscissasAndWeights(const INT32 n, const e_float& a = ef::zero());
+
+ virtual ~GuassLaguerreAbscissasAndWeights() { }
+
+ private:
+
+ class FindRootLaguerre : public Util::FindRootNewtonRaphson<e_float>
+ {
+ private:
+
+ static const e_float& my_tol(void)
+ {
+ static const e_float val = e_float("1E-" + Util::lexical_cast(ef::tol() - static_cast<INT64>(3)));
+ return val;
+ }
+
+ const INT32 N;
+ const e_float alpha;
+ mutable e_float my_d;
+ mutable e_float my_p2;
+
+ public:
+
+ FindRootLaguerre(const INT32 n,
+ const e_float& a,
+ const e_float& lo,
+ const e_float& hi) : Util::FindRootNewtonRaphson<e_float>(lo, hi, my_tol()),
+ N (n),
+ alpha(a),
+ my_d (static_cast<UINT32>(0u)),
+ my_p2(static_cast<UINT32>(0u)) { }
+
+ virtual ~FindRootLaguerre() { }
+
+ const e_float& poly_p2(void) const { return my_p2; }
+
+ private:
+
+ virtual void my_function_derivative(const e_float& x, e_float& f, e_float& d) const;
+
+ virtual e_float my_function(const e_float& x) const;
+ };
+
+ e_float my_laguerre(const e_float& x) const
+ {
+ static const FindRootLaguerre rl(N, alpha, ef::zero(), ef::million());
+ return rl.function(x);
+ }
+ };
+ }
+}
+
+examples::nr_008::GuassLaguerreAbscissasAndWeights::GuassLaguerreAbscissasAndWeights(const INT32 n,
+ const e_float& a) : GaussQuadAbscissasAndWeights(n),
+ alpha(a)
+{
+ // Walk through the function using a step-size of dynamic width
+ // to find the zero crossings of the Laguerre function.
+ // Store these zero crossings as bracketed root estimates.
+ std::vector<Util::point<e_float> > xlo_xhi_points;
+
+ // Estimate first Laguerre zero using the approximation of x1
+ // from Stroud and Secrest.
+ const double alf = ef::to_double(alpha);
+ const double x1 = ((1.0 + alf) * (3.0 + (0.92 * alf))) / (1.0 + (2.4 * static_cast<double>(N)) + (1.8 * alf));
+
+ e_float delta(x1 / 10.0);
+ e_float x = delta;
+
+ bool bo_is_neg = ef::isneg(my_laguerre(ef::zero()));
+
+ std::cout << "finding approximate roots..." << std::endl;
+
+ while(static_cast<INT32>(xlo_xhi_points.size()) < N)
+ {
+ x += delta;
+
+ if(ef::isneg(my_laguerre(x)) != bo_is_neg)
+ {
+ // Refine the approximate root with 16 bisection iteration steps.
+ e_float dx;
+ e_float rt;
+ e_float hi = x;
+ e_float lo = x - delta;
+
+ if(my_laguerre(lo) < ef::zero())
+ {
+ dx = hi - lo;
+ rt = lo;
+ }
+ else
+ {
+ dx = lo - hi;
+ rt = hi;
+ }
+
+ // Here are the 16 bisection iteration steps.
+ for(UINT32 j = static_cast<UINT32>(0u); j < static_cast<UINT32>(16u); j++)
+ {
+ dx /= static_cast<INT32>(2);
+
+ const e_float x_mid = rt + dx;
+ const e_float f_mid = my_laguerre(x_mid);
+
+ if(f_mid <= ef::zero())
+ {
+ rt = x_mid;
+ }
+ }
+
+ // Store the refined approximate root as a bracketed point.
+ xlo_xhi_points.push_back(Util::point<e_float>(rt - dx, rt + dx));
+
+ if(xlo_xhi_points.size() > static_cast<std::size_t>(1u))
+ {
+ const e_float new_delta2 = (xlo_xhi_points.end() - static_cast<std::size_t>(1u))->x
+ - (xlo_xhi_points.end() - static_cast<std::size_t>(2u))->x;
+
+ // Dynamically adjust and set the new step-size.
+ delta = new_delta2 / static_cast<INT32>(2);
+ }
+
+ bo_is_neg = !bo_is_neg;
+ }
+ }
+
+ // Calculate the abscissas and weights to full precision.
+ for(std::size_t i = static_cast<std::size_t>(0u); i < xlo_xhi_points.size(); i++)
+ {
+ std::cout << "calculating abscissa and weight for index: " << i << std::endl;
+
+ // Find the abscissas using Newton-Raphson iteration.
+ const FindRootLaguerre rooti(N, alpha, xlo_xhi_points[i].x, xlo_xhi_points[i].y);
+ const e_float ri = rooti.operation();
+
+ if(!rooti.success())
+ {
+ break;
+ }
+
+ // Calculate the weights.
+ e_float f;
+ e_float d;
+ rooti.function_derivative(ri, f, d);
+
+ // Store abscissas and weights.
+ xi.push_back(ri);
+
+ if(ef::iszero(alpha))
+ {
+ wi.push_back(-ef::one() / ((d * N) * rooti.poly_p2()));
+ }
+ else
+ {
+ const e_float norm_g = ef::gamma(alpha + N) / ef::factorial(static_cast<UINT32>(N - static_cast<INT32>(1)));
+ wi.push_back(-norm_g / ((d * N) * rooti.poly_p2()));
+ }
+ }
+}
+
+void examples::nr_008::GuassLaguerreAbscissasAndWeights::FindRootLaguerre::my_function_derivative(const e_float& x, e_float& f, e_float& d) const
+{
+ f = my_function(x);
+ d = my_d;
+}
+
+e_float examples::nr_008::GuassLaguerreAbscissasAndWeights::FindRootLaguerre::my_function(const e_float& x) const
+{
+ e_float p1(1);
+ e_float p2(0);
+
+ for(INT32 j = static_cast<INT32>(0); j < N; j++)
+ {
+ const e_float p3(p2);
+ p2 = p1;
+
+ const INT32 j_plus_one = static_cast<INT32>(j + static_cast<INT32>(1));
+ const INT32 two_j_plus_one = static_cast<INT32>(j + j_plus_one);
+
+ p1 = (((two_j_plus_one + (alpha - x)) * p2) - ((j + alpha) * p3)) / j_plus_one;
+ }
+
+ my_d = ((N * p1) - ((N + alpha) * p2)) / x;
+ my_p2 = p2;
+
+ return p1;
+}
+
+namespace examples
+{
+ namespace nr_008
+ {
+ struct GaussLaguerreAi : public Util::Function<e_float>
+ {
+ private:
+
+ const e_float x;
+ const e_float zeta;
+ const e_float one_over_zeta;
+ e_float factor;
+
+ public:
+
+ GaussLaguerreAi(const e_float& X) : x(X),
+ zeta(ef::two_third() * (ef::sqrt(x) * x)),
+ one_over_zeta(ef::one() / zeta),
+ factor(0u)
+ {
+ const e_float zeta_times_48_pow_sixth = ef::rootn(zeta * static_cast<INT32>(48), static_cast<INT32>(6));
+ factor = ef::one() / ((ef::sqrt_pi() * zeta_times_48_pow_sixth) * (ef::exp(zeta) * ef::gamma(ef::five() / static_cast<INT32>(6))));
+ }
+
+ virtual ~GaussLaguerreAi() { }
+
+ private:
+
+ virtual e_float my_function(const e_float& t) const
+ {
+ return factor / ef::rootn(ef::two() + (t * one_over_zeta), static_cast<INT32>(6));
+ }
+ };
+ }
+}
+
+e_float examples::nr_008::gauss_laguerre_airy_a(const e_float& x)
+{
+ static const INT32 n_terms = static_cast<INT32>(Util::DigitScale() * 400.0);
+
+ static const GuassLaguerreAbscissasAndWeights aw(n_terms, -ef::one() / static_cast<INT32>(6));
+
+ static const std::size_t sz = aw.weights().size();
+
+ const GaussLaguerreAi ai(x);
+
+ e_float sum(0);
+
+ for(std::size_t i = static_cast<std::size_t>(0u); i < sz; i++)
+ {
+ const e_float fi = ai.function(aw.abscissas()[i]);
+ const e_float wi = aw.weights()[i];
+
+ const e_float term = wi * fi;
+
+ const INT64 order_check = static_cast<INT64>(term.order() - sum.order());
+
+ if((i > static_cast<INT32>(20)) && (order_check < -ef::tol()))
+ {
+ break;
+ }
+
+ sum += term;
+ }
+
+ return sum;
+}
Added: sandbox/e_float/libs/e_float/example/examples.h
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/example/examples.h 2011-07-13 16:28:02 EDT (Wed, 13 Jul 2011)
@@ -0,0 +1,41 @@
+
+// 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)} © ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
+#ifndef _EXAMPLES_2010_01_02_H_
+ #define _EXAMPLES_2010_01_02_H_
+
+ #include <e_float/e_float.hpp>
+ #include <e_float/e_float_complex.hpp>
+
+ namespace examples
+ {
+ namespace nr_001
+ {
+ void basic_usage_real(void);
+ }
+
+ namespace nr_002
+ {
+ void basic_usage_imag(void);
+ }
+
+ namespace nr_005
+ {
+ e_float recursive_trapezoid_j0 (const e_float& x);
+ e_float recursive_trapezoid_j0_test(void);
+ }
+
+ namespace nr_008
+ {
+ e_float gauss_laguerre_airy_a(const e_float& x);
+ }
+ }
+
+#endif // _EXAMPLES_2010_01_02_H_