|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r78043 - in sandbox/math: boost/math/special_functions libs/math/test
From: john_at_[hidden]
Date: 2012-04-17 12:15:36
Author: johnmaddock
Date: 2012-04-17 12:15:34 EDT (Tue, 17 Apr 2012)
New Revision: 78043
URL: http://svn.boost.org/trac/boost/changeset/78043
Log:
General update for multiprecision support:
Added new accelerated version of T1.
Text files modified:
sandbox/math/boost/math/special_functions/owens_t.hpp | 221 +++++++++++++++++++++++++++------------
sandbox/math/libs/math/test/owens_t.ipp | 6
sandbox/math/libs/math/test/owens_t_large_data.ipp | 4
sandbox/math/libs/math/test/test_owens_t.cpp | 15 ++
4 files changed, 171 insertions(+), 75 deletions(-)
Modified: sandbox/math/boost/math/special_functions/owens_t.hpp
==============================================================================
--- sandbox/math/boost/math/special_functions/owens_t.hpp (original)
+++ sandbox/math/boost/math/special_functions/owens_t.hpp 2012-04-17 12:15:34 EDT (Tue, 17 Apr 2012)
@@ -22,6 +22,7 @@
#include <boost/throw_exception.hpp>
#include <boost/assert.hpp>
#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/big_constant.hpp>
#include <stdexcept>
@@ -274,37 +275,37 @@
static const RealType c2[] =
{
- 0.99999999999999999999999729978162447266851932041876728736094298092917625009873l,
- -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968l,
- 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536l,
- -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685l,
- 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147l,
- -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977l,
- 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267l,
- -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274l,
- 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048l,
- -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467l,
- 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967l,
- -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501l,
- 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716l,
- -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469l,
- 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483l,
- -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854l,
- 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567l,
- -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019l,
- 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673l,
- -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863l,
- 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755l,
- -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158l,
- 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263l,
- -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365l,
- 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689l,
- -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648l,
- 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115l,
- -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256l,
- 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142l,
- -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4l,
- 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6l
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4),
+ BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6)
};
const RealType as = a*a;
@@ -439,44 +440,48 @@
*/
const unsigned short m = 19;
- static const RealType pts[] = {0.0016634282895983227941l,
- 0.014904509242697054183l,
- 0.04103478879005817919l,
- 0.079359853513391511008l,
- 0.1288612130237615133l,
- 0.18822336642448518856l,
- 0.25586876186122962384l,
- 0.32999972011807857222l,
- 0.40864620815774761438l,
- 0.48971819306044782365l,
- 0.57106118513245543894l,
- 0.6505134942981533829l,
- 0.72596367859928091618l,
- 0.79540665919549865924l,
- 0.85699701386308739244l,
- 0.90909804422384697594l,
- 0.95032536436570154409l,
- 0.97958418733152273717l,
- 0.99610366384229088321l};
- static const RealType wts[] = {0.012975111395684900835l,
- 0.012888764187499150078l,
- 0.012716644398857307844l,
- 0.012459897461364705691l,
- 0.012120231988292330388l,
- 0.011699908404856841158l,
- 0.011201723906897224448l,
- 0.010628993848522759853l,
- 0.0099855296835573320047l,
- 0.0092756136096132857933l,
- 0.0085039700881139589055l,
- 0.0076757344408814561254l,
- 0.0067964187616556459109l,
- 0.005871875456524750363l,
- 0.0049082589542498110071l,
- 0.0039119870792519721409l,
- 0.0028897090921170700834l,
- 0.0018483371329504443947l,
- 0.00079623320100438873578l};
+ static const RealType pts[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321)
+ };
+ static const RealType wts[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578)
+ };
const RealType as = a*a;
const RealType hs = -h*h*boost::math::constants::half<RealType>();
@@ -528,6 +533,66 @@
return val;
} // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m)
+ template <class T, class Policy>
+ T owens_t_T1_accelerated(T h, T a, const Policy& pol)
+ {
+ //
+ // This is the same series as T1, but:
+ // * The Taylor series for atan has been combined with that for T1,
+ // reducing but not eliminating cancellation error.
+ // * The resulting alternating series is then accelerated using method 1
+ // from H. Cohen, F. Rodriguez Villegas, D. Zagier,
+ // "Convergence acceleration of alternating series", Bonn, (1991).
+ //
+ BOOST_MATH_STD_USING
+ static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)";
+ T half_h_h = h * h / 2;
+ T a_pow = a;
+ T aa = a * a;
+ T exp_term = exp(-h * h / 2);
+ T one_minus_dj_sum = exp_term;
+ T sum = a_pow * exp_term;
+ T dj_pow = exp_term;
+ T term = sum;
+ int j = 1;
+
+ int n = (boost::math::tools::digits<T>() * 392) / 1000;
+ n += itrunc(h*h*a*a / 2);
+ if(tools::log_max_value<T>() / n < 5.9)
+ policies::raise_evaluation_error(function, 0, T(0), pol);
+ T d = pow(3 + sqrt(T(8)), n);
+ d = (d + 1 / d) / 2;
+ T b = -1;
+ T c = -d;
+ c = b - c;
+ sum *= c;
+ b = -n * n * b * 2;
+
+ // We use this to check for gross cancellation error,
+ // if "sum" ever drops below this value
+ // then we would probably be better off using the long double
+ // precision routines. Note that this is a heuristic, it
+ // is not foolproof!
+ T lim = -sum;
+
+ while(j < n)
+ {
+ a_pow *= aa;
+ dj_pow *= half_h_h / j;
+ one_minus_dj_sum += dj_pow;
+ term = one_minus_dj_sum * a_pow / (2 * j + 1);
+ c = b - c;
+ sum += c * term;
+ if(sum < lim)
+ policies::raise_evaluation_error(function, 0, T(0), pol);
+ b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
+ ++j;
+ }
+ if(sum < 0) // sum must always be positive, if it's negative something really bad has happend:
+ policies::raise_evaluation_error(function, 0, T(0), pol);
+ return (sum / d) / boost::math::constants::two_pi<T>();
+ }
+
// This routine dispatches the call to one of six subroutines, depending on the values
// of h and a.
// preconditions: h >= 0, 0<=a<=1, ah=a*h
@@ -539,7 +604,23 @@
RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
- const unsigned short icode = owens_t_compute_code(h,a);
+ if(policies::digits<RealType, Policy>() > 64)
+ {
+ // Attempt arbitrary precision code, this will throw if it goes wrong:
+ try
+ {
+ typedef boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
+ val = owens_t_T1_accelerated(h, a, forwarding_policy());
+ return val;
+ }
+ catch(const boost::math::evaluation_error&)
+ {
+ // Do nothing, just fall through:
+ }
+ }
+
+
+ const unsigned short icode = owens_t_compute_code(h, a);
const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
// determine the appropriate method, T1 ... T6
Modified: sandbox/math/libs/math/test/owens_t.ipp
==============================================================================
--- sandbox/math/libs/math/test/owens_t.ipp (original)
+++ sandbox/math/libs/math/test/owens_t.ipp 2012-04-17 12:15:34 EDT (Tue, 17 Apr 2012)
@@ -1,4 +1,6 @@
-#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
+#ifndef SC_
+# define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
+#endif
static const boost::array<boost::array<T, 3>, 400> owens_t = {{
{{ SC_(1.9508080482482910156250000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), SC_(9.7540378570556640625000000000000000000000000000000000000000000000000000000000000000000000000000000000e-02), SC_(2.2942365980155351117754771526378292708740231264960539316842063696200884829628237039895484898529192823e-03) }},
{{ SC_(1.9508080482482910156250000000000000000000000000000000000000000000000000000000000000000000000000000000e+00), SC_(1.2698680162429809570312500000000000000000000000000000000000000000000000000000000000000000000000000000e-01), SC_(2.9680322689154555657058017082736661495213677947212098762543622722763825563827711568131850564334565506e-03) }},
@@ -401,5 +403,3 @@
{{ SC_(1.9929222106933593750000000000000000000000000000000000000000000000000000000000000000000000000000000000e+01), SC_(9.9288129806518554687500000000000000000000000000000000000000000000000000000000000000000000000000000000e-01), SC_(5.6765941572108102929182392406617202048399346202690008425365121190758925914149018427767218487038675738e-89) }},
{{ SC_(1.9929222106933593750000000000000000000000000000000000000000000000000000000000000000000000000000000000e+01), SC_(9.9646115303039550781250000000000000000000000000000000000000000000000000000000000000000000000000000000e-01), SC_(5.6765941572108102929182392406617202048399346202690008425365121190758925914149018427767301460804161903e-89) }}
}};
-#undef SC_
-
Modified: sandbox/math/libs/math/test/owens_t_large_data.ipp
==============================================================================
--- sandbox/math/libs/math/test/owens_t_large_data.ipp (original)
+++ sandbox/math/libs/math/test/owens_t_large_data.ipp 2012-04-17 12:15:34 EDT (Tue, 17 Apr 2012)
@@ -1,4 +1,6 @@
+#ifndef SC_
#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
+#endif
static const boost::array<boost::array<T, 3>, 657> owens_t_large_data = {{
{{ SC_(1.7306551853835117071866989135742187500000000000000000000000000000000000000000000000000000000000000000e-06), SC_(1.7306551853835117071866989135742187500000000000000000000000000000000000000000000000000000000000000000e-06), SC_(2.7544232754071888169814325928268239884486372691874342956812521082730648643186974640936808405777870468e-07) }},
{{ SC_(1.7306551853835117071866989135742187500000000000000000000000000000000000000000000000000000000000000000e-06), SC_(2.1657506295014172792434692382812500000000000000000000000000000000000000000000000000000000000000000000e-06), SC_(3.4468991818847944488358292096162674677254240401591803431829935684111943192188266678541525762196399584e-07) }},
@@ -658,5 +660,3 @@
{{ SC_(3.1387406250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+04), SC_(6.3875340856611728668212890625000000000000000000000000000000000000000000000000000000000000000000000000e-04), SC_(5.1627835650444864618495857270828390674008441557196974871479139803402308298850057953535052644800142655e-213926795) }},
{{ SC_(3.1387406250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e+04), SC_(1.0718167759478092193603515625000000000000000000000000000000000000000000000000000000000000000000000000e-03), SC_(5.1627835650444864618495857270828390674008441557196974871479139803402308298850057953535053712048715547e-213926795) }}
}};
-#undef SC_
-
Modified: sandbox/math/libs/math/test/test_owens_t.cpp
==============================================================================
--- sandbox/math/libs/math/test/test_owens_t.cpp (original)
+++ sandbox/math/libs/math/test/test_owens_t.cpp 2012-04-17 12:15:34 EDT (Tue, 17 Apr 2012)
@@ -41,6 +41,18 @@
#include "libs/math/test/table_type.hpp"
#include "libs/math/test/functor.hpp"
+//
+// Defining TEST_CPP_DEC_FLOAT enables testing of multiprecision support.
+// This requires the multiprecision library from sandbox/big_number.
+// Note that these tests *do not pass*, but they do give an idea of the
+// error rates that can be expected....
+//
+#ifdef TEST_CPP_DEC_FLOAT
+#include <boost/multiprecision/cpp_dec_float.hpp>
+using namespace boost;
+#define SC_(x) BOOST_MATH_BIG_CONSTANT(T, std::numeric_limits<T>::digits, x)
+#endif
+
#include "owens_t_T7.hpp"
#include <iostream>
@@ -260,6 +272,9 @@
test_owens_t(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
#endif
#endif
+#ifdef TEST_CPP_DEC_FLOAT
+ test_owens_t(boost::multiprecision::cpp_dec_float_100(0), "cpp_dec_float"); // Test real concept.
+#endif
return 0;
} // int test_main(int, char* [])
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk