voronoi_robust_fpt_test.cpp 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405
  1. // Boost.Polygon library voronoi_robust_fpt_test.cpp file
  2. // Copyright Andrii Sydorchuk 2010-2012.
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // (See accompanying file LICENSE_1_0.txt or copy at
  5. // http://www.boost.org/LICENSE_1_0.txt)
  6. // See http://www.boost.org for updates, documentation, and revision history.
  7. #include <boost/core/lightweight_test.hpp>
  8. #include <boost/polygon/detail/voronoi_ctypes.hpp>
  9. #include <boost/polygon/detail/voronoi_robust_fpt.hpp>
  10. #include <boost/random/mersenne_twister.hpp>
  11. #include <vector>
  12. #include <cmath>
  13. #include <ctime>
  14. using boost::polygon::detail::int32;
  15. using boost::polygon::detail::uint32;
  16. using boost::polygon::detail::int64;
  17. using boost::polygon::detail::fpt64;
  18. using boost::polygon::detail::efpt64;
  19. using boost::polygon::detail::extended_int;
  20. using boost::polygon::detail::extended_exponent_fpt;
  21. using boost::polygon::detail::robust_fpt;
  22. using boost::polygon::detail::robust_dif;
  23. using boost::polygon::detail::robust_sqrt_expr;
  24. using boost::polygon::detail::type_converter_fpt;
  25. using boost::polygon::detail::type_converter_efpt;
  26. using boost::polygon::detail::ulp_comparison;
  27. typedef robust_fpt<double> rfpt_type;
  28. typedef type_converter_fpt to_fpt_type;
  29. typedef type_converter_efpt to_efpt_type;
  30. type_converter_fpt to_fpt;
  31. void robust_fpt_constructors_test1()
  32. {
  33. rfpt_type a = rfpt_type();
  34. BOOST_TEST_EQ(a.fpv(), 0.0);
  35. BOOST_TEST_EQ(a.re(), 0.0);
  36. BOOST_TEST_EQ(a.ulp(), 0);
  37. }
  38. void robust_fpt_constructors_test2()
  39. {
  40. rfpt_type a(10.0, 1.0);
  41. BOOST_TEST_EQ(a.fpv(), 10.0);
  42. BOOST_TEST_EQ(a.re(), 1.0);
  43. BOOST_TEST_EQ(a.ulp(), 1.0);
  44. }
  45. void robust_fpt_constructors_test3()
  46. {
  47. rfpt_type a(10.0);
  48. BOOST_TEST_EQ(a.fpv(), 10.0);
  49. BOOST_TEST_EQ(a.re(), 0.0);
  50. BOOST_TEST_EQ(a.ulp(), 0.0);
  51. }
  52. void robust_fpt_constructors_test4()
  53. {
  54. rfpt_type a(10.0, 3.0);
  55. BOOST_TEST_EQ(a.fpv(), 10.0);
  56. BOOST_TEST_EQ(a.re(), 3.0);
  57. BOOST_TEST_EQ(a.ulp(), 3.0);
  58. rfpt_type b(10.0, 2.75);
  59. BOOST_TEST_EQ(b.fpv(), 10.0);
  60. BOOST_TEST_EQ(b.re(), 2.75);
  61. BOOST_TEST_EQ(b.ulp(), 2.75);
  62. }
  63. void robust_fpt_sum_test1()
  64. {
  65. rfpt_type a(2.0, 5.0);
  66. rfpt_type b(3.0, 4.0);
  67. rfpt_type c = a + b;
  68. BOOST_TEST_EQ(c.fpv(), 5.0);
  69. BOOST_TEST_EQ(c.re(), 6.0);
  70. BOOST_TEST_EQ(c.ulp(), 6.0);
  71. c += b;
  72. BOOST_TEST_EQ(c.fpv(), 8.0);
  73. BOOST_TEST_EQ(c.re(), 7.0);
  74. BOOST_TEST_EQ(c.ulp(), 7.0);
  75. }
  76. void robust_fpt_sum_test2()
  77. {
  78. rfpt_type a(3.0, 2.0);
  79. rfpt_type b(-2.0, 3.0);
  80. rfpt_type c = a + b;
  81. BOOST_TEST_EQ(c.fpv(), 1.0);
  82. BOOST_TEST_EQ(c.re(), 13.0);
  83. BOOST_TEST_EQ(c.ulp(), 13.0);
  84. c += b;
  85. BOOST_TEST_EQ(c.fpv(), -1.0);
  86. BOOST_TEST_EQ(c.re(), 20.0);
  87. BOOST_TEST_EQ(c.ulp(), 20.0);
  88. }
  89. void robust_fpt_dif_test1()
  90. {
  91. rfpt_type a(2.0, 5.0);
  92. rfpt_type b(-3.0, 4.0);
  93. rfpt_type c = a - b;
  94. BOOST_TEST_EQ(c.fpv(), 5.0);
  95. BOOST_TEST_EQ(c.re(), 6.0);
  96. BOOST_TEST_EQ(c.ulp(), 6.0);
  97. c -= b;
  98. BOOST_TEST_EQ(c.fpv(), 8.0);
  99. BOOST_TEST_EQ(c.re(), 7.0);
  100. BOOST_TEST_EQ(c.ulp(), 7.0);
  101. }
  102. void robust_fpt_dif_test2()
  103. {
  104. rfpt_type a(3.0, 2.0);
  105. rfpt_type b(2.0, 3.0);
  106. rfpt_type c = a - b;
  107. BOOST_TEST_EQ(c.fpv(), 1.0);
  108. BOOST_TEST_EQ(c.re(), 13.0);
  109. BOOST_TEST_EQ(c.ulp(), 13.0);
  110. c -= b;
  111. BOOST_TEST_EQ(c.fpv(), -1.0);
  112. BOOST_TEST_EQ(c.re(), 20.0);
  113. BOOST_TEST_EQ(c.ulp(), 20.0);
  114. }
  115. void robust_fpt_mult_test3()
  116. {
  117. rfpt_type a(2.0, 3.0);
  118. rfpt_type b(4.0, 1.0);
  119. rfpt_type c = a * b;
  120. BOOST_TEST_EQ(c.fpv(), 8.0);
  121. BOOST_TEST_EQ(c.re(), 5.0);
  122. BOOST_TEST_EQ(c.ulp(), 5.0);
  123. c *= b;
  124. BOOST_TEST_EQ(c.fpv(), 32.0);
  125. BOOST_TEST_EQ(c.re(), 7.0);
  126. BOOST_TEST_EQ(c.ulp(), 7.0);
  127. }
  128. void robust_fpt_div_test1()
  129. {
  130. rfpt_type a(2.0, 3.0);
  131. rfpt_type b(4.0, 1.0);
  132. rfpt_type c = a / b;
  133. BOOST_TEST_EQ(c.fpv(), 0.5);
  134. BOOST_TEST_EQ(c.re(), 5.0);
  135. BOOST_TEST_EQ(c.ulp(), 5.0);
  136. c /= b;
  137. BOOST_TEST_EQ(c.fpv(), 0.125);
  138. BOOST_TEST_EQ(c.re(), 7.0);
  139. BOOST_TEST_EQ(c.ulp(), 7.0);
  140. }
  141. void robust_dif_constructors_test()
  142. {
  143. robust_dif<int> rd1;
  144. BOOST_TEST_EQ(rd1.pos(), 0);
  145. BOOST_TEST_EQ(rd1.neg(), 0);
  146. BOOST_TEST_EQ(rd1.dif(), 0);
  147. robust_dif<int> rd2(1);
  148. BOOST_TEST_EQ(rd2.pos(), 1);
  149. BOOST_TEST_EQ(rd2.neg(), 0);
  150. BOOST_TEST_EQ(rd2.dif(), 1);
  151. robust_dif<int> rd3(-1);
  152. BOOST_TEST_EQ(rd3.pos(), 0);
  153. BOOST_TEST_EQ(rd3.neg(), 1);
  154. BOOST_TEST_EQ(rd3.dif(), -1);
  155. robust_dif<int> rd4(1, 2);
  156. BOOST_TEST_EQ(rd4.pos(), 1);
  157. BOOST_TEST_EQ(rd4.neg(), 2);
  158. BOOST_TEST_EQ(rd4.dif(), -1);
  159. }
  160. void robust_dif_operators_test1()
  161. {
  162. robust_dif<int> a(5, 2), b(1, 10);
  163. int dif_a = a.dif();
  164. int dif_b = b.dif();
  165. robust_dif<int> sum = a + b;
  166. robust_dif<int> dif = a - b;
  167. robust_dif<int> mult = a * b;
  168. robust_dif<int> umin = -a;
  169. BOOST_TEST_EQ(sum.dif(), dif_a + dif_b);
  170. BOOST_TEST_EQ(dif.dif(), dif_a - dif_b);
  171. BOOST_TEST_EQ(mult.dif(), dif_a * dif_b);
  172. BOOST_TEST_EQ(umin.dif(), -dif_a);
  173. }
  174. void robust_dif_operators_test2()
  175. {
  176. robust_dif<int> a(5, 2);
  177. for (int b = -3; b <= 3; b += 6) {
  178. int dif_a = a.dif();
  179. int dif_b = b;
  180. robust_dif<int> sum = a + b;
  181. robust_dif<int> dif = a - b;
  182. robust_dif<int> mult = a * b;
  183. robust_dif<int> div = a / b;
  184. BOOST_TEST_EQ(sum.dif(), dif_a + dif_b);
  185. BOOST_TEST_EQ(dif.dif(), dif_a - dif_b);
  186. BOOST_TEST_EQ(mult.dif(), dif_a * dif_b);
  187. BOOST_TEST_EQ(div.dif(), dif_a / dif_b);
  188. }
  189. }
  190. void robust_dif_operators_test3()
  191. {
  192. robust_dif<int> b(5, 2);
  193. for (int a = -3; a <= 3; a += 6) {
  194. int dif_a = a;
  195. int dif_b = b.dif();
  196. robust_dif<int> sum = a + b;
  197. robust_dif<int> dif = a - b;
  198. robust_dif<int> mult = a * b;
  199. BOOST_TEST_EQ(sum.dif(), dif_a + dif_b);
  200. BOOST_TEST_EQ(dif.dif(), dif_a - dif_b);
  201. BOOST_TEST_EQ(mult.dif(), dif_a * dif_b);
  202. }
  203. }
  204. void robust_dif_operators_test4()
  205. {
  206. std::vector< robust_dif<int> > a4(4, robust_dif<int>(5, 2));
  207. std::vector< robust_dif<int> > b4(4, robust_dif<int>(1, 2));
  208. std::vector< robust_dif<int> > c4 = a4;
  209. c4[0] += b4[0];
  210. c4[1] -= b4[1];
  211. c4[2] *= b4[2];
  212. BOOST_TEST_EQ(c4[0].dif(), a4[0].dif() + b4[0].dif());
  213. BOOST_TEST_EQ(c4[1].dif(), a4[1].dif() - b4[1].dif());
  214. BOOST_TEST_EQ(c4[2].dif(), a4[2].dif() * b4[2].dif());
  215. a4[0] += b4[0].dif();
  216. a4[1] -= b4[1].dif();
  217. a4[2] *= b4[2].dif();
  218. a4[3] /= b4[3].dif();
  219. BOOST_TEST_EQ(c4[0].dif(), a4[0].dif());
  220. BOOST_TEST_EQ(c4[1].dif(), a4[1].dif());
  221. BOOST_TEST_EQ(c4[2].dif(), a4[2].dif());
  222. BOOST_TEST_EQ(c4[3].dif() / b4[3].dif(), a4[3].dif());
  223. }
  224. void robust_sqrt_expr_test1()
  225. {
  226. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  227. int32 A[1] = {10};
  228. int32 B[1] = {100};
  229. BOOST_TEST_EQ(sqrt_expr.eval1(A, B), 100.0);
  230. }
  231. void robust_sqrt_expr_test2()
  232. {
  233. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  234. int32 A[2] = {10, 30};
  235. int32 B[2] = {400, 100};
  236. BOOST_TEST_EQ(sqrt_expr.eval2(A, B), 500.0);
  237. }
  238. void robust_sqrt_expr_test3()
  239. {
  240. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  241. int32 A[2] = {10, -30};
  242. int32 B[2] = {400, 100};
  243. BOOST_TEST_EQ(sqrt_expr.eval2(A, B), -100.0);
  244. }
  245. void robust_sqrt_expr_test4()
  246. {
  247. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  248. int32 A[3] = {10, 30, 20};
  249. int32 B[3] = {4, 1, 9};
  250. BOOST_TEST_EQ(sqrt_expr.eval3(A, B), 110.0);
  251. }
  252. void robust_sqrt_expr_test5()
  253. {
  254. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  255. int32 A[3] = {10, 30, -20};
  256. int32 B[3] = {4, 1, 9};
  257. BOOST_TEST_EQ(sqrt_expr.eval3(A, B), -10.0);
  258. }
  259. void robust_sqrt_expr_test6()
  260. {
  261. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  262. int32 A[4] = {10, 30, 20, 5};
  263. int32 B[4] = {4, 1, 9, 16};
  264. BOOST_TEST_EQ(sqrt_expr.eval4(A, B), 130.0);
  265. }
  266. void robust_sqrt_expr_test7()
  267. {
  268. robust_sqrt_expr<int32, fpt64, to_fpt_type> sqrt_expr;
  269. int32 A[4] = {10, 30, -20, -5};
  270. int32 B[4] = {4, 1, 9, 16};
  271. BOOST_TEST_EQ(sqrt_expr.eval4(A, B), -30.0);
  272. }
  273. void robust_sqrt_expr_test8()
  274. {
  275. typedef extended_int<16> eint512;
  276. robust_sqrt_expr<eint512, efpt64, to_efpt_type> sqrt_expr;
  277. int32 A[4] = {1000, 3000, -2000, -500};
  278. int32 B[4] = {400, 100, 900, 1600};
  279. eint512 AA[4], BB[4];
  280. for (std::size_t i = 0; i < 4; ++i) {
  281. AA[i] = A[i];
  282. BB[i] = B[i];
  283. }
  284. BOOST_TEST_EQ(to_fpt(sqrt_expr.eval4(AA, BB)), -30000.0);
  285. }
  286. template <typename _int, typename _fpt>
  287. class sqrt_expr_tester {
  288. public:
  289. static const std::size_t MX_SQRTS = 4;
  290. bool run() {
  291. static boost::mt19937 gen(static_cast<uint32>(time(NULL)));
  292. bool ret_val = true;
  293. for (std::size_t i = 0; i < MX_SQRTS; ++i) {
  294. a[i] = gen() & 1048575;
  295. int64 temp = gen() & 1048575;
  296. b[i] = temp * temp;
  297. }
  298. uint32 mask = (1 << MX_SQRTS);
  299. for (std::size_t i = 0; i < mask; i++) {
  300. fpt64 expected_val = 0.0;
  301. for (std::size_t j = 0; j < MX_SQRTS; j++) {
  302. if (i & (1 << j)) {
  303. A[j] = a[j];
  304. B[j] = b[j];
  305. expected_val += static_cast<fpt64>(a[j]) *
  306. std::sqrt(static_cast<fpt64>(b[j]));
  307. } else {
  308. A[j] = -a[j];
  309. B[j] = b[j];
  310. expected_val -= static_cast<fpt64>(a[j]) *
  311. std::sqrt(static_cast<fpt64>(b[j]));
  312. }
  313. }
  314. fpt64 received_val = to_fpt(sqrt_expr_.eval4(A, B));
  315. ret_val &= ulp_cmp(expected_val, received_val, 25) ==
  316. ulp_comparison<fpt64>::EQUAL;
  317. }
  318. return ret_val;
  319. }
  320. private:
  321. robust_sqrt_expr<_int, _fpt, to_efpt_type> sqrt_expr_;
  322. ulp_comparison<fpt64> ulp_cmp;
  323. _int A[MX_SQRTS];
  324. _int B[MX_SQRTS];
  325. int64 a[MX_SQRTS];
  326. int64 b[MX_SQRTS];
  327. };
  328. void mpz_sqrt_evaluator_test()
  329. {
  330. typedef extended_int<16> eint512;
  331. sqrt_expr_tester<eint512, efpt64> tester;
  332. for (int i = 0; i < 2000; ++i)
  333. BOOST_TEST(tester.run());
  334. }
  335. int main()
  336. {
  337. robust_fpt_constructors_test1();
  338. robust_fpt_constructors_test2();
  339. robust_fpt_constructors_test3();
  340. robust_fpt_constructors_test4();
  341. robust_fpt_sum_test1();
  342. robust_fpt_sum_test2();
  343. robust_fpt_dif_test1();
  344. robust_fpt_dif_test2();
  345. robust_fpt_mult_test3();
  346. robust_fpt_div_test1();
  347. robust_dif_constructors_test();
  348. robust_dif_operators_test1();
  349. robust_dif_operators_test2();
  350. robust_dif_operators_test3();
  351. robust_dif_operators_test4();
  352. robust_sqrt_expr_test1();
  353. robust_sqrt_expr_test2();
  354. robust_sqrt_expr_test3();
  355. robust_sqrt_expr_test4();
  356. robust_sqrt_expr_test5();
  357. robust_sqrt_expr_test6();
  358. robust_sqrt_expr_test7();
  359. robust_sqrt_expr_test8();
  360. mpz_sqrt_evaluator_test();
  361. return boost::report_errors();
  362. }