ceph/src/boost/libs/math/tools/erf_data.cpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #include <boost/math/special_functions/gamma.hpp>
   7 #include <boost/math/special_functions/erf.hpp> // for inverses
   8 #include <boost/math/constants/constants.hpp>
   9 #include <fstream>
  10 #include <boost/math/tools/test_data.hpp>
  11 #include "mp_t.hpp"
  12
  13 using namespace boost::math::tools;
  14 using namespace std;
  15
  16 float external_f;
  17 float force_truncate(const float* f)
  18 {
  19    external_f = *f;
  20    return external_f;
  21 }
  22
  23 float truncate_to_float(mp_t r)
  24 {
  25    float f = boost::math::tools::real_cast<float>(r);
  26    return force_truncate(&f);
  27 }
  28
  29 struct erf_data_generator
  30 {
  31    boost::math::tuple<mp_t, mp_t> operator()(mp_t z)
  32    {
  33       // very naively calculate spots using the gamma function at high precision:
  34       int sign = 1;
  35       if(z < 0)
  36       {
  37          sign = -1;
  38          z = -z;
  39       }
  40       mp_t g1, g2;
  41       g1 = boost::math::tgamma_lower(mp_t(0.5), z * z);
  42       g1 /= sqrt(boost::math::constants::pi<mp_t>());
  43       g1 *= sign;
  44
  45       if(z < 0.5)
  46       {
  47          g2 = 1 - (sign * g1);
  48       }
  49       else
  50       {
  51          g2 = boost::math::tgamma(mp_t(0.5), z * z);
  52          g2 /= sqrt(boost::math::constants::pi<mp_t>());
  53       }
  54       if(sign < 1)
  55          g2 = 2 - g2;
  56       return boost::math::make_tuple(g1, g2);
  57    }
  58 };
  59
  60 double double_factorial(int N)
  61 {
  62    double result = 1;
  63    while(N > 2)
  64    {
  65       N -= 2;
  66       result *= N;
  67    }
  68    return result;
  69 }
  70
  71 void asymptotic_limit(int Bits)
  72 {
  73    //
  74    // The following block of code estimates how large z has
  75    // to be before we can use the asymptotic expansion for
  76    // erf/erfc and still get convergence: the series becomes
  77    // divergent eventually so we have to be careful!
  78    //
  79    double result = (std::numeric_limits<double>::max)();
  80    int terms = 0;
  81    for(int n = 1; n < 15; ++n)
  82    {
  83       double lim = (Bits-n) * log(2.0) - log(sqrt(3.14)) + log(double_factorial(2*n+1));
  84       double x = 1;
  85       while(x*x + (2*n+1)*log(x) <= lim)
  86          x += 0.1;
  87       if(x < result)
  88       {
  89          result = x;
  90          terms = n;
  91       }
  92    }
  93
  94    std::cout << "Erf asymptotic limit for "
  95       << Bits << " bit numbers is "
  96       << result << " after approximately "
  97       << terms << " terms." << std::endl;
  98
  99    result = (std::numeric_limits<double>::max)();
 100    terms = 0;
 101    for(int n = 1; n < 30; ++n)
 102    {
 103       double x = pow(double_factorial(2*n+1)/pow(2.0, n-Bits), 1 / (2.0*n));
 104       if(x < result)
 105       {
 106          result = x;
 107          terms = n;
 108       }
 109    }
 110
 111    std::cout << "Erfc asymptotic limit for "
 112       << Bits << " bit numbers is "
 113       << result << " after approximately "
 114       << terms << " terms." << std::endl;
 115 }
 116
 117 boost::math::tuple<mp_t, mp_t> erfc_inv(mp_t r)
 118 {
 119    mp_t x = exp(-r * r);
 120    x = x.convert_to<double>();
 121    std::cout << x << "   ";
 122    mp_t result = boost::math::erfc_inv(x);
 123    std::cout << result << std::endl;
 124    return boost::math::make_tuple(x, result);
 125 }
 126
 127
 128 int main(int argc, char*argv [])
 129 {
 130    parameter_info<mp_t> arg1;
 131    test_data<mp_t> data;
 132
 133    bool cont;
 134    std::string line;
 135
 136    if(argc >= 2)
 137    {
 138       if(strcmp(argv[1], "--limits") == 0)
 139       {
 140          asymptotic_limit(24);
 141          asymptotic_limit(53);
 142          asymptotic_limit(64);
 143          asymptotic_limit(106);
 144          asymptotic_limit(113);
 145          return 0;
 146       }
 147       else if(strcmp(argv[1], "--erf_inv") == 0)
 148       {
 149          mp_t (*f)(mp_t);
 150          f = boost::math::erf_inv;
 151          std::cout << "Welcome.\n"
 152             "This program will generate spot tests for the inverse erf function:\n";
 153          std::cout << "Enter the number of data points: ";
 154          int points;
 155          std::cin >> points;
 156          data.insert(f, make_random_param(mp_t(-1), mp_t(1), points));
 157       }
 158       else if(strcmp(argv[1], "--erfc_inv") == 0)
 159       {
 160          boost::math::tuple<mp_t, mp_t> (*f)(mp_t);
 161          f = erfc_inv;
 162          std::cout << "Welcome.\n"
 163             "This program will generate spot tests for the inverse erfc function:\n";
 164          std::cout << "Enter the maximum *result* expected from erfc_inv: ";
 165          double max_val;
 166          std::cin >> max_val;
 167          std::cout << "Enter the number of data points: ";
 168          int points;
 169          std::cin >> points;
 170          parameter_info<mp_t> arg = make_random_param(mp_t(0), mp_t(max_val), points);
 171          arg.type |= dummy_param;
 172          data.insert(f, arg);
 173       }
 174    }
 175    else
 176    {
 177       std::cout << "Welcome.\n"
 178          "This program will generate spot tests for the erf and erfc functions:\n"
 179          "  erf(z) and erfc(z)\n\n";
 180
 181       do{
 182          if(0 == get_user_parameter_info(arg1, "a"))
 183             return 1;
 184          data.insert(erf_data_generator(), arg1);
 185
 186          std::cout << "Any more data [y/n]?";
 187          std::getline(std::cin, line);
 188          boost::algorithm::trim(line);
 189          cont = (line == "y");
 190       }while(cont);
 191    }
 192
 193    std::cout << "Enter name of test data file [default=erf_data.ipp]";
 194    std::getline(std::cin, line);
 195    boost::algorithm::trim(line);
 196    if(line == "")
 197       line = "erf_data.ipp";
 198    std::ofstream ofs(line.c_str());
 199    ofs << std::scientific << std::setprecision(40);
 200    write_code(ofs, data, "erf_data");
 201
 202    return 0;
 203 }
 204
 205 /* Output for asymptotic limits:
 206
 207 Erf asymptotic limit for 24 bit numbers is 2.8 after approximately 6 terms.
 208 Erfc asymptotic limit for 24 bit numbers is 4.12064 after approximately 17 terms.
 209 Erf asymptotic limit for 53 bit numbers is 4.3 after approximately 11 terms.
 210 Erfc asymptotic limit for 53 bit numbers is 6.19035 after approximately 29 terms.
 211 Erf asymptotic limit for 64 bit numbers is 4.8 after approximately 12 terms.
 212 Erfc asymptotic limit for 64 bit numbers is 7.06004 after approximately 29 terms.
 213 Erf asymptotic limit for 106 bit numbers is 6.5 after approximately 14 terms.
 214 Erfc asymptotic limit for 106 bit numbers is 11.6626 after approximately 29 terms.
 215 Erf asymptotic limit for 113 bit numbers is 6.8 after approximately 14 terms.
 216 Erfc asymptotic limit for 113 bit numbers is 12.6802 after approximately 29 terms.
 217 */
 218