libs/math/test/test_pFq.hpp

   1 // Copyright John Maddock 2006.
   2 // Copyright Paul A. Bristow 2007, 2009
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
   8
   9 #include <boost/math/concepts/real_concept.hpp>
  10 #include <boost/math/special_functions/math_fwd.hpp>
  11 #define BOOST_TEST_MAIN
  12 #include <boost/test/unit_test.hpp>
  13 #include <boost/test/tools/floating_point_comparison.hpp>
  14 #include <boost/math/tools/stats.hpp>
  15 #include <boost/math/tools/test.hpp>
  16 #include <boost/math/tools/big_constant.hpp>
  17 #include <boost/math/constants/constants.hpp>
  18 #include <boost/type_traits/is_floating_point.hpp>
  19 #include <boost/array.hpp>
  20 #include "functor.hpp"
  21
  22 #include "handle_test_result.hpp"
  23 #include "table_type.hpp"
  24
  25 #include <boost/math/special_functions/hypergeometric_pFq.hpp>
  26 #include <boost/math/special_functions/relative_difference.hpp>
  27
  28 #ifdef BOOST_MSVC
  29 #pragma warning(disable:4127)
  30 #endif
  31
  32 #ifndef SC_
  33 #define SC_(x) BOOST_MATH_BIG_CONSTANT(T, 1000000, x)
  34 #endif
  35
  36 template <class Seq>
  37 bool is_small_a(const Seq& a)
  38 {
  39    if (a.size() == 1)
  40    {
  41       auto v = *a.begin();
  42       if ((v > -14) && (v < 1))
  43          return true;
  44    }
  45    return false;
  46 }
  47
  48 template <class Seq>
  49 bool has_negative_ab(const Seq& a, const Seq& b)
  50 {
  51    for(auto p = a.begin(); p != a.end(); ++p)
  52    {
  53       if(*p < 0)
  54          return true;
  55    }
  56    for(auto p = b.begin(); p != b.end(); ++p)
  57    {
  58       if(*p < 0)
  59          return true;
  60    }
  61    return false;
  62 }
  63
  64 template <class T>
  65 void check_pFq_result(const T& result, const T& norm, const T& expect, const std::initializer_list<T>& a, const std::initializer_list<T>& b, const T& z)
  66 {
  67    //
  68    // Ideally the error rate we calculate from comparing norm to result
  69    // should be larger than the actual error.  However, in practice even
  70    // if all the terms are positive and norm == result there will still
  71    // be a small error from the actual summation (we could work out how
  72    // much from the number of terms summed, but that's overkill for this)
  73    // so we add a small fudge factor when comparing errors:
  74    //
  75    T err = boost::math::relative_difference(result, expect);
  76    T found_err = norm / fabs(result);
  77    T fudge_factor = 25;
  78    if (is_small_a(a))
  79       fudge_factor *= 4;  // not sure why??
  80    if ((has_negative_ab(a, b)) || ((a.size() == 2) && (b.size() == 1)) || (boost::math::tools::epsilon<T>() < boost::math::tools::epsilon<double>()))
  81    {
  82       T min_err = boost::math::tools::epsilon<T>() * 600 / found_err;
  83       fudge_factor = (std::max)(fudge_factor, min_err);
  84    }
  85    if ((((err > fudge_factor * found_err) && (found_err < 1)) || (boost::math::isnan)(found_err)) && (!(boost::math::isinf)(result)))
  86    {
  87       std::cout << "Found error = " << err << " error from norm = " << found_err << std::endl;
  88       std::cout << "Testing fudge factor = " << fudge_factor << std::endl;
  89       std::cout << "  a = ";
  90       for (auto pa = a.begin(); pa != a.end(); ++pa)
  91          std::cout << *pa << ",";
  92       std::cout << "\n  b = ";
  93       for (auto pb = b.begin(); pb != b.end(); ++pb)
  94          std::cout << *pb << ",";
  95       std::cout << "\n  z = " << z << std::endl;
  96       //
  97       // This will fail if we've got here:
  98       //
  99       BOOST_CHECK_LE(err, fudge_factor * found_err);
 100       BOOST_CHECK(!(boost::math::isnan)(found_err));
 101    }
 102 }
 103
 104 template <class T>
 105 void test_spots_1F0(T, const char*)
 106 {
 107    using std::pow;
 108
 109    T tolerance = boost::math::tools::epsilon<T>() * 1000;
 110
 111    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(2)), T(-1), tolerance);
 112    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(4)), T(-27), tolerance);
 113    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(0.5)), T(0.125), tolerance);
 114    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(0.5)), T(8), tolerance);
 115    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(2)), T(-1), tolerance);
 116    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(4)), T(T(-1) / 27), tolerance);
 117    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(-0.5)), pow(T(1.5), -3), tolerance);
 118    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(-2)), T(1 / T(27)), tolerance);
 119    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(3) }, {}, T(-4)), T(T(1) / 125), tolerance);
 120    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(-0.5)), pow(T(1.5), 3), tolerance);
 121    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(-2)), T(27), tolerance);
 122    BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(-4)), T(125), tolerance);
 123
 124    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(3) }, {}, T(1)), std::domain_error);
 125    //BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(-3) }, {}, T(1)), std::domain_error);
 126    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(3.25) }, {}, T(1)), std::domain_error);
 127    //BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(-3.25) }, {}, T(1)), std::domain_error);
 128    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(3.25) }, {}, T(2)), std::domain_error);
 129    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ T(-3.25) }, {}, T(2)), std::domain_error);
 130 }
 131
 132 template <class T>
 133 void test_spots_0F1(T, const char*)
 134 {
 135    T tolerance = boost::math::tools::epsilon<T>() * 50000;
 136
 137    BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({},  { T(3) }, T(0)), 1);
 138    BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({}, { T(-3) }, T(0)), 1);
 139    //BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({}, { T(0) }, T(0)), 1);
 140
 141    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({}, { T(0) }, T(-1)), std::domain_error);
 142    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({}, { T(-1) }, T(-1)), std::domain_error);
 143    BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({}, { T(-10) }, T(-5)), std::domain_error);
 144
 145    static const boost::array<boost::array<T, 3>, 35> hypergeometric_pFq_integer_data = { {
 146       { SC_(4.0), SC_(-20.0),  SC_(-0.012889714201783047561923257996127233830940165138385) },
 147       { SC_(8.0), SC_(-20.0),  SC_(0.046498609282365144223175012935939437508273248399881) },
 148       { SC_(12.0), SC_(-20.0),  SC_(0.16608847431869756642136191351311569335145459224622) },
 149       { SC_(16.0), SC_(-20.0),  SC_(0.27230484709157170329168048388841880599105216477631) },
 150       //{ SC_(20.0), SC_(-20.0),  SC_(0.35865872656868844615709101792040025253126126604266) },
 151       { SC_(4.0), SC_(-16.0),  SC_(-0.027293644412433023379286103818840667403690937153604) },
 152       { SC_(8.0), SC_(-16.0),  SC_(0.098618710511372349330666801041676087431136532039702) },
 153       { SC_(12.0), SC_(-16.0),  SC_(0.24360114226383905073379763460037817885919457531523) },
 154       //{ SC_(16.0), SC_(-16.0),  SC_(0.35635186318802906043824855864337727878754460163525) },
 155       //{ SC_(20.0), SC_(-16.0),  SC_(0.44218381382689101428948260613085371477815110358789) },
 156       { SC_(4.0), SC_(-12.0),  SC_(-0.021743572290699436419371120781513860006290363262907) },
 157       { SC_(8.0), SC_(-12.0),  SC_(0.19025625754362006866949730683824627505504067855043) },
 158       //{ SC_(12.0), SC_(-12.0),  SC_(0.35251228238278927379621049815222218665165551016489) },
 159       //{ SC_(16.0), SC_(-12.0),  SC_(0.46415411486674623230458980010115972932474705884865) },
 160       //{ SC_(20.0), SC_(-12.0),  SC_(0.54394918325286018927327004362535051310016558628741) },
 161       { SC_(4.0), SC_(-8.0),  SC_(0.056818744289274872033266550620647787396712125304880) },
 162       //{ SC_(8.0), SC_(-8.0),  SC_(0.34487371876996263249797701802458885718691612997456) },
 163       //{ SC_(12.0), SC_(-8.0),  SC_(0.50411654015891701804499796523449656998841355305043) },
 164       //{ SC_(16.0), SC_(-8.0),  SC_(0.60191459981670594041254437708158847428118361245442) },
 165       //{ SC_(20.0), SC_(-8.0),  SC_(0.66770752550930138035694866478078941681114294465418) },
 166       //{ SC_(4.0), SC_(-4.0),  SC_(0.32262860540671645526863760914000166725449779629143) },
 167       //{ SC_(8.0), SC_(-4.0),  SC_(0.59755773349355150397404772151441126513126998265958) },
 168       //{ SC_(12.0), SC_(-4.0),  SC_(0.71337465206009117934071859694314971137807212605147) },
 169       //{ SC_(16.0), SC_(-4.0),  SC_(0.77734333649378860739496954157535257278092349684783) },
 170       //{ SC_(20.0), SC_(-4.0),  SC_(0.81794177985447769150469288350369205683856312760890) },
 171
 172       { SC_(4.0), SC_(4.0),  SC_(2.5029568338152582758923890008139391395035041790831) },
 173       { SC_(8.0), SC_(4.0),  SC_(1.6273673128576761227855719910743734060605725722129) },
 174       { SC_(12.0), SC_(4.0),  SC_(1.3898419290864057799739567227851793491657442624207) },
 175       { SC_(16.0), SC_(4.0),  SC_(1.2817098157957427946677711269410726972209834860612) },
 176       { SC_(20.0), SC_(4.0),  SC_(1.2202539302152377230940386181201477276788392792437) },
 177       { SC_(4.0), SC_(8.0),  SC_(5.5616961007411965409200003309686924059253894118586) },
 178       { SC_(8.0), SC_(8.0),  SC_(2.5877053985451664722152913482683136948296873738479) },
 179       { SC_(12.0), SC_(8.0),  SC_(1.9166410733572697158003086323981583993970490592046) },
 180       { SC_(16.0), SC_(8.0),  SC_(1.6370675016890669952237854163997946987362497613701) },
 181       { SC_(20.0), SC_(8.0),  SC_(1.4862852701827990444915220582410007454379891584086) },
 182       { SC_(4.0), SC_(12.0),  SC_(11.419268276211177842169936131590385979116019595164) },
 183       { SC_(8.0), SC_(12.0),  SC_(4.0347215359576567066789638314925802225312840819037) },
 184       { SC_(12.0), SC_(12.0),  SC_(2.6242497527837800417573064942486918368886996538285) },
 185       { SC_(16.0), SC_(12.0),  SC_(2.0840468784170876805932772732753387258909164486511) },
 186       { SC_(20.0), SC_(12.0),  SC_(1.8071042457762091748544382847762106786633952487005) },
 187       { SC_(4.0), SC_(16.0),  SC_(22.132051970576036053853444648907108439504682530918) },
 188       { SC_(8.0), SC_(16.0),  SC_(6.1850485247748975008808779795786699492711191898792) },
 189       { SC_(12.0), SC_(16.0),  SC_(3.5694322843488018916484224923627864928705138154372) },
 190       { SC_(16.0), SC_(16.0),  SC_(2.6447371137201451261118187672029372265909501355722) },
 191       { SC_(20.0), SC_(16.0),  SC_(2.1934058398888071720297525592515838555602675797235) },
 192       { SC_(4.0), SC_(20.0),  SC_(41.021743268279206331672552645354782698296383424328) },
 193       { SC_(8.0), SC_(20.0),  SC_(9.3414225299809886395081381945971250426599939097753) },
 194       { SC_(12.0), SC_(20.0),  SC_(4.8253866205826406499959001774187695527272168375992) },
 195       { SC_(16.0), SC_(20.0),  SC_(3.3462305133519485784864062004430532216764447939942) },
 196       { SC_(20.0), SC_(20.0),  SC_(2.6578698872220394617444624241257799193518140676691) },
 197       } };
 198
 199    for (auto row = hypergeometric_pFq_integer_data.begin(); row != hypergeometric_pFq_integer_data.end(); ++row)
 200    {
 201       BOOST_CHECK_CLOSE(boost::math::hypergeometric_pFq({},  { (*row)[0] }, (*row)[1]), (*row)[2], tolerance);
 202    }
 203 }
 204
 205 template <class T>
 206 void test_spots_1F1(T, const char*)
 207 {
 208 #include "hypergeometric_1F1.ipp"
 209
 210    for (auto row = hypergeometric_1F1.begin(); row != hypergeometric_1F1.end(); ++row)
 211    {
 212       try {
 213          T norm;
 214          T result = boost::math::hypergeometric_pFq({ (*row)[0] }, { (*row)[1] }, (*row)[2], &norm);
 215          check_pFq_result(result, norm, (*row)[3], { (*row)[0] }, { (*row)[1] }, (*row)[2]);
 216       }
 217       catch (const boost::math::evaluation_error&) {}
 218    }
 219 }
 220
 221 template <class T>
 222 void test_spots_1F1_b(T, const char*)
 223 {
 224 #include "hypergeometric_1F1_big.ipp"
 225
 226    for (auto row = hypergeometric_1F1_big.begin(); row != hypergeometric_1F1_big.end(); ++row)
 227    {
 228       try {
 229          T norm;
 230          T result = boost::math::hypergeometric_pFq({ (*row)[0] }, { (*row)[1] }, (*row)[2], &norm);
 231          check_pFq_result(result, norm, (*row)[3], { (*row)[0] }, { (*row)[1] }, (*row)[2]);
 232       }
 233       catch (const boost::math::evaluation_error&) {}
 234    }
 235 }
 236
 237 template <class T>
 238 void test_spots_2F1(T, const char*)
 239 {
 240 #include "hypergeometric_2F1.ipp"
 241
 242    for (auto row = hypergeometric_2F1.begin(); row != hypergeometric_2F1.end(); ++row)
 243    {
 244       try {
 245          T norm;
 246          T result = boost::math::hypergeometric_pFq({ (*row)[0], (*row)[1] }, { (*row)[2] }, (*row)[3], &norm);
 247          check_pFq_result(result, norm, (*row)[4], { (*row)[0], (*row)[1] }, { (*row)[2] }, (*row)[3]);
 248       }
 249       catch (const boost::math::evaluation_error&) {}
 250    }
 251 }
 252
 253 template <class T>
 254 void test_spots_0F2(T, const char*)
 255 {
 256 #include "hypergeometric_0F2.ipp"
 257
 258    for (auto row = hypergeometric_0F2.begin(); row != hypergeometric_0F2.end(); ++row)
 259    {
 260       try {
 261          T norm;
 262          T result = boost::math::hypergeometric_pFq({}, { (*row)[0], (*row)[1] }, (*row)[2], &norm);
 263          check_pFq_result(result, norm, (*row)[3], {}, { (*row)[0], (*row)[1] }, (*row)[2]);
 264       }
 265       catch (const boost::math::evaluation_error&) {}
 266    }
 267 }
 268
 269 template <class T>
 270 void test_spots_1F2(T, const char*)
 271 {
 272 #include "hypergeometric_1F2.ipp"
 273
 274    for (auto row = hypergeometric_1F2.begin(); row != hypergeometric_1F2.end(); ++row)
 275    {
 276       try {
 277          T norm;
 278          T result = boost::math::hypergeometric_pFq({ (*row)[0] }, { (*row)[1], (*row)[2] }, (*row)[3], &norm);
 279          check_pFq_result(result, norm, (*row)[4], { (*row)[0] }, { (*row)[1], (*row)[2] }, (*row)[3]);
 280       }
 281       catch (const boost::math::evaluation_error&) {}
 282    }
 283 }
 284
 285 template <class T>
 286 void test_spots_2F2(T, const char*)
 287 {
 288 #include "hypergeometric_2F2.ipp"
 289
 290    for (auto row = hypergeometric_2F2.begin(); row != hypergeometric_2F2.end(); ++row)
 291    {
 292       try {
 293          T norm;
 294          T result = boost::math::hypergeometric_pFq({ (*row)[0], (*row)[1] }, { (*row)[2], (*row)[3] }, (*row)[4], &norm);
 295          check_pFq_result(result, norm, (*row)[5], { (*row)[0], (*row)[1] }, { (*row)[2], (*row)[3] }, (*row)[4]);
 296       }
 297       catch (const boost::math::evaluation_error&) {}
 298    }
 299 }
 300
 301 template <class T>
 302 void test_spots(T z, const char* type_name)
 303 {
 304    test_spots_1F0(z, type_name);
 305    test_spots_0F1(z, type_name);
 306    test_spots_1F1(z, type_name);
 307    test_spots_1F1_b(z, type_name);
 308    test_spots_0F2(z, type_name);
 309    test_spots_1F2(z, type_name);
 310    test_spots_2F2(z, type_name);
 311    test_spots_2F1(z, type_name);
 312 }