diff options
Diffstat (limited to 'macros/latex/required/l3kernel/l3fp-random.dtx')
-rw-r--r-- | macros/latex/required/l3kernel/l3fp-random.dtx | 629 |
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 |