Parallel Random Number Generation
=================================

There are five strategies implemented that can be used to produce
repeatable pseudo-random numbers across multiple processes (local
or distributed). The final method is suitable for parallel
applications where reproducibility is not a concern.


* :ref:`using-seed-sequence`
* :ref:`distinct-cryptographic-keys`
* :ref:`advancing`
* :ref:`jumping`
* :ref:`weyl-sequences`
* :ref:`non-reproducible-sequences`

.. _using-seed-sequence:

Using a SeedSequence
--------------------

Using a :class:`~numpy.random.SeedSequence` is a universal way to initialize many streams.
When using with a bit generator with a moderately large state (say 256-bits), the probability
of a collision is negligible. This said, some generators can have difficult to detect
correlation when producing many streams in corner cases.

.. code-block:: python

   import numpy as np
   from randomgen import SFC64

   NUM_STREAMS = 8196
   seed_seq = SeedSequence(603484028490308141)
   children = seed_seq.spawn(NUM_STREAMS)

   streams = [SFC64(child) for child in children]

.. _distinct-cryptographic-keys:

Distinct Keys in Cryptographic Generators
-----------------------------------------

The cryptographic pseudo-random number generators (PRNGs) support using distinct
keys to produce distinct sequence.
Generators that support this form of parallelization include :class:`~randomgen.aes.AESCounter`,
:class:`~randomgen.chacha.ChaCha`, :class:`~randomgen.hc128.HC128`,
:class:`~randomgen.threefry.ThreeFry`, :class:`~randomgen.philox.Philox`, and
:class:`~randomgen.speck128.SPECK128`.

This example shows how many streams can be created by passing in different index
values in the second input while using the same seed in the first.

.. code-block:: python

   import numpy as np
   from randomgen import AESCounter, random_entropy

   NUM_STREAMS = 8196
   keys = random_entropy(4 * NUM_STREAMS).view(np.uint64)
   keys = set([tuple(key.tolist()) for key in keys.reshape((-1, 2))])
   # Essentially 0 probability this is needed
   while len(keys) < NUM_STREAMS:
       new_keys = random_entropy(4 * (NUM_STREAMS - len(keys))).view(np.uint8)
       new_keys = set([tuple(key.tolist()) for key in new_keys.reshape((-1, 2))])
       keys.update(new_keys)

   # Distinct 128-bit numbers as a key encoded as 2 uint64 values
   streams = [AESCounter(key=np.array(key,dtype=np.uint64)) for key in keys]

.. _advancing:

Advancing the PRNG's state
--------------------------

Most of the cryptographic PRNGs are counter-based, and so support advancing
which increments the counter. Advancing a PRNG updates the underlying PRNG
state as if a given number of calls to the underlying PRNG have been made.
In general there is not a one-to-one relationship between the number output
random values from a particular distribution and the number of draws from
the core PRNG. This occurs for two reasons:

* The random values are simulated using a rejection-based method
  and so, on average, more than one value from the underlying
  PRNG is required to generate a single draw.
* The number of bits required to generate a simulated value
  differs from the number of bits generated by the underlying
  PRNG.  For example, two 16-bit integer values can be simulated
  from a single draw of a 32-bit PRNG.

Advancing the PRNG state resets any pre-computed random numbers. This is
required to ensure exact reproducibility.


.. code-block:: python

   import numpy as np
   from randomgen import SPECK128, SeedSequence

   PHI = (np.sqrt(5) - 1) / 2
   STEP = int(PHI * 2**96)
   NUM_STREAMS = 8196

   seed_seq = SeedSequence(603484028490308141)
   base = SPECK128(seed_seq)
   streams = [base]
   for i in range(1, NUM_STREAMS):
       next_gen = SPECK128(seed_seq)
       streams.append(next_gen.advance(i * STEP))

In addition to the cryptographic PRNGs,
the PCG-based generators also support ``advance``: :class:`~randomgen.pcg64.PCG64`,
:class:`~randomgen.pcg64.LCG128Mix`, and :class:`~randomgen.pcg32.PCG32`.
Note that :class:`~randomgen.hc128.HC128` is based on a stream cipher and so
does not support advancing a counter.

.. _jumping:

Jumping the PRNG state
----------------------

``jumped`` advances the state of the PRNG as if a large number of random
numbers have been drawn, and returns a new instance with this state.  Jumping
is more universal than ``advance`` since the multiplier needed to jump some PRNGs
is expensive to compute. However, this multiplier can be pre-computed for a fixed
step size when the PRNG does not support an arbitrary advance which enables the jump.
The specific number of as if draws varies by PRNG, and ranges from around :math:`2^{64}` to
:math:`2^{512}`.  Additionally, the as if draws also depend on the size of
the default random number produced by the specific PRNG.  The PRNGs that
support ``jumped``, along with the period of the PRNG, the size of the jump
and the bits in the default unsigned random are listed below.

+-----------------+-------------------------+-------------------------+-------------------------+
| PRNG            | Period                  |  Jump Size              | Bits                    |
+=================+=========================+=========================+=========================+
| AESCounter      | :math:`2^{128}`         | :math:`2^{64}`          | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| ChaCha          | :math:`2^{128}`         | :math:`2^{64}`          | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| LCG128Mix       | :math:`2^{128}`         | :math:`\phi`            | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| DSFMT           | :math:`2^{19937}`       | :math:`2^{128}`         | 53                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| LXM             | :math:`2^{256}`         | :math:`2^{128}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| MT19937         | :math:`2^{19937}`       | :math:`2^{128}`         | 32                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| PCG32           | :math:`2^{64}`          | :math:`\phi`            | 32                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| PCG64           | :math:`2^{128}`         | :math:`\phi`            | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| Philox          | :math:`2^{256}`         | :math:`2^{128}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| SFMT            | :math:`2^{19937}`       | :math:`2^{128}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| SPECK128        | :math:`2^{128}`         | :math:`2^{64}`          | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| ThreeFry        | :math:`2^{256}`         | :math:`2^{128}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| Xoroshiro128    | :math:`2^{128}`         | :math:`2^{64}`          | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| Xorshift1024    | :math:`2^{1024}`        | :math:`2^{512}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| Xoshiro256      | :math:`2^{256}`         | :math:`2^{128}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+
| Xoshiro512      | :math:`2^{512}`         | :math:`2^{256}`         | 64                      |
+-----------------+-------------------------+-------------------------+-------------------------+

``jumped`` can be used to produce long blocks that are long enough to not
overlap. A jump size of :math:`\phi` is the integer value of :math:`\sqrt{5}/2 - 1`, the
golden ratio, times the full period.

.. code-block:: python

   from randomgen.entropy import random_entropy
   from randomgen import Xoshiro512, SeedSequence

   NUM_STREAMS = 8196
   seed = SeedSequence()
   blocked_rng = []
   last_rng = rng = Xoshiro512(seed)
   for i in range(NUM_STREAMS):
       blocked_rng.append(last_rng)
       last_rng = last_rng.jumped()

.. _weyl-sequences:

Weyl Sequences
--------------

:class:`~randomgen.sfc.SFC64` uniquely supports generating streams using distinct
Weyl sequences. Distinct sequences are produced by setting unique values in ``k```
when initializing :class:`~randomgen.sfc.SFC64`. These increments must be odd.
There is some indication that sparse bit patterns in the Weyl sequence with 50%
or fewer non-zero bits perform best.

.. code-block:: python

   import numpy as np
   from randomgen.entropy import random_entropy
   from randomgen import SFC64, SeedSequence

   NUM_STREAMS = 8196


   seed_seq = SeedSequence()
   # Helper function to generate values with 32 non-zero bits
   # The number of non-zero bits is configurable
   weyl_inc = SFC64(seed_seq).weyl_increments(1000)
   streams = [SFC64(seed_seq, k=k) for k in weyl_inc]

.. _non-reproducible-sequences:

Non-reproducible Sequences
--------------------------

:class:`randomgen.rdrand.RDRAND` uses a hardware based random number generator to
produce non-reproducible sequences on random numbers. These can be used on any
system that supports the `RDRAND` instruction, subject to the key limitation that
results cannot be reproduced without storing the values generated.

.. code-block:: python

   from randomgen import RDRAND

   NUM_STREAMS = 8196
   streams = [RDRAND() for _ in range(NUM_STREAMS)]

.. end block