% \iffalse meta-comment
%
%% File: l3fp-random.dtx
%
% Copyright (C) 2016-2024 The LaTeX Project
%
% It may be distributed and/or modified under the conditions of the
% LaTeX Project Public License (LPPL), either version 1.3c of this
% license or (at your option) any later version.  The latest version
% of this license is in the file
%
%    https://www.latex-project.org/lppl.txt
%
% This file is part of the "l3kernel bundle" (The Work in LPPL)
% and all files in that bundle must be distributed together.
%
% -----------------------------------------------------------------------
%
% The development version of the bundle can be found at
%
%    https://github.com/latex3/latex3
%
% for those people who are interested.
%
%<*driver>
\documentclass[full,kernel]{l3doc}
\begin{document}
  \DocInput{\jobname.dtx}
\end{document}
%</driver>
% \fi
%
% \title{^^A
%   The \pkg{l3fp-random} module\\
%   Floating point random numbers
% }
% \author{^^A
%  The \LaTeX{} Project\thanks
%    {^^A
%      E-mail:
%        \href{mailto:latex-team@latex-project.org}
%          {latex-team@latex-project.org}^^A
%    }^^A
% }
% \date{Released 2024-12-25}
%
% \maketitle
%
% \begin{documentation}
%
% \end{documentation}
%
% \begin{implementation}
%
% \section{\pkg{l3fp-random} implementation}
%
%    \begin{macrocode}
%<*package>
%    \end{macrocode}
%
%    \begin{macrocode}
%<@@=fp>
%    \end{macrocode}
%
% \begin{macro}[EXP]{\@@_parse_word_rand:N , \@@_parse_word_randint:N}
%   Those functions may receive a variable number of arguments.  We
%   won't use the argument~|?|.
%    \begin{macrocode}
\cs_new:Npn \@@_parse_word_rand:N
  { \@@_parse_function:NNN \@@_rand_o:Nw ? }
\cs_new:Npn \@@_parse_word_randint:N
  { \@@_parse_function:NNN \@@_randint_o:Nw ? }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Engine support}
%
% Obviously, every word \enquote{random} below means
% \enquote{pseudo-random}, as we have no access to entropy (except a
% very unreliable source of entropy: the time it takes to run some
% code).
%
% The primitive random number generator (RNG) is provided as
% \cs{tex_uniformdeviate:D}.  Under the hood, it maintains an array of
% $55$ $28$-bit numbers, updated with a linear recursion relation
% (similar to Fibonacci numbers) modulo $2^{28}$.  When
% \cs{tex_uniformdeviate:D} \meta{integer} is called (for brevity denote
% by~$N$ the \meta{integer}), the next $28$-bit number is read from the
% array, scaled by $N/2^{28}$, and rounded.  To prevent $0$ and $N$ from
% appearing half as often as other numbers, they are both mapped to the
% result~$0$.
%
% This process means that \cs{tex_uniformdeviate:D} only gives a uniform
% distribution from $0$ to $N-1$ if $N$ is a divisor of $2^{28}$, so we
% will mostly call the RNG with such power of~$2$ arguments.  If $N$
% does not divide $2^{28}$, then the relative non-uniformity (difference
% between probabilities of getting different numbers) is about
% $N/2^{28}$.  This implies that detecting deviation from $1/N$ of the
% probability of a fixed value X requires about $2^{56}/N$ random
% trials. But collective patterns can reduce this to about
% $2^{56}/N^2$. For instance with $N=3\times 2^{k}$, the modulo~$3$
% repartition of such random numbers is biased with a non-uniformity
% about $2^k/2^{28}$ (which is much worse than the circa $3/2^{28}$
% non-uniformity from taking directly $N=3$).  This is detectable after
% about $2^{56}/2^{2k} = 9\cdot 2^{56}/N^2$ random numbers. For $k=15$,
% $N=98304$, this means roughly $2^{26}$ calls to the RNG
% (experimentally this takes at the very least 16 seconds on a 2 giga-hertz
% processor). While this bias is not quite problematic, it is
% uncomfortably close to being so, and it becomes worse as $N$ is
% increased.  In our code, we shall thus combine several results from
% the RNG\@.
%
% The RNG has three types of unexpected correlations.  First, everything
% is linear modulo~$2^{28}$, hence the lowest $k$ bits of the random
% numbers only depend on the lowest $k$ bits of the seed (and of course
% the number of times the RNG was called since setting the seed).  The
% recommended way to get a number from $0$ to $N-1$ is thus to scale the
% raw $28$-bit integer, as the engine's RNG does.  We will go further
% and in fact typically we discard some of the lowest bits.
%
% Second, suppose that we call the RNG with the same argument~$N$ to get
% a set of $K$ integers in $[0,N-1]$ (throwing away repeats), and
% suppose that $N>K^3$ and $K>55$.  The recursion used to construct more
% $28$-bit numbers from previous ones is linear:
% $x_n = x_{n-55} - x_{n-24}$ or $x_n = x_{n-55}-x_{n-24}+2^{28}$.
% After rescaling and rounding we find that the result $N_n\in[0,N-1]$
% is among $N_{n-55} - N_{n-24} + \{-1,0,1\}$ modulo~$N$ (a more
% detailed analysis shows that $0$ appears with frequency close
% to~$3/4$).  The resulting set thus has more triplets $(a,b,c)$ than
% expected obeying $a=b+c$ modulo~$N$.  Namely it will have of order
% $(K-55)\times 3/4$ such triplets, when one would expect $K^3/(6N)$.
% This starts to be detectable around $N=2^{18}>55^3$ (earlier if one
% keeps track of positions too, but this is more subtle than it looks
% because the array of $28$-bit integers is read backwards by the
% engine).  Hopefully the correlation is subtle enough to not affect
% realistic documents so we do not specifically mitigate against this.
% Since we typically use two calls to the RNG per \cs{int_rand:nn} we
% would need to investigate linear relations between the $x_{2n}$ on the
% one hand and between the $x_{2n+1}$ on the other hand.  Such relations
% will have more complicated coefficients than $\pm 1$, which alleviates
% the issue.
%
% Third, consider successive batches of $165$ calls to the RNG (with
% argument $2^{28}$ or with argument~$2$ for instance), then most
% batches have more odd than even numbers.  Note that this does not mean
% that there are more odd than even numbers overall.  Similar issues are
% discussed in Knuth's TAOCP volume 2 near exercise 3.3.2-31.  We do not
% have any mitigation strategy for this.
%
% Ideally, our algorithm should be:
% \begin{itemize}
% \item Uniform.  The result should be as uniform as possible assuming
%   that the RNG's underlying $28$-bit integers are uniform.
% \item Uncorrelated.  The result should not have detectable
%   correlations between different seeds, similar to the lowest-bit ones
%   mentioned earlier.
% \item Quick.  The algorithm should be fast in \TeX{}, so no
%   \enquote{bit twiddling}, but \enquote{digit twiddling} is ok.
% \item Simple.  The behaviour must be documentable precisely.
% \item Predictable.  The number of calls to the RNG should be the same
%   for any \cs{int_rand:nn}, because then the algorithm can be modified
%   later without changing the result of other uses of the RNG\@.
% \item Robust.  It should work even for \cs{int_rand:nn} |{| |-|
%   \cs{c_max_int} |}| |{| \cs{c_max_int} |}| where the range is not
%   representable as an integer.  In fact, we also provide later a
%   floating-point |randint| whose range can go all the way up to
%   $2\times 10^{16}-1$ possible values.
% \end{itemize}
% Some of these requirements conflict.  For instance, uniformity cannot
% be achieved with a fixed number of calls to the RNG\@.
%
% Denote by $\operatorname{random}(N)$ one call to
% \cs{tex_uniformdeviate:D} with argument~$N$, and by
% $\operatorname{ediv}(p,q)$ the \eTeX{} rounding division giving
% $\lfloor p/q+1/2\rfloor$.  Denote by $\meta{min}$, $\meta{max}$ and
% $R=\meta{max}-\meta{min}+1$ the arguments of \cs{int_min:nn} and the
% number of possible outcomes.  Note that $R\in [1,2^{32}-1]$ cannot
% necessarily be represented as an integer (however, $R-2^{31}$ can).
% Our strategy is to get two $28$-bit integers $X$ and $Y$ from the RNG,
% split each into $14$-bit integers, as $X=X_1\times 2^{14}+X_0$ and
% $Y=Y_1\times 2^{14}+Y_0$ then return essentially
% $\meta{min} + \lfloor R (X_1\times 2^{-14} + Y_1\times 2^{-28} +
% Y_0\times 2^{-42} + X_0\times 2^{-56})\rfloor$.  For small~$R$ the
% $X_0$ term has a tiny effect so we ignore it and we can compute
% $R\times Y/2^{28}$ much more directly by $\operatorname{random}(R)$.
% \begin{itemize}
% \item If $R \leq 2^{17}-1$ then return
%   $\operatorname{ediv}(R\operatorname{random}(2^{14}) +
%   \operatorname{random}(R) + 2^{13}, 2^{14}) - 1 + \meta{min}$.  The
%   shifts by $2^{13}$ and $-1$ convert \eTeX{} division to truncated
%   division.  The bound on $R$ ensures that the number obtained after
%   the shift is less than \cs{c_max_int}.  The non-uniformity is at
%   most of order $2^{17}/2^{42} = 2^{-25}$.
% \item Split $R=R_2\times 2^{28}+R_1\times 2^{14}+R_0$, where
%   $R_2\in [0,15]$.  Compute
%   $\meta{min} + R_2 X_1 2^{14} + (R_2 Y_1 + R_1 X_1) +
%   \operatorname{ediv}(R_2 Y_0 + R_1 Y_1 + R_0 X_1 +
%   \operatorname{ediv}(R_2 X_0 + R_0 Y_1 + \operatorname{ediv}((2^{14}
%   R_1 + R_0) (2^{14} Y_0 + X_0), 2^{28}), 2^{14}), 2^{14})$ then map a
%   result of $\meta{max}+1$ to $\meta{min}$.  Writing each
%   $\operatorname{ediv}$ in terms of truncated division with a shift,
%   and using
%   $\lfloor(p+\lfloor r/s\rfloor)/q\rfloor =
%   \lfloor(ps+r)/(sq)\rfloor$, what we compute is equal to
%   $\lfloor\meta{exact}+2^{-29}+2^{-15}+2^{-1}\rfloor$ with
%   $\meta{exact}=\meta{min} + R \times 0.X_1Y_1Y_0X_0$.  Given we map
%   $\meta{max}+1$ to $\meta{min}$, the shift has no effect on
%   uniformity.  The non-uniformity is bounded by $R/2^{56}<2^{-24}$.  It
%   may be possible to speed up the code by dropping tiny terms such as
%   $R_0 X_0$, but the analysis of non-uniformity proves too difficult.
%
%   To avoid the overflow when the computation yields $\meta{max}+1$
%   with $\meta{max}=2^{31}-1$ (note that $R$ is then arbitrary), we
%   compute the result in two pieces.  Compute
%   $\meta{first} = \meta{min} + R_2 X_1 2^{14}$ if $R_2<8$ or
%   $\meta{min} + 8 X_1 2^{14} + (R_2-8) X_1 2^{14}$ if $R_2\geq 8$, the
%   expressions being chosen to avoid overflow.  Compute
%   $\meta{second} = R_2 Y_1 + R_1 X_1 + \operatorname{ediv}({\dots})$,
%   at most
%   $R_2 2^{14} + R_1 2^{14} + R_0\leq 2^{28} + 15\times 2^{14} - 1$,
%   not at risk of overflowing.  We have
%   $\meta{first}+\meta{second}=\meta{max}+1=\meta{min}+R$ if and only
%   if $\meta{second} = R1 2^{14} + R_0 + R_2 2^{14}$ and
%   $2^{14} R_2 X_1 = 2^{28} R_2 - 2^{14} R_2$ (namely $R_2=0$ or
%   $X_1=2^{14}-1$).  In that case, return \meta{min}, otherwise return
%   $\meta{first}+\meta{second}$, which is safe because it is at most
%   \meta{max}.  Note that the decision of what to return does not need
%   \meta{first} explicitly so we don't actually compute it, just put it
%   in an integer expression in which \meta{second} is eventually added
%   (or not).
% \item To get a floating point number in $[0,1)$ just call the
%   $R=10000\leq 2^{17}-1$ procedure above to produce four blocks of four
%   digits.
% \item To get an integer floating point number in a range (whose size
%   can be up to $2\times 10^{16}-1$), work with fixed-point numbers:
%   get six times four digits to build a fixed point number, multiply by
%   $R$ and add $\meta{min}$.  This requires some care because
%   \pkg{l3fp-extended} only supports non-negative numbers.
% \end{itemize}
%
% \begin{variable}{\c__kernel_randint_max_int}
%   Constant equal to $2^{17}-1$, the maximal size of a range that
%   \cs{int_range:nn} can do with its \enquote{simple} algorithm.
%    \begin{macrocode}
\int_const:Nn \c__kernel_randint_max_int { 131071 }
%    \end{macrocode}
% \end{variable}
%
% \begin{macro}[EXP]{\__kernel_randint:n}
%   Used in an integer expression, \cs{__kernel_randint:n} |{|$R$|}|
%   gives a random number
%   $1+\lfloor(R\operatorname{random}(2^{14}) +
%   \operatorname{random}(R))/2^{14}\rfloor$ that is in $[1,R]$.
%   Previous code was computing $\lfloor p/2^{14}\rfloor$ as
%   $\operatorname{ediv}(p-2^{13},2^{14})$ but that wrongly gives $-1$
%   for $p=0$.
%    \begin{macrocode}
\cs_new:Npn \__kernel_randint:n #1
  {
    (#1 * \tex_uniformdeviate:D 16384
    + \tex_uniformdeviate:D #1 + 8192 ) / 16384
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[EXP]
%   {\@@_rand_myriads:n, \@@_rand_myriads_loop:w, \@@_rand_myriads_get:w}
%   Used as \cs{@@_rand_myriads:n} |{XXX}| with one letter |X|
%   (specifically) per block of four digit we want; it expands to |;|
%   followed by the requested number of brace groups, each containing
%   four (pseudo-random) digits.  Digits are produced as a random number
%   in $[10000,19999]$ for the usual reason of preserving leading zeros.
%    \begin{macrocode}
\cs_new:Npn \@@_rand_myriads:n #1
  { \@@_rand_myriads_loop:w #1 \prg_break: X \prg_break_point: ; }
\cs_new:Npn \@@_rand_myriads_loop:w #1 X
  {
    #1
    \exp_after:wN \@@_rand_myriads_get:w
    \int_value:w \@@_int_eval:w 9999 +
      \__kernel_randint:n { 10000 }
    \@@_rand_myriads_loop:w
  }
\cs_new:Npn \@@_rand_myriads_get:w 1 #1 ; { ; {#1} }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Random floating point}
%
% \begin{macro}[EXP]{\@@_rand_o:Nw, \@@_rand_o:w}
%   First we check that |random| was called without argument.  Then get
%   four blocks of four digits and convert that fixed point number to a
%   floating point number (this correctly sets the exponent).  This has
%   a minor bug: if all of the random numbers are zero then the result
%   is correctly~$0$ but it raises the \texttt{underflow} flag; it
%   should not do that.
%    \begin{macrocode}
\cs_new:Npn \@@_rand_o:Nw ? #1 @
  {
    \tl_if_empty:nTF {#1}
      {
        \exp_after:wN \@@_rand_o:w
        \exp:w \exp_end_continue_f:w
        \@@_rand_myriads:n { XXXX } { 0000 } { 0000 } ; 0
      }
      {
        \msg_expandable_error:nnnnn
          { fp } { num-args } { rand() } { 0 } { 0 }
        \exp_after:wN \c_nan_fp
      }
  }
\cs_new:Npn \@@_rand_o:w ;
  {
    \exp_after:wN \@@_sanitize:Nw
    \exp_after:wN 0
    \int_value:w \@@_int_eval:w \c_zero_int
      \@@_fixed_to_float_o:wN
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Random integer}
%
% \begin{macro}[EXP]{\@@_randint_o:Nw}
% \begin{macro}[EXP]
%   {
%     \@@_randint_default:w,
%     \@@_randint_badarg:w,
%     \@@_randint_o:w,
%     \@@_randint_auxi_o:ww,
%     \@@_randint_auxii:wn,
%     \@@_randint_auxiii_o:ww,
%     \@@_randint_auxiv_o:ww,
%     \@@_randint_auxv_o:w,
%   }
%     Enforce that there is one argument (then add first argument~$1$)
%     or two arguments.  Call \cs{@@_randint_badarg:w} on each; this
%     function inserts |1| \cs{exp_stop_f:} to end the \cs{if_case:w}
%     statement if either the argument is not an integer or if its
%     absolute value is $\geq 10^{16}$.  Also bail out if
%     \cs{@@_compare_back:ww} yields~|1|, meaning that the bounds are
%     not in the right order.  Otherwise an auxiliary converts each
%     argument times $10^{-16}$ (hence the shift in exponent) to a
%     $24$-digit fixed point number (see \pkg{l3fp-extended}).
%     Then compute the number of choices, $\meta{max}+1-\meta{min}$.
%     Create a random $24$-digit fixed-point number with
%     \cs{@@_rand_myriads:n}, then use a fused multiply-add instruction
%     to multiply the number of choices to that random number and add it
%     to \meta{min}.  Then truncate to $16$ digits (namely select the
%     integer part of $10^{16}$ times the result) before converting back
%     to a floating point number  (\cs{@@_sanitize:Nw} takes care of zero).
%     To avoid issues with negative numbers, add $1$ to all fixed point
%     numbers (namely $10^{16}$ to the integers they represent), except
%     of course when it is time to convert back to a float.
%    \begin{macrocode}
\cs_new:Npn \@@_randint_o:Nw ?
  {
    \@@_parse_function_one_two:nnw
      { randint }
      { \@@_randint_default:w \@@_randint_o:w }
  }
\cs_new:Npn \@@_randint_default:w #1 { \exp_after:wN #1 \c_one_fp }
\cs_new:Npn \@@_randint_badarg:w \s_@@ \@@_chk:w #1#2#3;
  {
    \@@_int:wTF \s_@@ \@@_chk:w #1#2#3;
      {
        \if_meaning:w 1 #1
          \if_int_compare:w
              \@@_use_i_until_s:nw #3 ; > \c_@@_prec_int
            \c_one_int
          \fi:
        \fi:
      }
      { \c_one_int }
  }
\cs_new:Npn \@@_randint_o:w #1; #2; @
  {
    \if_case:w
        \@@_randint_badarg:w #1;
        \@@_randint_badarg:w #2;
        \if:w 1 \@@_compare_back:ww #2; #1; \c_one_int \fi:
        \c_zero_int
      \@@_randint_auxi_o:ww #1; #2;
    \or:
      \@@_invalid_operation_tl_o:ff
        { randint } { \@@_array_to_clist:n { #1; #2; } }
      \exp:w
    \fi:
    \exp_after:wN \exp_end:
  }
\cs_new:Npn \@@_randint_auxi_o:ww #1 ; #2 ; #3 \exp_end:
  {
    \fi:
    \@@_randint_auxii:wn #2 ;
    { \@@_randint_auxii:wn #1 ; \@@_randint_auxiii_o:ww }
  }
\cs_new:Npn \@@_randint_auxii:wn \s_@@ \@@_chk:w #1#2#3#4 ;
  {
    \if_meaning:w 0 #1
      \exp_after:wN \use_i:nn
    \else:
      \exp_after:wN \use_ii:nn
    \fi:
    { \exp_after:wN \@@_fixed_continue:wn \c_@@_one_fixed_tl }
    {
      \exp_after:wN \@@_ep_to_fixed:wwn
      \int_value:w \@@_int_eval:w
        #3 - \c_@@_prec_int , #4 {0000} {0000} ;
      {
        \if_meaning:w 0 #2
          \exp_after:wN \use_i:nnnn
          \exp_after:wN \@@_fixed_add_one:wN
        \fi:
        \exp_after:wN \@@_fixed_sub:wwn \c_@@_one_fixed_tl
      }
      \@@_fixed_continue:wn
    }
  }
\cs_new:Npn \@@_randint_auxiii_o:ww #1 ; #2 ;
  {
    \@@_fixed_add:wwn #2 ;
      {0000} {0000} {0000} {0001} {0000} {0000} ;
    \@@_fixed_sub:wwn #1 ;
    {
      \exp_after:wN \use_i:nn
      \exp_after:wN \@@_fixed_mul_add:wwwn
      \exp:w \exp_end_continue_f:w \@@_rand_myriads:n { XXXXXX } ;
    }
    #1 ;
    \@@_randint_auxiv_o:ww
    #2 ;
    \@@_randint_auxv_o:w #1 ; @
  }
\cs_new:Npn \@@_randint_auxiv_o:ww #1#2#3#4#5 ; #6#7#8#9
  {
    \if_int_compare:w
      \if_int_compare:w #1#2 > #6#7 \exp_stop_f: 1 \else:
        \if_int_compare:w #1#2 < #6#7 \exp_stop_f: - \fi: \fi:
      #3#4 > #8#9 \exp_stop_f:
      \@@_use_i_until_s:nw
    \fi:
    \@@_randint_auxv_o:w {#1}{#2}{#3}{#4}#5
  }
\cs_new:Npn \@@_randint_auxv_o:w #1#2#3#4#5 ; #6 @
  {
    \exp_after:wN \@@_sanitize:Nw
    \int_value:w
    \if_int_compare:w #1 < 10000 \exp_stop_f:
      2
    \else:
      0
      \exp_after:wN \exp_after:wN
      \exp_after:wN \@@_reverse_args:Nww
    \fi:
    \exp_after:wN \@@_fixed_sub:wwn \c_@@_one_fixed_tl
    {#1} {#2} {#3} {#4} {0000} {0000} ;
    {
      \exp_after:wN \exp_stop_f:
      \int_value:w \@@_int_eval:w \c_@@_prec_int
        \@@_fixed_to_float_o:wN
    }
    0
    \exp:w \exp_after:wN \exp_end:
  }
%    \end{macrocode}
% \end{macro}
% \end{macro}
%
% \begin{macro}{\int_rand:nn, \@@_randint:ww}
%   Evaluate the argument and filter out the case where the lower
%   bound~|#1| is more than the upper bound~|#2|.  Then determine
%   whether the range is narrower than \cs{c__kernel_randint_max_int};
%   |#2-#1| may overflow for very large positive~|#2| and negative~|#1|.
%   If the range is narrow, call \cs{__kernel_randint:n} \Arg{choices}
%   where \meta{choices} is the number of possible outcomes.  If the
%   range is wide, use somewhat slower code.
%    \begin{macrocode}
\cs_new:Npn \int_rand:nn #1#2
  {
    \int_eval:n
      {
        \exp_after:wN \@@_randint:ww
        \int_value:w \int_eval:n {#1} \exp_after:wN ;
        \int_value:w \int_eval:n {#2} ;
      }
  }
\cs_new:Npn \@@_randint:ww #1; #2;
  {
    \if_int_compare:w #1 > #2 \exp_stop_f:
      \msg_expandable_error:nnnn
        { kernel } { randint-backward-range } {#1} {#2}
      \@@_randint:ww #2; #1;
    \else:
      \if_int_compare:w \@@_int_eval:w #2
          \if_int_compare:w #1 > \c_zero_int
            - #1 < \@@_int_eval:w
          \else:
            < \@@_int_eval:w #1 +
          \fi:
          \c__kernel_randint_max_int
          \@@_int_eval_end:
        \__kernel_randint:n
          { \@@_int_eval:w #2 - #1 + 1 \@@_int_eval_end: }
        - 1 + #1
      \else:
        \__kernel_randint:nn {#1} {#2}
      \fi:
    \fi:
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}
%   {
%     \__kernel_randint:nn, \@@_randint_split_o:Nw, \@@_randint_split_aux:w,
%     \@@_randinat_wide_aux:w, \@@_randinat_wide_auxii:w,
%   }
%   Any $n\in[-2^{31}+1,2^{31}-1]$ is uniquely written as
%   $2^{14}n_1+n_2$ with $n_1\in[-2^{17},2^{17}-1]$ and
%   $n_2\in[0,2^{14}-1]$.  Calling \cs{@@_randint_split_o:Nw} $n$ |;|
%   gives $n_1$|;| $n_2$|;| and expands the next token once.  We do this
%   for two random numbers and apply \cs{@@_randint_split_o:Nw} twice to
%   fully decompose the range~$R$.  One subtlety is that we compute
%   $R-2^{31}=\meta{max}-\meta{min}-(2^{31}-1)\in[-2^{31}+1,2^{31}-1]$
%   rather than $R$ to avoid overflow.
%
%   Then we have \cs{@@_randint_wide_aux:w} \meta{X_1}|;|\meta{X_0}|;|
%   \meta{Y_1}|;|\meta{Y_0}|;| \meta{R_2}|;|\meta{R_1}|;|\meta{R_0}|;.|
%   and we apply the algorithm described earlier.
%    \begin{macrocode}
\cs_new:Npn \__kernel_randint:nn #1#2
  {
    #1
    \exp_after:wN \@@_randint_wide_aux:w
    \int_value:w
      \exp_after:wN \@@_randint_split_o:Nw
      \tex_uniformdeviate:D 268435456 ;
    \int_value:w
      \exp_after:wN \@@_randint_split_o:Nw
      \tex_uniformdeviate:D 268435456 ;
    \int_value:w
      \exp_after:wN \@@_randint_split_o:Nw
      \int_value:w \@@_int_eval:w 131072 +
        \exp_after:wN \@@_randint_split_o:Nw
        \int_value:w
          \__kernel_int_add:nnn {#2} { -#1 } { -\c_max_int } ;
    .
  }
\cs_new:Npn \@@_randint_split_o:Nw #1#2 ;
  {
    \if_meaning:w 0 #1
      0 \exp_after:wN ; \int_value:w 0
    \else:
      \exp_after:wN \@@_randint_split_aux:w
      \int_value:w \@@_int_eval:w (#1#2 - 8192) / 16384 ;
      + #1#2
    \fi:
    \exp_after:wN ;
  }
\cs_new:Npn \@@_randint_split_aux:w #1 ;
  {
    #1 \exp_after:wN ;
    \int_value:w \@@_int_eval:w - #1 * 16384
  }
\cs_new:Npn \@@_randint_wide_aux:w #1;#2; #3;#4; #5;#6;#7; .
  {
    \exp_after:wN \@@_randint_wide_auxii:w
    \int_value:w \@@_int_eval:w #5 * #3 + #6 * #1 +
      (#5 * #4 + #6 * #3 + #7 * #1 +
        (#5 * #2 +           #7 * #3 +
          (16384 * #6 + #7) * (16384 * #4 + #2) / 268435456) / 16384
      ) / 16384 \exp_after:wN ;
    \int_value:w \@@_int_eval:w (#5 + #6) * 16384 + #7 ;
    #1 ; #5 ;
  }
\cs_new:Npn \@@_randint_wide_auxii:w #1; #2; #3; #4;
  {
    \if_int_odd:w 0
        \if_int_compare:w #1 = #2 \else: \exp_stop_f: \fi:
        \if_int_compare:w #4 = \c_zero_int 1 \fi:
        \if_int_compare:w #3 = 16383 ~ 1 \fi:
        \exp_stop_f:
      \exp_after:wN \prg_break:
    \fi:
    \if_int_compare:w #4 < 8 \exp_stop_f:
      + #4 * #3 * 16384
    \else:
      + 8 * #3 * 16384 + (#4 - 8) * #3 * 16384
    \fi:
    + #1
    \prg_break_point:
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\int_rand:n, \@@_randint:n}
%   Similar to \cs{int_rand:nn}, but needs fewer checks.
%    \begin{macrocode}
\cs_new:Npn \int_rand:n #1
  {
    \int_eval:n
      { \exp_args:Nf \@@_randint:n { \int_eval:n {#1} } }
  }
\cs_new:Npn \@@_randint:n #1
  {
    \if_int_compare:w #1 < \c_one_int
      \msg_expandable_error:nnnn
        { kernel } { randint-backward-range } { 1 } {#1}
      \@@_randint:ww #1; 1;
    \else:
      \if_int_compare:w #1 > \c__kernel_randint_max_int
        \__kernel_randint:nn { 1 } {#1}
      \else:
        \__kernel_randint:n {#1}
      \fi:
    \fi:
  }
%    \end{macrocode}
% \end{macro}
%
%    \begin{macrocode}
%</package>
%    \end{macrocode}
%
% \end{implementation}
%
% \PrintChanges
%
% \PrintIndex