@@ -363,207 +363,3 @@ void truncated_gauss_t::calculate_cdf()
363363}
364364
365365} // rng
366- #ifdef UNIT_TEST
367- // Code to test functionality and performance of our RNG implementations
368-
369- #include < random>
370- #include < tuple>
371-
372- #include " lib/fmt/format.h"
373- #include " util/generic.hpp"
374- #include " util/chrono.hpp"
375-
376- using test_clock = chrono::cpu_clock;
377-
378- template <typename Engine>
379- static void test_one ( rng::basic_rng_t <Engine>& rng, uint64_t n )
380- {
381- auto start_time = test_clock::now ();
382-
383- double average = 0 ;
384- for ( uint64_t i = 0 ; i < n; ++i )
385- {
386- double d = rng.real ();
387- // std::cout << d << "\n";
388- average += d;
389- }
390-
391- average /= n;
392- auto elapsed_cpu = chrono::elapsed_fp_seconds (start_time);
393-
394- fmt::print ( " {} calls to rng::{}::real(), average = {:.8f}, time = {}s, numbers/sec = {}\n\n " ,
395- n, rng.name (), average, elapsed_cpu, static_cast <uint64_t >( n * 1000.0 / elapsed_cpu ) );
396- }
397-
398- template <typename Engine>
399- static void test_seed ( rng::basic_rng_t <Engine>& rng, uint64_t n )
400- {
401- auto start_time = test_clock::now ();
402-
403- for ( uint64_t i = 0 ; i < n; ++i )
404- {
405- rng.seed ( uint64_t ( m_pi ) * n );
406- }
407-
408- auto elapsed_cpu = chrono::elapsed_fp_seconds (start_time);
409-
410- fmt::print ( " {} calls to rng::{}::seed(), time = {}s, numbers/sec = {}\n\n " ,
411- n, rng.name (), elapsed_cpu, static_cast <uint64_t >( n * 1000.0 / elapsed_cpu ) );
412- }
413-
414- // Monte-Carlo PI calculation.
415- template <typename Engine>
416- static void monte_carlo ( rng::basic_rng_t <Engine>& rng, uint64_t n )
417- {
418- auto start_time = test_clock::now ();
419-
420- uint64_t count_in_sphere = 0 ;
421- for ( uint64_t i = 0 ; i < n; ++i )
422- {
423- double x1 = rng.real ();
424- double x2 = rng.real ();
425- if ( x1 * x1 + x2 * x2 < 1.0 )
426- ++count_in_sphere;
427- }
428- auto elapsed_cpu = chrono::elapsed_fp_seconds (start_time);
429-
430- fmt::print ( " {} runs for pi-calculation with {}. count in sphere = {}, pi = {:.21f}, time = {}s, numbers/sec = {}\n\n " ,
431- n, rng.name (), count_in_sphere, static_cast <double >( count_in_sphere ) * 4 / n,
432- elapsed_cpu, static_cast <uint64_t >( n * 1000.0 / elapsed_cpu ) );
433- }
434-
435- template <typename Engine>
436- static void test_uniform_int ( rng::basic_rng_t <Engine>& rng, uint64_t n, unsigned num_buckets )
437- {
438- auto start_time = test_clock::now ();
439-
440- std::vector<unsigned > histogram (num_buckets);
441-
442- for ( uint64_t i = 0 ; i < n; ++i )
443- {
444- auto result = rng.range (histogram.size ());
445- histogram[ result ] += 1 ;
446- }
447-
448- auto elapsed_cpu = chrono::elapsed_fp_seconds (start_time);
449-
450- double expected_bucket_size = static_cast <double >(n) / histogram.size ();
451-
452- fmt::print (" {} call to rng::{}::range(size_t({}u))\n " , n, rng.name (), histogram.size ());
453- for (unsigned i = 0 ; i < histogram.size (); ++i)
454- {
455- double pct = static_cast <double >(histogram[ i ]) / n;
456- double diff = static_cast <double >(histogram[ i ]) / expected_bucket_size - 1.0 ;
457- fmt::print (" bucket {:2}: {:5.2f}% ({}) difference to expected: {:9.6f}%\n " , i, pct, histogram[ i ], diff);
458- }
459- fmt::print (" time = {} s\n\n " , elapsed_cpu);
460- }
461-
462- namespace detail {
463- template <typename Tuple, typename F, std::size_t ... I>
464- void for_each_impl (Tuple&& t, F&& f, std::index_sequence<I...>)
465- {
466- int _[] = { ( std::invoke ( std::forward<F>( f ), std::get<I>( std::forward<Tuple>(t) ) ), 0 )... };
467- (void )_;
468- }
469- }
470-
471- template <typename Tuple, typename F>
472- void for_each (Tuple&& t, F&& f)
473- {
474- return detail::for_each_impl ( std::forward<Tuple>(t), std::forward<F>(f),
475- std::make_index_sequence<std::tuple_size<std::remove_reference_t <Tuple>>::value>{});
476- }
477-
478- int main ( int /* argc*/ , char ** /* argv*/ )
479- {
480- std::tuple<
481- rng::basic_rng_t <rng::xoshiro256plus_t >,
482- rng::basic_rng_t <rng::xorshift128_t >,
483- rng::basic_rng_t <rng::xorshift1024_t >
484- > generators;
485-
486- std::random_device rd;
487- uint64_t seed = uint64_t (rd ()) | (uint64_t (rd ()) << 32 );
488- fmt::print ( " Seed: {}\n\n " , seed );
489-
490- for_each ( generators, [seed]( auto & g ) { g.seed ( seed ); } );
491-
492- for_each ( generators, []( auto & g ) { test_one ( g, 1'000'000'000 ); } );
493-
494- for_each ( generators, []( auto & g ) { monte_carlo ( g, 1'000'000'000 ); } );
495-
496- for_each ( generators, []( auto & g ) { test_seed ( g, 10'000'000 ); } );
497-
498- for_each ( generators, []( auto & g ) { test_uniform_int ( g, 10'000 , 10 ); } );
499-
500- fmt::print ( " random device: min={} max={}\n\n " , rd.min (), rd.max () );
501-
502- auto & rng = std::get<0 >( generators );
503- fmt::print ( " Testing {}\n\n " , rng.name () );
504-
505- const uint64_t n = 100'000'000 ;
506-
507- // double gauss
508- {
509- auto start_time = test_clock::now ();
510-
511- double average = 0 ;
512- for ( uint64_t i = 0 ; i < n; i++ )
513- average += rng.gauss ( 0 , 1 );
514- average /= n;
515-
516- fmt::print ( " {} calls to gauss(0, 1): average = {:.8f}, time = {}s\n\n " ,
517- n, average, chrono::elapsed_fp_seconds (start_time) );
518- }
519-
520- // double exgauss
521- {
522- auto start_time = test_clock::now ();
523-
524- double average = 0 ;
525- for ( uint64_t i = 0 ; i < n; i++ )
526- average += 0.1 + rng.exgauss ( 0.3 , 0.06 , 0.25 );
527- average /= n;
528-
529- fmt::print ( " {} calls to 0.1 + rng.exgauss( 0.3,0.06,0.25 );: average = {:.8f}, time = {}s\n\n " ,
530- n, average, chrono::elapsed_fp_seconds (start_time) );
531- }
532-
533- // exponential
534- {
535- auto start_time = test_clock::now ();
536-
537- double nu = 0.25 ;
538- double average = 0 ;
539- for ( uint64_t i = 0 ; i < n; i++ )
540- {
541- double result = rng.exponential ( nu );
542- average += result;
543- }
544- average /= n;
545-
546- fmt::print ( " {} calls exp nu=0.25: average = {}, time = {}s\n\n " ,
547- n, average, chrono::elapsed_fp_seconds (start_time) );
548- }
549-
550- fmt::print ( " \n real:\n " );
551- for ( unsigned i = 1 ; i <= 100 ; i++ )
552- fmt::print ( " {:>11.8f}{}" , rng.real (), i % 10 == 0 ? " \n " : " " );
553-
554- fmt::print ( " \n gauss mean=0, std_dev=1.0:\n " );
555- for ( unsigned i = 1 ; i <= 100 ; i++ )
556- fmt::print ( " {:>11.8f}{}" , rng.gauss ( 0.0 , 1.0 ), i % 10 == 0 ? " \n " : " " );
557-
558- fmt::print ( " \n roll 30%:\n " );
559- for ( unsigned i = 1 ; i <= 100 ; i++ )
560- fmt::print ( " {:d}{}" , rng.roll ( 0.30 ), i % 10 == 0 ? " \n " : " " );
561-
562- fmt::print ( " \n " );
563- fmt::print ( " m_pi={}\n " , m_pi );
564- fmt::print ( " calls to rng::stdnormal_inv( double x )\n " );
565- fmt::print ( " x=0.975: {:.7} should be equal to 1.959964\n " , rng::stdnormal_inv ( 0.975 ) );
566- fmt::print ( " x=0.995: {:.8} should be equal to 2.5758293\n " , rng::stdnormal_inv ( 0.995 ) );
567- }
568-
569- #endif // UNIT_TEST
0 commit comments