Lumiera  0.pre.03
»edit your freedom«
rational-test.cpp
Go to the documentation of this file.
1 /*
2  Rational(Test) - verify support for rational arithmetics
3 
4  Copyright (C)
5  2022, Hermann Vosseler <Ichthyostega@web.de>
6 
7   **Lumiera** is free software; you can redistribute it and/or modify it
8   under the terms of the GNU General Public License as published by the
9   Free Software Foundation; either version 2 of the License, or (at your
10   option) any later version. See the file COPYING for further details.
11 
12 * *****************************************************************/
13 
20 #include "lib/error.hpp"
21 #include "lib/test/run.hpp"
22 #include "lib/integral.hpp"
23 #include "lib/format-cout.hpp"
24 #include "lib/test/test-helper.hpp"
25 
26 #include "lib/rational.hpp"
27 
28 #include <chrono>
29 #include <array>
30 
31 using std::array;
32 
33 
34 namespace util {
35 namespace test {
36 
37 
38  /***************************************************************************/
46  class Rational_test : public Test
47  {
48 
49  virtual void
50  run (Arg)
51  {
52  demonstrate_basics();
53  verify_intLog2();
54  verify_limits();
55  verify_requant();
56  }
57 
58 
70  void
71  demonstrate_basics()
72  {
73  CHECK (Rat(10,3) == 10_r/3); // user-defined literal to construct a fraction
74  CHECK (Rat(10,3) == boost::rational<int64_t>(10,3)); // using Rat = boost::rational<int64_t>
75  CHECK (rational_cast<float> (10_r/3) == 3.3333333f); // rational_cast calculates division after type conversion
76 
77  CHECK (2_r/3 + 3_r/4 == 17_r/12);
78  CHECK (2_r/3 * 3_r/4 == 1_r/2);
79  CHECK (2_r/3 /(3_r/4)== 8_r/9);
80  CHECK (2_r/3 / 3 /4 == 1_r/18); // usual precedence and brace rules apply, yielding 2/36 here
81 
82  CHECK (util::toString(23_r/55) == "23/55sec"_expect);
83  CHECK (util::toString(24_r/56) == "3/7sec"_expect ); // rational numbers are normalised and reduced immediately
85 
86  CHECK (Rat(10,3).numerator() == int64_t(10));
87  CHECK (Rat(10,3).denominator() == int64_t(3));
88  CHECK (boost::rational<uint>(10,3).numerator() == uint(10));
89  CHECK (boost::rational<uint>(10,3).denominator() == uint(3));
90  CHECK (boost::rational<uint>(10,3) == rational_cast<boost::rational<uint>> (Rat(10,3)));
91  CHECK (boost::rational<uint>(11,3) != rational_cast<boost::rational<uint>> (Rat(10,3)));
92  }
93 
94 
101  void
102  verify_limits()
103  {
104  const Rat MAXI = Rat{std::numeric_limits<int64_t>::max()};
105  const Rat MINI = Rat{1, std::numeric_limits<int64_t>::max()};
106 
107  CHECK (rational_cast<int64_t>(MAXI) == std::numeric_limits<int64_t>::max());
108  CHECK (rational_cast<double> (MAXI) == 9.2233720368547758e+18);
109 
110  CHECK (MAXI > 0); // so this one still works
111  CHECK (MAXI+1 < 0); // but one more and we get a wrap-around
112  CHECK (MAXI+1 < -MAXI);
113  CHECK (util::toString(MAXI) == "9223372036854775807sec"_expect);
114  CHECK (util::toString(MAXI+1) == "-9223372036854775808sec"_expect);
115  CHECK (util::toString(-MAXI) == "-9223372036854775807sec"_expect);
116 
117  CHECK (MINI > 0); // smallest representable number above zero
118  CHECK (1-MINI < 1);
119  CHECK (0 < 1-MINI); // can be used below 1 just fine
120  CHECK (0 > 1+MINI); // but above we get a wrap-around in normalised numerator
121  CHECK (util::toString(MINI) == "1/9223372036854775807sec"_expect);
122  CHECK (util::toString(-MINI) == "-1/9223372036854775807sec"_expect);
123  CHECK (util::toString(1-MINI) == "9223372036854775806/9223372036854775807sec"_expect);
124  CHECK (util::toString(1+MINI) == "-9223372036854775808/9223372036854775807sec"_expect);
125 
126  CHECK ((MAXI-1)/MAXI == 1-MINI);
127  CHECK (MAXI/(MAXI-1) > 1); // as workaround we have to use a slightly larger ULP
128  CHECK (MAXI/(MAXI-1) - 1 > MINI); // ...this slightly larger one works without wrap-around
129  CHECK (1 - MAXI/(MAXI-1) < -MINI);
130  CHECK (util::toString(MAXI/(MAXI-1)) == "9223372036854775807/9223372036854775806sec"_expect);
131  CHECK (util::toString(MAXI/(MAXI-1) - 1) == "1/9223372036854775806sec"_expect);
132  CHECK (util::toString(1 - MAXI/(MAXI-1)) == "-1/9223372036854775806sec"_expect);
133 
134  // Now entering absolute danger territory....
135  const Rat MIMI = -MAXI-1; // this is the most extreme negative representable value
136  CHECK (MIMI < 0);
137  CHECK (util::toString(MIMI) == "-9223372036854775808sec"_expect);
138  CHECK (util::toString(1/MIMI) == "-1/-9223372036854775808sec"_expect);
139  try
140  {
141  -1-1/MIMI; // ...but it can't be used for any calculation without blowing up
142  NOTREACHED("expected boost::rational to flounder");
143  }
144  catch (std::exception& tilt)
145  {
146  CHECK (tilt.what() == string{"bad rational: non-zero singular denominator"});
147  }
148 
149 
150  // yet seemingly harmless values can be poisonous...
151  Rat poison = MAXI/49 / (MAXI/49-1);
152  Rat decoy = poison + 5;
153  CHECK (poison > 0);
154  CHECK (decoy > 6);
155  CHECK (rational_cast<double>(decoy) == 6); // looks innocuous...
156  CHECK (rational_cast<double>(decoy+5) == 11); // ...aaaand...
157  CHECK (rational_cast<double>(decoy+50) == -42); // ..ultimate answer..
158  CHECK (rational_cast<double>(decoy+500) == 15.999999999999996); // .dead in the water.
159 
160  // Heuristics to detect danger zone
161  CHECK ( can_represent_Sum(decoy,5));
162  CHECK (not can_represent_Sum(decoy,50));
163 
164  // alarm is given a bit too early
165  CHECK ( can_represent_Sum(decoy,15)); // ...because the check is based on bit positions
166  CHECK (not can_represent_Sum(decoy,16)); // ...and here the highest bit of one partner moved into danger zone
167  CHECK (decoy+16 > 0);
168  CHECK (decoy+43 > 0);
169  CHECK (decoy+44 < 0);
170 
171  // similar when increasing the denominator...
172  CHECK (poison + 1_r/10 > 0);
173  CHECK (poison + 1_r/90 > 0);
174  CHECK (poison + 1_r/98 < 0); // actually the flip already occurs at 1/91 but also causes an assertion failure
175  CHECK ( can_represent_Sum(poison,1_r/10));
176  CHECK ( can_represent_Sum(poison,1_r/15));
177  CHECK (not can_represent_Sum(poison,1_r/16));
178  CHECK (not can_represent_Sum(poison,1_r/91));
179  CHECK (not can_represent_Sum(poison,1_r/100));
180  }
181 
182 
190  void
191  verify_intLog2()
192  {
193  CHECK ( 5 == ilog2( int64_t(0b101010)));
194  CHECK ( 5 == ilog2(uint64_t(0b101010)));
195  CHECK ( 5 == ilog2( int32_t(0b101010)));
196  CHECK ( 5 == ilog2(uint32_t(0b101010)));
197  CHECK ( 5 == ilog2( int16_t(0b101010)));
198  CHECK ( 5 == ilog2(uint16_t(0b101010)));
199  CHECK ( 5 == ilog2( int8_t(0b101010)));
200  CHECK ( 5 == ilog2( uint8_t(0b101010)));
201  CHECK ( 5 == ilog2( int (0b101010)));
202  CHECK ( 5 == ilog2( uint (0b101010)));
203  CHECK ( 5 == ilog2( char (0b101010)));
204  CHECK ( 5 == ilog2( uchar (0b101010)));
205  CHECK ( 5 == ilog2( long (0b101010)));
206  CHECK ( 5 == ilog2( ulong (0b101010)));
207  CHECK ( 5 == ilog2( short (0b101010)));
208  CHECK ( 5 == ilog2( ushort (0b101010)));
209 
210  CHECK (63 == ilog2(std::numeric_limits<uint64_t>::max()));
211  CHECK (62 == ilog2(std::numeric_limits< int64_t>::max()));
212  CHECK (31 == ilog2(std::numeric_limits<uint32_t>::max()));
213  CHECK (30 == ilog2(std::numeric_limits< int32_t>::max()));
214  CHECK (15 == ilog2(std::numeric_limits<uint16_t>::max()));
215  CHECK (14 == ilog2(std::numeric_limits< int16_t>::max()));
216  CHECK ( 7 == ilog2(std::numeric_limits< uint8_t>::max()));
217  CHECK ( 6 == ilog2(std::numeric_limits< int8_t>::max()));
218 
219  CHECK ( 5 == ilog2(0b111111));
220  CHECK ( 5 == ilog2(0b101110));
221  CHECK ( 5 == ilog2(0b100100));
222  CHECK ( 5 == ilog2(0b100000));
223 
224  CHECK ( 2 == ilog2(4));
225  CHECK ( 1 == ilog2(2));
226  CHECK ( 0 == ilog2(1));
227  CHECK (-1 == ilog2(0));
228  CHECK (-1 == ilog2(-1));
229 
230  CHECK (-1 == ilog2(std::numeric_limits<uint64_t>::min()));
231  CHECK (-1 == ilog2(std::numeric_limits< int64_t>::min()));
232  CHECK (-1 == ilog2(std::numeric_limits<uint32_t>::min()));
233  CHECK (-1 == ilog2(std::numeric_limits< int32_t>::min()));
234  CHECK (-1 == ilog2(std::numeric_limits<uint16_t>::min()));
235  CHECK (-1 == ilog2(std::numeric_limits< int16_t>::min()));
236  CHECK (-1 == ilog2(std::numeric_limits< uint8_t>::min()));
237  CHECK (-1 == ilog2(std::numeric_limits< int8_t>::min()));
238 
239 
240  /* ==== compare with naive implementation ==== */
241 
242  auto floatLog = [](auto n)
243  {
244  return n <=0? -1 : ilogb(n);
245  };
246  auto bitshift = [](auto n)
247  {
248  if (n <= 0) return -1;
249  int logB = 0;
250  while (n >>= 1)
251  ++logB;
252  return logB;
253  };
254  auto do_nothing = [](auto n){ return n; };
255 
256  array<uint64_t, 1000> numz;
257  for (auto& n : numz)
258  {
259  n = rani() * uint64_t(1 << 31);
260  CHECK (ilog2(n) == floatLog(n));
261  CHECK (ilog2(n) == bitshift(n));
262  }
263 
264  int64_t dump{0}; // throw-away result to prevent optimisation
265  auto microbenchmark = [&numz,&dump](auto algo)
266  {
267  using std::chrono::system_clock;
268  using Dur = std::chrono::duration<double>;
269  const double SCALE = 1e9; // Results are in ns
270 
271  auto start = system_clock::now();
272  for (uint i=0; i<1000; ++i)
273  for (auto const& n : numz)
274  dump += algo(n);
275  Dur duration = system_clock::now () - start;
276  return duration.count()/(1000*1000) * SCALE;
277  };
278 
279 
280 
281  auto time_ilog2 = microbenchmark(ilog2<int64_t>);
282  auto time_float = microbenchmark(floatLog);
283  auto time_shift = microbenchmark(bitshift);
284  auto time_ident = microbenchmark(do_nothing);
285 
286  cout << "Microbenchmark integer-log2" <<endl
287  << "util::ilog2 :"<<time_ilog2<<"ns"<<endl
288  << "std::ilogb :"<<time_float<<"ns"<<endl
289  << "bit-shift :"<<time_shift<<"ns"<<endl
290  << "identity :"<<time_ident<<"ns"<<endl
291  << "(checksum="<<dump<<")" <<endl; // Warning: without outputting `dump`,
292  // the optimiser would eliminate most calls
293  // the following holds both with -O0 and -O3
294  CHECK (time_ilog2 < time_shift);
295  CHECK (time_ident < time_ilog2);
296 
297  /**** some numbers ****
298  *
299  * GCC-8, -O3, Debian-Buster, AMD FX83
300  *
301  * with uint64_t...
302  * - ilog2 : 5.6ns
303  * - ilogb : 5.0ns
304  * - shift : 44ns
305  * - ident : 0.6ns
306  *
307  * with uint8_t
308  * - ilog2 : 5.2ns
309  * - ilogb : 5.8ns
310  * - shift : 8.2ns
311  * - ident : 0.3ns
312  */
313  }
314 
315 
322  void
323  verify_requant()
324  {
325  const int64_t MAX{std::numeric_limits<int64_t>::max()};
326  const Rat MAXI = Rat{MAX};
327 
328  Rat poison = (MAXI-88)/(MAXI/7);
329 
330  auto approx = [](Rat rat){ return rational_cast<float> (rat); };
331  CHECK (poison > 0);
332  CHECK (poison+1 < 0); // wrap around!
333  CHECK (approx (poison ) == 6.99999952f); // wildly wrong results...
334  CHECK (approx (poison+1) == -6);
335  CHECK (approx (poison+7) == -6.83047369e-17f);
336  CHECK (approx (poison+9_r/5) == 0.400000006f);
337 
338  Rat sleazy = reQuant (poison, 1 << 24); // recast into multiples of an arbitrary other divisor,
339  CHECK (sleazy.denominator() == 1 << 24); // (here using a power of two as example)
340  // and now we can do all the slick stuff...
341  CHECK (sleazy > 0);
342  CHECK (sleazy+1 > 0);
343  CHECK (sleazy+7 > 0);
344  CHECK (approx (sleazy) == 7);
345  CHECK (approx (sleazy+1) == 8);
346  CHECK (approx (sleazy+7) == 14);
347  CHECK (approx (sleazy+9_r/5) == 8.80000019f);
348 
349  CHECK (util::toString (poison) == "9223372036854775719/1317624576693539401sec"_expect);
350  CHECK (util::toString (poison+1) =="-7905747460161236496/1317624576693539401sec"_expect);
351  CHECK (util::toString (sleazy) == "117440511/16777216sec"_expect);
352  CHECK (util::toString (sleazy+1) == "134217727/16777216sec"_expect);
353 
354  // also works towards larger denominator, or with negative numbers...
355  CHECK (reQuant (1/poison, MAX) == 1317624576693539413_r/9223372036854775807);
356  CHECK (reQuant (-poison, 7777) == -54438_r/ 7777);
357  CHECK (reQuant (poison, -7777) == -54438_r/-7777);
358 
359  CHECK (approx ( 1/poison ) == 0.142857149f);
360  CHECK (approx (reQuant (1/poison, MAX)) == 0.142857149f);
361  CHECK (approx (reQuant (poison, 7777)) == 6.99987125f);
362  }
363  };
364 
365  LAUNCHER (Rational_test, "unit common");
366 
367 
368 }} // namespace util::test
format-cout.hpp
Automatically use custom string conversion in C++ stream output.
util::test::Rational_test::verify_intLog2
void verify_intLog2()
Definition: rational-test.cpp:191
util::reQuant
int64_t reQuant(int64_t num, int64_t den, int64_t u)
Re-Quantise a number into a new grid, truncating to the next lower grid point.
Definition: rational.hpp:120
rational.hpp
Rational number support, based on boost::rational.
test
Definition: run.hpp:40
lib::rani
int rani(uint bound=_iBOUND())
Definition: random.hpp:135
util::test::Rational_test::verify_requant
void verify_requant()
Definition: rational-test.cpp:323
MAX
#define MAX(A, B)
the inevitable MAX macro, sometimes still necessary in template code
Definition: util.hpp:518
util::test::Rational_test::demonstrate_basics
void demonstrate_basics()
Definition: rational-test.cpp:71
run.hpp
Simplistic test class runner.
util::test::Rational_test
Definition: rational-test.cpp:46
util::toString
std::string toString(TY const &val) noexcept
get some string representation of any object, reliably.
Definition: format-obj.hpp:191
test-helper.hpp
A collection of frequently used helper functions to support unit testing.
error.hpp
Lumiera error handling (C++ interface).
util::test::Rational_test::verify_limits
void verify_limits()
Definition: rational-test.cpp:102
util::ilog2
constexpr int ilog2(I num)
Integral binary logarithm (disregarding fractional part)
Definition: util-quant.hpp:178