Showing pcg32.pl (raw)


  1. #!/usr/bin/env perl
  2.  
  3. ###############################################################################
  4. # Two implementations of PCG32. One in native Perl with no dependencies, and
  5. # one that uses Math::Int64. Surprisingly the native version is significantly
  6. # faster.
  7. #
  8. # A lot of work was done here to mimic how C handles overflow multiplication
  9. # on large uint64_t numbers. Perl converts scalars that are larger than 2^64-1
  10. # to floating point on the backend. We do *NOT* want that for PCG, because
  11. # PCG (and more PRNGs) rely on overflow math to do their magic. We utilize
  12. # 'use integer' to force Perl to do all math with regular 64bit values. When
  13. # overflow occurs Perl likes to convert those values to negative numbers. In
  14. # the original C all math is done with uint64_t, so we have to convert the
  15. # IV/negative numbers back into UV/unsigned (positive) values. PCG also uses
  16. # some uint32_t variables internally, so we mimic that by doing the math in
  17. # 64bit and then masking down to only the 32bit number.
  18. #
  19. ###############################################################################
  20. #
  21. # Original C code from: https://www.pcg-random.org/download.html
  22. #
  23. # typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;
  24. #
  25. # uint32_t pcg32_random_r(pcg32_random_t* rng) {
  26. #     uint64_t oldstate = rng->state;
  27. #     // Advance internal state
  28. #     rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
  29. #     // Calculate output function (XSH RR), uses old state for max ILP
  30. #     uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
  31. #     uint32_t rot = oldstate >> 59u;
  32. #     return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
  33. # }
  34. #
  35. ###############################################################################
  36.  
  37. use strict;
  38. use warnings;
  39. use v5.16;
  40. use Math::Int64 qw(uint64 uint64_to_number);
  41. use Getopt::Long;
  42. use Test::More;
  43.  
  44. ###############################################################################
  45. ###############################################################################
  46.  
  47. my $debug = 0;
  48. my $s1    = 15939250660798104135; # Default 64bit seed1
  49. my $s2    = 3988331200502121509;  # Default 64bit seed2
  50. my $seeds = [];
  51.  
  52. GetOptions(
  53.     'debug'      => \$debug,
  54.     'seed1=i'    => \$s1,
  55.     'seed2=i'    => \$s2,
  56.     'random'     => \&randomize_seeds,
  57.     'unit-tests' => \&run_unit_tests,
  58. );
  59.  
  60. my $num = $ARGV[0] || 8;
  61. my ($seed1, $seed2);
  62.  
  63. print color(83, "Seeding PRNG with: $s1 / $s2\n");
  64.  
  65. $seeds = [uint64($s1), uint64($s2)];
  66. my @x = ();
  67. for (my $i = 0; $i < $num; $i++) {
  68.     my $num = pcg32_math64($seeds);
  69.     push(@x, $num);
  70. }
  71. print color('yellow', "Math::Int64: ") . join(", ", @x);
  72. print "\n";
  73.  
  74. ##################################
  75.  
  76. $seeds = [$s1, $s2];
  77. my @y  = ();
  78. for (my $i = 0; $i < $num; $i++) {
  79.     my $num = pcg32_native($seeds);
  80.     push(@y, $num);
  81. }
  82. print color('yellow', "Native Perl: ") . join(", ", @y);
  83. print "\n";
  84.  
  85. ##################################
  86.  
  87. print "\n";
  88. $seeds = [$s1, $s2];
  89. my @z  = ();
  90. for (my $i = 0; $i < $num; $i++) {
  91.     my $num = pcg64_native($seeds);
  92.     push(@z, $num);
  93. }
  94. print color('yellow', "Native Perl 64bit: ") . join(", ", @z);
  95. print "\n";
  96.  
  97. ################################################################################
  98. ################################################################################
  99. ################################################################################
  100.  
  101. #my $seeds = [uint64(12), uint64(34)];
  102. #my $rand  = pcg32_math64($seeds);
  103. sub pcg32_math64 {
  104.     # state/inc are passed in by reference
  105.     my ($s) = @_;
  106.  
  107.     my $oldstate = $s->[0];
  108.     $s->[0]      = $oldstate * 6364136223846793005 + ($s->[1] | 1);
  109.  
  110.     my $xorshifted = (($oldstate >> 18) ^ $oldstate) >> 27;
  111.     $xorshifted    = $xorshifted & 0xFFFFFFFF; # Convert to uint32_t
  112.  
  113.     my $rot    = $oldstate >> 59;
  114.     my $invrot = 4294967296 - $rot;
  115.  
  116.     my $ret = ($xorshifted >> $rot) | ($xorshifted << ($invrot & 31));
  117.     $ret    = $ret & 0xFFFFFFFF; # Convert to uint32_t
  118.  
  119.     $ret = uint64_to_number($ret);
  120.  
  121.     if ($debug) {
  122.         # $oldstate is the state at the start of the function and $inc
  123.         # doesn't change so we can print out the initial values here
  124.         print color('orange', "State : $oldstate/$s->[1]\n");
  125.         print color('orange', "State2: $s->[0]\n");
  126.         print color('orange', "Xor   : $xorshifted\n");
  127.         print color('orange', "Rot   : $rot\n");
  128.     }
  129.  
  130.     return $ret;
  131. }
  132.  
  133. #my $seeds = [12, 34];
  134. #my $rand  = pcg32_native($seeds);
  135. sub pcg32_native {
  136.     # state/inc are passed in by reference
  137.     my ($s) = @_;
  138.  
  139.     my $oldstate = $s->[0]; # Save original state
  140.  
  141.     # We use interger math because Perl converts to floats any scalar
  142.     # larger than 2^64. PCG *requires* 64bit uint64_t math, with overflow,
  143.     # to calculate correctly. We have to unconvert the overflowed number
  144.     # from an IV to UV after the big math
  145.     use integer;
  146.     $s->[0] = $oldstate * 6364136223846793005 + ($s->[1] | 1);
  147.     $s->[0] = iv_2_uv($s->[0]);
  148.     no integer;
  149.  
  150.     my $xorshifted = (($oldstate >> 18) ^ $oldstate) >> 27;
  151.     $xorshifted    = $xorshifted & 0xFFFFFFFF; # Convert to uint32_t
  152.  
  153.     my $rot = ($oldstate >> 59);
  154.  
  155.     # -$rot on a uint32_t is the same as (2^32 - $rot)
  156.     my $invrot = 4294967296 - $rot;
  157.     my $ret    = ($xorshifted >> $rot) | ($xorshifted << ($invrot & 31));
  158.  
  159.     # Convert to uint32_t
  160.     $ret = $ret & 0xFFFFFFFF;
  161.  
  162.     if ($debug) {
  163.         # $oldstate is the state at the start of the function and $inc
  164.         # doesn't change so we can print out the initial values here
  165.         print color('orange', "State : $oldstate/$s->[1]\n");
  166.         print color('orange', "State2: $s->[0]\n");
  167.         print color('orange', "Xor   : $xorshifted\n");
  168.         print color('orange', "Rot   : $rot\n");
  169.     }
  170.  
  171.     return $ret;
  172. }
  173.  
  174. # During large integer math when a UV overflows and wraps back around
  175. # Perl casts it as a IV value. For the purposes of PCG we need that
  176. # wraparound math to stay in place. We need uint64_t all the time.
  177. sub iv_2_uv {
  178.     my $x = $_[0];
  179.  
  180.     # Flip it from a IV (signed) to a UV (unsigned)
  181.     # use Devel::Peek; Dump($var) # See the internal Perl type
  182.     if ($x < 0) {
  183.         no integer;
  184.         $x += 18446744073709551615;
  185.         $x += 1;
  186.     }
  187.  
  188.     return $x;
  189. }
  190.  
  191. # To get a 64bit number from PCG32 you create two different generators
  192. # and combine the results into a single 64bit value. All the examples
  193. # online show 1 for the inc/seed2 value. I'm not sure why that is, but
  194. # I copied it for my implementation.
  195. #
  196. #my $seeds = [12, 34];
  197. #my $rand  = pcg64_native($seeds);
  198. sub pcg64_native {
  199.     my ($s) = @_;
  200.  
  201.     # Build a new object to send to each pcg32 instance
  202.     my $inc = 1; # Can be any 64bit value
  203.     my $one = [$s->[0], $inc];
  204.     my $two = [$s->[1], $inc];
  205.  
  206.     # Get two 32bit ints
  207.     my $high = pcg32_native($one);
  208.     my $low  = pcg32_native($two);
  209.  
  210.     # We copy the data back into the original object
  211.     $s->[0] = $one->[0];
  212.     $s->[1] = $two->[0];
  213.  
  214.     # Combine the two 32bits into one 64bit int
  215.     my $ret = ($high << 32) | $low;
  216.  
  217.     return $ret;
  218. }
  219.  
  220. ###############################################################################
  221. ###############################################################################
  222.  
  223. # String format: '115', '165_bold', '10_on_140', 'reset', 'on_173', 'red', 'white_on_blue'
  224. sub color {
  225.     my ($str, $txt) = @_;
  226.  
  227.     # If we're NOT connected to a an interactive terminal don't do color
  228.     if (-t STDOUT == 0) { return $txt || ""; }
  229.  
  230.     # No string sent in, so we just reset
  231.     if (!length($str) || $str eq 'reset') { return "\e[0m"; }
  232.  
  233.     # Some predefined colors
  234.     my %color_map = qw(red 160 blue 27 green 34 yellow 226 orange 214 purple 93 white 15 black 0);
  235.     $str =~ s|([A-Za-z]+)|$color_map{$1} // $1|eg;
  236.  
  237.     # Get foreground/background and any commands
  238.     my ($fc,$cmd) = $str =~ /^(\d{1,3})?_?(\w+)?$/g;
  239.     my ($bc)      = $str =~ /on_(\d{1,3})$/g;
  240.  
  241.     if (defined($fc) && int($fc) > 255) { $fc = undef; } # above 255 is invalid
  242.  
  243.     # Some predefined commands
  244.     my %cmd_map = qw(bold 1 italic 3 underline 4 blink 5 inverse 7);
  245.     my $cmd_num = $cmd_map{$cmd // 0};
  246.  
  247.     my $ret = '';
  248.     if ($cmd_num)      { $ret .= "\e[${cmd_num}m"; }
  249.     if (defined($fc))  { $ret .= "\e[38;5;${fc}m"; }
  250.     if (defined($bc))  { $ret .= "\e[48;5;${bc}m"; }
  251.     if (defined($txt)) { $ret .= $txt . "\e[0m";   }
  252.  
  253.     return $ret;
  254. }
  255.  
  256. sub randomize_seeds {
  257.     print color(51, "Using random seeds\n\n");
  258.  
  259.     $s1 = int(rand() * (2**64 - 1));
  260.     $s2 = int(rand() * (2**64 - 1));
  261. }
  262.  
  263. # Creates methods k() and kd() to print, and print & die respectively
  264. BEGIN {
  265.     if (eval { require Data::Dump::Color }) {
  266.         *k = sub { Data::Dump::Color::dd(@_) };
  267.     } else {
  268.         require Data::Dumper;
  269.         *k = sub { print Data::Dumper::Dumper(\@_) };
  270.     }
  271.  
  272.     sub kd {
  273.         k(@_);
  274.  
  275.         printf("Died at %2\$s line #%3\$s\n",caller());
  276.         exit(15);
  277.     }
  278. }
  279.  
  280. # Run a test with a given seed and return a string of the results
  281. sub quick_test32 {
  282.     my $seed = $_[0];
  283.  
  284.     my @data = ();
  285.     for (my $i = 0; $i < 4; $i++) {
  286.         my $num = pcg32_native($seed);
  287.         push(@data, $num);
  288.     }
  289.  
  290.     my $ret = join(", ", @data);
  291.     return $ret;
  292. }
  293.  
  294. sub quick_test64 {
  295.     my $seed = $_[0];
  296.  
  297.     my @data = ();
  298.     for (my $i = 0; $i < 4; $i++) {
  299.         my $num = pcg64_native($seed);
  300.         push(@data, $num);
  301.     }
  302.  
  303.     my $ret = join(", ", @data);
  304.     return $ret;
  305. }
  306.  
  307. sub run_unit_tests {
  308.     # Seeds < 2**32
  309.     cmp_ok(quick_test32([11, 22])      , 'eq', '0, 1425092920, 3656087653, 1104107026');
  310.     cmp_ok(quick_test32([33, 44])      , 'eq', '0, 3850707138, 2930351490, 1110209703');
  311.     cmp_ok(quick_test32([55, 66])      , 'eq', '0, 1725101930, 224698313, 2870828486');
  312.     cmp_ok(quick_test32([12345, 67890]), 'eq', '0, 8251198, 44679150, 3046830521');
  313.     cmp_ok(quick_test32([9999, 9999])  , 'eq', '0, 521292032, 3698775557, 199399470');
  314.  
  315.     cmp_ok(quick_test64([11, 22])      , 'eq', '0, 6120727489207695446, 7904312005358798897, 14733674221366828425');
  316.     cmp_ok(quick_test64([33, 44])      , 'eq', '0, 16538661225628040268, 5269891931295187491, 5495286771333204711');
  317.     cmp_ok(quick_test64([55, 66])      , 'eq', '0, 7409256372025208996, 8212781881022671801, 8831782971077082788');
  318.     cmp_ok(quick_test64([12345, 67890]), 'eq', '0, 35438628484449140, 42862460907032573, 519456495312580246');
  319.     cmp_ok(quick_test64([9999, 9999])  , 'eq', '0, 2238932229626677504, 14236525402126437484, 10387246122801752400');
  320.  
  321.     # Seeds > 2**32
  322.     cmp_ok(quick_test32([42862460907032573, 519456495312580246])    , 'eq', '319349001, 562730850, 2229409754, 561058538');
  323.     cmp_ok(quick_test32([6120727489207695446, 7904312005358798897]) , 'eq', '635930912, 2099303707, 1638577555, 1426136496');
  324.     cmp_ok(quick_test32([4841811808465514507, 7141191103728083377]) , 'eq', '1986408540, 4264878569, 3066617590, 731859269');
  325.  
  326.     cmp_ok(quick_test64([42862460907032573, 519456495312580246])    , 'eq', '1371593519175525487, 17623029558467823369, 17850014000156247978, 768534907509427587');
  327.     cmp_ok(quick_test64([6120727489207695446, 7904312005358798897]) , 'eq', '2731302471965979098, 3465889473135782122, 4841811808465514507, 7141191103728083377');
  328.     cmp_ok(quick_test64([4841811808465514507, 7141191103728083377]) , 'eq', '8531559717926221063, 6031125200978744796, 3704366926003160989, 5594521440717127703');
  329.  
  330.     done_testing();
  331.     exit(0);
  332. }
  333.  
  334. # vim: tabstop=4 shiftwidth=4 noexpandtab autoindent softtabstop=4