#!/usr/bin/perl
# -*- mode: cperl; coding: utf-8; -*-

use strict;
use warnings;
use utf8;

use lib "/h/hamren/src/post/lib", "."; my $rval = do "common.pm" || die "$0: common.pm failed ($!) [$@]";  #--- Single-line common initializer

#--- End of header

sub description {
    my $arg1 = shift;
    p({style=>"padding-left:40pt; text-indent:-30pt;"}, "<b>$arg1</b>", @_);
}

sub file {
    my $arg1 = shift;
    p({style=>"padding-left:90pt; text-indent:-30pt;"}, "<span style='text-indent:-30pt; display: inline-block; min-width:15em'>««$arg1»»&nbsp;&nbsp;</span>", @_);
}

my $sp="&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;";

post(
    header(),

    p("In a ",
      href("post:marsaglia-random-numbers", "recent post"),
      " I wrote about five of George Marsaglia's random-number generators:"),

    table({class => "components-table"},
          trow(td("«shr3»"    ),  td("Three shifts, with a period of 2^32-1.")),
          trow(td("«mwc»"     ),  td("Multiply With Carry, with a period of about 2^60.")),
          trow(td("«kiss»"    ),  td("Keep It Simple Stupid, with a period of about 2^123")),
          trow(td("«lfib4»$sp"),  td("Lagged Fibonacci, with a period of about 2^287")),
          trow(td("«swb»"     ),  td("Subtract with borrow. Has a very long period: about 2^7578."))),

    p(" But how do they compare, speed-wise, to generators found in C, C++ and glibc?",
  " Here are five additional generators:"),

    table({class => "components-table"},
          trow(td("«rand»"     ),  td("The standard C routine. Returns  32-bit signed, but always non-negative, random number.",
                                      " Period is undocumented; assume that period is at most 2^31.",
                                      " It is not guaranteed to be thread-safe, but the GNU libc version is. ",
                                      " It uses a futex lock to protected its internal state.")),
          trow(td("«rand_r»"   ),  td("The reentrant version of «rand()». The state is a 32-bit unsigned integer, so the period can be at most 2^32.")),
          trow(td("«lrand48»"  ),  td("Defined by XOpen, now Open Group. Uses 48 bits of state. The return value is a 32-bit signed, but always non-negative, random number, just like for «rand()».",
                                      " The period is not documented, but the state is 48 bits.")),
          trow(td("«lrand48_r»"),  td("The reentrant version of «lrand48».")),
          trow(td("«mt19937»"  ),  td("A Mersenne Twister generator, from the C++ standard library «<random>».",
                                      " Its period is a whopping 2^19937."))
      ),

    p("For the purpose of the plot below, I have assumed that «rand()» and «rand_r()» have periods of 2^31,",
      " and that «lrand48()» and «lrand48_r()» have a periods of 2^48."),

    svg("plot-rng-times-1,000,000,000.svg"),

    p("The Marsaglia posting compared different ways of instantiating the generator.",
      " As it turned out, there was not significant difference between them.",
      " The times used here are those measured using a ⊂file-static object⊃."),

    p("It is not surprising that «lrand48()», using 48-bit arithmetic, is slower than the traditional «rand()».",
      " But why are the reentrant versions faster? This is strange, because of the overhead of passing state.",
      " The answer for ««rand()»» and ««rand_r()»» depends on the specifics of the glibc implementation. ",
      " While the standard does not require ««rand()»» to be reentrant, in glibc its state is protected by a futex lock, adding time.",
      " This does not apply to «lrand48()» and «lrand48_r()» .",
      " They both call the internal routine «__nrand48_r()», and nothing else, so exactly why «lrand48_r()» is consistently slightly faster than «lrand48()» is not obvious."),

    p("The generator mt19937, called a Mersenne Twister because 2^19927-1 is a Mersenne prime, is a lot more complicated, and probably much better, than any of the other.",
      " A good choice those who need a really godd RNG, and for everyone else who like the feeling of having a bit of a margin.",
      " It keeps a lot of state; an instance is exactly 5000 bytes is size with my setup."),

    h2("Measuring execution time"),
    p("It is strangely difficult to get consistent results when running the same test several times.",
      "On an multi-core CPU with no other load, the execution time should be virtually identical from one run to another",
      " at least when running for more than a few seconds.",
      " But that is not the case. ",
      " The variation probably has to do with the fact that Linux, on which I run these test, ",
      " places the stack and the heap at randomly addresses, that vary from run to run.",
      " This is a security feature, as it makes it more difficult to write malicious software."),
    p("To exemplify, here are all ten runs that the second plot above is based on."),
    svg("plot-rng-times-1,000,000,000-all-10.svg"),

    p("Each group of ten bars represent ten timing runs for one generator. ",
      " Within a group, the bars are ordered by decreasing execution time."),

    p("It is interesting to see how there are sometimes two levels, possibly depending on whether memory areas are ⊂close⊃ or ⊂distant⊃."),

    p("In the averaged plots, the two shortest and the two longest have been dropped, and the remaining six values have been averaged.",
      " While the shortest values show how fast the generator actually is, or at least sometimes can be, the ⊂middle average⊃ shows what can be expected."),

    h2("Download"),
    p("There is no separate download for this post.",
      " Use the one from the ",
      href("post:marsaglia-random-numbers", "previous post"),
      "."),

    h2("Final note"),

    p("I am not an expert in statistics or random numbers, and can neither judge nor make recommendations when it comes to random-number generators. ",
      " I have faith in George Marsaglia, the creator of the original diehard tests for RNGs.",
      " I also have faith in the the RNGs from GNU libc, because they get a lot of scrutiny. ",
      " Finally I have faith in the C++ «random» library, and the well-researched «mt19938» RNG."),
    footer()
);
__END__