www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.bugs - [Issue 3738] New: MinstdRand0 and MinstdRand very poor performance

reply d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738

           Summary: MinstdRand0 and MinstdRand very poor performance
           Product: D
           Version: 2.039
          Platform: x86
        OS/Version: Linux
            Status: NEW
          Severity: enhancement
          Priority: P2
         Component: Phobos
        AssignedTo: nobody puremagic.com
        ReportedBy: baryluk smp.if.uj.edu.pl



17:52:24 PST ---
MinstdRand0 is instantiation of CongruentialGenerator with constants
  a = 16807
  c = 0
  m = 2^^31 - 1

most important method of CongruetialGenerator popFront, doesn't take into the
account special structure of c,
which makes this generator quite slow especially compared to Mersen Twister
Mt19937.

This is important especially on the 32bit machines, because modulo operation
(%) is done using external (and notinlined) method __ULDIV__ which computes
unsigned long division and reminder. But it is also a problem on 64bit
machines, becuase even in hardware modulo operation is slow.

Currently compiler uses known trick for x % c, when c=2^^k is constant.
It changes it to binary and: x & (c - 1)

This is done to uint directly. For signed x it compiler emits 2 additional
'sub's to compensate for sign. For smaller than int it also works. If c=2^^k is
bigger than x.max it should be noop (and maybe small warning), but it isn't
really important. For ulong on 32bit it even works smarter by just anding with
only upper register (or lower and zeroing upper). This is all done even without
-O switch.

Not many programer know that there is also a trick for constant c = 2^^k - 1.

You can find some information about it (special case c= 2^^32 - 1) in this
article on page 6.

G. Marsaglia, 'On the randomness of Pi and other decimal expansions',

http://interstat.statjournals.net/YEAR/2005/articles/0510005.pdf


What is important for us is that currently MinstdRand{,0} is very slow. Much
slower than MersenTwistter,
which is more complex than MinstdRand, which is quite supprising.

I propose changing current popFront of LinearCongruentialEngine:

     void popFront()
     {
         static if (m)
-            _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
+            static if (is(UIntType == uint) && m == 4294967295uL) {
+                 const ulong x = (cast(ulong) a * _x + c);
+                 const ulong v = x >> 32;
+                 const ulong w = (x & 0xffff_ffffuL);
+                 const uint y = cast(uint)(v+w);
+                 _x = (y < v || y >= 4294967295u) ? (y+1) : y;
+            } else static if (is(UIntType == uint) && m == 2147483647u) {
+                 const ulong x = (cast(ulong) a * _x + c);
+                 const ulong v = x >> 31;
+                 const ulong w = (x & 0x7fff_ffffuL);
+                 const uint y = cast(uint)(v+w);
+                 _x = (y < v || y >= 2147483647u) ? ((y+1) & 0x7fff_ffffu) :
y;
+            } else {
+                _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
+            }
         else
             _x = a * _x + c;
     }

According to my test first 10_000_000_000 random numbers (few periods) are
exactly the same with and without this change. I provided also 2^^32-1 version,
I tested it in different place, but there are also some generator with such
modulus.

It can by also easly changed to 2^^63-1 and 2^^64-1, but it will need ucent
support,
which we lack currently. But it will allow implementing one very good
generators
(read bellow).

I hope it can be made general and puted directly into backends of DVD and LLVM,
GCC to make it simpler
for programer. Currently non of D or C compiler which I tested uses this trick.

Time mesurments for original code summarizes table.
All test performed under Debian Gnu/Linux, unstable, dmd 2.039, Linux
2.6.33-rc4-sredniczarny-00399-g24bc734-dirty.
Cpu:
model name    : Intel(R) Pentium(R) M processor 1.73GHz
stepping    : 8
cpu MHz        : 1733.000
cache size    : 2048 KB

user root, realtime RR scheduler, priority 50 (this gives deviations smaller
than 1%).

compiled in "-O -inline -release" mode, with all generators in the same file as
main function.

1_000_000_000 sample points.

MinstdRand0, realtime rr 50:
   standard: 44.62 s
   inline: 41.04
   standard with tricks: 17.03 s
   inline with tricks: 17.00 s

MinstdRand, realtime rr 50
  standard: 44.88 s
  inline: 41.06
  standard with tricks: 17.05 s
  inline with tricks: 17.02 s

Mt19937, realtime rr 50:
  standard: 33.93 s



"Inline" is manually inlined version (interesingly it isn't always faster),
probably some problems with registers, because dmd2 isn't inlining popFront and
front functions unfortunetly (even that they are very simple). Tricks is
version with patch above.

As you can see after changes this method is about 2.6 times faster. And it is
faster than Mersen Twister.

This is quite interesting, becuase patched code have 2 branches, 3 binary ands,
and 2 additions. But as maybe few know, DIV/MOD operation is done using quite
big external (notinlined) assembler routing __ULDIV__, which is just slow for
us. So it can be usefull to also create utility function for our special mod
operation.

Unittests still passes.






MinstdRand0 and MinstdRand are pretty good generators, for this speed.
But unfortunetly they fail lots of tests (see TestU01 documentation). Mosly
becuase of small period.

I know that they are in Phobos primarly becuase of simplicity (which maybe my
patch destroys, if so please incorporate it directly into compiler, so also
others can be beneficial of this).


There are only few which are better, so simple and very fast.

One is m=2^^61-1,a=2^^30-2^^19. But unfortunetly it is slow on 32bit machines
(128bit intermidiates). It is fast on 64bit machines (about 10% slower than
MinstdRand0), but have much better statistical properties (see Wu 1997).

One of very interesting and extra simple generator is Marsa-LFIB4 by Marsaglia.
It is just:
   x_i = (x_{i-55} + x_{i-119} + x_{i-179} + x_{i-256}) mod 2^^32
no multiplications and no complicated modulo operations. It is very fast
(especially on 64-bit machines), have period of 2^^287 and passes practically
all tests (actually all). 2 lines of code.

Other is KISS99 also by Marsaglia. it passes also all tests. Period 2^^123.
It is combination of 3 very simple generators. 4 Lines of code.

There also few LaggedFibbonaci generators which are ultra fast. But only on
64bit machines. They passes all tests. One example is LFib(2^^64, 1279, 861, *)
with period 2^^1340.

I also know one fast random generator derived from cryptographic primitieves,
notably from Rijandle/AES cipher. It is known as ISAAC by Jenkins and is fvery
ast, and passes all tests. Code have about 130 lines.

There is also Matlab randn function which is simple (xor combination of simple
LCG with m=2^^32 and simple shifting generator. both by Marsaglia) and have
period 2^^64, is very fast and passes most tests. It is used there for Normal
varietes, but as integer generator it is very good.


All above presented generators are faster and better (according to tests) than
Mt19937.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 23 2010
next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738




18:29:30 PST ---
Ok, sorry. Not all are better than Mt19937. But definitly better than
MinstdRand{,0}. And few are acutally better (empirically).

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 23 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738


Andrei Alexandrescu <andrei metalanguage.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |ASSIGNED
                 CC|                            |andrei metalanguage.com



23:09:46 PST ---
Thanks for the fix. I'll incorporate the patch into Phobos. Walter will decide
whether he wants to pursue a compiler change.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 23 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738


Witold Baryluk <baryluk smp.if.uj.edu.pl> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |NEW



08:35:08 PST ---
One can observe that ?: line (branch2)
        _x = (y < v || y >= 2147483647) ? ((y+1) & 0x7fff_ffffu) : y;
can be changed , first to not contain not needed y<v, to  (branch1)
        _x = (y >= 2147483647) ? ((y+1) & 0x7fff_ffffu) : y;
or even remove branch completly (branchlessA)
        _x = (y & 0x7fff_ffff) + (y >> 31);
and possibly rearange terms (branchlessB)
        _x = (y + (y >> 31)) & 0x7fff_ffff;


MinstdRnd0
  standard: 44.65 s
  inline: 41.11 s
  standard with tricks (2 branches): 17.08 s
  inline with tricks (2 branches): 17.04 s
  standard with tricks (1 branch): 15.24 s
  inline with tricks (1 branch): 13.59 s
  standard with tricks (branchlesA): 16.50 s
  inline with tricks (branchelesA): 15.82 s
  standard with tricks (branchlesB):  17.74 s
  inline with tricks (branchelesB): 15.30 s

MinstdRnd
  standard: 44.95 s
  inline: 41.17 s
  standard with tricks (2 branch): 17.04 s
  inline with tricks (2 branche): 17.03 s
  standard with tricks (1 branch): 15.29 s
  inline with tricks (1 branche): 13.54 s
  standard with tricks (branchlesA): 16.46 s
  inline with tricks (branchelessA): 15.88 s
  standard with tricks (branchlesB): 18.10 s
  inline with tricks (branchelessB): 15.31 s

So fastest is currently branch1 version. 

Updated patch:

     void popFront()
     {
         static if (m)
-            _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
+            static if (is(UIntType == uint) && m == 4294967295uL) {
+                 const ulong x = (cast(ulong) a * _x + c);
+                 const ulong v = x >> 32;
+                 const ulong w = (x & 0xffff_ffffuL);
+                 const uint y = cast(uint)(v+w);
+                 _x = (y < v || y >= 4294967295u) ? (y+1) : y;
+            } else static if (is(UIntType == uint) && m == 2147483647u) {
+                 const ulong x = (cast(ulong) a * _x + c);
+                 const ulong v = x >> 31;
+                 const ulong w = (x & 0x7fff_ffffuL);
+                 const uint y = cast(uint)(v+w);
+                 _x = (y >= 2147483647u) ? ((y+1) & 0x7fff_ffffu) : y;
+            } else {
+                _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
+            }
         else
             _x = a * _x + c;
     }

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 24 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738


Andrei Alexandrescu <andrei metalanguage.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |ASSIGNED
         AssignedTo|nobody puremagic.com        |andrei metalanguage.com



09:19:04 PST ---
Here's the modified version as I plan to commit:

    void popFront()
    {
        static if (m)
        {
            static if (is(UIntType == uint) && m == uint.max)
            {
                immutable ulong
                    x = (cast(ulong) a * _x + c),
                    v = x >> 32,
                    w = x & uint.max;
                immutable y = cast(uint)(v + w);
                _x = (y < v || y == uint.max) ? (y + 1) : y;
            }
            else static if (is(UIntType == uint) && m == int.max)
            {
                immutable ulong
                    x = (cast(ulong) a * _x + c),
                    v = x >> 31,
                    w = x & int.max;
                const uint y = cast(uint)(v + w);
                _x = (y >= int.max) ? ((y + 1) & int.max) : y;
            }
            else
            {
                _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
            }
        }
        else
        {
            _x = a * _x + c;
        }
    }

I changed formatting a bit and literals with symbolic constants. Unfortunately
unittests don't fail when I slightly change some constants :o|.

Looks good?

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 24 2010
prev sibling next sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738




12:03:20 PST ---
All my tests pases on your version.
10_000_000_000 (more than all possibilities) passes (for MinstdRand0 and
MinstdRand, it is possible that for other a it will fail, but i tested and done
some calculation.  and it should work for all a,c<m.
I also tested versions with slightly changed constants (like 30 instad of 31,
or int.max-1, int.max-2, int.max-3, int.max+1 in multiple comibinations an
allone), and all them fail.


Shouldn't const uint y by immucatable y, just like in first static if?


I think that our pretty optimised line, can be further optimalised: (1 branch
sub)

-             _x = (y >= int.max) ? ((y + 1) & int.max) : y;
+             _x = (y >= int.max) ? (y - int.max) : y;

:)

This gives:
MinstdRand0:
  standard with tricks (1 branch sub): 15.25
  inline with tricks (1 branch sub): 13.53 s
MinstdRand:
  standard with tricks (1 branch sub): 15.26
  inline with tricks (1 branch sub): 13.52 s

It doesn't give much difference (maybe 2%). But i thinkg it is better to have
one substraction, than one add and one binary and.



I found few interesting materials about this m=2^^31-1 operation.
One is http://www.firstpr.com.au/dsp/rand31/

I tested this Carta's versions: But they are only relevant for a=16807 (out
Park-Miller's StdminRand0 generator), so pretty limited:

Carta1 (original):
      uint lo = 16807*(_x & 0xFFFFu);
      immutable uint hi = 16807*(_x >> 16);
      lo += (hi & 0x7FFF) << 16;
      lo += (hi >> 15);
      _x = (lo > int.max) ? (lo - int.max) : lo;

Carta3:
      uint lo = 16807*(_x & 0xFFFFu);
      immutable uint hi = 16807*(_x >> 16);
      lo += ((hi & 0x7FFF) << 16) + (hi >> 15);
      _x = (lo > int.max) ? (lo - int.max) : lo;

Carta2: 
        immutable uint
           hi = 16807*(_x >> 16),
           lo 16807*(_x & 0xFFFFu) + ((hi & 0x7FFF) << 16) + (hi >> 15);
           _x = (lo > int.max) ? (lo - int.max) : lo;

Last lines can also be changed to branchles:
           _x = (lo > int.max) ? (lo - int.max) : lo;

Timings:
Carta code1: 22.30 s
Carta code1 branchless: 18.42
Carta code3: 21.84 s
Carta code3 branchless: 18.19 s
Carta code2: 23.84 s
Carta code2 branchless: 19.84 s

So it is slower in all cases, and it is only for speciall cases (c=0, a <
2^^15). Of course answers are still correct.


Performing uint*uint+uint -> ulong is just faster than their tricks, and more
general, and can be reused in many other places.

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
Jan 24 2010
prev sibling parent d-bugmail puremagic.com writes:
http://d.puremagic.com/issues/show_bug.cgi?id=3738


Andrei Alexandrescu <andrei metalanguage.com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED



15:32:46 PST ---

 All my tests pases on your version.
[snip] Thanks! Committed with revision 1406. Any chance you send us your unittests so I can integrate them within the module? -- Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email ------- You are receiving this mail because: -------
Jan 24 2010