|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r48737 - in sandbox/math_toolkit: boost/math/bindings boost/math/bindings/detail boost/math/concepts boost/math/distributions boost/math/distributions/detail boost/math/special_functions boost/math/special_functions/detail boost/math/tools libs/math/test
From: john_at_[hidden]
Date: 2008-09-11 13:58:15
Author: johnmaddock
Date: 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
New Revision: 48737
URL: http://svn.boost.org/trac/boost/changeset/48737
Log:
Added mpfr support, and updated NTL::RR support.
Added concept checks for those.
Added type casts where required to get mpfr_class support working OK.
Added:
sandbox/math_toolkit/boost/math/bindings/detail/
sandbox/math_toolkit/boost/math/bindings/detail/big_digamma.hpp (contents, props changed)
sandbox/math_toolkit/boost/math/bindings/detail/big_lanczos.hpp (contents, props changed)
sandbox/math_toolkit/boost/math/bindings/mpfr.hpp (contents, props changed)
sandbox/math_toolkit/boost/math/concepts/real_type_concept.hpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/mpfr_concept_check.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/ntl_concept_check.cpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/bindings/rr.hpp | 162 +++++++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/boost/math/concepts/distributions.hpp | 2
sandbox/math_toolkit/boost/math/distributions/binomial.hpp | 4
sandbox/math_toolkit/boost/math/distributions/cauchy.hpp | 2
sandbox/math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp | 10 +-
sandbox/math_toolkit/boost/math/distributions/negative_binomial.hpp | 8
sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp | 6
sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp | 4
sandbox/math_toolkit/boost/math/distributions/poisson.hpp | 4
sandbox/math_toolkit/boost/math/distributions/triangular.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/bessel.hpp | 8
sandbox/math_toolkit/boost/math/special_functions/beta.hpp | 30 +++---
sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp | 6
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp | 18 ++--
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp | 16 +-
sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp | 6
sandbox/math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp | 6
sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp | 93 ++++++++++++++--------
sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp | 6
sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/erf.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/expint.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/expm1.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/factorials.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/gamma.hpp | 23 +++-
sandbox/math_toolkit/boost/math/special_functions/log1p.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/next.hpp | 10 +-
sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp | 8
sandbox/math_toolkit/boost/math/special_functions/trunc.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/zeta.hpp | 2
sandbox/math_toolkit/boost/math/tools/minima.hpp | 4
sandbox/math_toolkit/boost/math/tools/rational.hpp | 12 +-
sandbox/math_toolkit/boost/math/tools/roots.hpp | 6
sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp | 12 +-
39 files changed, 347 insertions(+), 159 deletions(-)
Added: sandbox/math_toolkit/boost/math/bindings/detail/big_digamma.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/bindings/detail/big_digamma.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,294 @@
+// (C) Copyright John Maddock 2006-8.
+// 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)
+
+#ifndef BOOST_MATH_NTL_DIGAMMA
+#define BOOST_MATH_NTL_DIGAMMA
+
+#include <boost/math/tools/rational.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/constants/constants.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T>
+T big_digamma_helper(T x)
+{
+ static const T P[61] = {
+ boost::lexical_cast<T>("0.6660133691143982067148122682345055274952e81"),
+ boost::lexical_cast<T>("0.6365271516829242456324234577164675383137e81"),
+ boost::lexical_cast<T>("0.2991038873096202943405966144203628966976e81"),
+ boost::lexical_cast<T>("0.9211116495503170498076013367421231351115e80"),
+ boost::lexical_cast<T>("0.2090792764676090716286400360584443891749e80"),
+ boost::lexical_cast<T>("0.3730037777359591428226035156377978092809e79"),
+ boost::lexical_cast<T>("0.5446396536956682043376492370432031543834e78"),
+ boost::lexical_cast<T>("0.6692523966335177847425047827449069256345e77"),
+ boost::lexical_cast<T>("0.7062543624100864681625612653756619116848e76"),
+ boost::lexical_cast<T>("0.6499914905966283735005256964443226879158e75"),
+ boost::lexical_cast<T>("0.5280364564853225211197557708655426736091e74"),
+ boost::lexical_cast<T>("0.3823205608981176913075543599005095206953e73"),
+ boost::lexical_cast<T>("0.2486733714214237704739129972671154532415e72"),
+ boost::lexical_cast<T>("0.1462562139602039577983434547171318011675e71"),
+ boost::lexical_cast<T>("0.7821169065036815012381267259559910324285e69"),
+ boost::lexical_cast<T>("0.3820552182348155468636157988764435365078e68"),
+ boost::lexical_cast<T>("0.1711618296983598244658239925535632505062e67"),
+ boost::lexical_cast<T>("0.7056661618357643731419080738521475204245e65"),
+ boost::lexical_cast<T>("0.2685246896473614017356264531791459936036e64"),
+ boost::lexical_cast<T>("0.9455168125599643085283071944864977592391e62"),
+ boost::lexical_cast<T>("0.3087541626972538362237309145177486236219e61"),
+ boost::lexical_cast<T>("0.9367928873352980208052601301625005737407e59"),
+ boost::lexical_cast<T>("0.2645306130689794942883818547314327466007e58"),
+ boost::lexical_cast<T>("0.6961815141171454309161007351079576190079e56"),
+ boost::lexical_cast<T>("0.1709637824471794552313802669803885946843e55"),
+ boost::lexical_cast<T>("0.3921553258481531526663112728778759311158e53"),
+ boost::lexical_cast<T>("0.8409006354449988687714450897575728228696e51"),
+ boost::lexical_cast<T>("0.1686755204461325935742097669030363344927e50"),
+ boost::lexical_cast<T>("0.3166653542877070999007425197729038754254e48"),
+ boost::lexical_cast<T>("0.5566029092358215049069560272835654229637e46"),
+ boost::lexical_cast<T>("0.9161766287916328133080586672953875116242e44"),
+ boost::lexical_cast<T>("1412317772330871298317974693514430627922000"),
+ boost::lexical_cast<T>("20387991986727877473732570146112459874790"),
+ boost::lexical_cast<T>("275557928713904105182512535678580359839.3"),
+ boost::lexical_cast<T>("3485719851040516559072031256589598330.723"),
+ boost::lexical_cast<T>("41247046743564028399938106707656877.40859"),
+ boost::lexical_cast<T>("456274078125709314602601667471879.0147312"),
+ boost::lexical_cast<T>("4714450683242899367025707077155.310613012"),
+ boost::lexical_cast<T>("45453933537925041680009544258.75073849996"),
+ boost::lexical_cast<T>("408437900487067278846361972.302331241052"),
+ boost::lexical_cast<T>("3415719344386166273085838.705771571751035"),
+ boost::lexical_cast<T>("26541502879185876562320.93134691487351145"),
+ boost::lexical_cast<T>("191261415065918713661.1571433274648417668"),
+ boost::lexical_cast<T>("1275349770108718421.645275944284937551702"),
+ boost::lexical_cast<T>("7849171120971773.318910987434906905704272"),
+ boost::lexical_cast<T>("44455946386549.80866460312682983576538056"),
+ boost::lexical_cast<T>("230920362395.3198137186361608905136598046"),
+ boost::lexical_cast<T>("1095700096.240863858624279930600654130254"),
+ boost::lexical_cast<T>("4727085.467506050153744334085516289728134"),
+ boost::lexical_cast<T>("18440.75118859447173303252421991479005424"),
+ boost::lexical_cast<T>("64.62515887799460295677071749181651317052"),
+ boost::lexical_cast<T>("0.201851568864688406206528472883512147547"),
+ boost::lexical_cast<T>("0.0005565091674187978029138500039504078098143"),
+ boost::lexical_cast<T>("0.1338097668312907986354698683493366559613e-5"),
+ boost::lexical_cast<T>("0.276308225077464312820179030238305271638e-8"),
+ boost::lexical_cast<T>("0.4801582970473168520375942100071070575043e-11"),
+ boost::lexical_cast<T>("0.6829184144212920949740376186058541800175e-14"),
+ boost::lexical_cast<T>("0.7634080076358511276617829524639455399182e-17"),
+ boost::lexical_cast<T>("0.6290035083727140966418512608156646142409e-20"),
+ boost::lexical_cast<T>("0.339652245667538733044036638506893821352e-23"),
+ boost::lexical_cast<T>("0.9017518064256388530773585529891677854909e-27")
+ };
+ static const T Q[61] = {
+ boost::lexical_cast<T>("0"),
+ boost::lexical_cast<T>("0.1386831185456898357379390197203894063459e81"),
+ boost::lexical_cast<T>("0.6467076379487574703291056110838151259438e81"),
+ boost::lexical_cast<T>("0.1394967823848615838336194279565285465161e82"),
+ boost::lexical_cast<T>("0.1872927317344192945218570366455046340458e82"),
+ boost::lexical_cast<T>("0.1772461045338946243584650759986310355937e82"),
+ boost::lexical_cast<T>("0.1267294892200258648315971144069595555118e82"),
+ boost::lexical_cast<T>("0.7157764838362416821508872117623058626589e81"),
+ boost::lexical_cast<T>("0.329447266909948668265277828268378274513e81"),
+ boost::lexical_cast<T>("0.1264376077317689779509250183194342571207e81"),
+ boost::lexical_cast<T>("0.4118230304191980787640446056583623228873e80"),
+ boost::lexical_cast<T>("0.1154393529762694616405952270558316515261e80"),
+ boost::lexical_cast<T>("0.281655612889423906125295485693696744275e79"),
+ boost::lexical_cast<T>("0.6037483524928743102724159846414025482077e78"),
+ boost::lexical_cast<T>("0.1145927995397835468123576831800276999614e78"),
+ boost::lexical_cast<T>("0.1938624296151985600348534009382865995154e77"),
+ boost::lexical_cast<T>("0.293980925856227626211879961219188406675e76"),
+ boost::lexical_cast<T>("0.4015574518216966910319562902099567437832e75"),
+ boost::lexical_cast<T>("0.4961475457509727343545565970423431880907e74"),
+ boost::lexical_cast<T>("0.5565482348278933960215521991000378896338e73"),
+ boost::lexical_cast<T>("0.5686112924615820754631098622770303094938e72"),
+ boost::lexical_cast<T>("0.5305988545844796293285410303747469932856e71"),
+ boost::lexical_cast<T>("0.4533363413802585060568537458067343491358e70"),
+ boost::lexical_cast<T>("0.3553932059473516064068322757331575565718e69"),
+ boost::lexical_cast<T>("0.2561198565218704414618802902533972354203e68"),
+ boost::lexical_cast<T>("0.1699519313292900324098102065697454295572e67"),
+ boost::lexical_cast<T>("0.1039830160862334505389615281373574959236e66"),
+ boost::lexical_cast<T>("0.5873082967977428281000961954715372504986e64"),
+ boost::lexical_cast<T>("0.3065255179030575882202133042549783442446e63"),
+ boost::lexical_cast<T>("0.1479494813481364701208655943688307245459e62"),
+ boost::lexical_cast<T>("0.6608150467921598615495180659808895663164e60"),
+ boost::lexical_cast<T>("0.2732535313770902021791888953487787496976e59"),
+ boost::lexical_cast<T>("0.1046402297662493314531194338414508049069e58"),
+ boost::lexical_cast<T>("0.3711375077192882936085049147920021549622e56"),
+ boost::lexical_cast<T>("0.1219154482883895482637944309702972234576e55"),
+ boost::lexical_cast<T>("0.3708359374149458741391374452286837880162e53"),
+ boost::lexical_cast<T>("0.1044095509971707189716913168889769471468e52"),
+ boost::lexical_cast<T>("0.271951506225063286130946773813524945052e50"),
+ boost::lexical_cast<T>("0.6548016291215163843464133978454065823866e48"),
+ boost::lexical_cast<T>("0.1456062447610542135403751730809295219344e47"),
+ boost::lexical_cast<T>("0.2986690175077969760978388356833006028929e45"),
+ boost::lexical_cast<T>("5643149706574013350061247429006443326844000"),
+ boost::lexical_cast<T>("98047545414467090421964387960743688053480"),
+ boost::lexical_cast<T>("1563378767746846395507385099301468978550"),
+ boost::lexical_cast<T>("22823360528584500077862274918382796495"),
+ boost::lexical_cast<T>("304215527004115213046601295970388750"),
+ boost::lexical_cast<T>("3690289075895685793844344966820325"),
+ boost::lexical_cast<T>("40584512015702371433911456606050"),
+ boost::lexical_cast<T>("402834190897282802772754873905"),
+ boost::lexical_cast<T>("3589522158493606918146495750"),
+ boost::lexical_cast<T>("28530557707503483723634725"),
+ boost::lexical_cast<T>("200714561335055753000730"),
+ boost::lexical_cast<T>("1237953783437761888641"),
+ boost::lexical_cast<T>("6614698701445762950"),
+ boost::lexical_cast<T>("30155495647727505"),
+ boost::lexical_cast<T>("114953256021450"),
+ boost::lexical_cast<T>("356398020013"),
+ boost::lexical_cast<T>("863113950"),
+ boost::lexical_cast<T>("1531345"),
+ boost::lexical_cast<T>("1770"),
+ boost::lexical_cast<T>("1")
+ };
+ static const T PD[60] = {
+ boost::lexical_cast<T>("0.6365271516829242456324234577164675383137e81"),
+ 2*boost::lexical_cast<T>("0.2991038873096202943405966144203628966976e81"),
+ 3*boost::lexical_cast<T>("0.9211116495503170498076013367421231351115e80"),
+ 4*boost::lexical_cast<T>("0.2090792764676090716286400360584443891749e80"),
+ 5*boost::lexical_cast<T>("0.3730037777359591428226035156377978092809e79"),
+ 6*boost::lexical_cast<T>("0.5446396536956682043376492370432031543834e78"),
+ 7*boost::lexical_cast<T>("0.6692523966335177847425047827449069256345e77"),
+ 8*boost::lexical_cast<T>("0.7062543624100864681625612653756619116848e76"),
+ 9*boost::lexical_cast<T>("0.6499914905966283735005256964443226879158e75"),
+ 10*boost::lexical_cast<T>("0.5280364564853225211197557708655426736091e74"),
+ 11*boost::lexical_cast<T>("0.3823205608981176913075543599005095206953e73"),
+ 12*boost::lexical_cast<T>("0.2486733714214237704739129972671154532415e72"),
+ 13*boost::lexical_cast<T>("0.1462562139602039577983434547171318011675e71"),
+ 14*boost::lexical_cast<T>("0.7821169065036815012381267259559910324285e69"),
+ 15*boost::lexical_cast<T>("0.3820552182348155468636157988764435365078e68"),
+ 16*boost::lexical_cast<T>("0.1711618296983598244658239925535632505062e67"),
+ 17*boost::lexical_cast<T>("0.7056661618357643731419080738521475204245e65"),
+ 18*boost::lexical_cast<T>("0.2685246896473614017356264531791459936036e64"),
+ 19*boost::lexical_cast<T>("0.9455168125599643085283071944864977592391e62"),
+ 20*boost::lexical_cast<T>("0.3087541626972538362237309145177486236219e61"),
+ 21*boost::lexical_cast<T>("0.9367928873352980208052601301625005737407e59"),
+ 22*boost::lexical_cast<T>("0.2645306130689794942883818547314327466007e58"),
+ 23*boost::lexical_cast<T>("0.6961815141171454309161007351079576190079e56"),
+ 24*boost::lexical_cast<T>("0.1709637824471794552313802669803885946843e55"),
+ 25*boost::lexical_cast<T>("0.3921553258481531526663112728778759311158e53"),
+ 26*boost::lexical_cast<T>("0.8409006354449988687714450897575728228696e51"),
+ 27*boost::lexical_cast<T>("0.1686755204461325935742097669030363344927e50"),
+ 28*boost::lexical_cast<T>("0.3166653542877070999007425197729038754254e48"),
+ 29*boost::lexical_cast<T>("0.5566029092358215049069560272835654229637e46"),
+ 30*boost::lexical_cast<T>("0.9161766287916328133080586672953875116242e44"),
+ 31*boost::lexical_cast<T>("1412317772330871298317974693514430627922000"),
+ 32*boost::lexical_cast<T>("20387991986727877473732570146112459874790"),
+ 33*boost::lexical_cast<T>("275557928713904105182512535678580359839.3"),
+ 34*boost::lexical_cast<T>("3485719851040516559072031256589598330.723"),
+ 35*boost::lexical_cast<T>("41247046743564028399938106707656877.40859"),
+ 36*boost::lexical_cast<T>("456274078125709314602601667471879.0147312"),
+ 37*boost::lexical_cast<T>("4714450683242899367025707077155.310613012"),
+ 38*boost::lexical_cast<T>("45453933537925041680009544258.75073849996"),
+ 39*boost::lexical_cast<T>("408437900487067278846361972.302331241052"),
+ 40*boost::lexical_cast<T>("3415719344386166273085838.705771571751035"),
+ 41*boost::lexical_cast<T>("26541502879185876562320.93134691487351145"),
+ 42*boost::lexical_cast<T>("191261415065918713661.1571433274648417668"),
+ 43*boost::lexical_cast<T>("1275349770108718421.645275944284937551702"),
+ 44*boost::lexical_cast<T>("7849171120971773.318910987434906905704272"),
+ 45*boost::lexical_cast<T>("44455946386549.80866460312682983576538056"),
+ 46*boost::lexical_cast<T>("230920362395.3198137186361608905136598046"),
+ 47*boost::lexical_cast<T>("1095700096.240863858624279930600654130254"),
+ 48*boost::lexical_cast<T>("4727085.467506050153744334085516289728134"),
+ 49*boost::lexical_cast<T>("18440.75118859447173303252421991479005424"),
+ 50*boost::lexical_cast<T>("64.62515887799460295677071749181651317052"),
+ 51*boost::lexical_cast<T>("0.201851568864688406206528472883512147547"),
+ 52*boost::lexical_cast<T>("0.0005565091674187978029138500039504078098143"),
+ 53*boost::lexical_cast<T>("0.1338097668312907986354698683493366559613e-5"),
+ 54*boost::lexical_cast<T>("0.276308225077464312820179030238305271638e-8"),
+ 55*boost::lexical_cast<T>("0.4801582970473168520375942100071070575043e-11"),
+ 56*boost::lexical_cast<T>("0.6829184144212920949740376186058541800175e-14"),
+ 57*boost::lexical_cast<T>("0.7634080076358511276617829524639455399182e-17"),
+ 58*boost::lexical_cast<T>("0.6290035083727140966418512608156646142409e-20"),
+ 59*boost::lexical_cast<T>("0.339652245667538733044036638506893821352e-23"),
+ 60*boost::lexical_cast<T>("0.9017518064256388530773585529891677854909e-27")
+ };
+ static const T QD[60] = {
+ boost::lexical_cast<T>("0.1386831185456898357379390197203894063459e81"),
+ 2*boost::lexical_cast<T>("0.6467076379487574703291056110838151259438e81"),
+ 3*boost::lexical_cast<T>("0.1394967823848615838336194279565285465161e82"),
+ 4*boost::lexical_cast<T>("0.1872927317344192945218570366455046340458e82"),
+ 5*boost::lexical_cast<T>("0.1772461045338946243584650759986310355937e82"),
+ 6*boost::lexical_cast<T>("0.1267294892200258648315971144069595555118e82"),
+ 7*boost::lexical_cast<T>("0.7157764838362416821508872117623058626589e81"),
+ 8*boost::lexical_cast<T>("0.329447266909948668265277828268378274513e81"),
+ 9*boost::lexical_cast<T>("0.1264376077317689779509250183194342571207e81"),
+ 10*boost::lexical_cast<T>("0.4118230304191980787640446056583623228873e80"),
+ 11*boost::lexical_cast<T>("0.1154393529762694616405952270558316515261e80"),
+ 12*boost::lexical_cast<T>("0.281655612889423906125295485693696744275e79"),
+ 13*boost::lexical_cast<T>("0.6037483524928743102724159846414025482077e78"),
+ 14*boost::lexical_cast<T>("0.1145927995397835468123576831800276999614e78"),
+ 15*boost::lexical_cast<T>("0.1938624296151985600348534009382865995154e77"),
+ 16*boost::lexical_cast<T>("0.293980925856227626211879961219188406675e76"),
+ 17*boost::lexical_cast<T>("0.4015574518216966910319562902099567437832e75"),
+ 18*boost::lexical_cast<T>("0.4961475457509727343545565970423431880907e74"),
+ 19*boost::lexical_cast<T>("0.5565482348278933960215521991000378896338e73"),
+ 20*boost::lexical_cast<T>("0.5686112924615820754631098622770303094938e72"),
+ 21*boost::lexical_cast<T>("0.5305988545844796293285410303747469932856e71"),
+ 22*boost::lexical_cast<T>("0.4533363413802585060568537458067343491358e70"),
+ 23*boost::lexical_cast<T>("0.3553932059473516064068322757331575565718e69"),
+ 24*boost::lexical_cast<T>("0.2561198565218704414618802902533972354203e68"),
+ 25*boost::lexical_cast<T>("0.1699519313292900324098102065697454295572e67"),
+ 26*boost::lexical_cast<T>("0.1039830160862334505389615281373574959236e66"),
+ 27*boost::lexical_cast<T>("0.5873082967977428281000961954715372504986e64"),
+ 28*boost::lexical_cast<T>("0.3065255179030575882202133042549783442446e63"),
+ 29*boost::lexical_cast<T>("0.1479494813481364701208655943688307245459e62"),
+ 30*boost::lexical_cast<T>("0.6608150467921598615495180659808895663164e60"),
+ 31*boost::lexical_cast<T>("0.2732535313770902021791888953487787496976e59"),
+ 32*boost::lexical_cast<T>("0.1046402297662493314531194338414508049069e58"),
+ 33*boost::lexical_cast<T>("0.3711375077192882936085049147920021549622e56"),
+ 34*boost::lexical_cast<T>("0.1219154482883895482637944309702972234576e55"),
+ 35*boost::lexical_cast<T>("0.3708359374149458741391374452286837880162e53"),
+ 36*boost::lexical_cast<T>("0.1044095509971707189716913168889769471468e52"),
+ 37*boost::lexical_cast<T>("0.271951506225063286130946773813524945052e50"),
+ 38*boost::lexical_cast<T>("0.6548016291215163843464133978454065823866e48"),
+ 39*boost::lexical_cast<T>("0.1456062447610542135403751730809295219344e47"),
+ 40*boost::lexical_cast<T>("0.2986690175077969760978388356833006028929e45"),
+ 41*boost::lexical_cast<T>("5643149706574013350061247429006443326844000"),
+ 42*boost::lexical_cast<T>("98047545414467090421964387960743688053480"),
+ 43*boost::lexical_cast<T>("1563378767746846395507385099301468978550"),
+ 44*boost::lexical_cast<T>("22823360528584500077862274918382796495"),
+ 45*boost::lexical_cast<T>("304215527004115213046601295970388750"),
+ 46*boost::lexical_cast<T>("3690289075895685793844344966820325"),
+ 47*boost::lexical_cast<T>("40584512015702371433911456606050"),
+ 48*boost::lexical_cast<T>("402834190897282802772754873905"),
+ 49*boost::lexical_cast<T>("3589522158493606918146495750"),
+ 50*boost::lexical_cast<T>("28530557707503483723634725"),
+ 51*boost::lexical_cast<T>("200714561335055753000730"),
+ 52*boost::lexical_cast<T>("1237953783437761888641"),
+ 53*boost::lexical_cast<T>("6614698701445762950"),
+ 54*boost::lexical_cast<T>("30155495647727505"),
+ 55*boost::lexical_cast<T>("114953256021450"),
+ 56*boost::lexical_cast<T>("356398020013"),
+ 57*boost::lexical_cast<T>("863113950"),
+ 58*boost::lexical_cast<T>("1531345"),
+ 59*boost::lexical_cast<T>("1770"),
+ 60*boost::lexical_cast<T>("1")
+ };
+ static const double g = 63.192152;
+
+ T zgh = x + g - 0.5;
+
+ T result = (x - 0.5) / zgh;
+ result += log(zgh);
+ result += tools::evaluate_polynomial(PD, x) / tools::evaluate_polynomial(P, x);
+ result -= tools::evaluate_polynomial(QD, x) / tools::evaluate_polynomial(Q, x);
+ result -= 1;
+
+ return result;
+}
+
+template <class T>
+T big_digamma(T x)
+{
+ BOOST_MATH_STD_USING
+ if(x < 0)
+ {
+ return big_digamma_helper(static_cast<T>(1-x)) + constants::pi<T>() / tan(constants::pi<T>() * (1-x));
+ }
+ return big_digamma_helper(x);
+}
+
+}}}
+
+#endif // include guard
Added: sandbox/math_toolkit/boost/math/bindings/detail/big_lanczos.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/bindings/detail/big_lanczos.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,887 @@
+// (C) Copyright John Maddock 2006-8.
+// 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)
+
+#ifndef BOOST_BIG_LANCZOS_HPP
+#define BOOST_BIG_LANCZOS_HPP
+
+#include <boost/math/special_functions/lanczos.hpp>
+#include <boost/lexical_cast.hpp>
+
+namespace boost{ namespace math{ namespace lanczos{
+
+struct lanczos13UDT
+{
+ template <class T>
+ static T lanczos_sum(const T& z)
+ {
+ static const T num[13] = {
+ boost::lexical_cast<T>("44012138428004.60895436261759919070125699"),
+ boost::lexical_cast<T>("41590453358593.20051581730723108131357995"),
+ boost::lexical_cast<T>("18013842787117.99677796276038389462742949"),
+ boost::lexical_cast<T>("4728736263475.388896889723995205703970787"),
+ boost::lexical_cast<T>("837910083628.4046470415724300225777912264"),
+ boost::lexical_cast<T>("105583707273.4299344907359855510105321192"),
+ boost::lexical_cast<T>("9701363618.494999493386608345339104922694"),
+ boost::lexical_cast<T>("654914397.5482052641016767125048538245644"),
+ boost::lexical_cast<T>("32238322.94213356530668889463945849409184"),
+ boost::lexical_cast<T>("1128514.219497091438040721811544858643121"),
+ boost::lexical_cast<T>("26665.79378459858944762533958798805525125"),
+ boost::lexical_cast<T>("381.8801248632926870394389468349331394196"),
+ boost::lexical_cast<T>("2.506628274631000502415763426076722427007"),
+ };
+ static const boost::uint32_t denom[13] = {
+ static_cast<boost::uint32_t>(0u),
+ static_cast<boost::uint32_t>(39916800u),
+ static_cast<boost::uint32_t>(120543840u),
+ static_cast<boost::uint32_t>(150917976u),
+ static_cast<boost::uint32_t>(105258076u),
+ static_cast<boost::uint32_t>(45995730u),
+ static_cast<boost::uint32_t>(13339535u),
+ static_cast<boost::uint32_t>(2637558u),
+ static_cast<boost::uint32_t>(357423u),
+ static_cast<boost::uint32_t>(32670u),
+ static_cast<boost::uint32_t>(1925u),
+ static_cast<boost::uint32_t>(66u),
+ static_cast<boost::uint32_t>(1u),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 13);
+ }
+
+ template <class T>
+ static T lanczos_sum_expG_scaled(const T& z)
+ {
+ static const T num[13] = {
+ boost::lexical_cast<T>("86091529.53418537217994842267760536134841"),
+ boost::lexical_cast<T>("81354505.17858011242874285785316135398567"),
+ boost::lexical_cast<T>("35236626.38815461910817650960734605416521"),
+ boost::lexical_cast<T>("9249814.988024471294683815872977672237195"),
+ boost::lexical_cast<T>("1639024.216687146960253839656643518985826"),
+ boost::lexical_cast<T>("206530.8157641225032631778026076868855623"),
+ boost::lexical_cast<T>("18976.70193530288915698282139308582105936"),
+ boost::lexical_cast<T>("1281.068909912559479885759622791374106059"),
+ boost::lexical_cast<T>("63.06093343420234536146194868906771599354"),
+ boost::lexical_cast<T>("2.207470909792527638222674678171050209691"),
+ boost::lexical_cast<T>("0.05216058694613505427476207805814960742102"),
+ boost::lexical_cast<T>("0.0007469903808915448316510079585999893674101"),
+ boost::lexical_cast<T>("0.4903180573459871862552197089738373164184e-5"),
+ };
+ static const boost::uint32_t denom[13] = {
+ static_cast<boost::uint32_t>(0u),
+ static_cast<boost::uint32_t>(39916800u),
+ static_cast<boost::uint32_t>(120543840u),
+ static_cast<boost::uint32_t>(150917976u),
+ static_cast<boost::uint32_t>(105258076u),
+ static_cast<boost::uint32_t>(45995730u),
+ static_cast<boost::uint32_t>(13339535u),
+ static_cast<boost::uint32_t>(2637558u),
+ static_cast<boost::uint32_t>(357423u),
+ static_cast<boost::uint32_t>(32670u),
+ static_cast<boost::uint32_t>(1925u),
+ static_cast<boost::uint32_t>(66u),
+ static_cast<boost::uint32_t>(1u),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 13);
+ }
+
+
+ template<class T>
+ static T lanczos_sum_near_1(const T& dz)
+ {
+ static const T d[12] = {
+ boost::lexical_cast<T>("4.832115561461656947793029596285626840312"),
+ boost::lexical_cast<T>("-19.86441536140337740383120735104359034688"),
+ boost::lexical_cast<T>("33.9927422807443239927197864963170585331"),
+ boost::lexical_cast<T>("-31.41520692249765980987427413991250886138"),
+ boost::lexical_cast<T>("17.0270866009599345679868972409543597821"),
+ boost::lexical_cast<T>("-5.5077216950865501362506920516723682167"),
+ boost::lexical_cast<T>("1.037811741948214855286817963800439373362"),
+ boost::lexical_cast<T>("-0.106640468537356182313660880481398642811"),
+ boost::lexical_cast<T>("0.005276450526660653288757565778182586742831"),
+ boost::lexical_cast<T>("-0.0001000935625597121545867453746252064770029"),
+ boost::lexical_cast<T>("0.462590910138598083940803704521211569234e-6"),
+ boost::lexical_cast<T>("-0.1735307814426389420248044907765671743012e-9"),
+ };
+ T result = 0;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(k*dz + k*k);
+ }
+ return result;
+ }
+
+ template<class T>
+ static T lanczos_sum_near_2(const T& dz)
+ {
+ static const T d[12] = {
+ boost::lexical_cast<T>("26.96979819614830698367887026728396466395"),
+ boost::lexical_cast<T>("-110.8705424709385114023884328797900204863"),
+ boost::lexical_cast<T>("189.7258846119231466417015694690434770085"),
+ boost::lexical_cast<T>("-175.3397202971107486383321670769397356553"),
+ boost::lexical_cast<T>("95.03437648691551457087250340903980824948"),
+ boost::lexical_cast<T>("-30.7406022781665264273675797983497141978"),
+ boost::lexical_cast<T>("5.792405601630517993355102578874590410552"),
+ boost::lexical_cast<T>("-0.5951993240669148697377539518639997795831"),
+ boost::lexical_cast<T>("0.02944979359164017509944724739946255067671"),
+ boost::lexical_cast<T>("-0.0005586586555377030921194246330399163602684"),
+ boost::lexical_cast<T>("0.2581888478270733025288922038673392636029e-5"),
+ boost::lexical_cast<T>("-0.9685385411006641478305219367315965391289e-9"),
+ };
+ T result = 0;
+ T z = dz + 2;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
+ }
+ return result;
+ }
+
+ static double g(){ return 13.1445650000000000545696821063756942749; }
+};
+
+
+//
+// Lanczos Coefficients for N=22 G=22.61891
+// Max experimental error (with arbitary precision arithmetic) 2.9524e-38
+// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006
+//
+struct lanczos22UDT
+{
+ template <class T>
+ static T lanczos_sum(const T& z)
+ {
+ static const T num[22] = {
+ boost::lexical_cast<T>("46198410803245094237463011094.12173081986"),
+ boost::lexical_cast<T>("43735859291852324413622037436.321513777"),
+ boost::lexical_cast<T>("19716607234435171720534556386.97481377748"),
+ boost::lexical_cast<T>("5629401471315018442177955161.245623932129"),
+ boost::lexical_cast<T>("1142024910634417138386281569.245580222392"),
+ boost::lexical_cast<T>("175048529315951173131586747.695329230778"),
+ boost::lexical_cast<T>("21044290245653709191654675.41581372963167"),
+ boost::lexical_cast<T>("2033001410561031998451380.335553678782601"),
+ boost::lexical_cast<T>("160394318862140953773928.8736211601848891"),
+ boost::lexical_cast<T>("10444944438396359705707.48957290388740896"),
+ boost::lexical_cast<T>("565075825801617290121.1466393747967538948"),
+ boost::lexical_cast<T>("25475874292116227538.99448534450411942597"),
+ boost::lexical_cast<T>("957135055846602154.6720835535232270205725"),
+ boost::lexical_cast<T>("29874506304047462.23662392445173880821515"),
+ boost::lexical_cast<T>("769651310384737.2749087590725764959689181"),
+ boost::lexical_cast<T>("16193289100889.15989633624378404096011797"),
+ boost::lexical_cast<T>("273781151680.6807433264462376754578933261"),
+ boost::lexical_cast<T>("3630485900.32917021712188739762161583295"),
+ boost::lexical_cast<T>("36374352.05577334277856865691538582936484"),
+ boost::lexical_cast<T>("258945.7742115532455441786924971194951043"),
+ boost::lexical_cast<T>("1167.501919472435718934219997431551246996"),
+ boost::lexical_cast<T>("2.50662827463100050241576528481104525333"),
+ };
+ static const T denom[22] = {
+ boost::lexical_cast<T>("0"),
+ boost::lexical_cast<T>("2432902008176640000"),
+ boost::lexical_cast<T>("8752948036761600000"),
+ boost::lexical_cast<T>("13803759753640704000"),
+ boost::lexical_cast<T>("12870931245150988800"),
+ boost::lexical_cast<T>("8037811822645051776"),
+ boost::lexical_cast<T>("3599979517947607200"),
+ boost::lexical_cast<T>("1206647803780373360"),
+ boost::lexical_cast<T>("311333643161390640"),
+ boost::lexical_cast<T>("63030812099294896"),
+ boost::lexical_cast<T>("10142299865511450"),
+ boost::lexical_cast<T>("1307535010540395"),
+ boost::lexical_cast<T>("135585182899530"),
+ boost::lexical_cast<T>("11310276995381"),
+ boost::lexical_cast<T>("756111184500"),
+ boost::lexical_cast<T>("40171771630"),
+ boost::lexical_cast<T>("1672280820"),
+ boost::lexical_cast<T>("53327946"),
+ boost::lexical_cast<T>("1256850"),
+ boost::lexical_cast<T>("20615"),
+ boost::lexical_cast<T>("210"),
+ boost::lexical_cast<T>("1"),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 22);
+ }
+
+ template <class T>
+ static T lanczos_sum_expG_scaled(const T& z)
+ {
+ static const T num[22] = {
+ boost::lexical_cast<T>("6939996264376682180.277485395074954356211"),
+ boost::lexical_cast<T>("6570067992110214451.87201438870245659384"),
+ boost::lexical_cast<T>("2961859037444440551.986724631496417064121"),
+ boost::lexical_cast<T>("845657339772791245.3541226499766163431651"),
+ boost::lexical_cast<T>("171556737035449095.2475716923888737881837"),
+ boost::lexical_cast<T>("26296059072490867.7822441885603400926007"),
+ boost::lexical_cast<T>("3161305619652108.433798300149816829198706"),
+ boost::lexical_cast<T>("305400596026022.4774396904484542582526472"),
+ boost::lexical_cast<T>("24094681058862.55120507202622377623528108"),
+ boost::lexical_cast<T>("1569055604375.919477574824168939428328839"),
+ boost::lexical_cast<T>("84886558909.02047889339710230696942513159"),
+ boost::lexical_cast<T>("3827024985.166751989686050643579753162298"),
+ boost::lexical_cast<T>("143782298.9273215199098728674282885500522"),
+ boost::lexical_cast<T>("4487794.24541641841336786238909171265944"),
+ boost::lexical_cast<T>("115618.2025760830513505888216285273541959"),
+ boost::lexical_cast<T>("2432.580773108508276957461757328744780439"),
+ boost::lexical_cast<T>("41.12782532742893597168530008461874360191"),
+ boost::lexical_cast<T>("0.5453771709477689805460179187388702295792"),
+ boost::lexical_cast<T>("0.005464211062612080347167337964166505282809"),
+ boost::lexical_cast<T>("0.388992321263586767037090706042788910953e-4"),
+ boost::lexical_cast<T>("0.1753839324538447655939518484052327068859e-6"),
+ boost::lexical_cast<T>("0.3765495513732730583386223384116545391759e-9"),
+ };
+ static const T denom[22] = {
+ boost::lexical_cast<T>("0"),
+ boost::lexical_cast<T>("2432902008176640000"),
+ boost::lexical_cast<T>("8752948036761600000"),
+ boost::lexical_cast<T>("13803759753640704000"),
+ boost::lexical_cast<T>("12870931245150988800"),
+ boost::lexical_cast<T>("8037811822645051776"),
+ boost::lexical_cast<T>("3599979517947607200"),
+ boost::lexical_cast<T>("1206647803780373360"),
+ boost::lexical_cast<T>("311333643161390640"),
+ boost::lexical_cast<T>("63030812099294896"),
+ boost::lexical_cast<T>("10142299865511450"),
+ boost::lexical_cast<T>("1307535010540395"),
+ boost::lexical_cast<T>("135585182899530"),
+ boost::lexical_cast<T>("11310276995381"),
+ boost::lexical_cast<T>("756111184500"),
+ boost::lexical_cast<T>("40171771630"),
+ boost::lexical_cast<T>("1672280820"),
+ boost::lexical_cast<T>("53327946"),
+ boost::lexical_cast<T>("1256850"),
+ boost::lexical_cast<T>("20615"),
+ boost::lexical_cast<T>("210"),
+ boost::lexical_cast<T>("1"),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 22);
+ }
+
+
+ template<class T>
+ static T lanczos_sum_near_1(const T& dz)
+ {
+ static const T d[21] = {
+ boost::lexical_cast<T>("8.318998691953337183034781139546384476554"),
+ boost::lexical_cast<T>("-63.15415991415959158214140353299240638675"),
+ boost::lexical_cast<T>("217.3108224383632868591462242669081540163"),
+ boost::lexical_cast<T>("-448.5134281386108366899784093610397354889"),
+ boost::lexical_cast<T>("619.2903759363285456927248474593012711346"),
+ boost::lexical_cast<T>("-604.1630177420625418522025080080444177046"),
+ boost::lexical_cast<T>("428.8166750424646119935047118287362193314"),
+ boost::lexical_cast<T>("-224.6988753721310913866347429589434550302"),
+ boost::lexical_cast<T>("87.32181627555510833499451817622786940961"),
+ boost::lexical_cast<T>("-25.07866854821128965662498003029199058098"),
+ boost::lexical_cast<T>("5.264398125689025351448861011657789005392"),
+ boost::lexical_cast<T>("-0.792518936256495243383586076579921559914"),
+ boost::lexical_cast<T>("0.08317448364744713773350272460937904691566"),
+ boost::lexical_cast<T>("-0.005845345166274053157781068150827567998882"),
+ boost::lexical_cast<T>("0.0002599412126352082483326238522490030412391"),
+ boost::lexical_cast<T>("-0.6748102079670763884917431338234783496303e-5"),
+ boost::lexical_cast<T>("0.908824383434109002762325095643458603605e-7"),
+ boost::lexical_cast<T>("-0.5299325929309389890892469299969669579725e-9"),
+ boost::lexical_cast<T>("0.994306085859549890267983602248532869362e-12"),
+ boost::lexical_cast<T>("-0.3499893692975262747371544905820891835298e-15"),
+ boost::lexical_cast<T>("0.7260746353663365145454867069182884694961e-20"),
+ };
+ T result = 0;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(k*dz + k*k);
+ }
+ return result;
+ }
+
+ template<class T>
+ static T lanczos_sum_near_2(const T& dz)
+ {
+ static const T d[21] = {
+ boost::lexical_cast<T>("75.39272007105208086018421070699575462226"),
+ boost::lexical_cast<T>("-572.3481967049935412452681346759966390319"),
+ boost::lexical_cast<T>("1969.426202741555335078065370698955484358"),
+ boost::lexical_cast<T>("-4064.74968778032030891520063865996757519"),
+ boost::lexical_cast<T>("5612.452614138013929794736248384309574814"),
+ boost::lexical_cast<T>("-5475.357667500026172903620177988213902339"),
+ boost::lexical_cast<T>("3886.243614216111328329547926490398103492"),
+ boost::lexical_cast<T>("-2036.382026072125407192448069428134470564"),
+ boost::lexical_cast<T>("791.3727954936062108045551843636692287652"),
+ boost::lexical_cast<T>("-227.2808432388436552794021219198885223122"),
+ boost::lexical_cast<T>("47.70974355562144229897637024320739257284"),
+ boost::lexical_cast<T>("-7.182373807798293545187073539819697141572"),
+ boost::lexical_cast<T>("0.7537866989631514559601547530490976100468"),
+ boost::lexical_cast<T>("-0.05297470142240154822658739758236594717787"),
+ boost::lexical_cast<T>("0.00235577330936380542539812701472320434133"),
+ boost::lexical_cast<T>("-0.6115613067659273118098229498679502138802e-4"),
+ boost::lexical_cast<T>("0.8236417010170941915758315020695551724181e-6"),
+ boost::lexical_cast<T>("-0.4802628430993048190311242611330072198089e-8"),
+ boost::lexical_cast<T>("0.9011113376981524418952720279739624707342e-11"),
+ boost::lexical_cast<T>("-0.3171854152689711198382455703658589996796e-14"),
+ boost::lexical_cast<T>("0.6580207998808093935798753964580596673177e-19"),
+ };
+ T result = 0;
+ T z = dz + 2;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
+ }
+ return result;
+ }
+
+ static double g(){ return 22.61890999999999962710717227309942245483; }
+};
+
+//
+// Lanczos Coefficients for N=31 G=32.08067
+// Max experimental error (with arbitary precision arithmetic) 0.162e-52
+// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at May 9 2006
+//
+struct lanczos31UDT
+{
+ template <class T>
+ static T lanczos_sum(const T& z)
+ {
+ static const T num[31] = {
+ boost::lexical_cast<T>("0.2579646553333513328235723061836959833277e46"),
+ boost::lexical_cast<T>("0.2444796504337453845497419271639377138264e46"),
+ boost::lexical_cast<T>("0.1119885499016017172212179730662673475329e46"),
+ boost::lexical_cast<T>("0.3301983829072723658949204487793889113715e45"),
+ boost::lexical_cast<T>("0.7041171040503851585152895336505379417066e44"),
+ boost::lexical_cast<T>("0.1156687509001223855125097826246939403504e44"),
+ boost::lexical_cast<T>("1522559363393940883866575697565974893306000"),
+ boost::lexical_cast<T>("164914363507650839510801418717701057005700"),
+ boost::lexical_cast<T>("14978522943127593263654178827041568394060"),
+ boost::lexical_cast<T>("1156707153701375383907746879648168666774"),
+ boost::lexical_cast<T>("76739431129980851159755403434593664173.2"),
+ boost::lexical_cast<T>("4407916278928188620282281495575981079.306"),
+ boost::lexical_cast<T>("220487883931812802092792125175269667.3004"),
+ boost::lexical_cast<T>("9644828280794966468052381443992828.433924"),
+ boost::lexical_cast<T>("369996467042247229310044531282837.6549068"),
+ boost::lexical_cast<T>("12468380890717344610932904378961.13494291"),
+ boost::lexical_cast<T>("369289245210898235894444657859.0529720075"),
+ boost::lexical_cast<T>("9607992460262594951559461829.34885209022"),
+ boost::lexical_cast<T>("219225935074853412540086410.981421315799"),
+ boost::lexical_cast<T>("4374309943598658046326340.720767382079549"),
+ boost::lexical_cast<T>("76008779092264509404014.10530947173485581"),
+ boost::lexical_cast<T>("1143503533822162444712.335663112617754987"),
+ boost::lexical_cast<T>("14779233719977576920.37884890049671578409"),
+ boost::lexical_cast<T>("162409028440678302.9992838032166348069916"),
+ boost::lexical_cast<T>("1496561553388385.733407609544964535634135"),
+ boost::lexical_cast<T>("11347624460661.81008311053190661436107043"),
+ boost::lexical_cast<T>("68944915931.32004991941950530448472223832"),
+ boost::lexical_cast<T>("322701221.6391432296123937035480931903651"),
+ boost::lexical_cast<T>("1092364.213992634267819050120261755371294"),
+ boost::lexical_cast<T>("2380.151399852411512711176940867823024864"),
+ boost::lexical_cast<T>("2.506628274631000502415765284811045253007"),
+ };
+ static const T denom[31] = {
+ boost::lexical_cast<T>("0"),
+ boost::lexical_cast<T>("0.8841761993739701954543616e31"),
+ boost::lexical_cast<T>("0.3502799997985980526649278464e32"),
+ boost::lexical_cast<T>("0.622621928420356134910574592e32"),
+ boost::lexical_cast<T>("66951000306085302338993639424000"),
+ boost::lexical_cast<T>("49361465831621147825759587123200"),
+ boost::lexical_cast<T>("26751280755793398822580822142976"),
+ boost::lexical_cast<T>("11139316913434780466101123891200"),
+ boost::lexical_cast<T>("3674201658710345201899117607040"),
+ boost::lexical_cast<T>("981347603630155088295475765440"),
+ boost::lexical_cast<T>("215760462268683520394805979744"),
+ boost::lexical_cast<T>("39539238727270799376544542000"),
+ boost::lexical_cast<T>("6097272817323042122728617800"),
+ boost::lexical_cast<T>("796974693974455191377937300"),
+ boost::lexical_cast<T>("88776380550648116217781890"),
+ boost::lexical_cast<T>("8459574446076318147830625"),
+ boost::lexical_cast<T>("691254538651580660999025"),
+ boost::lexical_cast<T>("48487623689430693038025"),
+ boost::lexical_cast<T>("2918939500751087661105"),
+ boost::lexical_cast<T>("150566737512021319125"),
+ boost::lexical_cast<T>("6634460278534540725"),
+ boost::lexical_cast<T>("248526574856284725"),
+ boost::lexical_cast<T>("7860403394108265"),
+ boost::lexical_cast<T>("207912996295875"),
+ boost::lexical_cast<T>("4539323721075"),
+ boost::lexical_cast<T>("80328850875"),
+ boost::lexical_cast<T>("1122686019"),
+ boost::lexical_cast<T>("11921175"),
+ boost::lexical_cast<T>("90335"),
+ boost::lexical_cast<T>("435"),
+ boost::lexical_cast<T>("1"),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 31);
+ }
+
+ template <class T>
+ static T lanczos_sum_expG_scaled(const T& z)
+ {
+ static const T num[31] = {
+ boost::lexical_cast<T>("30137154810677525966583148469478.52374216"),
+ boost::lexical_cast<T>("28561746428637727032849890123131.36314653"),
+ boost::lexical_cast<T>("13083250730789213354063781611435.74046294"),
+ boost::lexical_cast<T>("3857598154697777600846539129354.783647"),
+ boost::lexical_cast<T>("822596651552555685068015316144.0952185852"),
+ boost::lexical_cast<T>("135131964033213842052904200372.039133532"),
+ boost::lexical_cast<T>("17787555889683709693655685146.19771358863"),
+ boost::lexical_cast<T>("1926639793777927562221423874.149673297196"),
+ boost::lexical_cast<T>("174989113988888477076973808.6991839697774"),
+ boost::lexical_cast<T>("13513425905835560387095425.01158383184045"),
+ boost::lexical_cast<T>("896521313378762433091075.1446749283094845"),
+ boost::lexical_cast<T>("51496223433749515758124.71524415105430686"),
+ boost::lexical_cast<T>("2575886794780078381228.37205955912263407"),
+ boost::lexical_cast<T>("112677328855422964200.4155776009524490958"),
+ boost::lexical_cast<T>("4322545967487943330.625233358130724324796"),
+ boost::lexical_cast<T>("145663957202380774.0362027607207590519723"),
+ boost::lexical_cast<T>("4314283729473470.686566233465428332496534"),
+ boost::lexical_cast<T>("112246988185485.8877916434026906290603878"),
+ boost::lexical_cast<T>("2561143864972.040563435178307062626388193"),
+ boost::lexical_cast<T>("51103611767.9626550674442537989885239605"),
+ boost::lexical_cast<T>("887985348.0369447209508500133077232094491"),
+ boost::lexical_cast<T>("13359172.3954672607019822025834072685839"),
+ boost::lexical_cast<T>("172660.8841147568768783928167105965064459"),
+ boost::lexical_cast<T>("1897.370795407433013556725714874693719617"),
+ boost::lexical_cast<T>("17.48383210090980598861217644749573257178"),
+ boost::lexical_cast<T>("0.1325705316732132940835251054350153028901"),
+ boost::lexical_cast<T>("0.0008054605783673449641889260501816356090452"),
+ boost::lexical_cast<T>("0.377001130700104515644336869896819162464e-5"),
+ boost::lexical_cast<T>("0.1276172868883867038813825443204454996531e-7"),
+ boost::lexical_cast<T>("0.2780651912081116274907381023821492811093e-10"),
+ boost::lexical_cast<T>("0.2928410648650955854121639682890739211234e-13"),
+ };
+ static const T denom[31] = {
+ boost::lexical_cast<T>("0"),
+ boost::lexical_cast<T>("0.8841761993739701954543616e31"),
+ boost::lexical_cast<T>("0.3502799997985980526649278464e32"),
+ boost::lexical_cast<T>("0.622621928420356134910574592e32"),
+ boost::lexical_cast<T>("66951000306085302338993639424000"),
+ boost::lexical_cast<T>("49361465831621147825759587123200"),
+ boost::lexical_cast<T>("26751280755793398822580822142976"),
+ boost::lexical_cast<T>("11139316913434780466101123891200"),
+ boost::lexical_cast<T>("3674201658710345201899117607040"),
+ boost::lexical_cast<T>("981347603630155088295475765440"),
+ boost::lexical_cast<T>("215760462268683520394805979744"),
+ boost::lexical_cast<T>("39539238727270799376544542000"),
+ boost::lexical_cast<T>("6097272817323042122728617800"),
+ boost::lexical_cast<T>("796974693974455191377937300"),
+ boost::lexical_cast<T>("88776380550648116217781890"),
+ boost::lexical_cast<T>("8459574446076318147830625"),
+ boost::lexical_cast<T>("691254538651580660999025"),
+ boost::lexical_cast<T>("48487623689430693038025"),
+ boost::lexical_cast<T>("2918939500751087661105"),
+ boost::lexical_cast<T>("150566737512021319125"),
+ boost::lexical_cast<T>("6634460278534540725"),
+ boost::lexical_cast<T>("248526574856284725"),
+ boost::lexical_cast<T>("7860403394108265"),
+ boost::lexical_cast<T>("207912996295875"),
+ boost::lexical_cast<T>("4539323721075"),
+ boost::lexical_cast<T>("80328850875"),
+ boost::lexical_cast<T>("1122686019"),
+ boost::lexical_cast<T>("11921175"),
+ boost::lexical_cast<T>("90335"),
+ boost::lexical_cast<T>("435"),
+ boost::lexical_cast<T>("1"),
+ };
+ return boost::math::tools::evaluate_rational(num, denom, z, 31);
+ }
+
+
+ template<class T>
+ static T lanczos_sum_near_1(const T& dz)
+ {
+ static const T d[30] = {
+ boost::lexical_cast<T>("11.80038544942943603508206880307972596807"),
+ boost::lexical_cast<T>("-130.6355975335626214564236363322099481079"),
+ boost::lexical_cast<T>("676.2177719145993049893392276809256538927"),
+ boost::lexical_cast<T>("-2174.724497783850503069990936574060452057"),
+ boost::lexical_cast<T>("4869.877180638131076410069103742986502022"),
+ boost::lexical_cast<T>("-8065.744271864238179992762265472478229172"),
+ boost::lexical_cast<T>("10245.03825618572106228191509520638651539"),
+ boost::lexical_cast<T>("-10212.83902362683215459850403668669647192"),
+ boost::lexical_cast<T>("8110.289185383288952562767679576754140336"),
+ boost::lexical_cast<T>("-5179.310892558291062401828964000776095156"),
+ boost::lexical_cast<T>("2673.987492589052370230989109591011091273"),
+ boost::lexical_cast<T>("-1118.342574651205183051884250033505609141"),
+ boost::lexical_cast<T>("378.5812742511620662650096436471920295596"),
+ boost::lexical_cast<T>("-103.3725999812126067084828735543906768961"),
+ boost::lexical_cast<T>("22.62913974335996321848099677797888917792"),
+ boost::lexical_cast<T>("-3.936414819950859548507275533569588041446"),
+ boost::lexical_cast<T>("0.5376818198843817355682124535902641644854"),
+ boost::lexical_cast<T>("-0.0567827903603478957483409124122554243201"),
+ boost::lexical_cast<T>("0.004545544993648879420352693271088478106482"),
+ boost::lexical_cast<T>("-0.0002689795568951033950042375135970897959935"),
+ boost::lexical_cast<T>("0.1139493459006846530734617710847103572122e-4"),
+ boost::lexical_cast<T>("-0.3316581197839213921885210451302820192794e-6"),
+ boost::lexical_cast<T>("0.6285613334898374028443777562554713906213e-8"),
+ boost::lexical_cast<T>("-0.7222145115734409070310317999856424167091e-10"),
+ boost::lexical_cast<T>("0.4562976983547274766890241815002584238219e-12"),
+ boost::lexical_cast<T>("-0.1380593023819058919640038942493212141072e-14"),
+ boost::lexical_cast<T>("0.1629663871586410129307496385264268190679e-17"),
+ boost::lexical_cast<T>("-0.5429994291916548849493889660077076739993e-21"),
+ boost::lexical_cast<T>("0.2922682842441892106795386303084661338957e-25"),
+ boost::lexical_cast<T>("-0.8456967065309046044689041041336866118459e-31"),
+ };
+ T result = 0;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(k*dz + k*k);
+ }
+ return result;
+ }
+
+ template<class T>
+ static T lanczos_sum_near_2(const T& dz)
+ {
+ static const T d[30] = {
+ boost::lexical_cast<T>("147.9979641587472136175636384176549713358"),
+ boost::lexical_cast<T>("-1638.404318611773924210055619836375434296"),
+ boost::lexical_cast<T>("8480.981744216135641122944743711911653273"),
+ boost::lexical_cast<T>("-27274.93942104458448200467097634494071176"),
+ boost::lexical_cast<T>("61076.98019918759324489193232276937262854"),
+ boost::lexical_cast<T>("-101158.8762737154296509560513952101409264"),
+ boost::lexical_cast<T>("128491.1252383947174824913796141607174379"),
+ boost::lexical_cast<T>("-128087.2892038336581928787480535905496026"),
+ boost::lexical_cast<T>("101717.5492545853663296795562084430123258"),
+ boost::lexical_cast<T>("-64957.8330410311808907869707511362206858"),
+ boost::lexical_cast<T>("33536.59139229792478811870738772305570317"),
+ boost::lexical_cast<T>("-14026.01847115365926835980820243003785821"),
+ boost::lexical_cast<T>("4748.087094096186515212209389240715050212"),
+ boost::lexical_cast<T>("-1296.477510211815971152205100242259733245"),
+ boost::lexical_cast<T>("283.8099337545793198947620951499958085157"),
+ boost::lexical_cast<T>("-49.36969067101255103452092297769364837753"),
+ boost::lexical_cast<T>("6.743492833270653628580811118017061581404"),
+ boost::lexical_cast<T>("-0.7121578704864048548351804794951487823626"),
+ boost::lexical_cast<T>("0.0570092738016915476694118877057948681298"),
+ boost::lexical_cast<T>("-0.003373485297696102660302960722607722438643"),
+ boost::lexical_cast<T>("0.0001429128843527532859999752593761934089751"),
+ boost::lexical_cast<T>("-0.41595867130858508233493767243236888636e-5"),
+ boost::lexical_cast<T>("0.7883284669307241040059778207492255409785e-7"),
+ boost::lexical_cast<T>("-0.905786322462384670803148223703187214379e-9"),
+ boost::lexical_cast<T>("0.5722790216999820323272452464661250331451e-11"),
+ boost::lexical_cast<T>("-0.1731510870832349779315841757234562309727e-13"),
+ boost::lexical_cast<T>("0.2043890314358438601429048378015983874378e-16"),
+ boost::lexical_cast<T>("-0.6810185176079344204740000170500311171065e-20"),
+ boost::lexical_cast<T>("0.3665567641131713928114853776588342403919e-24"),
+ boost::lexical_cast<T>("-0.1060655106553177007425710511436497259484e-29"),
+ };
+ T result = 0;
+ T z = dz + 2;
+ for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
+ }
+ return result;
+ }
+
+ static double g(){ return 32.08066999999999779902282170951366424561; }
+};
+
+//
+// Lanczos Coefficients for N=61 G=63.192152
+// Max experimental error (with 1000-bit precision arithmetic) 3.740e-113
+// Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 12 2006
+//
+struct lanczos61UDT
+{
+ template <class T>
+ static T lanczos_sum(const T& z)
+ {
+ using namespace boost;
+ static const T d[61] = {
+ boost::lexical_cast<T>("2.50662827463100050241576528481104525300698674060993831662992357634229365460784197494659584"),
+ boost::lexical_cast<T>("13349415823254323512107320481.3495396037261649201426994438803767191136434970492309775123879"),
+ boost::lexical_cast<T>("-300542621510568204264185787475.230003734889859348050696493467253861933279360152095861484548"),
+ boost::lexical_cast<T>("3273919938390136737194044982676.40271056035622723775417608127544182097346526115858803376474"),
+ boost::lexical_cast<T>("-22989594065095806099337396006399.5874206181563663855129141706748733174902067950115092492439"),
+ boost::lexical_cast<T>("116970582893952893160414263796102.542775878583510989850142808618916073286745084692189044738"),
+ boost::lexical_cast<T>("-459561969036479455224850813196807.283291532483532558959779434457349912822256480548436066098"),
+ boost::lexical_cast<T>("1450959909778264914956547227964788.89797179379520834974601372820249784303794436366366810477"),
+ boost::lexical_cast<T>("-3782846865486775046285288437885921.41537699732805465141128848354901016102326190612528503251"),
+ boost::lexical_cast<T>("8305043213936355459145388670886540.09976337905520168067329932809302445437208115570138102767"),
+ boost::lexical_cast<T>("-15580988484396722546934484726970745.4927787160273626078250810989811865283255762028143642311"),
+ boost::lexical_cast<T>("25262722284076250779006793435537600.0822901485517345545978818780090308947301031347345640449"),
+ boost::lexical_cast<T>("-35714428027687018805443603728757116.5304655170478705341887572982734901197345415291580897698"),
+ boost::lexical_cast<T>("44334726194692443174715432419157343.2294160783772787096321009453791271387235388689346602833"),
+ boost::lexical_cast<T>("-48599573547617297831555162417695106.187829304963846482633791012658974681648157963911491985"),
+ boost::lexical_cast<T>("47258466493366798944386359199482189.0753349196625125615316002614813737880755896979754845101"),
+ boost::lexical_cast<T>("-40913448215392412059728312039233342.142914753896559359297977982314043378636755884088383226"),
+ boost::lexical_cast<T>("31626312914486892948769164616982902.7262756989418188077611392594232674722318027323102462687"),
+ boost::lexical_cast<T>("-21878079174441332123064991795834438.4699982361692990285700077902601657354101259411789722708"),
+ boost::lexical_cast<T>("13567268503974326527361474986354265.3136632133935430378937191911532112778452274286122946396"),
+ boost::lexical_cast<T>("-7551494211746723529747611556474669.62996644923557605747803028485900789337467673523741066527"),
+ boost::lexical_cast<T>("3775516572689476384052312341432597.70584966904950490541958869730702790312581801585742038997"),
+ boost::lexical_cast<T>("-1696271471453637244930364711513292.79902955514107737992185368006225264329876265486853482449"),
+ boost::lexical_cast<T>("684857608019352767999083000986166.20765273693720041519286231015176745354062413008561259139"),
+ boost::lexical_cast<T>("-248397566275708464679881624417990.410438107634139924805871051723444048539177890346227250473"),
+ boost::lexical_cast<T>("80880368999557992138783568858556.1512378233079327986518410244522800950609595592170022878937"),
+ boost::lexical_cast<T>("-23618197945394013802495450485616.9025005749893350650829964098117490779655546610665927669729"),
+ boost::lexical_cast<T>("6176884636893816103087134481332.06708966653493024119556843727320635285468825056891248447124"),
+ boost::lexical_cast<T>("-1444348683723439589948246285262.64080678953468490544615312565485170860503207005915261691108"),
+ boost::lexical_cast<T>("301342031656979076702313946827.961658905182634508217626783081241074250132289461989777865387"),
+ boost::lexical_cast<T>("-55959656587719766738301589651.3940625826610668990368881615587469329021742236397809951765678"),
+ boost::lexical_cast<T>("9223339169004064297247180402.36227016155682738556103138196079389248843082157924368301293963"),
+ boost::lexical_cast<T>("-1344882881571942601385730283.42710150182526891377514071881534880944872423492272147871101373"),
+ boost::lexical_cast<T>("172841913316760599352601139.54409257740173055624405575900164401527761357324625574736896079"),
+ boost::lexical_cast<T>("-19496120443876233531343952.3802212016691702737346568192063937387825469602063310488794471653"),
+ boost::lexical_cast<T>("1920907372583710284097959.44121420322495784420169085871802458519363819782779653621724219067"),
+ boost::lexical_cast<T>("-164429314798240461613359.399597503536657962383155875723527581699785846599051112719962464604"),
+ boost::lexical_cast<T>("12154026644351189572525.1452249886865981747374191977801688548318519692423556934568426042152"),
+ boost::lexical_cast<T>("-770443988366210815096.519382051977221101156336663806705367929328924137169970381042234329058"),
+ boost::lexical_cast<T>("41558909851418707920.4696085656889424895313728719601503526476333404973280596225722152966128"),
+ boost::lexical_cast<T>("-1890879946549708819.24562220042687554209318172044783707920086716716717574156283898330017796"),
+ boost::lexical_cast<T>("71844996557297623.9583461685535340440524052492427928388171299145330229958643439878608673403"),
+ boost::lexical_cast<T>("-2253785109518255.55600197759875781765803818232939130127735487613049577235879610065545755637"),
+ boost::lexical_cast<T>("57616883849355.997562563968344493719962252675875692642406455612671495250543228005045106721"),
+ boost::lexical_cast<T>("-1182495730353.08218118278997948852215670614084013289033222774171548915344541249351599628436"),
+ boost::lexical_cast<T>("19148649358.6196967288062261380599423925174178776792840639099120170800869284300966978300613"),
+ boost::lexical_cast<T>("-239779605.891370259668403359614360511661030470269478602533200704639655585967442891496784613"),
+ boost::lexical_cast<T>("2267583.00284368310957842936892685032434505866445291643236133553754152047677944820353796872"),
+ boost::lexical_cast<T>("-15749.490806784673108773558070497383604733010677027764233749920147549999361110299643477893"),
+ boost::lexical_cast<T>("77.7059495149052727171505425584459982871343274332635726864135949842508025564999785370162956"),
+ boost::lexical_cast<T>("-0.261619987273930331397625130282851629108569607193781378836014468617185550622160348688297247"),
+ boost::lexical_cast<T>("0.000572252321659691600529444769356185993188551770817110673186068921175991328434642504616377475"),
+ boost::lexical_cast<T>("-0.765167220661540041663007112207436426423746402583423562585653954743978584117929356523307954e-6"),
+ boost::lexical_cast<T>("0.579179571056209077507916813937971472839851499147582627425979879366849876944438724610663401e-9"),
+ boost::lexical_cast<T>("-0.224804733043915149719206760378355636826808754704148660354494460792713189958510735070096991e-12"),
+ boost::lexical_cast<T>("0.392711975389579343321746945135488409914483448652884894759297584020979857734289645857714768e-16"),
+ boost::lexical_cast<T>("-0.258603588346412049542768766878162221817684639789901440429511261589010049357907537684380983e-20"),
+ boost::lexical_cast<T>("0.499992460848751668441190360024540741752242879565548017176883304716370989218399797418478685e-25"),
+ boost::lexical_cast<T>("-0.196211614533318174187346267877390498735734213905206562766348625767919122511096089367364025e-30"),
+ boost::lexical_cast<T>("0.874722648949676363732094858062907290148733370978226751462386623191111439121706262759209573e-37"),
+ boost::lexical_cast<T>("-0.163907874717737848669759890242660846846105433791283903651924563157080252845636658802930428e-44"),
+ };
+ T result = d[0];
+ for(int k = 1; k < sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += d[k]/(z+(k-1));
+ }
+ return result;
+ }
+
+ template <class T>
+ static T lanczos_sum_expG_scaled(const T& z)
+ {
+ using namespace boost;
+ static const T d[61] = {
+ boost::lexical_cast<T>("0.901751806425638853077358552989167785490911341809902155556127108480303870921448984935411583e-27"),
+ boost::lexical_cast<T>("4.80241125306810017699523302110401965428995345115391817406006361151407344955277298373661032"),
+ boost::lexical_cast<T>("-108.119283021710869401330097315436214587270846871451487282117128515476478251641970487922552"),
+ boost::lexical_cast<T>("1177.78262074811362219818923738088833932279000985161077740440010901595132448469513438139561"),
+ boost::lexical_cast<T>("-8270.43570321334374279057432173814835581983913167617217749736484999327758232081395983082867"),
+ boost::lexical_cast<T>("42079.807161997077661752306902088979258826568702655699995911391774839958572703348502730394"),
+ boost::lexical_cast<T>("-165326.003834443330215001219988296482004968548294447320869281647211603153902436231468280089"),
+ boost::lexical_cast<T>("521978.361504895300685499370463597042533432134369277742485307843747923127933979566742421213"),
+ boost::lexical_cast<T>("-1360867.51629992863544553419296636395576666570468519805449755596254912681418267100657262281"),
+ boost::lexical_cast<T>("2987713.73338656161102517003716335104823650191612448011720936412226357385029800040631754755"),
+ boost::lexical_cast<T>("-5605212.64915921452169919008770165304171481697939254152852673009005154549681704553438450709"),
+ boost::lexical_cast<T>("9088186.58332916818449459635132673652700922052988327069535513580836143146727832380184335474"),
+ boost::lexical_cast<T>("-12848155.5543636394746355365819800465364996596310467415907815393346205151090486190421959769"),
+ boost::lexical_cast<T>("15949281.2867656960575878805158849857756293807220033618942867694361569866468996967761600898"),
+ boost::lexical_cast<T>("-17483546.9948295433308250581770557182576171673272450149400973735206019559576269174369907171"),
+ boost::lexical_cast<T>("17001087.8599749419660906448951424280111036786456594738278573653160553115681287326064596964"),
+ boost::lexical_cast<T>("-14718487.0643665950346574802384331125115747311674609017568623694516187494204567579595827859"),
+ boost::lexical_cast<T>("11377468.7255609717716845971105161298889777425898291183866813269222281486121330837791392732"),
+ boost::lexical_cast<T>("-7870571.64253038587947746661946939286858490057774448573157856145556080330153403858747655263"),
+ boost::lexical_cast<T>("4880783.08440908743641723492059912671377140680710345996273343885045364205895751515063844239"),
+ boost::lexical_cast<T>("-2716626.7992639625103140035635916455652302249897918893040695025407382316653674141983775542"),
+ boost::lexical_cast<T>("1358230.46602865696544327299659410214201327791319846880787515156343361311278133805428800255"),
+ boost::lexical_cast<T>("-610228.440751458395860905749312275043435828322076830117123636938979942213530882048883969802"),
+ boost::lexical_cast<T>("246375.416501158654327780901087115642493055617468601787093268312234390446439555559050129729"),
+ boost::lexical_cast<T>("-89360.2599028475206119333931211015869139511677735549267100272095432140508089207221096740632"),
+ boost::lexical_cast<T>("29096.4637987498328341260960356772198979319790332957125131055960448759586930781530063775634"),
+ boost::lexical_cast<T>("-8496.57401431514433694413130585404918350686834939156759654375188338156288564260152505382438"),
+ boost::lexical_cast<T>("2222.11523574301594407443285016240908726891841242444092960094015874546135316534057865883047"),
+ boost::lexical_cast<T>("-519.599993280949289705514822058693289933461756514489674453254304308040888101533569480646682"),
+ boost::lexical_cast<T>("108.406868361306987817730701109400305482972790224573776407166683184990131682003417239181112"),
+ boost::lexical_cast<T>("-20.1313142142558596796857948064047373605439974799116521459977609253211918146595346493447238"),
+ boost::lexical_cast<T>("3.31806787671783168020012913552384112429614503798293169239082032849759155847394955909684383"),
+ boost::lexical_cast<T>("-0.483817477111537877685595212919784447924875428848331771524426361483392903320495411973587861"),
+ boost::lexical_cast<T>("0.0621793463102927384924303843912913542297042029136293808338022462765755471002366206711862637"),
+ boost::lexical_cast<T>("-0.00701366932085103924241526535768453911099671087892444015581511551813219720807206445462785293"),
+ boost::lexical_cast<T>("0.000691040514756294308758606917671220770856269030526647010461217455799229645004351524024364997"),
+ boost::lexical_cast<T>("-0.591529398871361458428147660822525365922497109038495896497692806150033516658042357799869656e-4"),
+ boost::lexical_cast<T>("0.437237367535177689875119370170434437030440227275583289093139147244747901678407875809020739e-5"),
+ boost::lexical_cast<T>("-0.277164853397051135996651958345647824709602266382721185838782221179129726200661453504250697e-6"),
+ boost::lexical_cast<T>("0.149506899012035980148891401548317536032574502641368034781671941165064546410613201579653674e-7"),
+ boost::lexical_cast<T>("-0.68023824066463262779882895193964639544061678698791279217407325790147925675797085217462974e-9"),
+ boost::lexical_cast<T>("0.258460163734186329938721529982859244969655253624066115559707985878606277800329299821882688e-10"),
+ boost::lexical_cast<T>("-0.810792256024669306744649981276512583535251727474303382740940985102669076169168931092026491e-12"),
+ boost::lexical_cast<T>("0.207274966207031327521921078048021807442500113231320959236850963529304158700096495799022922e-13"),
+ boost::lexical_cast<T>("-0.425399199286327802950259994834798737777721414442095221716122926637623478450472871269742436e-15"),
+ boost::lexical_cast<T>("0.688866766744198529064607574117697940084548375790020728788313274612845280173338912495478431e-17"),
+ boost::lexical_cast<T>("-0.862599751805643281578607291655858333628582704771553874199104377131082877406179933909898885e-19"),
+ boost::lexical_cast<T>("0.815756005678735657200275584442908437977926312650210429668619446123450972547018343768177988e-21"),
+ boost::lexical_cast<T>("-0.566583084099007858124915716926967268295318152203932871370429534546567151650626184750291695e-23"),
+ boost::lexical_cast<T>("0.279544761599725082805446124351997692260093135930731230328454667675190245860598623539891708e-25"),
+ boost::lexical_cast<T>("-0.941169851584987983984201821679114408126982142904386301937192011680047611188837432096199601e-28"),
+ boost::lexical_cast<T>("0.205866011331040736302780507155525142187875191518455173304638008169488993406425201915370746e-30"),
+ boost::lexical_cast<T>("-0.27526655245712584371295491216289353976964567057707464008951584303679019796521332324352501e-33"),
+ boost::lexical_cast<T>("0.208358067979444305082929004102609366169534624348056112144990933897581971394396210379638792e-36"),
+ boost::lexical_cast<T>("-0.808728107661779323263133007119729988596844663194254976820030366188579170595441991680169012e-40"),
+ boost::lexical_cast<T>("0.141276924383478964519776436955079978012672985961918248012931336621229652792338950573694356e-43"),
+ boost::lexical_cast<T>("-0.930318449287651389310440021745842417218125582685428432576258687100661462527604238849332053e-48"),
+ boost::lexical_cast<T>("0.179870748819321661641184169834635133045197146966203370650788171790610563029431722308057539e-52"),
+ boost::lexical_cast<T>("-0.705865243912790337263229413370018392321238639017433365017168104310561824133229343750737083e-58"),
+ boost::lexical_cast<T>("0.3146787805734405996448268100558028857930560442377698646099945108125281507396722995748398e-64"),
+ boost::lexical_cast<T>("-0.589653534231618730406843260601322236697428143603814281282790370329151249078338470962782338e-72"),
+ };
+ T result = d[0];
+ for(int k = 1; k < sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += d[k]/(z+(k-1));
+ }
+ return result;
+ }
+
+ template<class T>
+ static T lanczos_sum_near_1(const T& dz)
+ {
+ using namespace boost;
+ static const T d[60] = {
+ boost::lexical_cast<T>("23.2463658527729692390378860713647146932236940604550445351214987229819352880524561852919518"),
+ boost::lexical_cast<T>("-523.358012551815715084547614110229469295755088686612838322817729744722233637819564673967396"),
+ boost::lexical_cast<T>("5701.12892340421080714956066268650092612647620400476183901625272640935853188559347587495571"),
+ boost::lexical_cast<T>("-40033.5506451901904954336453419007623117537868026994808919201793803506999271787018654246699"),
+ boost::lexical_cast<T>("203689.884259074923009439144410340256983393397995558814367995938668111650624499963153485034"),
+ boost::lexical_cast<T>("-800270.648969745331278757692597096167418585957703057412758177038340791380011708874081291202"),
+ boost::lexical_cast<T>("2526668.23380061659863999395867315313385499515711742092815402701950519696944982260718031476"),
+ boost::lexical_cast<T>("-6587362.57559198722630391278043503867973853468105110382293763174847657538179665571836023631"),
+ boost::lexical_cast<T>("14462211.3454541602975917764900442754186801975533106565506542322063393991552960595701762805"),
+ boost::lexical_cast<T>("-27132375.1879227404375395522940895789625516798992585980380939378508607160857820002128106898"),
+ boost::lexical_cast<T>("43991923.8735251977856804364757478459275087361742168436524951824945035673768875988985478116"),
+ boost::lexical_cast<T>("-62192284.0030124679010201921841372967696262036115679150017649233887633598058364494608060812"),
+ boost::lexical_cast<T>("77203473.0770033513405070667417251568915937590689150831268228886281254637715669678358204991"),
+ boost::lexical_cast<T>("-84630180.2217173903516348977915150565994784278120192219937728967986198118628659866582594874"),
+ boost::lexical_cast<T>("82294807.2253549409847505891112074804416229757832871133388349982640444405470371147991706317"),
+ boost::lexical_cast<T>("-71245738.2484649177928765605893043553453557808557887270209768163561363857395639001251515788"),
+ boost::lexical_cast<T>("55073334.3180266913441333534260714059077572215147571872597651029894142803987981342430068594"),
+ boost::lexical_cast<T>("-38097984.1648990787690036742690550656061009857688125101191167768314179751258568264424911474"),
+ boost::lexical_cast<T>("23625729.5032184580395130592017474282828236643586203914515183078852982915252442161768527976"),
+ boost::lexical_cast<T>("-13149998.4348054726172055622442157883429575511528431835657668083088902165366619827169829685"),
+ boost::lexical_cast<T>("6574597.77221556423683199818131482663205682902023554728024972161230111356285973623550338976"),
+ boost::lexical_cast<T>("-2953848.1483469149918109110050192571921018042012905892107136410603990336401921230407043408"),
+ boost::lexical_cast<T>("1192595.29584357246380113611351829515963605337523874715861849584306265512819543347806085356"),
+ boost::lexical_cast<T>("-432553.812019608638388918135375154289816441900572658692369491476137741687213006403648722272"),
+ boost::lexical_cast<T>("140843.215385933866391177743292449477205328233960902455943995092958295858485718885800427129"),
+ boost::lexical_cast<T>("-41128.186992679630058614841985110676526199977321524879849001760603476646382839182691529968"),
+ boost::lexical_cast<T>("10756.2849191854701631989789887757784944313743544315113894758328432005767448056040879337769"),
+ boost::lexical_cast<T>("-2515.15559672041299884426826962296210458052543246529646213159169885444118227871246315458787"),
+ boost::lexical_cast<T>("524.750087004805200600237632074908875763734578390662349666321453103782638818305404274166951"),
+ boost::lexical_cast<T>("-97.4468596421732493988298219295878130651986131393383646674144877163795143930682205035917965"),
+ boost::lexical_cast<T>("16.0613108128210806736384551896802799172445298357754547684100294231532127326987175444453058"),
+ boost::lexical_cast<T>("-2.34194813526540240672426202485306862230641838409943369059203455578340880900483887447559874"),
+ boost::lexical_cast<T>("0.300982934748016059399829007219431333744032924923669397318820178988611410275964499475465815"),
+ boost::lexical_cast<T>("-0.033950095985367909789000959795708551814461844488183964432565731809399824963680858993718525"),
+ boost::lexical_cast<T>("0.00334502394288921146242772614150438404658527112198421937945605441697314216921393987758378122"),
+ boost::lexical_cast<T>("-0.000286333429067523984607730553301991502191011265745476190940771685897687956262049750683013485"),
+ boost::lexical_cast<T>("0.211647426149364947402896718485536530479491688838087899435991994237067890628274492042231115e-4"),
+ boost::lexical_cast<T>("-0.134163345121302758109675190598169832775248626443483098532368562186356128620805552609175683e-5"),
+ boost::lexical_cast<T>("0.723697303042715085329782938306424498336642078597508935450663080894255773653328980495047891e-7"),
+ boost::lexical_cast<T>("-0.329273487343139063533251321553223583999676337945788660475231347828772272134656322947906888e-8"),
+ boost::lexical_cast<T>("0.12510922551028971731767784013117088894558604838820475961392154031378891971216135267744134e-9"),
+ boost::lexical_cast<T>("-0.392468958215589939603666430583400537413757786057185505426804034745840192914621891690369903e-11"),
+ boost::lexical_cast<T>("0.100332717101049934370760667782927946803279422001380028508200697081188326364078428184546051e-12"),
+ boost::lexical_cast<T>("-0.205917088291197705194762747639836655808855850989058813560983717575008725393428497910009756e-14"),
+ boost::lexical_cast<T>("0.333450178247893143608439314203175490705783992567107481617660357577257627854979230819461489e-16"),
+ boost::lexical_cast<T>("-0.417546693906616047110563550428133589051498362676394888715581845170969319500638944065594319e-18"),
+ boost::lexical_cast<T>("0.394871691642184410859178529844325390739857256666676534513661579365702353214518478078730801e-20"),
+ boost::lexical_cast<T>("-0.274258012587811199757875927548699264063511843669070634471054184977334027224611843434000922e-22"),
+ boost::lexical_cast<T>("0.135315354265459854889496635967343027244391821142592599244505313738163473730636430399785048e-24"),
+ boost::lexical_cast<T>("-0.455579032003288910408487905303223613382276173706542364543918076752861628464036586507967767e-27"),
+ boost::lexical_cast<T>("0.99650703862462739161520123768147312466695159780582506041370833824093136783202694548427718e-30"),
+ boost::lexical_cast<T>("-0.1332444609228706921659395801935919548447859029572115502899861345555006827214220195650058e-32"),
+ boost::lexical_cast<T>("0.100856999148765307000182397631280249632761913433008379786888200467467364474581430670889392e-35"),
+ boost::lexical_cast<T>("-0.39146979455613683472384690509165395074425354524713697428673406058157887065953366609738731e-39"),
+ boost::lexical_cast<T>("0.683859606707931248105140296850112494069265272540298100341919970496564103098283709868586478e-43"),
+ boost::lexical_cast<T>("-0.450326344248604222735147147805963966503893913752040066400766411031387063854141246972061792e-47"),
+ boost::lexical_cast<T>("0.870675378039492904184581895322153006461319724931909799151743284769901586333730037761678891e-52"),
+ boost::lexical_cast<T>("-0.341678395249272265744518787745356400350877656459401143889000625280131819505857966769964401e-57"),
+ boost::lexical_cast<T>("0.152322191370871666358069530949353871960316638394428595988162174042653299702098929238880862e-63"),
+ boost::lexical_cast<T>("-0.285425405297633795767452984791738825078111150078605004958179057245980222485147999495352632e-71"),
+ };
+ T result = 0;
+ for(int k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(k*dz + k*k);
+ }
+ return result;
+ }
+
+ template<class T>
+ static T lanczos_sum_near_2(const T& dz)
+ {
+ using namespace boost;
+ static const T d[60] = {
+ boost::lexical_cast<T>("557.56438192770795764344217888434355281097193198928944200046501607026919782564033547346298"),
+ boost::lexical_cast<T>("-12552.748616427645475141433405567201788681683808077269330800392600825597799119572762385222"),
+ boost::lexical_cast<T>("136741.650054039199076788077149441364242294724343897779563222338447737802381279007988884806"),
+ boost::lexical_cast<T>("-960205.223613240309942047656967301131022760634321049075674684679438471862998829007639437133"),
+ boost::lexical_cast<T>("4885504.47588736223774859617054275229642041717942140469884121916073195308537421162982679289"),
+ boost::lexical_cast<T>("-19194501.738192166918904824982935279260356575935661514109550613809352009246483412530545583"),
+ boost::lexical_cast<T>("60602169.8633537742937457094837494059849674261357199218329545854990149896822944801504450843"),
+ boost::lexical_cast<T>("-157997975.522506767297528502540724511908584668874487506510120462561270100749019845014382885"),
+ boost::lexical_cast<T>("346876323.86092543685419723290495817330608574729543216092477261152247521712190505658568876"),
+ boost::lexical_cast<T>("-650770365.471136883718747607976242475416651908858429752332176373467422603953536408709972919"),
+ boost::lexical_cast<T>("1055146856.05909309330903130910708372830487315684258450293308627289334336871273080305128138"),
+ boost::lexical_cast<T>("-1491682726.25614447429071368736790697283307005456720192465860871846879804207692411263924608"),
+ boost::lexical_cast<T>("1851726287.94866167094858600116562210167031458934987154557042242638980748286188183300900268"),
+ boost::lexical_cast<T>("-2029855953.68334371445800569238095379629407314338521720558391277508374519526853827142679839"),
+ boost::lexical_cast<T>("1973842002.53354868177824629525448788555435194808657489238517523691040148611221295436287925"),
+ boost::lexical_cast<T>("-1708829941.98209573247426625323314413060108441455314880934710595647408841619484540679859098"),
+ boost::lexical_cast<T>("1320934627.12433688699625456833930317624783222321555050330381730035733198249283009357314954"),
+ boost::lexical_cast<T>("-913780636.858542526294419197161614811332299249415125108737474024007693329922089123296358727"),
+ boost::lexical_cast<T>("566663423.929632170286007468016419798879660054391183200464733820209439185545886930103546787"),
+ boost::lexical_cast<T>("-315402880.436816230388857961460509181823167373029384218959199936902955049832392362044305869"),
+ boost::lexical_cast<T>("157691811.550465734461741500275930418786875005567018699867961482249002625886064187146134966"),
+ boost::lexical_cast<T>("-70848085.5705405970640618473551954585013808128304384354476488268600720054598122945113512731"),
+ boost::lexical_cast<T>("28604413.4050137708444142264980840059788755325900041515286382001704951527733220637586013815"),
+ boost::lexical_cast<T>("-10374808.7067303054787164054055989420809074792801592763124972441820101840292558840131568633"),
+ boost::lexical_cast<T>("3378126.32016207486657791623723515804931231041318964254116390764473281291389374196880720069"),
+ boost::lexical_cast<T>("-986460.090390653140964189383080344920103075349535669020623874668558777188889544398718979744"),
+ boost::lexical_cast<T>("257989.631187387317948158483575125380011593209850756066176921901006833159795100137743395985"),
+ boost::lexical_cast<T>("-60326.0391159227288325790327830741260824763549807922845004854215660451399733578621565837087"),
+ boost::lexical_cast<T>("12586.1375399649496159880821645216260841794563919652590583420570326276086323953958907053394"),
+ boost::lexical_cast<T>("-2337.26417330316922535871922886167444795452055677161063205953141105726549966801856628447293"),
+ boost::lexical_cast<T>("385.230745012608736644117458716226876976056390433401632749144285378123105481506733917763829"),
+ boost::lexical_cast<T>("-56.1716559403731491675970177460841997333796694857076749852739159067307309470690838101179615"),
+ boost::lexical_cast<T>("7.21907953468550196348585224042498727840087634483369357697580053424523903859773769748599575"),
+ boost::lexical_cast<T>("-0.814293485887386870805786409956942772883474224091975496298369877683530509729332902182019049"),
+ boost::lexical_cast<T>("0.0802304419995150047616460464220884371214157889148846405799324851793571580868840034085001373"),
+ boost::lexical_cast<T>("-0.00686771095380619535195996193943858680694970000948742557733102777115482017857981277171196115"),
+ boost::lexical_cast<T>("0.000507636621977556438232617777542864427109623356049335590894564220687567763620803789858345916"),
+ boost::lexical_cast<T>("-0.32179095465362720747836116655088181481893063531138957363431280817392443948706754917605911e-4"),
+ boost::lexical_cast<T>("0.173578890579848508947329833426585354230744194615295570820295052665075101971588563893718407e-5"),
+ boost::lexical_cast<T>("-0.789762880006288893888161070734302768702358633525134582027140158619195373770299678322596835e-7"),
+ boost::lexical_cast<T>("0.300074637200885066788470310738617992259402710843493097610337134266720909870967550606601658e-8"),
+ boost::lexical_cast<T>("-0.941337297619721713151110244242536308266701344868601679868536153775533330272973088246835359e-10"),
+ boost::lexical_cast<T>("0.24064815013182536657310186836079333949814111498828401548170442715552017773994482539471456e-11"),
+ boost::lexical_cast<T>("-0.493892399304811910466345686492277504628763169549384435563232052965821874553923373100791477e-13"),
+ boost::lexical_cast<T>("0.799780678476644196901221989475355609743387528732994566453160178199295104357319723742820593e-15"),
+ boost::lexical_cast<T>("-0.100148627870893347527249092785257443532967736956154251497569191947184705954310833302770086e-16"),
+ boost::lexical_cast<T>("0.947100256812658897084619699699028861352615460106539259289295071616221848196411749449858071e-19"),
+ boost::lexical_cast<T>("-0.657808193528898116367845405906343884364280888644748907819280236995018351085371701094007759e-21"),
+ boost::lexical_cast<T>("0.324554050057463845012469010247790763753999056976705084126950591088538742509983426730851614e-23"),
+ boost::lexical_cast<T>("-0.10927068902162908990029309141242256163212535730975970310918370355165185052827948996110107e-25"),
+ boost::lexical_cast<T>("0.239012340507870646690121104637467232366271566488184795459093215303237974655782634371712486e-28"),
+ boost::lexical_cast<T>("-0.31958700972990573259359660326375143524864710953063781494908314884519046349402409667329667e-31"),
+ boost::lexical_cast<T>("0.241905641292988284384362036555782113275737930713192053073501265726048089991747342247551645e-34"),
+ boost::lexical_cast<T>("-0.93894080230619233745797029179332447129464915420290457418429337322820997038069119047864035e-38"),
+ boost::lexical_cast<T>("0.164023814025085488413251990798690797467351995518990067783355251949198292596815470576539877e-41"),
+ boost::lexical_cast<T>("-0.108010831192689925518484618970761942019888832176355541674171850211917230280206410356465451e-45"),
+ boost::lexical_cast<T>("0.208831600642796805563854019033577205240227465154130766898180386564934443551840379116390645e-50"),
+ boost::lexical_cast<T>("-0.819516067465171848863933747691434138146789031214932473898084756489529673230665363014007306e-56"),
+ boost::lexical_cast<T>("0.365344970579318347488211604761724311582675708113250505307342682118101409913523622073678179e-62"),
+ boost::lexical_cast<T>("-0.684593199208628857931267904308244537968349564351534581234005234847904343404822808648361291e-70"),
+ };
+ T result = 0;
+ T z = dz + 2;
+ for(int k = 1; k <= sizeof(d)/sizeof(d[0]); ++k)
+ {
+ result += (-d[k-1]*dz)/(z + k*z + k*k - 1);
+ }
+ return result;
+ }
+
+ static double g(){ return 63.19215200000000010049916454590857028961181640625; }
+};
+
+}}} // namespaces
+
+#endif
+
+
Added: sandbox/math_toolkit/boost/math/bindings/mpfr.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/bindings/mpfr.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,850 @@
+// Copyright John Maddock 2008.
+// 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)
+//
+// Wrapper that works with mpfr_class defined in gmpfrxx.h
+// See http://math.berkeley.edu/~wilken/code/gmpfrxx/
+// Also requires the gmp and mpfr libraries.
+//
+
+#ifndef BOOST_MATH_MPLFR_BINDINGS_HPP
+#define BOOST_MATH_MPLFR_BINDINGS_HPP
+
+#include <gmpfrxx.h>
+#include <boost/math/tools/precision.hpp>
+#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/bindings/detail/big_digamma.hpp>
+#include <boost/math/bindings/detail/big_lanczos.hpp>
+
+inline mpfr_class fabs(const mpfr_class& v)
+{
+ return abs(v);
+}
+
+inline mpfr_class pow(const mpfr_class& b, const mpfr_class e)
+{
+ mpfr_class result;
+ mpfr_pow(result.__get_mp(), b.__get_mp(), e.__get_mp(), GMP_RNDN);
+ return result;
+}
+
+inline mpfr_class ldexp(const mpfr_class& v, int e)
+{
+ //int e = mpfr_get_exp(*v.__get_mp());
+ mpfr_class result(v);
+ mpfr_set_exp(result.__get_mp(), e);
+ return result;
+}
+
+inline mpfr_class frexp(const mpfr_class& v, int* expon)
+{
+ int e = mpfr_get_exp(v.__get_mp());
+ mpfr_class result(v);
+ mpfr_set_exp(result.__get_mp(), 0);
+ *expon = e;
+ return result;
+}
+
+mpfr_class fmod(const mpfr_class& v1, const mpfr_class& v2)
+{
+ mpfr_class n;
+ if(v1 < 0)
+ n = ceil(v1 / v2);
+ else
+ n = floor(v1 / v2);
+ return v1 - n * v2;
+}
+
+template <class Policy>
+inline mpfr_class modf(const mpfr_class& v, long long* ipart, const Policy& pol)
+{
+ *ipart = lltrunc(v, pol);
+ return v - boost::math::tools::real_cast<mpfr_class>(*ipart);
+}
+template <class Policy>
+inline int iround(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<int>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline long lround(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline long long llround(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long long>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline int itrunc(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<int>(boost::math::trunc(x, pol));
+}
+
+template <class Policy>
+inline long ltrunc(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long>(boost::math::trunc(x, pol));
+}
+
+template <class Policy>
+inline long long lltrunc(mpfr_class const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long long>(boost::math::trunc(x, pol));
+}
+
+namespace boost{ namespace math{
+
+#if defined(__GNUC__) && (__GNUC__ < 4)
+ using ::iround;
+ using ::lround;
+ using ::llround;
+ using ::itrunc;
+ using ::ltrunc;
+ using ::lltrunc;
+ using ::modf;
+#endif
+
+namespace lanczos{
+
+struct mpfr_lanczos
+{
+ static mpfr_class lanczos_sum(const mpfr_class& z)
+ {
+ unsigned long p = z.get_dprec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum(z);
+ }
+ static mpfr_class lanczos_sum_expG_scaled(const mpfr_class& z)
+ {
+ unsigned long p = z.get_dprec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_expG_scaled(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_expG_scaled(z);
+ }
+ static mpfr_class lanczos_sum_near_1(const mpfr_class& z)
+ {
+ unsigned long p = z.get_dprec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_1(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_1(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_1(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_1(z);
+ }
+ static mpfr_class lanczos_sum_near_2(const mpfr_class& z)
+ {
+ unsigned long p = z.get_dprec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_2(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_2(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_2(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_2(z);
+ }
+ static mpfr_class g()
+ {
+ unsigned long p = mpfr_class::get_dprec();
+ if(p <= 72)
+ return lanczos13UDT::g();
+ else if(p <= 120)
+ return lanczos22UDT::g();
+ else if(p <= 170)
+ return lanczos31UDT::g();
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::g();
+ }
+};
+
+template<class Policy>
+struct lanczos<mpfr_class, Policy>
+{
+ typedef mpfr_lanczos type;
+};
+
+} // namespace lanczos
+
+namespace tools
+{
+
+template <class T, class U>
+struct promote_arg<__gmp_expr<T,U> >
+{ // If T is integral type, then promote to double.
+ typedef mpfr_class type;
+};
+
+template<>
+inline int digits<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ return mpfr_class::get_dprec();
+}
+
+namespace detail{
+
+template<class I>
+void convert_to_long_result(mpfr_class const& r, I& result)
+{
+ result = 0;
+ I last_result(0);
+ mpfr_class t(r);
+ double term;
+ do
+ {
+ term = real_cast<double>(t);
+ last_result = result;
+ result += static_cast<I>(term);
+ t -= term;
+ }while(result != last_result);
+}
+
+}
+
+template <>
+inline mpfr_class real_cast<mpfr_class, long long>(long long t)
+{
+ mpfr_class result;
+ int expon = 0;
+ int sign = 1;
+ if(t < 0)
+ {
+ sign = -1;
+ t = -t;
+ }
+ while(t)
+ {
+ result += ldexp((double)(t & 0xffffL), expon);
+ expon += 32;
+ t >>= 32;
+ }
+ return result * sign;
+}
+template <>
+inline unsigned real_cast<unsigned, mpfr_class>(mpfr_class t)
+{
+ return t.get_ui();
+}
+template <>
+inline int real_cast<int, mpfr_class>(mpfr_class t)
+{
+ return t.get_si();
+}
+template <>
+inline double real_cast<double, mpfr_class>(mpfr_class t)
+{
+ return t.get_d();
+}
+template <>
+inline float real_cast<float, mpfr_class>(mpfr_class t)
+{
+ return static_cast<float>(t.get_d());
+}
+template <>
+inline long real_cast<long, mpfr_class>(mpfr_class t)
+{
+ long result;
+ detail::convert_to_long_result(t, result);
+ return result;
+}
+template <>
+inline long long real_cast<long long, mpfr_class>(mpfr_class t)
+{
+ long long result;
+ detail::convert_to_long_result(t, result);
+ return result;
+}
+
+template <>
+inline mpfr_class max_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ static bool has_init = false;
+ static mpfr_class val;
+ if(!has_init)
+ {
+ val = 0.5;
+ mpfr_set_exp(val.__get_mp(), mpfr_get_emax());
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr_class min_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ static bool has_init = false;
+ static mpfr_class val;
+ if(!has_init)
+ {
+ val = 0.5;
+ mpfr_set_exp(val.__get_mp(), mpfr_get_emin());
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr_class log_max_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ static bool has_init = false;
+ static mpfr_class val = max_value<mpfr_class>();
+ if(!has_init)
+ {
+ val = log(val);
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr_class log_min_value<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ static bool has_init = false;
+ static mpfr_class val = max_value<mpfr_class>();
+ if(!has_init)
+ {
+ val = log(val);
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr_class epsilon<mpfr_class>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr_class))
+{
+ return ldexp(mpfr_class(1), 1-boost::math::policies::digits<mpfr_class, boost::math::policies::policy<> >());
+}
+
+} // namespace tools
+
+namespace policies{
+
+template <class T, class U, class Policy>
+struct evaluation<__gmp_expr<T, U>, Policy>
+{
+ typedef mpfr_class type;
+};
+
+}
+
+template <class Policy>
+inline mpfr_class skewness(const extreme_value_distribution<mpfr_class, Policy>& /*dist*/)
+{
+ //
+ // This is 12 * sqrt(6) * zeta(3) / pi^3:
+ // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
+ //
+ return boost::lexical_cast<mpfr_class>("1.1395470994046486574927930193898461120875997958366");
+}
+
+template <class Policy>
+inline mpfr_class skewness(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)
+{
+ // using namespace boost::math::constants;
+ return boost::lexical_cast<mpfr_class>("0.63111065781893713819189935154422777984404221106391");
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
+}
+
+template <class Policy>
+inline mpfr_class kurtosis(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)
+{
+ // using namespace boost::math::constants;
+ return boost::lexical_cast<mpfr_class>("3.2450893006876380628486604106197544154170667057995");
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
+ // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
+}
+
+template <class Policy>
+inline mpfr_class kurtosis_excess(const rayleigh_distribution<mpfr_class, Policy>& /*dist*/)
+{
+ //using namespace boost::math::constants;
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ return boost::lexical_cast<mpfr_class>("0.2450893006876380628486604106197544154170667057995");
+ // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
+ // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
+} // kurtosis
+
+namespace detail{
+
+//
+// Version of Digamma accurate to ~100 decimal digits.
+//
+template <class Policy>
+mpfr_class digamma_imp(mpfr_class x, const mpl::int_<0>* , const Policy& pol)
+{
+ //
+ // This handles reflection of negative arguments, and all our
+ // empfr_classor handling, then forwards to the T-specific approximation.
+ //
+ BOOST_MATH_STD_USING // ADL of std functions.
+
+ mpfr_class result = 0;
+ //
+ // Check for negative arguments and use reflection:
+ //
+ if(x < 0)
+ {
+ // Reflect:
+ x = 1 - x;
+ // Argument reduction for tan:
+ mpfr_class remainder = x - floor(x);
+ // Shift to negative if > 0.5:
+ if(remainder > 0.5)
+ {
+ remainder -= 1;
+ }
+ //
+ // check for evaluation at a negative pole:
+ //
+ if(remainder == 0)
+ {
+ return policies::raise_pole_error<mpfr_class>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
+ }
+ result = constants::pi<mpfr_class>() / tan(constants::pi<mpfr_class>() * remainder);
+ }
+ result += big_digamma(x);
+ return result;
+}
+//
+// Specialisations of this function provides the initial
+// starting guess for Halley iteration:
+//
+template <class Policy>
+mpfr_class erf_inv_imp(const mpfr_class& p, const mpfr_class& q, const Policy&, const boost::mpl::int_<64>*)
+{
+ BOOST_MATH_STD_USING // for ADL of std names.
+
+ mpfr_class result = 0;
+
+ if(p <= 0.5)
+ {
+ //
+ // Evaluate inverse erf using the rational approximation:
+ //
+ // x = p(p+10)(Y+R(p))
+ //
+ // Where Y is a constant, and R(p) is optimised for a low
+ // absolute empfr_classor compared to |Y|.
+ //
+ // double: Max empfr_classor found: 2.001849e-18
+ // long double: Max empfr_classor found: 1.017064e-20
+ // Maximum Deviation Found (actual empfr_classor term at infinite precision) 8.030e-21
+ //
+ static const float Y = 0.0891314744949340820313f;
+ static const mpfr_class P[] = {
+ -0.000508781949658280665617,
+ -0.00836874819741736770379,
+ 0.0334806625409744615033,
+ -0.0126926147662974029034,
+ -0.0365637971411762664006,
+ 0.0219878681111168899165,
+ 0.00822687874676915743155,
+ -0.00538772965071242932965
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ -0.970005043303290640362,
+ -1.56574558234175846809,
+ 1.56221558398423026363,
+ 0.662328840472002992063,
+ -0.71228902341542847553,
+ -0.0527396382340099713954,
+ 0.0795283687341571680018,
+ -0.00233393759374190016776,
+ 0.000886216390456424707504
+ };
+ mpfr_class g = p * (p + 10);
+ mpfr_class r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
+ result = g * Y + g * r;
+ }
+ else if(q >= 0.25)
+ {
+ //
+ // Rational approximation for 0.5 > q >= 0.25
+ //
+ // x = sqrt(-2*log(q)) / (Y + R(q))
+ //
+ // Where Y is a constant, and R(q) is optimised for a low
+ // absolute empfr_classor compared to Y.
+ //
+ // double : Max empfr_classor found: 7.403372e-17
+ // long double : Max empfr_classor found: 6.084616e-20
+ // Maximum Deviation Found (empfr_classor term) 4.811e-20
+ //
+ static const float Y = 2.249481201171875f;
+ static const mpfr_class P[] = {
+ -0.202433508355938759655,
+ 0.105264680699391713268,
+ 8.37050328343119927838,
+ 17.6447298408374015486,
+ -18.8510648058714251895,
+ -44.6382324441786960818,
+ 17.445385985570866523,
+ 21.1294655448340526258,
+ -3.67192254707729348546
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 6.24264124854247537712,
+ 3.9713437953343869095,
+ -28.6608180499800029974,
+ -20.1432634680485188801,
+ 48.5609213108739935468,
+ 10.8268667355460159008,
+ -22.6436933413139721736,
+ 1.72114765761200282724
+ };
+ mpfr_class g = sqrt(-2 * log(q));
+ mpfr_class xs = q - 0.25;
+ mpfr_class r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = g / (Y + r);
+ }
+ else
+ {
+ //
+ // For q < 0.25 we have a series of rational approximations all
+ // of the general form:
+ //
+ // let: x = sqrt(-log(q))
+ //
+ // Then the result is given by:
+ //
+ // x(Y+R(x-B))
+ //
+ // where Y is a constant, B is the lowest value of x for which
+ // the approximation is valid, and R(x-B) is optimised for a low
+ // absolute empfr_classor compared to Y.
+ //
+ // Note that almost all code will really go through the first
+ // or maybe second approximation. After than we're dealing with very
+ // small input values indeed: 80 and 128 bit long double's go all the
+ // way down to ~ 1e-5000 so the "tail" is rather long...
+ //
+ mpfr_class x = sqrt(-log(q));
+ if(x < 3)
+ {
+ // Max empfr_classor found: 1.089051e-20
+ static const float Y = 0.807220458984375f;
+ static const mpfr_class P[] = {
+ -0.131102781679951906451,
+ -0.163794047193317060787,
+ 0.117030156341995252019,
+ 0.387079738972604337464,
+ 0.337785538912035898924,
+ 0.142869534408157156766,
+ 0.0290157910005329060432,
+ 0.00214558995388805277169,
+ -0.679465575181126350155e-6,
+ 0.285225331782217055858e-7,
+ -0.681149956853776992068e-9
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 3.46625407242567245975,
+ 5.38168345707006855425,
+ 4.77846592945843778382,
+ 2.59301921623620271374,
+ 0.848854343457902036425,
+ 0.152264338295331783612,
+ 0.01105924229346489121
+ };
+ mpfr_class xs = x - 1.125;
+ mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 6)
+ {
+ // Max empfr_classor found: 8.389174e-21
+ static const float Y = 0.93995571136474609375f;
+ static const mpfr_class P[] = {
+ -0.0350353787183177984712,
+ -0.00222426529213447927281,
+ 0.0185573306514231072324,
+ 0.00950804701325919603619,
+ 0.00187123492819559223345,
+ 0.000157544617424960554631,
+ 0.460469890584317994083e-5,
+ -0.230404776911882601748e-9,
+ 0.266339227425782031962e-11
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 1.3653349817554063097,
+ 0.762059164553623404043,
+ 0.220091105764131249824,
+ 0.0341589143670947727934,
+ 0.00263861676657015992959,
+ 0.764675292302794483503e-4
+ };
+ mpfr_class xs = x - 3;
+ mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 18)
+ {
+ // Max empfr_classor found: 1.481312e-19
+ static const float Y = 0.98362827301025390625f;
+ static const mpfr_class P[] = {
+ -0.0167431005076633737133,
+ -0.00112951438745580278863,
+ 0.00105628862152492910091,
+ 0.000209386317487588078668,
+ 0.149624783758342370182e-4,
+ 0.449696789927706453732e-6,
+ 0.462596163522878599135e-8,
+ -0.281128735628831791805e-13,
+ 0.99055709973310326855e-16
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 0.591429344886417493481,
+ 0.138151865749083321638,
+ 0.0160746087093676504695,
+ 0.000964011807005165528527,
+ 0.275335474764726041141e-4,
+ 0.282243172016108031869e-6
+ };
+ mpfr_class xs = x - 6;
+ mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 44)
+ {
+ // Max empfr_classor found: 5.697761e-20
+ static const float Y = 0.99714565277099609375f;
+ static const mpfr_class P[] = {
+ -0.0024978212791898131227,
+ -0.779190719229053954292e-5,
+ 0.254723037413027451751e-4,
+ 0.162397777342510920873e-5,
+ 0.396341011304801168516e-7,
+ 0.411632831190944208473e-9,
+ 0.145596286718675035587e-11,
+ -0.116765012397184275695e-17
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 0.207123112214422517181,
+ 0.0169410838120975906478,
+ 0.000690538265622684595676,
+ 0.145007359818232637924e-4,
+ 0.144437756628144157666e-6,
+ 0.509761276599778486139e-9
+ };
+ mpfr_class xs = x - 18;
+ mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else
+ {
+ // Max empfr_classor found: 1.279746e-20
+ static const float Y = 0.99941349029541015625f;
+ static const mpfr_class P[] = {
+ -0.000539042911019078575891,
+ -0.28398759004727721098e-6,
+ 0.899465114892291446442e-6,
+ 0.229345859265920864296e-7,
+ 0.225561444863500149219e-9,
+ 0.947846627503022684216e-12,
+ 0.135880130108924861008e-14,
+ -0.348890393399948882918e-21
+ };
+ static const mpfr_class Q[] = {
+ 1,
+ 0.0845746234001899436914,
+ 0.00282092984726264681981,
+ 0.468292921940894236786e-4,
+ 0.399968812193862100054e-6,
+ 0.161809290887904476097e-8,
+ 0.231558608310259605225e-11
+ };
+ mpfr_class xs = x - 44;
+ mpfr_class R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ }
+ return result;
+}
+
+mpfr_class bessel_i0(mpfr_class x)
+{
+ static const mpfr_class P1[] = {
+ boost::lexical_cast<mpfr_class>("-2.2335582639474375249e+15"),
+ boost::lexical_cast<mpfr_class>("-5.5050369673018427753e+14"),
+ boost::lexical_cast<mpfr_class>("-3.2940087627407749166e+13"),
+ boost::lexical_cast<mpfr_class>("-8.4925101247114157499e+11"),
+ boost::lexical_cast<mpfr_class>("-1.1912746104985237192e+10"),
+ boost::lexical_cast<mpfr_class>("-1.0313066708737980747e+08"),
+ boost::lexical_cast<mpfr_class>("-5.9545626019847898221e+05"),
+ boost::lexical_cast<mpfr_class>("-2.4125195876041896775e+03"),
+ boost::lexical_cast<mpfr_class>("-7.0935347449210549190e+00"),
+ boost::lexical_cast<mpfr_class>("-1.5453977791786851041e-02"),
+ boost::lexical_cast<mpfr_class>("-2.5172644670688975051e-05"),
+ boost::lexical_cast<mpfr_class>("-3.0517226450451067446e-08"),
+ boost::lexical_cast<mpfr_class>("-2.6843448573468483278e-11"),
+ boost::lexical_cast<mpfr_class>("-1.5982226675653184646e-14"),
+ boost::lexical_cast<mpfr_class>("-5.2487866627945699800e-18"),
+ };
+ static const mpfr_class Q1[] = {
+ boost::lexical_cast<mpfr_class>("-2.2335582639474375245e+15"),
+ boost::lexical_cast<mpfr_class>("7.8858692566751002988e+12"),
+ boost::lexical_cast<mpfr_class>("-1.2207067397808979846e+10"),
+ boost::lexical_cast<mpfr_class>("1.0377081058062166144e+07"),
+ boost::lexical_cast<mpfr_class>("-4.8527560179962773045e+03"),
+ boost::lexical_cast<mpfr_class>("1.0"),
+ };
+ static const mpfr_class P2[] = {
+ boost::lexical_cast<mpfr_class>("-2.2210262233306573296e-04"),
+ boost::lexical_cast<mpfr_class>("1.3067392038106924055e-02"),
+ boost::lexical_cast<mpfr_class>("-4.4700805721174453923e-01"),
+ boost::lexical_cast<mpfr_class>("5.5674518371240761397e+00"),
+ boost::lexical_cast<mpfr_class>("-2.3517945679239481621e+01"),
+ boost::lexical_cast<mpfr_class>("3.1611322818701131207e+01"),
+ boost::lexical_cast<mpfr_class>("-9.6090021968656180000e+00"),
+ };
+ static const mpfr_class Q2[] = {
+ boost::lexical_cast<mpfr_class>("-5.5194330231005480228e-04"),
+ boost::lexical_cast<mpfr_class>("3.2547697594819615062e-02"),
+ boost::lexical_cast<mpfr_class>("-1.1151759188741312645e+00"),
+ boost::lexical_cast<mpfr_class>("1.3982595353892851542e+01"),
+ boost::lexical_cast<mpfr_class>("-6.0228002066743340583e+01"),
+ boost::lexical_cast<mpfr_class>("8.5539563258012929600e+01"),
+ boost::lexical_cast<mpfr_class>("-3.1446690275135491500e+01"),
+ boost::lexical_cast<mpfr_class>("1.0"),
+ };
+ mpfr_class value, factor, r;
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+
+ if (x < 0)
+ {
+ x = -x; // even function
+ }
+ if (x == 0)
+ {
+ return static_cast<mpfr_class>(1);
+ }
+ if (x <= 15) // x in (0, 15]
+ {
+ mpfr_class y = x * x;
+ value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
+ }
+ else // x in (15, \infty)
+ {
+ mpfr_class y = 1 / x - 1 / 15;
+ r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
+ factor = exp(x) / sqrt(x);
+ value = factor * r;
+ }
+
+ return value;
+}
+
+mpfr_class bessel_i1(mpfr_class x)
+{
+ static const mpfr_class P1[] = {
+ static_cast<mpfr_class>("-1.4577180278143463643e+15"),
+ static_cast<mpfr_class>("-1.7732037840791591320e+14"),
+ static_cast<mpfr_class>("-6.9876779648010090070e+12"),
+ static_cast<mpfr_class>("-1.3357437682275493024e+11"),
+ static_cast<mpfr_class>("-1.4828267606612366099e+09"),
+ static_cast<mpfr_class>("-1.0588550724769347106e+07"),
+ static_cast<mpfr_class>("-5.1894091982308017540e+04"),
+ static_cast<mpfr_class>("-1.8225946631657315931e+02"),
+ static_cast<mpfr_class>("-4.7207090827310162436e-01"),
+ static_cast<mpfr_class>("-9.1746443287817501309e-04"),
+ static_cast<mpfr_class>("-1.3466829827635152875e-06"),
+ static_cast<mpfr_class>("-1.4831904935994647675e-09"),
+ static_cast<mpfr_class>("-1.1928788903603238754e-12"),
+ static_cast<mpfr_class>("-6.5245515583151902910e-16"),
+ static_cast<mpfr_class>("-1.9705291802535139930e-19"),
+ };
+ static const mpfr_class Q1[] = {
+ static_cast<mpfr_class>("-2.9154360556286927285e+15"),
+ static_cast<mpfr_class>("9.7887501377547640438e+12"),
+ static_cast<mpfr_class>("-1.4386907088588283434e+10"),
+ static_cast<mpfr_class>("1.1594225856856884006e+07"),
+ static_cast<mpfr_class>("-5.1326864679904189920e+03"),
+ static_cast<mpfr_class>("1.0"),
+ };
+ static const mpfr_class P2[] = {
+ static_cast<mpfr_class>("1.4582087408985668208e-05"),
+ static_cast<mpfr_class>("-8.9359825138577646443e-04"),
+ static_cast<mpfr_class>("2.9204895411257790122e-02"),
+ static_cast<mpfr_class>("-3.4198728018058047439e-01"),
+ static_cast<mpfr_class>("1.3960118277609544334e+00"),
+ static_cast<mpfr_class>("-1.9746376087200685843e+00"),
+ static_cast<mpfr_class>("8.5591872901933459000e-01"),
+ static_cast<mpfr_class>("-6.0437159056137599999e-02"),
+ };
+ static const mpfr_class Q2[] = {
+ static_cast<mpfr_class>("3.7510433111922824643e-05"),
+ static_cast<mpfr_class>("-2.2835624489492512649e-03"),
+ static_cast<mpfr_class>("7.4212010813186530069e-02"),
+ static_cast<mpfr_class>("-8.5017476463217924408e-01"),
+ static_cast<mpfr_class>("3.2593714889036996297e+00"),
+ static_cast<mpfr_class>("-3.8806586721556593450e+00"),
+ static_cast<mpfr_class>("1.0"),
+ };
+ mpfr_class value, factor, r, w;
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+
+ w = abs(x);
+ if (x == 0)
+ {
+ return static_cast<mpfr_class>(0);
+ }
+ if (w <= 15) // w in (0, 15]
+ {
+ mpfr_class y = x * x;
+ r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
+ factor = w;
+ value = factor * r;
+ }
+ else // w in (15, \infty)
+ {
+ mpfr_class y = 1 / w - mpfr_class(1) / 15;
+ r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
+ factor = exp(w) / sqrt(w);
+ value = factor * r;
+ }
+
+ if (x < 0)
+ {
+ value *= -value; // odd function
+ }
+ return value;
+}
+
+} // namespace detail
+
+}}
+
+#endif // BOOST_MATH_MPLFR_BINDINGS_HPP
+
Modified: sandbox/math_toolkit/boost/math/bindings/rr.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/bindings/rr.hpp (original)
+++ sandbox/math_toolkit/boost/math/bindings/rr.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -3,6 +3,8 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#ifndef BOOST_MATH_NTL_RR_HPP
+#define BOOST_MATH_NTL_RR_HPP
#include <boost/config.hpp>
#include <boost/limits.hpp>
@@ -11,15 +13,14 @@
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/bindings/detail/big_digamma.hpp>
+#include <boost/math/bindings/detail/big_lanczos.hpp>
#include <ostream>
#include <istream>
#include <boost/config/no_tr1/cmath.hpp>
#include <NTL/RR.h>
-#ifndef BOOST_MATH_NTL_RR_HPP
-#define BOOST_MATH_NTL_RR_HPP
-
namespace boost{ namespace math{
namespace ntl
@@ -421,6 +422,80 @@
} // namespace ntl
+namespace lanczos{
+
+struct ntl_lanczos
+{
+ static ntl::RR lanczos_sum(const ntl::RR& z)
+ {
+ unsigned long p = ntl::RR::precision();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum(z);
+ }
+ static ntl::RR lanczos_sum_expG_scaled(const ntl::RR& z)
+ {
+ unsigned long p = ntl::RR::precision();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_expG_scaled(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_expG_scaled(z);
+ }
+ static ntl::RR lanczos_sum_near_1(const ntl::RR& z)
+ {
+ unsigned long p = ntl::RR::precision();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_1(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_1(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_1(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_1(z);
+ }
+ static ntl::RR lanczos_sum_near_2(const ntl::RR& z)
+ {
+ unsigned long p = ntl::RR::precision();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_2(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_2(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_2(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_2(z);
+ }
+ static ntl::RR g()
+ {
+ unsigned long p = ntl::RR::precision();
+ if(p <= 72)
+ return lanczos13UDT::g();
+ else if(p <= 120)
+ return lanczos22UDT::g();
+ else if(p <= 170)
+ return lanczos31UDT::g();
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::g();
+ }
+};
+
+template<class Policy>
+struct lanczos<ntl::RR, Policy>
+{
+ typedef ntl_lanczos type;
+};
+
+} // namespace lanczos
+
namespace tools
{
@@ -487,7 +562,21 @@
template <>
inline int real_cast<int, boost::math::ntl::RR>(boost::math::ntl::RR t)
{
- unsigned result;
+ int result;
+ detail::convert_to_long_result(t.value(), result);
+ return result;
+}
+template <>
+inline long real_cast<long, boost::math::ntl::RR>(boost::math::ntl::RR t)
+{
+ long result;
+ detail::convert_to_long_result(t.value(), result);
+ return result;
+}
+template <>
+inline long long real_cast<long long, boost::math::ntl::RR>(boost::math::ntl::RR t)
+{
+ long long result;
detail::convert_to_long_result(t.value(), result);
return result;
}
@@ -704,13 +793,78 @@
}
template <class Policy>
+ inline long lround(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<long>(round(x, pol));
+ }
+
+ template <class Policy>
+ inline long long llround(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<long long>(round(x, pol));
+ }
+
+ template <class Policy>
inline int itrunc(RR const& x, const Policy& pol)
{
return tools::real_cast<int>(trunc(x, pol));
}
+ template <class Policy>
+ inline long ltrunc(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<long>(trunc(x, pol));
+ }
+
+ template <class Policy>
+ inline long long lltrunc(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<long long>(trunc(x, pol));
+ }
+
} // namespace ntl
+namespace detail{
+
+template <class Policy>
+ntl::RR digamma_imp(ntl::RR x, const mpl::int_<0>* t, const Policy& pol)
+{
+ //
+ // This handles reflection of negative arguments, and all our
+ // error handling, then forwards to the T-specific approximation.
+ //
+ BOOST_MATH_STD_USING // ADL of std functions.
+
+ ntl::RR result = 0;
+ //
+ // Check for negative arguments and use reflection:
+ //
+ if(x < 0)
+ {
+ // Reflect:
+ x = 1 - x;
+ // Argument reduction for tan:
+ ntl::RR remainder = x - floor(x);
+ // Shift to negative if > 0.5:
+ if(remainder > 0.5)
+ {
+ remainder -= 1;
+ }
+ //
+ // check for evaluation at a negative pole:
+ //
+ if(remainder == 0)
+ {
+ return policies::raise_pole_error<ntl::RR>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
+ }
+ result = constants::pi<ntl::RR>() / tan(constants::pi<ntl::RR>() * remainder);
+ }
+ result += big_digamma(x);
+ return result;
+}
+
+} // namespace detail
+
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/boost/math/concepts/distributions.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/distributions.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/distributions.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -162,6 +162,7 @@
v = quantile(complement(dist, d));
v = hazard(dist, d);
v = chf(dist, d);
+#ifndef TEST_MPFR
long double ld = 1;
v = cdf(dist, ld);
v = cdf(complement(dist, ld));
@@ -170,6 +171,7 @@
v = quantile(complement(dist, ld));
v = hazard(dist, ld);
v = chf(dist, ld);
+#endif
int i = 1;
v = cdf(dist, i);
v = cdf(complement(dist, i));
Added: sandbox/math_toolkit/boost/math/concepts/real_type_concept.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/concepts/real_type_concept.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,109 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+#ifndef BOOST_MATH_REAL_TYPE_CONCEPT_HPP
+#define BOOST_MATH_REAL_TYPE_CONCEPT_HPP
+
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable: 4100)
+#pragma warning(disable: 4510)
+#pragma warning(disable: 4610)
+#endif
+#include <boost/concept_check.hpp>
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+#include <boost/math/tools/config.hpp>
+#include <boost/math/tools/precision.hpp>
+
+
+namespace boost{ namespace math{ namespace concepts{
+
+template <class RealType>
+struct RealTypeConcept
+{
+ template <class Other>
+ void check_binary_ops(Other o)
+ {
+ RealType r(o);
+ r = o;
+ r -= o;
+ r += o;
+ r *= o;
+ r /= o;
+ r = r - o;
+ r = o - r;
+ r = r + o;
+ r = o + r;
+ r = o * r;
+ r = r * o;
+ r = r / o;
+ r = o / r;
+ bool b;
+ b = r == o;
+ b = o == r;
+ b = r != o;
+ b = o != r;
+ b = r <= o;
+ b = o <= r;
+ b = r >= o;
+ b = o >= r;
+ b = r < o;
+ b = o < r;
+ b = r > o;
+ b = o > r;
+ }
+
+ void constraints()
+ {
+ BOOST_MATH_STD_USING
+
+ RealType r;
+ check_binary_ops(r);
+ check_binary_ops(0.5f);
+ check_binary_ops(0.5);
+ //check_binary_ops(0.5L);
+ check_binary_ops(1);
+ //check_binary_ops(1u);
+ check_binary_ops(1L);
+ //check_binary_ops(1uL);
+#ifndef BOOST_HAS_LONG_LONG
+ check_binary_ops(1LL);
+#endif
+ RealType r2 = +r;
+ r2 = -r;
+
+ r2 = fabs(r);
+ r2 = abs(r);
+ r2 = ceil(r);
+ r2 = floor(r);
+ r2 = exp(r);
+ r2 = pow(r, r2);
+ r2 = sqrt(r);
+ r2 = log(r);
+ r2 = cos(r);
+ r2 = sin(r);
+ r2 = tan(r);
+ r2 = asin(r);
+ r2 = acos(r);
+ r2 = atan(r);
+ int i;
+ r2 = ldexp(r, i);
+ r2 = frexp(r, &i);
+ i = boost::math::tools::digits<RealType>();
+ r2 = boost::math::tools::max_value<RealType>();
+ r2 = boost::math::tools::min_value<RealType>();
+ r2 = boost::math::tools::log_max_value<RealType>();
+ r2 = boost::math::tools::log_min_value<RealType>();
+ r2 = boost::math::tools::epsilon<RealType>();
+ }
+}; // struct DistributionConcept
+
+
+}}} // namespaces
+
+#endif
+
Modified: sandbox/math_toolkit/boost/math/distributions/binomial.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/binomial.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/binomial.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -649,13 +649,13 @@
template <class RealType, class Policy>
inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
{
- return binomial_detail::quantile_imp(dist, p, 1-p);
+ return binomial_detail::quantile_imp(dist, p, RealType(1-p));
} // quantile
template <class RealType, class Policy>
RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
{
- return binomial_detail::quantile_imp(c.dist, 1-c.param, c.param);
+ return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param);
} // quantile
template <class RealType, class Policy>
Modified: sandbox/math_toolkit/boost/math/distributions/cauchy.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/cauchy.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/cauchy.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -140,7 +140,7 @@
return location;
}
result = -scale / tan(constants::pi<RealType>() * P);
- return complement ? location - result : location + result;
+ return complement ? RealType(location - result) : RealType(location + result);
} // quantile
} // namespace detail
Modified: sandbox/math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/detail/inv_discrete_quantile.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -128,7 +128,7 @@
else
{
b = a;
- a = (std::max)(b - 1, value_type(0));
+ a = (std::max)(value_type(b - 1), value_type(0));
if(a < min_bound)
a = min_bound;
fa = f(a);
@@ -158,7 +158,7 @@
}
else
{
- b = (std::max)(a - adder, value_type(0));
+ b = (std::max)(value_type(a - adder), value_type(0));
if(b < min_bound)
b = min_bound;
}
@@ -182,7 +182,7 @@
}
else
{
- b = (std::max)(a - adder, value_type(0));
+ b = (std::max)(value_type(a - adder), value_type(0));
if(b < min_bound)
b = min_bound;
}
@@ -330,7 +330,7 @@
dist,
p,
q,
- (guess < 1 ? value_type(1) : floor(guess)),
+ (guess < 1 ? value_type(1) : (value_type)floor(guess)),
multiplier,
adder,
tools::equal_floor(),
@@ -340,7 +340,7 @@
dist,
p,
q,
- ceil(guess),
+ (value_type)ceil(guess),
multiplier,
adder,
tools::equal_ceil(),
Modified: sandbox/math_toolkit/boost/math/distributions/negative_binomial.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/negative_binomial.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/negative_binomial.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -468,14 +468,14 @@
RealType guess = 0;
RealType factor = 5;
if(r * r * r * P * p > 0.005)
- guess = detail::inverse_negative_binomial_cornish_fisher(r, p, 1-p, P, 1-P, Policy());
+ guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), P, RealType(1-P), Policy());
if(guess < 10)
{
//
// Cornish-Fisher Negative binomial approximation not accurate in this area:
//
- guess = (std::min)(r * 2, RealType(10));
+ guess = (std::min)(RealType(r * 2), RealType(10));
}
else
factor = (1-P < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
@@ -545,14 +545,14 @@
RealType guess = 0;
RealType factor = 5;
if(r * r * r * (1-Q) * p > 0.005)
- guess = detail::inverse_negative_binomial_cornish_fisher(r, p, 1-p, 1-Q, Q, Policy());
+ guess = detail::inverse_negative_binomial_cornish_fisher(r, p, RealType(1-p), RealType(1-Q), Q, Policy());
if(guess < 10)
{
//
// Cornish-Fisher Negative binomial approximation not accurate in this area:
//
- guess = (std::min)(r * 2, RealType(10));
+ guess = (std::min)(RealType(r * 2), RealType(10));
}
else
factor = (Q < sqrt(tools::epsilon<RealType>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -597,7 +597,7 @@
if(l == 0)
return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
- non_central_beta_pdf(a, b, l, static_cast<value_type>(x), 1 - static_cast<value_type>(x), forwarding_policy()),
+ non_central_beta_pdf(a, b, l, static_cast<value_type>(x), value_type(1 - static_cast<value_type>(x)), forwarding_policy()),
"function");
}
@@ -781,7 +781,7 @@
if(l == 0)
return cdf(beta_distribution<RealType, Policy>(a, b), x);
- return detail::non_central_beta_cdf(x, 1 - x, a, b, l, false, Policy());
+ return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, false, Policy());
} // cdf
template <class RealType, class Policy>
@@ -818,7 +818,7 @@
if(l == 0)
return cdf(complement(beta_distribution<RealType, Policy>(a, b), x));
- return detail::non_central_beta_cdf(x, 1 - x, a, b, l, true, Policy());
+ return detail::non_central_beta_cdf(x, RealType(1 - x), a, b, l, true, Policy());
} // ccdf
template <class RealType, class Policy>
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -1042,7 +1042,7 @@
{ // Quantile (or Percent Point) function.
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
- return detail::non_central_t_quantile(v, l, p, 1-p, Policy());
+ return detail::non_central_t_quantile(v, l, p, RealType(1-p), Policy());
} // quantile
template <class RealType, class Policy>
@@ -1052,7 +1052,7 @@
RealType q = c.param;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
- return detail::non_central_t_quantile(v, l, 1-q, q, Policy());
+ return detail::non_central_t_quantile(v, l, RealType(1-q), q, Policy());
} // quantile complement.
} // namespace math
Modified: sandbox/math_toolkit/boost/math/distributions/poisson.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/poisson.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/poisson.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -479,7 +479,7 @@
if(z < 1)
guess = z;
else
- guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, 1-p, Policy());
+ guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, RealType(1-p), Policy());
if(z > 5)
{
if(z > 1000)
@@ -547,7 +547,7 @@
if(z < 1)
guess = z;
else
- guess = boost::math::detail::inverse_poisson_cornish_fisher(z, 1-q, q, Policy());
+ guess = boost::math::detail::inverse_poisson_cornish_fisher(z, RealType(1-q), q, Policy());
if(z > 5)
{
if(z > 1000)
Modified: sandbox/math_toolkit/boost/math/distributions/triangular.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/triangular.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/triangular.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -214,11 +214,11 @@
}
if (x == lower)
{ // (mode - lower) == 0 which would lead to divide by zero!
- return (mode == lower) ? 2 / (upper - lower) : 0;
+ return (mode == lower) ? 2 / (upper - lower) : RealType(0);
}
else if (x == upper)
{
- return (mode == upper) ? 2 / (upper - lower) : 0;
+ return (mode == upper) ? 2 / (upper - lower) : RealType(0);
}
else if (x <= mode)
{
Modified: sandbox/math_toolkit/boost/math/special_functions/bessel.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/bessel.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/bessel.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -128,7 +128,7 @@
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_j_imp(v, -x, t, pol);
+ T r = cyl_bessel_j_imp(v, T(-x), t, pol);
if(iround(v, pol) & 1)
r = -r;
return r;
@@ -209,7 +209,7 @@
// Default case is just a naive evaluation of the definition:
//
return sqrt(constants::pi<T>() / (2 * x))
- * cyl_bessel_j_imp(T(n)+T(0.5f), x, bessel_no_int_tag(), pol);
+ * cyl_bessel_j_imp(T(T(n)+T(0.5f)), x, bessel_no_int_tag(), pol);
}
template <class T, class Policy>
@@ -227,7 +227,7 @@
// better have integer v:
if(floor(v) == v)
{
- T r = cyl_bessel_i_imp(v, -x, pol);
+ T r = cyl_bessel_i_imp(v, T(-x), pol);
if(iround(v, pol) & 1)
r = -r;
return r;
@@ -380,7 +380,7 @@
if(x < 2 * tools::min_value<T>())
return -policies::raise_overflow_error<T>(function, 0, pol);
- T result = cyl_neumann_imp(T(v)+0.5f, x, bessel_no_int_tag(), pol);
+ T result = cyl_neumann_imp(T(T(v)+0.5f), x, bessel_no_int_tag(), pol);
T tx = sqrt(constants::pi<T>() / (2 * x));
if((tx > 1) && (tools::max_value<T>() / tx < result))
Modified: sandbox/math_toolkit/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/beta.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -822,7 +822,7 @@
BOOST_MATH_STD_USING // ADL of std names
T result = pow(x, n);
T term = result;
- for(unsigned i = itrunc(n - 1); i > k; --i)
+ for(unsigned i = itrunc(T(n - 1)); i > k; --i)
{
term *= ((i + 1) * y) / ((n - i) * x) ;
result += term;
@@ -874,15 +874,15 @@
{
if(p_derivative)
{
- *p_derivative = (a == 1) ? 1 : (a < 1) ? tools::max_value<T>() / 2 : tools::min_value<T>() * 2;
+ *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
}
- return (invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
+ return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
}
if(x == 1)
{
if(p_derivative)
{
- *p_derivative = (b == 1) ? 1 : (b < 1) ? tools::max_value<T>() / 2 : tools::min_value<T>() * 2;
+ *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
}
return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
}
@@ -931,7 +931,7 @@
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+b, a, 20);
+ prefix = rising_factorial_ratio(T(a+b), a, 20);
}
else
{
@@ -939,12 +939,12 @@
}
fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
if(!invert)
- fract = beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
else
{
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
- fract = -beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
}
}
}
@@ -997,7 +997,7 @@
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+b, a, 20);
+ prefix = rising_factorial_ratio(T(a+b), a, 20);
}
else
{
@@ -1005,12 +1005,12 @@
}
fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative);
if(!invert)
- fract = beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
else
{
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
invert = false;
- fract = -beta_small_b_large_a_series(a + 20, b, x, y, fract, prefix, pol, normalised);
+ fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised);
}
}
}
@@ -1060,14 +1060,14 @@
else if(a > 15)
{
// sidestep so we can use the series representation:
- int n = itrunc(floor(b), pol);
+ int n = itrunc(T(floor(b)), pol);
if(n == b)
--n;
T bbar = b - n;
T prefix;
if(!normalised)
{
- prefix = rising_factorial_ratio(a+bbar, bbar, n);
+ prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
}
else
{
@@ -1082,7 +1082,7 @@
// the formula here for the non-normalised case is tricky to figure
// out (for me!!), and requires two pochhammer calculations rather
// than one, so leave it for now....
- int n = itrunc(floor(b), pol);
+ int n = itrunc(T(floor(b)), pol);
T bbar = b - n;
if(bbar <= 0)
{
@@ -1094,7 +1094,7 @@
if(invert)
fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
//fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y);
- fract = beta_small_b_large_a_series(a+20, bbar, x, y, fract, T(1), pol, normalised);
+ fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised);
if(invert)
{
fract = -fract;
@@ -1167,7 +1167,7 @@
// Now the regular cases:
//
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
- T f1 = ibeta_power_terms(a, b, x, 1 - x, lanczos_type(), true, pol);
+ T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol);
T y = (1 - x) * x;
if(f1 == 0)
Modified: sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -44,7 +44,7 @@
return 0;
rem = cos(constants::pi<T>() * rem);
- return invert ? -rem : rem;
+ return invert ? T(-rem) : rem;
}
} // namespace detail
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -49,11 +49,11 @@
b = exp(v * a);
sigma = -a * v;
c = abs(v) < tools::epsilon<T>() ?
- 1 : boost::math::sin_pi(v) / (v * pi<T>());
+ T(1) : T(boost::math::sin_pi(v) / (v * pi<T>()));
d = abs(sigma) < tools::epsilon<T>() ?
- 1 : sinh(sigma) / sigma;
+ T(1) : T(sinh(sigma) / sigma);
gamma1 = abs(v) < tools::epsilon<T>() ?
- -euler<T>() : (0.5f / v) * (gp - gm) * c;
+ T(-euler<T>()) : T((0.5f / v) * (gp - gm) * c);
gamma2 = (2 + gp + gm) * c / 2;
// initial values
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -55,12 +55,12 @@
sigma = -a * v;
d = abs(sigma) < tools::epsilon<T>() ?
T(1) : sinh(sigma) / sigma;
- e = abs(v) < tools::epsilon<T>() ? v*pi<T>()*pi<T>() / 2
- : 2 * spv2 * spv2 / v;
+ e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
+ : T(2 * spv2 * spv2 / v);
- T g1 = (v == 0) ? -euler<T>() : (gp - gm) / ((1 + gp) * (1 + gm) * 2 * v);
+ T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
- T vspv = (fabs(v) < tools::epsilon<T>()) ? 1/constants::pi<T>() : v / spv;
+ T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
p = vspv / (xp * (1 + gm));
@@ -118,7 +118,7 @@
tolerance = 2 * tools::epsilon<T>();
tiny = sqrt(tools::min_value<T>());
C = f = tiny; // b0 = 0, replace with tiny
- D = 0.0L;
+ D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
{
a = -1;
@@ -131,7 +131,7 @@
delta = C * D;
f *= delta;
if (D < 0) { s = -s; }
- if (abs(delta - 1.0L) < tolerance)
+ if (abs(delta - 1) < tolerance)
{ break; }
}
policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
@@ -158,7 +158,7 @@
typedef typename complex_trait<T>::type complex_type;
complex_type C, D, f, a, b, delta, one(1);
- T tiny, zero(0.0L);
+ T tiny, zero(0);
unsigned long k;
// |x| >= |v|, CF2_jy converges rapidly
@@ -169,7 +169,7 @@
// Lentz, Applied Optics, vol 15, 668 (1976)
T tolerance = 2 * tools::epsilon<T>();
tiny = sqrt(tools::min_value<T>());
- C = f = complex_type(-0.5f/x, 1.0L);
+ C = f = complex_type(-0.5f/x, 1);
D = 0;
for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
{
@@ -289,7 +289,7 @@
if(kind&need_y)
{
Yu = asymptotic_bessel_y_large_x_2(u, x);
- Yu1 = asymptotic_bessel_y_large_x_2(u + 1, x);
+ Yu1 = asymptotic_bessel_y_large_x_2(T(u + 1), x);
}
else
Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy_asym.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -208,7 +208,7 @@
{
// default case:
BOOST_MATH_STD_USING
- T v2 = (std::max)(T(3), v * v);
+ T v2 = (std::max)(T(3), T(v * v));
return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
}
template <class T>
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/gamma_inva.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -176,7 +176,7 @@
detail::gamma_inva_imp(
static_cast<value_type>(x),
static_cast<value_type>(p),
- 1 - static_cast<value_type>(p),
+ static_cast<value_type>(1 - static_cast<value_type>(p)),
pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");
}
@@ -205,7 +205,7 @@
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
detail::gamma_inva_imp(
static_cast<value_type>(x),
- 1 - static_cast<value_type>(q),
+ static_cast<value_type>(1 - static_cast<value_type>(q)),
static_cast<value_type>(q),
pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");
}
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inv_ab.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -120,11 +120,11 @@
//
if((p < q) != swap_ab)
{
- guess = (std::min)(b * 2, T(1));
+ guess = (std::min)(T(b * 2), T(1));
}
else
{
- guess = (std::min)(b / 2, T(1));
+ guess = (std::min)(T(b / 2), T(1));
}
}
if(n * n * n * u * sf > 0.005)
@@ -137,11 +137,11 @@
//
if((p < q) != swap_ab)
{
- guess = (std::min)(b * 2, T(10));
+ guess = (std::min)(T(b * 2), T(10));
}
else
{
- guess = (std::min)(b / 2, T(10));
+ guess = (std::min)(T(b / 2), T(10));
}
}
else
@@ -186,7 +186,7 @@
static_cast<value_type>(b),
static_cast<value_type>(x),
static_cast<value_type>(p),
- 1 - static_cast<value_type>(p),
+ static_cast<value_type>(1 - static_cast<value_type>(p)),
false, pol),
"boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
}
@@ -217,7 +217,7 @@
detail::ibeta_inv_ab_imp(
static_cast<value_type>(b),
static_cast<value_type>(x),
- 1 - static_cast<value_type>(q),
+ static_cast<value_type>(1 - static_cast<value_type>(q)),
static_cast<value_type>(q),
false, pol),
"boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
@@ -250,7 +250,7 @@
static_cast<value_type>(a),
static_cast<value_type>(x),
static_cast<value_type>(p),
- 1 - static_cast<value_type>(p),
+ static_cast<value_type>(1 - static_cast<value_type>(p)),
true, pol),
"boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
}
@@ -281,7 +281,7 @@
detail::ibeta_inv_ab_imp(
static_cast<value_type>(a),
static_cast<value_type>(x),
- 1 - static_cast<value_type>(q),
+ static_cast<value_type>(1 - static_cast<value_type>(q)),
static_cast<value_type>(q),
true, pol),
"boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/ibeta_inverse.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -106,7 +106,7 @@
//
// Bring them together to get a final estimate for eta:
//
- T eta = tools::evaluate_polynomial(terms, 1/a, 4);
+ T eta = tools::evaluate_polynomial(terms, T(1/a), 4);
//
// now we need to convert eta to x, by solving the appropriate
// quadratic equation:
@@ -207,7 +207,7 @@
// Bring the correction terms together to evaluate eta,
// this is the last equation on page 151:
//
- T eta = tools::evaluate_polynomial(terms, 1/r, 4);
+ T eta = tools::evaluate_polynomial(terms, T(1/r), 4);
//
// Now that we have eta we need to back solve for x,
// we seek the value of x that gives eta in Eq 3.2.
@@ -660,7 +660,7 @@
//
T lx = log(p * a * boost::math::beta(a, b, pol)) / a;
x = exp(lx);
- y = x < 0.9 ? 1 - x : -boost::math::expm1(lx, pol);
+ y = x < 0.9 ? T(1 - x) : (T)(-boost::math::expm1(lx, pol));
if((b < a) && (x < 0.2))
{
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/igamma_inverse.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -209,7 +209,7 @@
}
else
{
- T D = (std::max)(T(2), a * (a - 1));
+ T D = (std::max)(T(2), T(a * (a - 1)));
T lg = boost::math::lgamma(a, pol);
T lb = log(q) + lg;
if(lb < -D * 2.3)
@@ -357,7 +357,7 @@
return tools::max_value<T>();
if(p == 0)
return 0;
- T guess = detail::find_inverse_gamma(a, p, 1 - p, pol);
+ T guess = detail::find_inverse_gamma<T>(a, p, 1 - p, pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
@@ -399,7 +399,7 @@
return tools::max_value<T>();
if(q == 1)
return 0;
- T guess = detail::find_inverse_gamma(a, 1 - q, q, pol);
+ T guess = detail::find_inverse_gamma<T>(a, 1 - q, q, pol);
T lower = tools::min_value<T>();
if(guess <= lower)
guess = tools::min_value<T>();
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -142,35 +142,60 @@
// Figure out what the coefficients are, note these depend
// only on the degrees of freedom (Eq 57 of Shaw):
//
- c[2] = T(1) / 6 + T(1) / (6 * df);
+ c[2] = 0.16666666666666666667 + 0.16666666666666666667 / df;
T in = 1 / df;
- c[3] = (((T(1) / 120) * in) + (T(1) / 15)) * in + (T(7) / 120);
- c[4] = ((((T(1) / 5040) * in + (T(1) / 560)) * in + (T(3) / 112)) * in + T(127) / 5040);
- c[5] = ((((T(1) / 362880) * in + (T(17) / 45360)) * in + (T(67) / 60480)) * in + (T(479) / 45360)) * in + (T(4369) / 362880);
- c[6] = ((((((T(1) / 39916800) * in + (T(2503) / 39916800)) * in + (T(11867) / 19958400)) * in + (T(1285) / 798336)) * in + (T(153161) / 39916800)) * in + (T(34807) / 5702400));
- c[7] = (((((((T(1) / 6227020800LL) * in + (T(37) / 2402400)) * in +
- (T(339929) / 2075673600LL)) * in + (T(67217) / 97297200)) * in +
- (T(870341) / 691891200LL)) * in + (T(70691) / 64864800LL)) * in +
- (T(20036983LL) / 6227020800LL));
- c[8] = (((((((T(1) / 1307674368000LL) * in + (T(1042243LL) / 261534873600LL)) * in +
- (T(21470159) / 435891456000LL)) * in + (T(326228899LL) / 1307674368000LL)) * in +
- (T(843620579) / 1307674368000LL)) * in + (T(332346031LL) / 435891456000LL)) * in +
- (T(43847599) / 1307674368000LL)) * in + (T(2280356863LL) / 1307674368000LL);
- c[9] = (((((((((T(1) / 355687428096000LL)) * in + (T(24262727LL) / 22230464256000LL)) * in +
- (T(123706507LL) / 8083805184000LL)) * in + (T(404003599LL) / 4446092851200LL)) * in +
- (T(51811946317LL) / 177843714048000LL)) * in + (T(91423417LL) / 177843714048LL)) * in +
- (T(32285445833LL) / 88921857024000LL)) * in + (T(531839683LL) / 1710035712000LL)) * in +
- (T(49020204823LL) / 50812489728000LL);
- c[10] = (((((((((T(1) / 121645100408832000LL) * in +
- (T(4222378423LL) / 13516122267648000LL)) * in +
- (T(49573465457LL) / 10137091700736000LL)) * in +
- (T(176126809LL) / 5304600576000LL)) * in +
- (T(44978231873LL) / 355687428096000LL)) * in +
- (T(5816850595639LL) / 20274183401472000LL)) * in +
- (T(73989712601LL) / 206879422464000LL)) * in +
- (T(26591354017LL) / 259925428224000LL)) * in +
- (T(14979648446341LL) / 40548366802944000LL)) * in +
- (T(65967241200001LL) / 121645100408832000LL);
+ c[3] = (0.0083333333333333333333 * in
+ + 0.066666666666666666667) * in
+ + 0.058333333333333333333;
+ c[4] = ((0.00019841269841269841270 * in
+ + 0.0017857142857142857143) * in
+ + 0.026785714285714285714) * in
+ + 0.025198412698412698413;
+ c[5] = (((2.7557319223985890653e10-6 * in
+ + 0.00037477954144620811287) * in
+ - 0.0011078042328042328042) * in
+ + 0.010559964726631393298) * in
+ + 0.012039792768959435626;
+ c[6] = ((((2.5052108385441718775e-8 * in
+ - 0.000062705427288760622094) * in
+ + 0.00059458674042007375341) * in
+ - 0.0016095979637646304313) * in
+ + 0.0061039211560044893378) * in
+ + 0.0038370059724226390893;
+ c[7] = (((((1.6059043836821614599e-10 * in
+ + 0.000015401265401265401265) * in
+ - 0.00016376804137220803887) * in
+ + 0.00069084207973096861986) * in
+ - 0.0012579159844784844785) * in
+ + 0.0010898206731540064873) * in
+ + 0.0032177478835464946576;
+ c[8] = ((((((7.6471637318198164759e-13 * in
+ - 3.9851014346715404916e-6) * in
+ + 0.000049255746366361445727) * in
+ - 0.00024947258047043099953) * in
+ + 0.00064513046951456342991) * in
+ - 0.00076245135440323932387) * in
+ + 0.000033530976880017885309) * in
+ + 0.0017438262298340009980;
+ c[9] = (((((((2.8114572543455207632e-15 * in
+ + 1.0914179173496789432e-6) * in
+ - 0.000015303004486655377567) * in
+ + 0.000090867107935219902229) * in
+ - 0.00029133414466938067350) * in
+ + 0.00051406605788341121363) * in
+ - 0.00036307660358786885787) * in
+ - 0.00031101086326318780412) * in
+ + 0.00096472747321388644237;
+ c[10] = ((((((((8.2206352466243297170e-18 * in
+ - 3.1239569599829868045e-7) * in
+ + 4.8903045291975346210e-6) * in
+ - 0.000033202652391372058698) * in
+ + 0.00012645437628698076975) * in
+ - 0.00028690924218514613987) * in
+ + 0.00035764655430568632777) * in
+ - 0.00010230378073700412687) * in
+ - 0.00036942667800009661203) * in
+ + 0.00054229262813129686486;
//
// The result is then a polynomial in v (see Eq 56 of Shaw):
//
@@ -240,7 +265,7 @@
T root_alpha = sqrt(alpha);
T r = 4 * cos(acos(root_alpha) / 3) / root_alpha;
T x = sqrt(r - 4);
- result = u - 0.5f < 0 ? -x : x;
+ result = u - 0.5f < 0 ? (T)-x : x;
if(pexact)
*pexact = true;
break;
@@ -259,7 +284,7 @@
//
T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
T b = boost::math::cbrt(a);
- static const T c = 0.85498797333834849467655443627193L;
+ static const T c = 0.85498797333834849467655443627193;
T p = 6 * (1 + c * (1 / b - 1));
T p0;
do{
@@ -274,7 +299,7 @@
// Use Eq 45 to extract the result:
//
p = sqrt(p - df);
- result = (u - 0.5f) < 0 ? -p : p;
+ result = (u - 0.5f) < 0 ? (T)-p : p;
break;
}
#if 0
@@ -369,7 +394,7 @@
// where we use Shaw's tail series.
// The crossover point is roughly exponential in -df:
//
- T crossover = ldexp(1.0f, iround(df / -0.654f, pol));
+ T crossover = ldexp(1.0f, iround(T(df / -0.654f), pol));
if(u > crossover)
{
result = boost::math::detail::inverse_students_t_hill(df, u, pol);
@@ -380,13 +405,13 @@
}
}
}
- return invert ? -result : result;
+ return invert ? (T)-result : result;
}
template <class T, class Policy>
inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol)
{
- T u = (p > q) ? 0.5f - q / 2 : p / 2;
+ T u = (p > q) ? T(0.5f - q) / T(2) : T(p / 2);
T v = 1 - u; // u < 0.5 so no cancellation error
T df = a * 2;
T t = boost::math::detail::inverse_students_t(df, u, v, pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_1.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -88,7 +88,7 @@
// so rewritten to use fmod instead:
//
BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi<T>() / 2);
- T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(phi, T(constants::pi<T>() / 2));
BOOST_MATH_INSTRUMENT_VARIABLE(rphi);
T m = 2 * (phi - rphi) / constants::pi<T>();
BOOST_MATH_INSTRUMENT_VARIABLE(m);
@@ -104,7 +104,7 @@
T cosp = cos(rphi);
BOOST_MATH_INSTRUMENT_VARIABLE(sinp);
BOOST_MATH_INSTRUMENT_VARIABLE(cosp);
- result = s * sinp * ellint_rf_imp(cosp * cosp, 1 - k * k * sinp * sinp, T(1), pol);
+ result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(m != 0)
{
@@ -112,7 +112,7 @@
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
}
- return invert ? -result : result;
+ return invert ? T(-result) : result;
}
// Complete elliptic integral (Legendre form) of the first kind
Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_2.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -74,7 +74,7 @@
// but that fails if T has more digits than a long long,
// so rewritten to use fmod instead:
//
- T rphi = boost::math::tools::fmod_workaround(phi, constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(phi, T(constants::pi<T>() / 2));
T m = 2 * (phi - rphi) / constants::pi<T>();
int s = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
@@ -93,7 +93,7 @@
if(m != 0)
result += m * ellint_e_imp(k, pol);
}
- return invert ? -result : result;
+ return invert ? T(-result) : result;
}
// Complete elliptic integral (Legendre form) of the second kind
Modified: sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/ellint_3.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -182,7 +182,7 @@
}
else
{
- T rphi = boost::math::tools::fmod_workaround(fabs(phi), constants::pi<T>() / 2);
+ T rphi = boost::math::tools::fmod_workaround(fabs(phi), T(constants::pi<T>() / 2));
T m = 2 * (fabs(phi) - rphi) / constants::pi<T>();
int sign = 1;
if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5)
Modified: sandbox/math_toolkit/boost/math/special_functions/erf.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/erf.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/erf.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -122,9 +122,9 @@
if(z < 0)
{
if(!invert)
- return -erf_imp(-z, invert, pol, t);
+ return -erf_imp(T(-z), invert, pol, t);
else
- return 1 + erf_imp(-z, false, pol, t);
+ return 1 + erf_imp(T(-z), false, pol, t);
}
T result;
Modified: sandbox/math_toolkit/boost/math/special_functions/expint.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/expint.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/expint.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -436,7 +436,7 @@
if(z < 0)
return policies::raise_domain_error<T>(function, "Function requires z >= 0 but got %1%.", z, pol);
if(z == 0)
- return n == 1 ? policies::raise_overflow_error<T>(function, 0, pol) : 1 / (static_cast<T>(n - 1));
+ return n == 1 ? policies::raise_overflow_error<T>(function, 0, pol) : T(1 / (static_cast<T>(n - 1)));
T result;
@@ -504,7 +504,7 @@
{
static const char* function = "boost::math::expint<%1%>(%1%)";
if(z < 0)
- return -expint_imp(1, -z, pol, tag);
+ return -expint_imp(1, T(-z), pol, tag);
if(z == 0)
return -policies::raise_overflow_error<T>(function, 0, pol);
return expint_i_as_series(z, pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/expm1.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/expm1.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/expm1.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -75,7 +75,7 @@
BOOST_MATH_STD_USING
T a = fabs(x);
- if(a > T(0.5L))
+ if(a > T(0.5f))
return exp(x) - T(1);
if(a < tools::epsilon<T>())
return x;
Modified: sandbox/math_toolkit/boost/math/special_functions/factorials.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/factorials.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/factorials.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -163,7 +163,7 @@
// handle it, split the product up into three parts:
//
T xp1 = x + 1;
- unsigned n2 = itrunc(floor(xp1), pol);
+ unsigned n2 = itrunc((T)floor(xp1), pol);
if(n2 == xp1)
return 0;
T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/gamma.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -44,6 +44,7 @@
#include <boost/assert.hpp>
#include <boost/mpl/greater.hpp>
#include <boost/mpl/equal_to.hpp>
+#include <boost/mpl/greater.hpp>
#include <boost/config/no_tr1/cmath.hpp>
#include <algorithm>
@@ -148,7 +149,7 @@
return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
if(z <= -20)
{
- result = gamma_imp(-z, pol, l) * sinpx(z);
+ result = gamma_imp(T(-z), pol, l) * sinpx(z);
if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
result = -boost::math::constants::pi<T>() / result;
@@ -236,13 +237,19 @@
{
typedef typename policies::precision<T, Policy>::type precision_type;
typedef typename mpl::if_<
- mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::and_<
+ mpl::less_equal<precision_type, mpl::int_<64> >,
+ mpl::greater<precision_type, mpl::int_<0> >
+ >,
mpl::int_<64>,
typename mpl::if_<
- mpl::less_equal<precision_type, mpl::int_<113> >,
+ mpl::and_<
+ mpl::less_equal<precision_type, mpl::int_<113> >,
+ mpl::greater<precision_type, mpl::int_<0> >
+ >,
mpl::int_<113>, mpl::int_<0> >::type
>::type tag_type;
- result = lgamma_small_imp(z, z - 1, z - 2, tag_type(), pol, l);
+ result = lgamma_small_imp<T>(z, z - 1, z - 2, tag_type(), pol, l);
}
else if((z >= 3) && (z < 100))
{
@@ -465,7 +472,7 @@
{
// Use expm1 on lgamma:
result = boost::math::expm1(-boost::math::log1p(dz, pol)
- + lgamma_small_imp(dz+2, dz + 1, dz, tag_type(), pol, l));
+ + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
BOOST_MATH_INSTRUMENT_CODE(result);
}
}
@@ -474,7 +481,7 @@
if(dz < 2)
{
// Use expm1 on lgamma:
- result = boost::math::expm1(lgamma_small_imp(dz+1, dz, dz-1, tag_type(), pol, l), pol);
+ result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
BOOST_MATH_INSTRUMENT_CODE(result);
}
else
@@ -1097,7 +1104,7 @@
//
if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
{
- return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(z + delta) - 1);
+ return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta)) - 1);
}
}
if(fabs(delta) < 20)
@@ -1428,7 +1435,7 @@
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
- return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b) - static_cast<value_type>(a), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
+ return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type
Modified: sandbox/math_toolkit/boost/math/special_functions/log1p.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/log1p.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/log1p.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -86,7 +86,7 @@
function, 0, pol);
result_type a = abs(result_type(x));
- if(a > result_type(0.5L))
+ if(a > result_type(0.5f))
return log(1 + result_type(x));
// Note that without numeric_limits specialisation support,
// epsilon just returns zero, and our "optimisation" will always fail:
@@ -432,7 +432,7 @@
function, 0, pol);
result_type a = abs(result_type(x));
- if(a > result_type(0.95L))
+ if(a > result_type(0.95f))
return log(1 + result_type(x)) - result_type(x);
// Note that without numeric_limits specialisation support,
// epsilon just returns zero, and our "optimisation" will always fail:
Modified: sandbox/math_toolkit/boost/math/special_functions/next.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/next.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/next.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -183,18 +183,18 @@
if(a == b)
return 0;
if(a == 0)
- return 1 + fabs(float_distance(boost::math::sign(b) * detail::get_smallest_value<T>(), b, pol));
+ return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol));
if(b == 0)
- return 1 + fabs(float_distance(boost::math::sign(a) * detail::get_smallest_value<T>(), a, pol));
+ return 1 + fabs(float_distance(static_cast<T>(boost::math::sign(a) * detail::get_smallest_value<T>()), a, pol));
if(boost::math::sign(a) != boost::math::sign(b))
- return 2 + fabs(float_distance(boost::math::sign(b) * detail::get_smallest_value<T>(), b, pol))
- + fabs(float_distance(boost::math::sign(a) * detail::get_smallest_value<T>(), a, pol));
+ return 2 + fabs(float_distance(static_cast<T>(boost::math::sign(b) * detail::get_smallest_value<T>()), b, pol))
+ + fabs(float_distance(static_cast<T>(boost::math::sign(a) * detail::get_smallest_value<T>()), a, pol));
//
// By the time we get here, both a and b must have the same sign, we want
// b > a and both postive for the following logic:
//
if(a < 0)
- return float_distance(-b, -a);
+ return float_distance(static_cast<T>(-b), static_cast<T>(-a));
BOOST_ASSERT(a >= 0);
BOOST_ASSERT(b >= a);
Modified: sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -46,7 +46,7 @@
return static_cast<T>(invert ? -1 : 1);
rem = sin(constants::pi<T>() * rem);
- return invert ? -rem : rem;
+ return invert ? T(-rem) : rem;
}
} // namespace detail
Modified: sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/spherical_harmonic.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -61,7 +61,7 @@
if(m&1)
{
// Check phase if theta is outside [0, PI]:
- T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
if(mod < 0)
mod += 2 * constants::pi<T>();
if(mod > constants::pi<T>())
@@ -70,7 +70,7 @@
// Get the value and adjust sign as required:
T prefix = spherical_harmonic_prefix(n, m, theta, pol);
prefix *= cos(m * phi);
- return sign ? -prefix : prefix;
+ return sign ? T(-prefix) : prefix;
}
template <class T, class Policy>
@@ -88,7 +88,7 @@
if(m&1)
{
// Check phase if theta is outside [0, PI]:
- T mod = boost::math::tools::fmod_workaround(theta, 2 * constants::pi<T>());
+ T mod = boost::math::tools::fmod_workaround(theta, T(2 * constants::pi<T>()));
if(mod < 0)
mod += 2 * constants::pi<T>();
if(mod > constants::pi<T>())
@@ -97,7 +97,7 @@
// Get the value and adjust sign as required:
T prefix = spherical_harmonic_prefix(n, m, theta, pol);
prefix *= sin(m * phi);
- return sign ? -prefix : prefix;
+ return sign ? T(-prefix) : prefix;
}
template <class T, class U, class Policy>
Modified: sandbox/math_toolkit/boost/math/special_functions/trunc.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/trunc.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/trunc.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -22,7 +22,7 @@
BOOST_MATH_STD_USING
if(!(boost::math::isfinite)(v))
return policies::raise_rounding_error("boost::math::trunc<%1%>(%1%)", 0, v, pol);
- return (v >= 0) ? floor(v) : ceil(v);
+ return (v >= 0) ? static_cast<T>(floor(v)) : static_cast<T>(ceil(v));
}
template <class T>
inline T trunc(const T& v)
Modified: sandbox/math_toolkit/boost/math/special_functions/zeta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/zeta.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/zeta.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -884,7 +884,7 @@
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
static_cast<value_type>(s),
- 1 - static_cast<value_type>(s),
+ static_cast<value_type>(1 - static_cast<value_type>(s)),
forwarding_policy(),
tag_type()), "boost::math::zeta<%1%>(%1%)");
}
Modified: sandbox/math_toolkit/boost/math/tools/minima.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/minima.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/minima.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -77,7 +77,7 @@
delta = p / q;
u = x + delta;
if(((u - min) < fract2) || ((max- u) < fract2))
- delta = (mid - x) < 0 ? -fabs(fract1) : fabs(fract1);
+ delta = (mid - x) < 0 ? (T)-fabs(fract1) : (T)fabs(fract1);
}
}
else
@@ -87,7 +87,7 @@
delta = golden * delta2;
}
// update current position:
- u = (fabs(delta) >= fract1) ? x + delta : (delta > 0 ? x + fabs(fract1) : x - fabs(fract1));
+ u = (fabs(delta) >= fract1) ? T(x + delta) : (delta > 0 ? T(x + fabs(fract1)) : T(x - fabs(fract1)));
fu = f(u);
if(fu <= fx)
{
Modified: sandbox/math_toolkit/boost/math/tools/rational.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/rational.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/rational.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -220,19 +220,19 @@
template <class T, class U>
inline U evaluate_even_polynomial(const T* poly, U z, std::size_t count)
{
- return evaluate_polynomial(poly, z*z, count);
+ return evaluate_polynomial(poly, U(z*z), count);
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const T(&a)[N], const V& z)
{
- return evaluate_polynomial(a, z*z);
+ return evaluate_polynomial(a, V(z*z));
}
template <std::size_t N, class T, class V>
inline V evaluate_even_polynomial(const boost::array<T,N>& a, const V& z)
{
- return evaluate_polynomial(a, z*z);
+ return evaluate_polynomial(a, V(z*z));
}
//
// Odd polynomials come next:
@@ -240,21 +240,21 @@
template <class T, class U>
inline U evaluate_odd_polynomial(const T* poly, U z, std::size_t count)
{
- return poly[0] + z * evaluate_polynomial(poly+1, z*z, count-1);
+ return poly[0] + z * evaluate_polynomial(poly+1, U(z*z), count-1);
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const T(&a)[N], const V& z)
{
typedef mpl::int_<N-1> tag_type;
- return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, z*z, static_cast<tag_type const*>(0));
+ return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a) + 1, V(z*z), static_cast<tag_type const*>(0));
}
template <std::size_t N, class T, class V>
inline V evaluate_odd_polynomial(const boost::array<T,N>& a, const V& z)
{
typedef mpl::int_<N-1> tag_type;
- return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, z*z, static_cast<tag_type const*>(0));
+ return a[0] + z * detail::evaluate_polynomial_c_imp(static_cast<const T*>(a.data()) + 1, V(z*z), static_cast<tag_type const*>(0));
}
template <class T, class U, class V>
Modified: sandbox/math_toolkit/boost/math/tools/roots.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/roots.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/roots.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -279,7 +279,7 @@
T result = guess;
T factor = static_cast<T>(ldexp(1.0, 1 - digits));
- T delta = (std::max)(10000000 * guess, T(10000000)); // arbitarily large delta
+ T delta = (std::max)(T(10000000 * guess), T(10000000)); // arbitarily large delta
T last_f0 = 0;
T delta1 = delta;
T delta2 = delta;
@@ -347,7 +347,7 @@
// check for out of bounds step:
if(result < min)
{
- T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? 1000 : result / min;
+ T diff = ((fabs(min) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(min))) ? T(1000) : T(result / min);
if(fabs(diff) < 1)
diff = 1 / diff;
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
@@ -368,7 +368,7 @@
}
else if(result > max)
{
- T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? 1000 : result / max;
+ T diff = ((fabs(max) < 1) && (fabs(result) > 1) && (tools::max_value<T>() / fabs(result) < fabs(max))) ? T(1000) : T(result / max);
if(fabs(diff) < 1)
diff = 1 / diff;
if(!out_of_bounds_sentry && (diff > 0) && (diff < 3))
Modified: sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -26,7 +26,7 @@
eps_tolerance(unsigned bits)
{
BOOST_MATH_STD_USING
- eps = (std::max)(T(ldexp(1.0F, 1-bits)), 2 * tools::epsilon<T>());
+ eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(2 * tools::epsilon<T>()));
}
bool operator()(const T& a, const T& b)
{
@@ -193,9 +193,9 @@
//
// Start by obtaining the coefficients of the quadratic polynomial:
//
- T B = safe_div(fb - fa, b - a, tools::max_value<T>());
- T A = safe_div(fd - fb, d - b, tools::max_value<T>());
- A = safe_div(A - B, d - a, T(0));
+ T B = safe_div(T(fb - fa), T(b - a), tools::max_value<T>());
+ T A = safe_div(T(fd - fb), T(d - b), tools::max_value<T>());
+ A = safe_div(T(A - B), T(d - a), T(0));
if(a == 0)
{
@@ -220,7 +220,7 @@
for(unsigned i = 1; i <= count; ++i)
{
//c -= safe_div(B * c, (B + A * (2 * c - a - b)), 1 + c - a);
- c -= safe_div(fa+(B+A*(c-b))*(c-a), (B + A * (2 * c - a - b)), 1 + c - a);
+ c -= safe_div(T(fa+(B+A*(c-b))*(c-a)), T(B + A * (2 * c - a - b)), T(1 + c - a));
}
if((c <= a) || (c >= b))
{
@@ -436,7 +436,7 @@
//
e = d;
fe = fd;
- detail::bracket(f, a, b, a + (b - a) / 2, fa, fb, d, fd);
+ detail::bracket(f, a, b, T(a + (b - a) / 2), fa, fb, d, fd);
--count;
BOOST_MATH_INSTRUMENT_CODE("Not converging: Taking a bisection!!!!");
BOOST_MATH_INSTRUMENT_CODE(" a = " << a << " b = " << b);
Added: sandbox/math_toolkit/libs/math/test/mpfr_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/mpfr_concept_check.cpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,32 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that mpfr_class meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+#define TEST_MPFR
+
+#pragma warning(disable:4800)
+#pragma warning(disable:4512)
+#pragma warning(disable:4127)
+
+#include <boost/math/bindings/mpfr.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+void foo()
+{
+ instantiate(mpfr_class());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<mpfr_class>));
+}
+
+
Added: sandbox/math_toolkit/libs/math/test/ntl_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/ntl_concept_check.cpp 2008-09-11 13:58:11 EDT (Thu, 11 Sep 2008)
@@ -0,0 +1,26 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that boost::math::ntl::RR meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+
+#include <boost/math/bindings/rr.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+void foo()
+{
+ instantiate(boost::math::ntl::RR());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<boost::math::ntl::RR>));
+}
+
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