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/constants/constants.hpp>
8 #include <boost/lexical_cast.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
15 //
16 // Force truncation to float precision of input values:
17 // we must ensure that the input values are exactly representable
18 // in whatever type we are testing, or the output values will all
19 // be thrown off:
20 //
21 float external_f;
force_truncate(const float * f)22 float force_truncate(const float* f)
23 {
24 external_f = *f;
25 return external_f;
26 }
27
truncate_to_float(mp_t r)28 float truncate_to_float(mp_t r)
29 {
30 float f = boost::math::tools::real_cast<float>(r);
31 return force_truncate(&f);
32 }
33
34 //
35 // Our generator takes two arguments, but the second is interpreted
36 // as an instruction not a value, the second argument won't be placed
37 // in the output matrix by class test_data, so we add our algorithmically
38 // derived second argument to the output.
39 //
40 struct igamma_data_generator
41 {
operator ()igamma_data_generator42 boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t a, mp_t x)
43 {
44 // very naively calculate spots:
45 mp_t z;
46 switch((int)real_cast<float>(x))
47 {
48 case 1:
49 z = truncate_to_float((std::min)(mp_t(1), a/100));
50 break;
51 case 2:
52 z = truncate_to_float(a / 2);
53 break;
54 case 3:
55 z = truncate_to_float((std::max)(0.9*a, a - 2));
56 break;
57 case 4:
58 z = a;
59 break;
60 case 5:
61 z = truncate_to_float((std::min)(1.1*a, a + 2));
62 break;
63 case 6:
64 z = truncate_to_float(a * 2);
65 break;
66 case 7:
67 z = truncate_to_float((std::max)(mp_t(100), a*100));
68 break;
69 default:
70 BOOST_ASSERT(0 == "Can't get here!!");
71 }
72
73 //mp_t g = boost::math::tgamma(a);
74 mp_t lg = boost::math::tgamma_lower(a, z);
75 mp_t ug = boost::math::tgamma(a, z);
76 mp_t rlg = boost::math::gamma_p(a, z);
77 mp_t rug = boost::math::gamma_q(a, z);
78
79 return boost::math::make_tuple(z, ug, rug, lg, rlg);
80 }
81 };
82
83 struct gamma_inverse_generator
84 {
operator ()gamma_inverse_generator85 boost::math::tuple<mp_t, mp_t> operator()(const mp_t a, const mp_t p)
86 {
87 mp_t x1 = boost::math::gamma_p_inv(a, p);
88 mp_t x2 = boost::math::gamma_q_inv(a, p);
89 std::cout << "Inverse for " << a << " " << p << std::endl;
90 return boost::math::make_tuple(x1, x2);
91 }
92 };
93
94 struct gamma_inverse_generator_a
95 {
operator ()gamma_inverse_generator_a96 boost::math::tuple<mp_t, mp_t> operator()(const mp_t x, const mp_t p)
97 {
98 mp_t x1 = boost::math::gamma_p_inva(x, p);
99 mp_t x2 = boost::math::gamma_q_inva(x, p);
100 std::cout << "Inverse for " << x << " " << p << std::endl;
101 return boost::math::make_tuple(x1, x2);
102 }
103 };
104
105
main(int argc,char * argv[])106 int main(int argc, char*argv [])
107 {
108 bool cont;
109 std::string line;
110
111 parameter_info<mp_t> arg1, arg2;
112 test_data<mp_t> data;
113
114 if((argc >= 2) && (std::strcmp(argv[1], "-inverse") == 0))
115 {
116 std::cout << "Welcome.\n"
117 "This program will generate spot tests for the inverse incomplete gamma function:\n"
118 " gamma_p_inv(a, p) and gamma_q_inv(a, q)\n\n";
119 do{
120 get_user_parameter_info(arg1, "a");
121 get_user_parameter_info(arg2, "p");
122 data.insert(gamma_inverse_generator(), arg1, arg2);
123
124 std::cout << "Any more data [y/n]?";
125 std::getline(std::cin, line);
126 boost::algorithm::trim(line);
127 cont = (line == "y");
128 }while(cont);
129 }
130 else if((argc >= 2) && (std::strcmp(argv[1], "-inverse_a") == 0))
131 {
132 std::cout << "Welcome.\n"
133 "This program will generate spot tests for the inverse incomplete gamma function:\n"
134 " gamma_p_inva(a, p) and gamma_q_inva(a, q)\n\n";
135 do{
136 get_user_parameter_info(arg1, "x");
137 get_user_parameter_info(arg2, "p");
138 data.insert(gamma_inverse_generator_a(), arg1, arg2);
139
140 std::cout << "Any more data [y/n]?";
141 std::getline(std::cin, line);
142 boost::algorithm::trim(line);
143 cont = (line == "y");
144 }while(cont);
145 }
146 else
147 {
148 arg2 = make_periodic_param(mp_t(1), mp_t(8), 7);
149 arg2.type |= boost::math::tools::dummy_param;
150
151 std::cout << "Welcome.\n"
152 "This program will generate spot tests for the incomplete gamma function:\n"
153 " gamma(a, z)\n\n";
154
155 do{
156 get_user_parameter_info(arg1, "a");
157 data.insert(igamma_data_generator(), arg1, arg2);
158
159 std::cout << "Any more data [y/n]?";
160 std::getline(std::cin, line);
161 boost::algorithm::trim(line);
162 cont = (line == "y");
163 }while(cont);
164 }
165
166 std::cout << "Enter name of test data file [default=igamma_data.ipp]";
167 std::getline(std::cin, line);
168 boost::algorithm::trim(line);
169 if(line == "")
170 line = "igamma_data.ipp";
171 std::ofstream ofs(line.c_str());
172 ofs << std::scientific << std::setprecision(40);
173 write_code(ofs, data, "igamma_data");
174
175 return 0;
176 }
177
178