]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. | |
4 | // (See accompanying file LICENSE_1_0.txt | |
5 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | ||
20effc67 | 7 | #include "mp_t.hpp" |
7c673cae FG |
8 | #include <boost/math/tools/test_data.hpp> |
9 | #include <boost/test/included/prg_exec_monitor.hpp> | |
10 | #include <boost/math/special_functions/ellint_rj.hpp> | |
11 | #include <boost/math/special_functions/ellint_rd.hpp> | |
12 | #include <fstream> | |
13 | #include <boost/math/tools/test_data.hpp> | |
14 | #include <boost/random.hpp> | |
7c673cae FG |
15 | |
16 | float extern_val; | |
17 | // confuse the compilers optimiser, and force a truncation to float precision: | |
18 | float truncate_to_float(float const * pf) | |
19 | { | |
20 | extern_val = *pf; | |
21 | return *pf; | |
22 | } | |
23 | ||
24 | // | |
25 | // Archived here is the original implementation of this | |
26 | // function by Xiaogang Zhang, we can use this to | |
27 | // generate special test cases for the new version: | |
28 | // | |
29 | template <typename T, typename Policy> | |
30 | T ellint_rj_old(T x, T y, T z, T p, const Policy& pol) | |
31 | { | |
32 | T value, u, lambda, alpha, beta, sigma, factor, tolerance; | |
33 | T X, Y, Z, P, EA, EB, EC, E2, E3, S1, S2, S3; | |
34 | unsigned long k; | |
35 | ||
36 | BOOST_MATH_STD_USING | |
37 | using namespace boost::math; | |
38 | ||
39 | static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)"; | |
40 | ||
41 | if(x < 0) | |
42 | { | |
43 | return policies::raise_domain_error<T>(function, | |
44 | "Argument x must be non-negative, but got x = %1%", x, pol); | |
45 | } | |
46 | if(y < 0) | |
47 | { | |
48 | return policies::raise_domain_error<T>(function, | |
49 | "Argument y must be non-negative, but got y = %1%", y, pol); | |
50 | } | |
51 | if(z < 0) | |
52 | { | |
53 | return policies::raise_domain_error<T>(function, | |
54 | "Argument z must be non-negative, but got z = %1%", z, pol); | |
55 | } | |
56 | if(p == 0) | |
57 | { | |
58 | return policies::raise_domain_error<T>(function, | |
59 | "Argument p must not be zero, but got p = %1%", p, pol); | |
60 | } | |
61 | if(x + y == 0 || y + z == 0 || z + x == 0) | |
62 | { | |
63 | return policies::raise_domain_error<T>(function, | |
64 | "At most one argument can be zero, " | |
65 | "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); | |
66 | } | |
67 | ||
68 | // error scales as the 6th power of tolerance | |
69 | tolerance = pow(T(1) * tools::epsilon<T>() / 3, T(1) / 6); | |
70 | ||
71 | // for p < 0, the integral is singular, return Cauchy principal value | |
72 | if(p < 0) | |
73 | { | |
74 | // | |
75 | // We must ensure that (z - y) * (y - x) is positive. | |
76 | // Since the integral is symmetrical in x, y and z | |
77 | // we can just permute the values: | |
78 | // | |
79 | if(x > y) | |
80 | std::swap(x, y); | |
81 | if(y > z) | |
82 | std::swap(y, z); | |
83 | if(x > y) | |
84 | std::swap(x, y); | |
85 | ||
86 | T q = -p; | |
87 | T pmy = (z - y) * (y - x) / (y + q); // p - y | |
88 | ||
89 | BOOST_ASSERT(pmy >= 0); | |
90 | ||
91 | p = pmy + y; | |
92 | value = ellint_rj_old(x, y, z, p, pol); | |
93 | value *= pmy; | |
94 | value -= 3 * boost::math::ellint_rf(x, y, z, pol); | |
95 | value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(x * z + p * q, p * q, pol); | |
96 | value /= (y + q); | |
97 | return value; | |
98 | } | |
99 | ||
100 | // duplication | |
101 | sigma = 0; | |
102 | factor = 1; | |
103 | k = 1; | |
104 | do | |
105 | { | |
106 | u = (x + y + z + p + p) / 5; | |
107 | X = (u - x) / u; | |
108 | Y = (u - y) / u; | |
109 | Z = (u - z) / u; | |
110 | P = (u - p) / u; | |
111 | ||
112 | if((tools::max)(abs(X), abs(Y), abs(Z), abs(P)) < tolerance) | |
113 | break; | |
114 | ||
115 | T sx = sqrt(x); | |
116 | T sy = sqrt(y); | |
117 | T sz = sqrt(z); | |
118 | ||
119 | lambda = sy * (sx + sz) + sz * sx; | |
120 | alpha = p * (sx + sy + sz) + sx * sy * sz; | |
121 | alpha *= alpha; | |
122 | beta = p * (p + lambda) * (p + lambda); | |
123 | sigma += factor * boost::math::ellint_rc(alpha, beta, pol); | |
124 | factor /= 4; | |
125 | x = (x + lambda) / 4; | |
126 | y = (y + lambda) / 4; | |
127 | z = (z + lambda) / 4; | |
128 | p = (p + lambda) / 4; | |
129 | ++k; | |
130 | } while(k < policies::get_max_series_iterations<Policy>()); | |
131 | ||
132 | // Check to see if we gave up too soon: | |
133 | policies::check_series_iterations<T>(function, k, pol); | |
134 | ||
135 | // Taylor series expansion to the 5th order | |
136 | EA = X * Y + Y * Z + Z * X; | |
137 | EB = X * Y * Z; | |
138 | EC = P * P; | |
139 | E2 = EA - 3 * EC; | |
140 | E3 = EB + 2 * P * (EA - EC); | |
141 | S1 = 1 + E2 * (E2 * T(9) / 88 - E3 * T(9) / 52 - T(3) / 14); | |
142 | S2 = EB * (T(1) / 6 + P * (T(-6) / 22 + P * T(3) / 26)); | |
143 | S3 = P * ((EA - EC) / 3 - P * EA * T(3) / 22); | |
144 | value = 3 * sigma + factor * (S1 + S2 + S3) / (u * sqrt(u)); | |
145 | ||
146 | return value; | |
147 | } | |
148 | ||
149 | template <typename T, typename Policy> | |
150 | T ellint_rd_imp_old(T x, T y, T z, const Policy& pol) | |
151 | { | |
152 | T value, u, lambda, sigma, factor, tolerance; | |
153 | T X, Y, Z, EA, EB, EC, ED, EE, S1, S2; | |
154 | unsigned long k; | |
155 | ||
156 | BOOST_MATH_STD_USING | |
157 | using namespace boost::math; | |
158 | ||
159 | static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)"; | |
160 | ||
161 | if(x < 0) | |
162 | { | |
163 | return policies::raise_domain_error<T>(function, | |
164 | "Argument x must be >= 0, but got %1%", x, pol); | |
165 | } | |
166 | if(y < 0) | |
167 | { | |
168 | return policies::raise_domain_error<T>(function, | |
169 | "Argument y must be >= 0, but got %1%", y, pol); | |
170 | } | |
171 | if(z <= 0) | |
172 | { | |
173 | return policies::raise_domain_error<T>(function, | |
174 | "Argument z must be > 0, but got %1%", z, pol); | |
175 | } | |
176 | if(x + y == 0) | |
177 | { | |
178 | return policies::raise_domain_error<T>(function, | |
179 | "At most one argument can be zero, but got, x + y = %1%", x + y, pol); | |
180 | } | |
181 | ||
182 | // error scales as the 6th power of tolerance | |
183 | tolerance = pow(tools::epsilon<T>() / 3, T(1) / 6); | |
184 | ||
185 | // duplication | |
186 | sigma = 0; | |
187 | factor = 1; | |
188 | k = 1; | |
189 | do | |
190 | { | |
191 | u = (x + y + z + z + z) / 5; | |
192 | X = (u - x) / u; | |
193 | Y = (u - y) / u; | |
194 | Z = (u - z) / u; | |
195 | if((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) | |
196 | break; | |
197 | T sx = sqrt(x); | |
198 | T sy = sqrt(y); | |
199 | T sz = sqrt(z); | |
200 | lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x); | |
201 | sigma += factor / (sz * (z + lambda)); | |
202 | factor /= 4; | |
203 | x = (x + lambda) / 4; | |
204 | y = (y + lambda) / 4; | |
205 | z = (z + lambda) / 4; | |
206 | ++k; | |
207 | } while(k < policies::get_max_series_iterations<Policy>()); | |
208 | ||
209 | // Check to see if we gave up too soon: | |
210 | policies::check_series_iterations<T>(function, k, pol); | |
211 | ||
212 | // Taylor series expansion to the 5th order | |
213 | EA = X * Y; | |
214 | EB = Z * Z; | |
215 | EC = EA - EB; | |
216 | ED = EA - 6 * EB; | |
217 | EE = ED + EC + EC; | |
218 | S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14); | |
219 | S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26)); | |
220 | value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u)); | |
221 | ||
222 | return value; | |
223 | } | |
224 | ||
225 | template <typename T, typename Policy> | |
226 | T ellint_rf_imp_old(T x, T y, T z, const Policy& pol) | |
227 | { | |
228 | T value, X, Y, Z, E2, E3, u, lambda, tolerance; | |
229 | unsigned long k; | |
230 | BOOST_MATH_STD_USING | |
231 | using namespace boost::math; | |
232 | static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; | |
233 | if(x < 0 || y < 0 || z < 0) | |
234 | { | |
235 | return policies::raise_domain_error<T>(function, | |
236 | "domain error, all arguments must be non-negative, " | |
237 | "only sensible result is %1%.", | |
238 | std::numeric_limits<T>::quiet_NaN(), pol); | |
239 | } | |
240 | if(x + y == 0 || y + z == 0 || z + x == 0) | |
241 | { | |
242 | return policies::raise_domain_error<T>(function, | |
243 | "domain error, at most one argument can be zero, " | |
244 | "only sensible result is %1%.", | |
245 | std::numeric_limits<T>::quiet_NaN(), pol); | |
246 | } | |
247 | // Carlson scales error as the 6th power of tolerance, | |
248 | // but this seems not to work for types larger than | |
249 | // 80-bit reals, this heuristic seems to work OK: | |
250 | if(policies::digits<T, Policy>() > 64) | |
251 | { | |
252 | tolerance = pow(tools::epsilon<T>(), T(1) / 4.25f); | |
253 | BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
254 | } | |
255 | else | |
256 | { | |
257 | tolerance = pow(4 * tools::epsilon<T>(), T(1) / 6); | |
258 | BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); | |
259 | } | |
260 | // duplication | |
261 | k = 1; | |
262 | do | |
263 | { | |
264 | u = (x + y + z) / 3; | |
265 | X = (u - x) / u; | |
266 | Y = (u - y) / u; | |
267 | Z = (u - z) / u; | |
268 | // Termination condition: | |
269 | if((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) | |
270 | break; | |
271 | T sx = sqrt(x); | |
272 | T sy = sqrt(y); | |
273 | T sz = sqrt(z); | |
274 | lambda = sy * (sx + sz) + sz * sx; | |
275 | x = (x + lambda) / 4; | |
276 | y = (y + lambda) / 4; | |
277 | z = (z + lambda) / 4; | |
278 | ++k; | |
279 | } while(k < policies::get_max_series_iterations<Policy>()); | |
280 | // Check to see if we gave up too soon: | |
281 | policies::check_series_iterations<T>(function, k, pol); | |
282 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
283 | // Taylor series expansion to the 5th order | |
284 | E2 = X * Y - Z * Z; | |
285 | E3 = X * Y * Z; | |
286 | value = (1 + E2*(E2 / 24 - E3*T(3) / 44 - T(0.1)) + E3 / 14) / sqrt(u); | |
287 | BOOST_MATH_INSTRUMENT_VARIABLE(value); | |
288 | return value; | |
289 | } | |
290 | ||
291 | ||
292 | ||
293 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rj_data_4e(mp_t n) | |
294 | { | |
295 | mp_t result = ellint_rj_old(n, n, n, n, boost::math::policies::policy<>()); | |
296 | return boost::math::make_tuple(n, n, n, result); | |
297 | } | |
298 | ||
299 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_3e(mp_t x, mp_t p) | |
300 | { | |
301 | mp_t r = ellint_rj_old(x, x, x, p, boost::math::policies::policy<>()); | |
302 | return boost::math::make_tuple(x, x, x, p, r); | |
303 | } | |
304 | ||
305 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_1(mp_t x, mp_t y, mp_t p) | |
306 | { | |
307 | mp_t r = ellint_rj_old(x, x, y, p, boost::math::policies::policy<>()); | |
308 | return boost::math::make_tuple(x, x, y, p, r); | |
309 | } | |
310 | ||
311 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_2(mp_t x, mp_t y, mp_t p) | |
312 | { | |
313 | mp_t r = ellint_rj_old(x, y, x, p, boost::math::policies::policy<>()); | |
314 | return boost::math::make_tuple(x, y, x, p, r); | |
315 | } | |
316 | ||
317 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_3(mp_t x, mp_t y, mp_t p) | |
318 | { | |
319 | mp_t r = ellint_rj_old(y, x, x, p, boost::math::policies::policy<>()); | |
320 | return boost::math::make_tuple(y, x, x, p, r); | |
321 | } | |
322 | ||
323 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data_2e_4(mp_t x, mp_t y, mp_t p) | |
324 | { | |
325 | mp_t r = ellint_rj_old(x, y, p, p, boost::math::policies::policy<>()); | |
326 | return boost::math::make_tuple(x, y, p, p, r); | |
327 | } | |
328 | ||
329 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_1(mp_t x, mp_t y) | |
330 | { | |
331 | mp_t r = ellint_rd_imp_old(x, y, y, boost::math::policies::policy<>()); | |
332 | return boost::math::make_tuple(x, y, y, r); | |
333 | } | |
334 | ||
335 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_2(mp_t x, mp_t y) | |
336 | { | |
337 | mp_t r = ellint_rd_imp_old(x, x, y, boost::math::policies::policy<>()); | |
338 | return boost::math::make_tuple(x, x, y, r); | |
339 | } | |
340 | ||
341 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_2e_3(mp_t x) | |
342 | { | |
343 | mp_t r = ellint_rd_imp_old(mp_t(0), x, x, boost::math::policies::policy<>()); | |
344 | return boost::math::make_tuple(0, x, x, r); | |
345 | } | |
346 | ||
347 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_3e(mp_t x) | |
348 | { | |
349 | mp_t r = ellint_rd_imp_old(x, x, x, boost::math::policies::policy<>()); | |
350 | return boost::math::make_tuple(x, x, x, r); | |
351 | } | |
352 | ||
353 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data_0xy(mp_t x, mp_t y) | |
354 | { | |
355 | mp_t r = ellint_rd_imp_old(mp_t(0), x, y, boost::math::policies::policy<>()); | |
356 | return boost::math::make_tuple(mp_t(0), x, y, r); | |
357 | } | |
358 | ||
359 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xxx(mp_t x) | |
360 | { | |
361 | mp_t r = ellint_rf_imp_old(x, x, x, boost::math::policies::policy<>()); | |
362 | return boost::math::make_tuple(x, x, x, r); | |
363 | } | |
364 | ||
365 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xyy(mp_t x, mp_t y) | |
366 | { | |
367 | mp_t r = ellint_rf_imp_old(x, y, y, boost::math::policies::policy<>()); | |
368 | return boost::math::make_tuple(x, y, y, r); | |
369 | } | |
370 | ||
371 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xxy(mp_t x, mp_t y) | |
372 | { | |
373 | mp_t r = ellint_rf_imp_old(x, x, y, boost::math::policies::policy<>()); | |
374 | return boost::math::make_tuple(x, x, y, r); | |
375 | } | |
376 | ||
377 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xyx(mp_t x, mp_t y) | |
378 | { | |
379 | mp_t r = ellint_rf_imp_old(x, y, x, boost::math::policies::policy<>()); | |
380 | return boost::math::make_tuple(x, y, x, r); | |
381 | } | |
382 | ||
383 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_0yy(mp_t y) | |
384 | { | |
385 | mp_t r = ellint_rf_imp_old(mp_t(0), y, y, boost::math::policies::policy<>()); | |
386 | return boost::math::make_tuple(mp_t(0), y, y, r); | |
387 | } | |
388 | ||
389 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data_xy0(mp_t x, mp_t y) | |
390 | { | |
391 | mp_t r = ellint_rf_imp_old(x, y, mp_t(0), boost::math::policies::policy<>()); | |
392 | return boost::math::make_tuple(x, y, mp_t(0), r); | |
393 | } | |
394 | ||
395 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rf_data(mp_t n) | |
396 | { | |
397 | static boost::mt19937 r; | |
398 | boost::uniform_real<float> ur(0, 1); | |
399 | boost::uniform_int<int> ui(-100, 100); | |
400 | float x = ur(r); | |
401 | x = ldexp(x, ui(r)); | |
402 | mp_t xr(truncate_to_float(&x)); | |
403 | float y = ur(r); | |
404 | y = ldexp(y, ui(r)); | |
405 | mp_t yr(truncate_to_float(&y)); | |
406 | float z = ur(r); | |
407 | z = ldexp(z, ui(r)); | |
408 | mp_t zr(truncate_to_float(&z)); | |
409 | ||
410 | mp_t result = boost::math::ellint_rf(xr, yr, zr); | |
411 | return boost::math::make_tuple(xr, yr, zr, result); | |
412 | } | |
413 | ||
414 | boost::math::tuple<mp_t, mp_t, mp_t> generate_rc_data(mp_t n) | |
415 | { | |
416 | static boost::mt19937 r; | |
417 | boost::uniform_real<float> ur(0, 1); | |
418 | boost::uniform_int<int> ui(-100, 100); | |
419 | float x = ur(r); | |
420 | x = ldexp(x, ui(r)); | |
421 | mp_t xr(truncate_to_float(&x)); | |
422 | float y = ur(r); | |
423 | y = ldexp(y, ui(r)); | |
424 | mp_t yr(truncate_to_float(&y)); | |
425 | ||
426 | mp_t result = boost::math::ellint_rc(xr, yr); | |
427 | return boost::math::make_tuple(xr, yr, result); | |
428 | } | |
429 | ||
430 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t> generate_rj_data(mp_t n) | |
431 | { | |
432 | static boost::mt19937 r; | |
433 | boost::uniform_real<float> ur(0, 1); | |
434 | boost::uniform_real<float> nur(-1, 1); | |
435 | boost::uniform_int<int> ui(-100, 100); | |
436 | float x = ur(r); | |
437 | x = ldexp(x, ui(r)); | |
438 | mp_t xr(truncate_to_float(&x)); | |
439 | float y = ur(r); | |
440 | y = ldexp(y, ui(r)); | |
441 | mp_t yr(truncate_to_float(&y)); | |
442 | float z = ur(r); | |
443 | z = ldexp(z, ui(r)); | |
444 | mp_t zr(truncate_to_float(&z)); | |
445 | float p = nur(r); | |
446 | p = ldexp(p, ui(r)); | |
447 | mp_t pr(truncate_to_float(&p)); | |
448 | ||
449 | boost::math::ellint_rj(x, y, z, p); | |
450 | ||
451 | mp_t result = boost::math::ellint_rj(xr, yr, zr, pr); | |
452 | return boost::math::make_tuple(xr, yr, zr, pr, result); | |
453 | } | |
454 | ||
455 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rd_data(mp_t n) | |
456 | { | |
457 | static boost::mt19937 r; | |
458 | boost::uniform_real<float> ur(0, 1); | |
459 | boost::uniform_int<int> ui(-100, 100); | |
460 | float x = ur(r); | |
461 | x = ldexp(x, ui(r)); | |
462 | mp_t xr(truncate_to_float(&x)); | |
463 | float y = ur(r); | |
464 | y = ldexp(y, ui(r)); | |
465 | mp_t yr(truncate_to_float(&y)); | |
466 | float z = ur(r); | |
467 | z = ldexp(z, ui(r)); | |
468 | mp_t zr(truncate_to_float(&z)); | |
469 | ||
470 | mp_t result = boost::math::ellint_rd(xr, yr, zr); | |
471 | return boost::math::make_tuple(xr, yr, zr, result); | |
472 | } | |
473 | ||
474 | mp_t rg_imp(mp_t x, mp_t y, mp_t z) | |
475 | { | |
476 | using std::swap; | |
477 | // If z is zero permute so the call to RD is valid: | |
478 | if(z == 0) | |
479 | swap(x, z); | |
480 | return (z * ellint_rf_imp_old(x, y, z, boost::math::policies::policy<>()) | |
481 | - (x - z) * (y - z) * ellint_rd_imp_old(x, y, z, boost::math::policies::policy<>()) / 3 | |
482 | + sqrt(x * y / z)) / 2; | |
483 | } | |
484 | ||
485 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_data(mp_t n) | |
486 | { | |
487 | static boost::mt19937 r; | |
488 | boost::uniform_real<float> ur(0, 1); | |
489 | boost::uniform_int<int> ui(-100, 100); | |
490 | float x = ur(r); | |
491 | x = ldexp(x, ui(r)); | |
492 | mp_t xr(truncate_to_float(&x)); | |
493 | float y = ur(r); | |
494 | y = ldexp(y, ui(r)); | |
495 | mp_t yr(truncate_to_float(&y)); | |
496 | float z = ur(r); | |
497 | z = ldexp(z, ui(r)); | |
498 | mp_t zr(truncate_to_float(&z)); | |
499 | ||
500 | mp_t result = rg_imp(xr, yr, zr); | |
501 | return boost::math::make_tuple(xr, yr, zr, result); | |
502 | } | |
503 | ||
504 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xxx(mp_t x) | |
505 | { | |
506 | mp_t result = rg_imp(x, x, x); | |
507 | return boost::math::make_tuple(x, x, x, result); | |
508 | } | |
509 | ||
510 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xyy(mp_t x, mp_t y) | |
511 | { | |
512 | mp_t result = rg_imp(x, y, y); | |
513 | return boost::math::make_tuple(x, y, y, result); | |
514 | } | |
515 | ||
516 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xxy(mp_t x, mp_t y) | |
517 | { | |
518 | mp_t result = rg_imp(x, x, y); | |
519 | return boost::math::make_tuple(x, x, y, result); | |
520 | } | |
521 | ||
522 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xyx(mp_t x, mp_t y) | |
523 | { | |
524 | mp_t result = rg_imp(x, y, x); | |
525 | return boost::math::make_tuple(x, y, x, result); | |
526 | } | |
527 | ||
528 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_0xx(mp_t x) | |
529 | { | |
530 | mp_t result = rg_imp(mp_t(0), x, x); | |
531 | return boost::math::make_tuple(mp_t(0), x, x, result); | |
532 | } | |
533 | ||
534 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_x0x(mp_t x) | |
535 | { | |
536 | mp_t result = rg_imp(x, mp_t(0), x); | |
537 | return boost::math::make_tuple(x, mp_t(0), x, result); | |
538 | } | |
539 | ||
540 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xx0(mp_t x) | |
541 | { | |
542 | mp_t result = rg_imp(x, x, mp_t(0)); | |
543 | return boost::math::make_tuple(x, x, mp_t(0), result); | |
544 | } | |
545 | ||
546 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_00x(mp_t x) | |
547 | { | |
548 | mp_t result = sqrt(x) / 2; | |
549 | return boost::math::make_tuple(mp_t(0), mp_t(0), x, result); | |
550 | } | |
551 | ||
552 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_0x0(mp_t x) | |
553 | { | |
554 | mp_t result = sqrt(x) / 2; | |
555 | return boost::math::make_tuple(mp_t(0), x, mp_t(0), result); | |
556 | } | |
557 | ||
558 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_x00(mp_t x) | |
559 | { | |
560 | mp_t result = sqrt(x) / 2; | |
561 | return boost::math::make_tuple(x, mp_t(0), mp_t(0), result); | |
562 | } | |
563 | ||
564 | boost::math::tuple<mp_t, mp_t, mp_t, mp_t> generate_rg_xy0(mp_t x, mp_t y) | |
565 | { | |
566 | mp_t result = rg_imp(x, y, mp_t(0)); | |
567 | return boost::math::make_tuple(x, y, mp_t(0), result); | |
568 | } | |
569 | ||
570 | int cpp_main(int argc, char*argv[]) | |
571 | { | |
572 | using namespace boost::math::tools; | |
573 | ||
574 | parameter_info<mp_t> arg1, arg2, arg3; | |
575 | test_data<mp_t> data; | |
576 | ||
577 | bool cont; | |
578 | std::string line; | |
579 | ||
580 | if(argc < 1) | |
581 | return 1; | |
582 | ||
583 | do{ | |
584 | #if 0 | |
585 | int count; | |
586 | std::cout << "Number of points: "; | |
587 | std::cin >> count; | |
588 | ||
589 | arg1 = make_periodic_param(mp_t(0), mp_t(1), count); | |
590 | arg1.type |= dummy_param; | |
591 | ||
592 | // | |
593 | // Change this next line to get the R variant you want: | |
594 | // | |
595 | data.insert(&generate_rd_data, arg1); | |
596 | ||
597 | std::cout << "Any more data [y/n]?"; | |
598 | std::getline(std::cin, line); | |
599 | boost::algorithm::trim(line); | |
600 | cont = (line == "y"); | |
601 | #else | |
602 | get_user_parameter_info(arg1, "x"); | |
603 | get_user_parameter_info(arg2, "y"); | |
604 | //get_user_parameter_info(arg3, "p"); | |
605 | arg1.type |= dummy_param; | |
606 | arg2.type |= dummy_param; | |
607 | //arg3.type |= dummy_param; | |
608 | data.insert(generate_rd_data_0xy, arg1, arg2); | |
609 | ||
610 | std::cout << "Any more data [y/n]?"; | |
611 | std::getline(std::cin, line); | |
612 | boost::algorithm::trim(line); | |
613 | cont = (line == "y"); | |
614 | #endif | |
615 | }while(cont); | |
616 | ||
617 | std::cout << "Enter name of test data file [default=ellint_rf_data.ipp]"; | |
618 | std::getline(std::cin, line); | |
619 | boost::algorithm::trim(line); | |
620 | if(line == "") | |
621 | line = "ellint_rf_data.ipp"; | |
622 | std::ofstream ofs(line.c_str()); | |
623 | line.erase(line.find('.')); | |
624 | ofs << std::scientific << std::setprecision(40); | |
625 | write_code(ofs, data, line.c_str()); | |
626 | ||
627 | return 0; | |
628 | } | |
629 | ||
630 |