summaryrefslogtreecommitdiff
path: root/macros/latex/contrib/l3kernel/l3fp-random.dtx
blob: 17bd04d9eed15084ffc82c5cf561fd424e340c63 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
% \iffalse meta-comment
%
%% File: l3fp-random.dtx
%
% Copyright (C) 2016-2020 The LaTeX3 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{The \textsf{l3fp-random} package\\
%   Floating point random numbers}
% \author{^^A
%  The \LaTeX3 Project\thanks
%    {^^A
%      E-mail:
%        \href{mailto:latex-team@latex-project.org}
%          {latex-team@latex-project.org}^^A
%    }^^A
% }
% \date{Released 2020-03-06}
%
% \maketitle
%
% \begin{documentation}
%
% \end{documentation}
%
% \begin{implementation}
%
% \section{\pkg{l3fp-random} Implementation}
%
%    \begin{macrocode}
%<*initex|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}
%
% Most engines provide random numbers, but not all.  We write the test
% twice simply in order to write the \texttt{false} branch first.
%    \begin{macrocode}
\sys_if_rand_exist:F
  {
    \__kernel_msg_new:nnn { kernel } { fp-no-random }
      { Random~numbers~unavailable~for~#1 }
    \cs_new:Npn \@@_rand_o:Nw ? #1 @
      {
        \__kernel_msg_expandable_error:nnn { kernel } { fp-no-random }
          { fp~rand }
        \exp_after:wN \c_nan_fp
      }
    \cs_new_eq:NN \@@_randint_o:Nw \@@_rand_o:Nw
    \cs_new:Npn \int_rand:nn #1#2
      {
        \__kernel_msg_expandable_error:nnn { kernel } { fp-no-random }
          { \int_rand:nn {#1} {#2} }
        \int_eval:n {#1}
      }
    \cs_new:Npn \int_rand:n #1
      {
        \__kernel_msg_expandable_error:nnn { kernel } { fp-no-random }
          { \int_rand:n {#1} }
        1
      }
  }
\sys_if_rand_exist:T
  {
%    \end{macrocode}
%
% 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
          }
          {
            \__kernel_msg_expandable_error:nnnnn
              { kernel } { 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
                1 \exp_stop_f:
              \fi:
            \fi:
          }
          { 1 \exp_stop_f: }
      }
    \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; 1 \exp_stop_f: \fi:
            0 \exp_stop_f:
          \@@_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:
          \__kernel_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 < 1 \exp_stop_f:
          \__kernel_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}
%
% End the initial conditional that ensures these commands are only
% defined in engines that support random numbers.
%    \begin{macrocode}
  }
%    \end{macrocode}
%
%    \begin{macrocode}
%</initex|package>
%    \end{macrocode}
%
% \end{implementation}
%
% \PrintChanges
%
% \PrintIndex