in reply to PRNG/TRNG Cesaro's theorem

Of course you will want to do it differently, but here is something I wrote to play with this. I used the ntheory module because it has different RNGs as well as gcd and Pi.

#!/usr/bin/env perl
use warnings;
use strict;
use ntheory ":all";

cesaro("drand48", sub { int(rand(1<<32)) } );
cesaro("ChaCha20", sub { irand64 } );
cesaro("/dev/urandom", sub { unpack("Q",random_bytes(8)) } );

sub cesaro {
  my($name, $rng) = @_;

  print "Using $name:\n";
  for my $e (1..7) {
    my $n = 10**$e;
    my $t = 0;
    for (1..$n) { $t++ if gcd( $rng->(), $rng->() ) == 1; }
    printf "%8d  %10.8f\n", $n, sqrt(6*$n/$t);
  }
  print "       Pi ", Pi(), "\n\n";
}

[download]

I saved cut-n-paste by passing in the RNG code. One for Perl's default rand (which is drand48 on modern Perl), one for ntheory's CSPRNG, and one that gets data from /dev/urandom. You could replace that with something that got data from random.org. There are even a couple modules that do it automatically (Math::RandomOrg and Data::Entropy). Since we expect t/n ~ 6/Pi^2, a little algebra gives us Pi ~ sqrt(6n/t). I don't believe we can read anything into what the results might mean for the different RNG methods. The results will be different every time unless we seed (and the last can't be seeded since it is O/S collected entropy (long technical discussion of /dev/random vs. /dev/urandom vs. hardware could be had here)).

Using drand48:
      10  2.73861279
     100  3.06186218
    1000  3.15964572
   10000  3.14632429
  100000  3.13988122
 1000000  3.14035598
10000000  3.14143144
       Pi 3.14159265358979

Using ChaCha20:
      10  3.16227766
     100  3.18896402
    1000  3.17287158
   10000  3.15675816
  100000  3.14329189
 1000000  3.13978836
10000000  3.14138726
       Pi 3.14159265358979

Using /dev/urandom:
      10  3.46410162
     100  3.03821810
    1000  3.11588476
   10000  3.14217962
  100000  3.13643021
 1000000  3.14020372
10000000  3.14130744
       Pi 3.14159265358979

[download]

Replies are listed 'Best First'.
Re^2: PRNG/TRNG Cesaro's theorem
by marioroy (Prior) on Oct 09, 2017 at 01:30 UTC

    Hello brothers and sisters of the Monastery,

    Seeing danaj's post made me think of passing arguments to workers for a parallel demonstration. But first, I need to check if random numbers are unique between workers. They are not for non-threaded workers, irand64 and random_bytes.

    Here is the test code. Calling MCE::relay is a way to have workers display output orderly, starting with worker 1, 2, ..., 8. The init_relay value isn't used here, but the option tells MCE to load MCE::Relay and enable relay capability. Workers persist between each run.

    use strict;
use warnings;

use ntheory ":all";
use MCE::Flow;

my ( $name, $rng );

my %rand = (
  "drand48" => sub { int(rand(1 << 32)) },
  "ChaCha20" => sub { irand64() },
  "/dev/urandom" => sub { unpack("Q", random_bytes(8)) }
);

MCE::Flow::init(
  max_workers => 8,
  init_relay  => 0,
  user_begin  => sub {
    $name = MCE->user_args()->[0];
    $rng  = $rand{ $name };
  }
);

sub func {
  MCE::relay {
    print MCE->wid(), ": ", $rng->(), "\n";
  };
}

for my $name ( "drand48", "ChaCha20", "/dev/urandom" ) {
  print "Usage $name:\n";
  mce_flow { user_args => [$name] }, \&func;
  print "\n";
}

    [download]

    Output.

    Usage drand48:
1: 3498494761
2: 2506930441
3: 1515366121
4: 523801801
5: 3827204777
6: 2835640457
7: 1844076137
8: 852511817

Usage ChaCha20:
1: 4471005142860083063
2: 4471005142860083063
3: 4471005142860083063
4: 4471005142860083063
5: 4471005142860083063
6: 4471005142860083063
7: 4471005142860083063
8: 4471005142860083063

Usage /dev/urandom:
1: 15746895497052787399
2: 15746895497052787399
3: 15746895497052787399
4: 15746895497052787399
5: 15746895497052787399
6: 15746895497052787399
7: 15746895497052787399
8: 15746895497052787399

    [download]

    Later today will release MCE 1.831 and MCE::Shared 1.832 containing the fix.

    Usage drand48:
1: 600074529
2: 3903477505
3: 2911913185
4: 1920348865
5: 928784545
6: 4232187521
7: 3240623201
8: 2249058881

Usage ChaCha20:
1: 8740887910466299010
2: 12948789762855324085
3: 7574729187958724006
4: 14608687740989597345
5: 10145950054018120246
6: 11767641523694169551
7: 5811941457879652367
8: 2397613489984096139

Usage /dev/urandom:
1: 14391656456731294109
2: 2750708286643159769
3: 10844675827853246458
4: 7920672879166021322
5: 16939013845838223421
6: 9482848646152826462
7: 11535629003375292447
8: 6903845178044907896

    [download]

    Once the release hits CPAN, I'd come back and post a parallel demonstration.

    Regards, Mario

[reply]
[d/l]
[select]

      MCE 1.831 and MCE::Shared 1.832 have been released containing the fix. What follows is the parallel demonstration for danaj's example. For sequence of numbers, the MCE bounds_only option is handy. Workers receive the begin and end values only.

      use strict;
use warnings;

use ntheory 0.67 ":all";
use MCE::Flow 1.831;

my ( $name, $rng );

my %rand = (
  "drand48" => sub { int(rand(1 << 32)) },
  "ChaCha20" => sub { irand64() },
  "/dev/urandom" => sub { unpack("Q", random_bytes(8)) }
);

# Workers receive [ begin, end ] values.

MCE::Flow::init(
  max_workers => MCE::Util::get_ncpu(),
  chunk_size  => 10000,
  bounds_only => 1,
  user_begin  => sub {
    $name = MCE->user_args()->[0];
    $rng  = $rand{ $name };
  }
);

sub func {
  my ( $beg_seq, $end_seq ) = @{ $_ };
  my ( $t ) = ( 0 );

  for ( $beg_seq .. $end_seq ) {
    $t++ if gcd( $rng->(), $rng->() ) == 1;
  }

  MCE->gather($t);
}

# The user_args option is how to pass arguments.
# Workers persist between each run.

sub cesaro {
  my ( $name ) = @_;
  print "Usage $name:\n";

  for my $e ( 1..7 ) {
    my $n   = 10 ** $e;
    my @ret = mce_flow_s { user_args => [$name] }, \&func, 0, $n - 1;
    my $t   = 0;   $t += $_ for @ret;

    printf "%8d  %0.8f\n", $n, sqrt(6 * $n / $t);
  }

  printf "%9s %s\n\n", "Pi", Pi();
}

for ( "drand48", "ChaCha20", "/dev/urandom" ) {
  cesaro($_);
}

      [download]

      Output.

      Usage drand48:
      10  3.46410162
     100  3.11085508
    1000  3.12347524
   10000  3.15335577
  100000  3.14122349
 1000000  3.14092649
10000000  3.14209456
       Pi 3.14159265358979

Usage ChaCha20:
      10  3.46410162
     100  3.08606700
    1000  3.11588476
   10000  3.14606477
  100000  3.14461263
 1000000  3.14453748
10000000  3.14269499
       Pi 3.14159265358979

Usage /dev/urandom:
      10  3.16227766
     100  3.13625024
    1000  3.24727816
   10000  3.13959750
  100000  3.14238646
 1000000  3.14247180
10000000  3.14023908
       Pi 3.14159265358979

      [download]

      The parallel code scales linearly. It runs about 4x faster on a machine with 4 "real" cores. A little faster with extra hyper-threads.

      Regards, Mario

[reply]
[d/l]
[select]

        The OP mentioned (here) on marking down the number pairs whenever gcd(x, y) is equal to 1.

        Here, workers write the number pairs to a shared STDERR file handle. Locking is handled automatically by the shared-manager. For best performance, do not have workers write each pair individually to the shared file handle. Instead, have workers append to a local variable. Writing once per chunk/segment relieves pressure on the shared-manager process, which must keep up with many workers.

        To run, redirect STDERR to a file: perl demo.pl 2> pairs.txt

        use strict;
use warnings;

use ntheory 0.67 ":all";
use MCE::Flow 1.831;
use MCE::Shared;

mce_open my $err_fh, '>>', \*STDERR;

my ( $name, $rng );

my %rand = (
  "drand48" => sub { int(rand(1 << 32)) },
  "ChaCha20" => sub { irand64() },
  "/dev/urandom" => sub { unpack("Q", random_bytes(8)) }
);

# Workers receive [ begin, end ] values.

MCE::Flow::init(
  max_workers => MCE::Util::get_ncpu(),
  chunk_size  => 10000,
  bounds_only => 1,
  user_begin  => sub {
    $name = MCE->user_args()->[0];
    $rng  = $rand{ $name };
  }
);

sub func {
  my ( $beg_seq, $end_seq ) = @{ $_ };
  my ( $t, $pairs, $r1, $r2 ) = ( 0, '' );

  for ( $beg_seq .. $end_seq ) {
    ( $r1, $r2 ) = ( $rng->(), $rng->() );
    if ( gcd($r1, $r2) == 1 ) {
      $pairs .= sprintf("%20s %20s\n", $r1, $r2);
      $t++;
    }
  }

  print $err_fh $pairs;

  MCE->gather($t);
}

# The user_args option is how to pass arguments.

sub cesaro {
  my ( $name ) = @_;
  print "Usage $name:\n";

  my $n   = 10 ** 5;
  my @ret = mce_flow_s { user_args => [$name] }, \&func, 0, $n - 1;
  my $t   = 0;   $t += $_ for @ret;

  printf "%8d  %0.8f\n", $n, sqrt(6 * $n / $t);
  printf "%9s %s\n\n", "Pi", Pi();
}

cesaro("/dev/urandom");

        [download]

        stdout

        Usage /dev/urandom:
  100000  3.14168852
       Pi 3.14159265358979

        [download]

        stderr

         3287478503828099458 17020534165422626509
12874128583739175743 15858105624010204897
12114767169191945777  1381808520787800067
 6217235583185614040 13880267056100467759
 6095107227201871385 14484880859615821612
 6023474899938562107  6188664627301334099
 2843055836738994412  3086747352859159329
10461621510749606095   302140866635797242
11160050491040241465  8409335412227119072
 3718439936437935345  5121716120103428116

...

        [download]

        Regards, Mario

[reply]
[d/l]
[select]
Re^2: PRNG/TRNG Cesaro's theorem
by Anonymous Monk on Oct 08, 2017 at 13:45 UTC
    Notice that the number of correct digits of π is roughly half of the number of digits in $n, because this is basically a random walk process.
[reply]

In Section Seekers of Perl Wisdom