summaryrefslogtreecommitdiff
path: root/macros/latex/required/l3kernel/l3fp-random.dtx
diff options
context:
space:
mode:
Diffstat (limited to 'macros/latex/required/l3kernel/l3fp-random.dtx')
-rw-r--r--macros/latex/required/l3kernel/l3fp-random.dtx629
1 files changed, 629 insertions, 0 deletions
diff --git a/macros/latex/required/l3kernel/l3fp-random.dtx b/macros/latex/required/l3kernel/l3fp-random.dtx
new file mode 100644
index 0000000000..ce5a0c48f7
--- /dev/null
+++ b/macros/latex/required/l3kernel/l3fp-random.dtx
@@ -0,0 +1,629 @@
+% \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-04-11}
+%
+% \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