******************************************************
:f:mod:`spral_random` - Pseudo-random number generator
******************************************************
.. f:module:: spral_random
   :synopsis: Random number generator

=======
Purpose
=======

This package generates pseudo-random numbers using a linear congruential
generator. It should generate the same random numbers using any standards
compliant Fortran compiler on any architecture so long as the default
integer and real kinds are the same.

The seed can optionally be observed or specified by the user. Otherwise
a default seed of 486502 is used.

Version history
---------------

2016-09-08 Version 1.1.0
   Add support for long integers

2014-04-07 Version 1.0.0
   Initial release

========
Routines
========

Random Number Generation
------------------------

.. f:function:: random_real(state[, positive])

   Return a real uniformly at random from the interval :math:`(-1,1)`
   (positive =.true.) or :math:`(0,1)` (positive =.false.).

   :p random_state state [inout]: current state of RNG.
   :o logical positive [in,default=.false.]:
      if .true., sample from :math:`(0,1)`;
      otherwise, sample from :math:`(-1,1)`.
   :r random_real: Sampled value.
   :rtype random_real: real

.. f:function:: random_integer(state, n)

   Return an integer uniformly at random from the interval :math:`[1,n]`.

   :p random_state state [inout]: current state of the RNG.
   :p integer(kind) n [in]: largest value in range to be sampled. `kind` may be
      either default or long integer (return type will match).
   :r random_integer: Sampled value.
   :rtype random_integer: integer(kind)

.. f:function:: random_logical(state)

   Return a logical with equal probability of being .true. or .false..

   :p random_state state [inout]: current state of the RNG.
   :r random_logical: Sampled value.
   :rtype random_logical: logical

Get/Set Random Seed
-------------------

.. f:function:: random_get_seed(state)
   
   Return the current random seed stored in state.

   The stream of random numbers generated after this call can be reproduced
   through the same sequence of calls after seed has been passed to
   :f:subr:`random_set_seed`.

   :p random_state state [in]: state variable to extract state from
   :r random_get_seed: current random seed
   :rtype: integer

.. f:subroutine:: random_set_seed(state, seed)

   Set the random seed stored in the state variable

   :p random_state state [inout]: state variable to set seed for.
   :p integer seed [in]: new seed.

==========
Data Types
==========

.. f:type:: random_state

   State of the random number generator.
   Compontents are not available to the user, but may be examined and
   altered through calls to :f:func:`random_get_seed` and
   :f:func:`random_set_seed` respectively.

=======
Example
=======
The following code:

.. literalinclude:: ../../examples/Fortran/random.f90
   :language: Fortran

Produces the following output::

   Some random values
   Sample Unif(-1,1)               =   0.951878630556
   Sample Unif(0,1)                =   0.395779648796
   Sample Unif(1, ..., 20)         =                3
   Sample Unif(1, ..., 20*huge(0)) =      33572664025
   Sample B(1,0.5)                 =                F

   The same random values again
   Sample Unif(-1,1)               =   0.951878630556
   Sample Unif(0,1)                =   0.395779648796
   Sample Unif(1, ..., 20)         =                3
   Sample Unif(1, ..., 20*huge(0)) =      33572664025
   Sample B(1,0.5)                 =                F

======
Method
======
We use a linear congruential generator of the following form:

.. math::

   X_{n+1} = (aX_n + c)\quad \mathrm{mod}\; m

with the following constants

.. math::

   a = 1103515245, \\
   c = 12345, \\
   m = 2^{31}.

According to Wikipedia, this is the same as used in glibc.

The LCG is evolved before each sample is taken, and the sample is based on the
new value.

The routines :f:func:`random_get_seed` and :f:func:`random_set_seed` allow the
user to get and set the current value of :math:`X_n`. The default seed is
:math:`X_0 = 486502`.

In :f:func:`random_real`
------------------------

Samples from :math:`\mathrm{Unif}(0,1)` are generated as

.. math::

   \frac{\mathrm{real}(X_n)}{\mathrm{real}(m)},

and samples from :math:`\mathrm{Unif}(-1,1)` are generated as

.. math::

   1−\frac{\mathrm{real}(2X_n)}{\mathrm{real}(m)}.

In :f:func:`random_integer`
---------------------------

Samples from :math:`\mathrm{Unif}(1,\ldots,n)` are generated as

.. math::

   \mathrm{int}\left(X_n\frac{\mathrm{real}(n)}{\mathrm{real}(m)}\right) + 1

In :f:func:`random_logical`
----------------------------

Returns the value of the Fortran expression

.. code-block:: Fortran

   (1 .eq. random_integer(state,2))