-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathexporand.html
executable file
·3479 lines (2910 loc) · 326 KB
/
exporand.html
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
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html><html xmlns:dc="http://purl.org/dc/terms/" itemscope itemtype="http://schema.org/Article"><head><meta http-equiv=Content-Type content="text/html; charset=utf-8"><meta name="citation_pdf_url" content="https://peteroupc.github.io/exporand.pdf"><meta name="citation_url" content="https://peteroupc.github.io/exporand.html"><meta name="citation_date" content="2025/01/17"><meta name="citation_online_date" content="2025/01/17"><meta name="og:type" content="article"><meta name="og:url" content="https://peteroupc.github.io/exporand.html"><meta name="og:site_name" content="peteroupc.github.io"><meta name="author" content="Peter Occil"/><meta name="citation_author" content="Peter Occil"/><meta name="viewport" content="width=device-width"><link rel=stylesheet type="text/css" href="/style.css">
<script type="text/x-mathjax-config"> MathJax.Hub.Config({"HTML-CSS": { availableFonts: ["STIX","TeX"], linebreaks: { automatic:true }, preferredFont: "TeX" },
tex2jax: { displayMath: [ ["$$","$$"], ["\\[", "\\]"] ], inlineMath: [ ["$", "$"], ["\\\\(","\\\\)"] ], processEscapes: true } });
</script><script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-AMS_HTML-full"></script></head><body> <div class="header">
<nav><p><a href="#navigation">Menu</a> - <a href="#top">Top</a> - <a href="/">Home</a></nav></div>
<div class="mainarea" id="top">
<h1 id="partially-sampled-random-numbers-for-accurate-sampling-of-continuous-distributions">Partially-Sampled Random Numbers for Accurate Sampling of Continuous Distributions</h1>
<p><a href="mailto:poccil14@gmail.com"><strong>Peter Occil</strong></a></p>
<p><strong>2020 Mathematics Subject Classification:</strong> 68W20, 60-08, 60-04.</p>
<p><a id="Introduction"></a></p>
<h2 id="introduction">Introduction</h2>
<p>This page introduces a implementation of <em>partially-sampled random numbers</em> (PSRNs) in the Python programming language. Although structures for PSRNs were largely described before this work, this document unifies the concepts for these kinds of numbers from prior works and shows how they can be used to sample the beta distribution (for most sets of parameters), the exponential distribution (with an arbitrary rate parameter), and many other continuous distributions—</p>
<ul>
<li>while avoiding floating-point arithmetic, and</li>
<li>to an arbitrary precision and with user-specified error bounds (and thus in an “exact” manner in the sense defined in (Karney 2016)<sup id="fnref:1" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>).</li>
</ul>
<p>For instance, these two points distinguish the beta sampler in this document from any other specially-designed beta sampler I am aware of. As for the exponential distribution, there are papers that discuss generating exponential random variates using random bits (Flajolet and Saheb 1982)<sup id="fnref:2" role="doc-noteref"><a href="#fn:2" class="footnote" rel="footnote">2</a></sup>, (Karney 2016)<sup id="fnref:1:1" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>, (Devroye and Gravel 2020)<sup id="fnref:3" role="doc-noteref"><a href="#fn:3" class="footnote" rel="footnote">3</a></sup>, (Thomas and Luk 2008)<sup id="fnref:4" role="doc-noteref"><a href="#fn:4" class="footnote" rel="footnote">4</a></sup>, but most if not all of them don’t deal with generating exponential PSRNs using an arbitrary rate, not just 1. (Habibizad Navin et al., 2010)<sup id="fnref:5" role="doc-noteref"><a href="#fn:5" class="footnote" rel="footnote">5</a></sup> is perhaps an exception; however the approach appears to involve pregenerated tables of digit probabilities.</p>
<p>The samplers discussed here also draw on work dealing with a construct called the <em>Bernoulli factory</em> (Keane and O’Brien 1994)<sup id="fnref:6" role="doc-noteref"><a href="#fn:6" class="footnote" rel="footnote">6</a></sup> (Flajolet et al., 2010)<sup id="fnref:7" role="doc-noteref"><a href="#fn:7" class="footnote" rel="footnote">7</a></sup>, which can simulate a new probability given a coin that shows heads with an unknown probability. One important feature of Bernoulli factories is that they can simulate a given probability <em>exactly</em>, without having to calculate that probability manually, which is important if the probability can be an irrational number whose value no computer can compute exactly (such as <code>sqrt(p)</code> or <code>exp(-2)</code>).</p>
<p>This page shows <a href="#Sampler_Code"><strong>Python code</strong></a> for these samplers.</p>
<p><a id="About_This_Document"></a></p>
<h3 id="about-this-document">About This Document</h3>
<p><strong>This is an open-source document; for an updated version, see the</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/raw/master/exporand.md"><strong>source code</strong></a> <strong>or its</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/blob/master/exporand.md"><strong>rendering on GitHub</strong></a><strong>. You can send comments on this document either on</strong> <a href="https://www.codeproject.com/Articles/5272482/Partially-Sampled-Random-Numbers-for-Accurate-Samp"><strong>CodeProject</strong></a> <strong>or on the</strong> <a href="https://github.com/peteroupc/peteroupc.github.io/issues"><strong>GitHub issues page</strong></a><strong>.</strong></p>
<p>My audience for this article is <strong>computer programmers with mathematics knowledge, but little or no familiarity with calculus</strong>.</p>
<p>I encourage readers to implement any of the algorithms given in this page, and report their implementation experiences. In particular, <a href="https://github.com/peteroupc/peteroupc.github.io/issues/18"><strong>I seek comments on the following aspects</strong></a>:</p>
<ul>
<li>Are the algorithms in this article easy to implement? Is each algorithm written so that someone could write code for that algorithm after reading the article?</li>
<li>Does this article have errors that should be corrected?</li>
<li>Are there ways to make this article more useful to the target audience?</li>
</ul>
<p>Comments on other aspects of this document are welcome.</p>
<p><a id="Contents"></a></p>
<h2 id="contents">Contents</h2>
<ul>
<li><a href="#Introduction"><strong>Introduction</strong></a>
<ul>
<li><a href="#About_This_Document"><strong>About This Document</strong></a></li>
</ul>
</li>
<li><a href="#Contents"><strong>Contents</strong></a></li>
<li><a href="#About_the_Beta_Distribution"><strong>About the Beta Distribution</strong></a></li>
<li><a href="#About_the_Exponential_Distribution"><strong>About the Exponential Distribution</strong></a></li>
<li><a href="#About_Partially_Sampled_Random_Numbers"><strong>About Partially-Sampled Random Numbers</strong></a>
<ul>
<li><a href="#Uniform_Partially_Sampled_Random_Numbers"><strong>Uniform Partially-Sampled Random Numbers</strong></a></li>
<li><a href="#Exponential_Partially_Sampled_Random_Numbers"><strong>Exponential Partially-Sampled Random Numbers</strong></a></li>
<li><a href="#Other_Distributions"><strong>Other Distributions</strong></a></li>
<li><a href="#Properties"><strong>Properties</strong></a></li>
<li><a href="#Limitations"><strong>Limitations</strong></a></li>
<li><a href="#Relation_to_Constructive_Reals"><strong>Relation to Constructive Reals</strong></a></li>
</ul>
</li>
<li><a href="#Sampling_Uniform_and_Exponential_PSRNs"><strong>Sampling Uniform and Exponential PSRNs</strong></a>
<ul>
<li><a href="#Sampling_Uniform_PSRNs"><strong>Sampling Uniform PSRNs</strong></a></li>
<li><a href="#Sampling_E_rands"><strong>Sampling E-rands</strong></a></li>
</ul>
</li>
<li><a href="#Arithmetic_and_Comparisons_with_PSRNs"><strong>Arithmetic and Comparisons with PSRNs</strong></a>
<ul>
<li><a href="#Addition_and_Subtraction"><strong>Addition and Subtraction</strong></a></li>
<li><a href="#Multiplication"><strong>Multiplication</strong></a></li>
<li><a href="#Reciprocal_and_Division"><strong>Reciprocal and Division</strong></a></li>
<li><a href="#Using_the_Arithmetic_Algorithms"><strong>Using the Arithmetic Algorithms</strong></a></li>
<li><a href="#Comparisons"><strong>Comparisons</strong></a></li>
<li><a href="#Discussion"><strong>Discussion</strong></a></li>
</ul>
</li>
<li><a href="#Building_Blocks"><strong>Building Blocks</strong></a>
<ul>
<li><a href="#SampleGeometricBag"><strong>SampleGeometricBag</strong></a></li>
<li><a href="#FillGeometricBag"><strong>FillGeometricBag</strong></a></li>
<li><a href="#kthsmallest"><strong>kthsmallest</strong></a></li>
<li><a href="#Power_of_Uniform_Subalgorithm"><strong>Power-of-Uniform Subalgorithm</strong></a></li>
</ul>
</li>
<li><a href="#Algorithms_for_the_Beta_and_Exponential_Distributions"><strong>Algorithms for the Beta and Exponential Distributions</strong></a>
<ul>
<li><a href="#Beta_Distribution"><strong>Beta Distribution</strong></a></li>
<li><a href="#Exponential_Distribution"><strong>Exponential Distribution</strong></a></li>
</ul>
</li>
<li><a href="#Sampler_Code"><strong>Sampler Code</strong></a></li>
<li><a href="#Correctness_Testing"><strong>Correctness Testing</strong></a>
<ul>
<li><a href="#Beta_Sampler"><strong>Beta Sampler</strong></a></li>
<li><a href="#ExpRandFill"><strong>ExpRandFill</strong></a></li>
<li><a href="#ExpRandLess"><strong>ExpRandLess</strong></a></li>
</ul>
</li>
<li><a href="#Accurate_Simulation_of_Continuous_Distributions"><strong>Accurate Simulation of Continuous Distributions</strong></a>
<ul>
<li><a href="#General_Arbitrary_Precision_Samplers"><strong>General Arbitrary-Precision Samplers</strong></a>
<ul>
<li><a href="#Uniform_Distribution_Inside_N_Dimensional_Shapes"><strong>Uniform Distribution Inside N-Dimensional Shapes</strong></a></li>
<li><a href="#Building_an_Arbitrary_Precision_Sampler"><strong>Building an Arbitrary-Precision Sampler</strong></a></li>
<li><a href="#Continuous_Distributions_Supported_on_0_to_1"><strong>Continuous Distributions Supported on 0 to 1</strong></a></li>
</ul>
</li>
<li><a href="#Specific_Arbitrary_Precision_Samplers"><strong>Specific Arbitrary-Precision Samplers</strong></a>
<ul>
<li><a href="#Rayleigh_Distribution"><strong>Rayleigh Distribution</strong></a></li>
<li><a href="#Hyperbolic_Secant_Distribution"><strong>Hyperbolic Secant Distribution</strong></a></li>
<li><a href="#Sum_of_Uniform_Random_Variates"><strong>Sum of Uniform Random Variates</strong></a></li>
<li><a href="#Ratio_of_Two_Uniform_Random_Variates"><strong>Ratio of Two Uniform Random Variates</strong></a></li>
<li><a href="#Reciprocal_of_Uniform_Random_Variate"><strong>Reciprocal of Uniform Random Variate</strong></a></li>
<li><a href="#Reciprocal_of_Power_of_Uniform_Random_Variate"><strong>Reciprocal of Power of Uniform Random Variate</strong></a></li>
<li><a href="#Distribution_of__U__1_minus__U"><strong>Distribution of <em>U</em>/(1−<em>U</em>)</strong></a></li>
<li><a href="#Arc_Cosine_Distribution"><strong>Arc-Cosine Distribution</strong></a></li>
<li><a href="#Logistic_Distribution"><strong>Logistic Distribution</strong></a></li>
<li><a href="#Cauchy_Distribution"><strong>Cauchy Distribution</strong></a></li>
<li><a href="#Exponential_Distribution_with_Unknown_Small_Rate"><strong>Exponential Distribution with Unknown Small Rate</strong></a></li>
<li><a href="#Exponential_Distribution_with_Rate_ln__x"><strong>Exponential Distribution with Rate ln(<em>x</em>)</strong></a></li>
<li><a href="#Lindley_Distribution_and_Lindley_Like_Mixtures"><strong>Lindley Distribution and Lindley-Like Mixtures</strong></a></li>
<li><a href="#Gamma_Distribution"><strong>Gamma Distribution</strong></a></li>
<li><a href="#One_Dimensional_Epanechnikov_Kernel"><strong>One-Dimensional Epanechnikov Kernel</strong></a></li>
<li><a href="#Uniform_Distribution_Inside_Rectellipse"><strong>Uniform Distribution Inside Rectellipse</strong></a></li>
<li><a href="#Tulap_distribution"><strong>Tulap distribution</strong></a></li>
<li><a href="#Continuous_Bernoulli_Distribution"><strong>Continuous Bernoulli Distribution</strong></a></li>
<li><a href="#Exchangeable_Farlie_ndash_Gumbel_ndash_Morgenstern_copula"><strong>Exchangeable Farlie–Gumbel–Morgenstern copula</strong></a></li>
</ul>
</li>
</ul>
</li>
<li><a href="#Complexity"><strong>Complexity</strong></a>
<ul>
<li><a href="#General_Principles"><strong>General Principles</strong></a></li>
<li><a href="#Complexity_of_Specific_Algorithms"><strong>Complexity of Specific Algorithms</strong></a></li>
</ul>
</li>
<li><a href="#Application_to_Weighted_Reservoir_Sampling"><strong>Application to Weighted Reservoir Sampling</strong></a></li>
<li><a href="#Acknowledgments"><strong>Acknowledgments</strong></a></li>
<li><a href="#Other_Documents"><strong>Other Documents</strong></a></li>
<li><a href="#Notes"><strong>Notes</strong></a></li>
<li><a href="#Appendix"><strong>Appendix</strong></a>
<ul>
<li><a href="#Equivalence_of_SampleGeometricBag_Algorithms"><strong>Equivalence of SampleGeometricBag Algorithms</strong></a></li>
<li><a href="#UniformMultiply_Algorithm"><strong>UniformMultiply Algorithm</strong></a></li>
<li><a href="#Uniform_of_Uniforms_Produces_a_Product_of_Uniforms"><strong>Uniform of Uniforms Produces a Product of Uniforms</strong></a></li>
<li><a href="#Oberhoff_s_Exact_Rejection_Sampling_Method"><strong>Oberhoff’s “Exact Rejection Sampling” Method</strong></a></li>
<li><a href="#Probability_Transformations"><strong>Probability Transformations</strong></a></li>
<li><a href="#Ratio_of_Uniforms"><strong>Ratio of Uniforms</strong></a></li>
<li><a href="#Setting_Digits_by_Digit_Probabilities"><strong>Setting Digits by Digit Probabilities</strong></a></li>
<li><a href="#Security_Considerations"><strong>Security Considerations</strong></a></li>
</ul>
</li>
<li><a href="#License"><strong>License</strong></a></li>
</ul>
<p><a id="About_the_Beta_Distribution"></a></p>
<h2 id="about-the-beta-distribution">About the Beta Distribution</h2>
<p>The <a href="https://en.wikipedia.org/wiki/Beta_distribution"><strong>beta distribution</strong></a> is a bounded-domain probability distribution; its two parameters, <code>alpha</code> and <code>beta</code>, are both greater than 0 and describe the distribution’s shape. Depending on <code>alpha</code> and <code>beta</code>, the shape can be a smooth peak or a smooth valley. The beta distribution can take on values in the interval [0, 1]. Any value in this interval (<code>x</code>) can occur with a probability proportional to—</p>
<pre>pow(x, alpha - 1) * pow(1 - x, beta - 1). (1)
</pre>
<p>Although <code>alpha</code> and <code>beta</code> can each be greater than 0, the sampler presented in this document only works if—</p>
<ul>
<li>both parameters are 1 or greater, or</li>
<li>in the case of base-2 numbers, one parameter equals 1 and the other is greater than 0.</li>
</ul>
<p><a id="About_the_Exponential_Distribution"></a></p>
<h2 id="about-the-exponential-distribution">About the Exponential Distribution</h2>
<p>The <em>exponential distribution</em> takes a parameter <em>λ</em>. Informally speaking, a random variate that follows an exponential distribution is the number of units of time between one event and the next, and <em>λ</em> is the expected average number of events per unit of time. Usually, <em>λ</em> is equal to 1.</p>
<p>An exponential random variate is commonly generated as follows: <code>-ln(1 - X) / lamda</code>, where <code>X</code> is a uniformly-distributed random real number in the interval (0, 1). (This particular algorithm, however, is not robust in practice, for reasons that are outside the scope of this document, but see (Pedersen 2018)<sup id="fnref:8" role="doc-noteref"><a href="#fn:8" class="footnote" rel="footnote">8</a></sup>.) This page presents an alternative way to sample exponential random variates.</p>
<p><a id="About_Partially_Sampled_Random_Numbers"></a></p>
<h2 id="about-partially-sampled-random-numbers">About Partially-Sampled Random Numbers</h2>
<p>In this document, a <em>partially-sampled random number</em> (PSRN) is a data structure that stores a real number of unlimited precision, but whose contents are sampled only when necessary. PSRNs open the door to algorithms that sample a random variate that “exactly” follows a probability distribution, <em>with arbitrary precision</em>, and <em>without floating-point arithmetic</em> (see “<a href="#Properties"><strong>Properties</strong></a>” later in this section).</p>
<p>PSRNs specified here consist of the following three things:</p>
<ul>
<li>
<p>A <em>fractional part</em> with an arbitrary number of digits. This can be implemented as an array of digits or as a packed integer containing all the digits. Some algorithms care whether those digits were <em>sampled</em> or <em>unsampled</em>; in that case, if a digit is unsampled, its unsampled status can be noted in a way that distinguishes it from sampled digits (for example, by using the <code>None</code> keyword in the Python programming language, or the number −1, or by storing a separate bit array indicating which bits are sampled and unsampled). The base in which all the digits are stored (such as base 10 for decimal or base 2 for binary) is arbitrary. The fractional part’s digits form a so-called <em>digit expansion</em> (for example, <em>binary expansion</em> in the case of binary or base-2 digits). Digits beyond those stored in the fractional part are unsampled.</p>
<p>For example, if the fractional part stores the base-10 digits [1, 3, 5], in that order, then it represents a random variate in the interval [0.135, 0.136], reflecting the fact that the digits between 0.135 and 0.136 are unknown.</p>
</li>
<li>An optional <em>integer part</em> (more specifically, the integer part of the number’s absolute value, that is, <code>floor(abs(x))</code>).</li>
<li>An optional <em>sign</em> (positive or negative).</li>
</ul>
<p>If the integer part is not stored, it’s assumed to be 0. If the sign is not stored, it’s assumed to be positive. For example, an implementation can care only about PSRNs in the interval [0, 1] by storing only a fractional part.</p>
<p>PSRNs ultimately represent a random variate between two numbers; one of the variate’s two bounds has the following form: sign * (integer part + fractional part), which is a lower bound if the PSRN is positive, or an upper bound if it’s negative. For example, if the PSRN stores a positive sign, the integer 3, and the fractional part [3, 5, 6] (in base 10), then the PSRN represents a random variate in the interval [3.356, 3.357]. Here, one of the bounds is built using the PSRN’s sign, integer part, and fractional part, and because the PSRN is positive, this is a lower bound.</p>
<p>This section specifies two kinds of PSRNs: uniform and exponential.</p>
<p><a id="Uniform_Partially_Sampled_Random_Numbers"></a></p>
<h3 id="uniform-partially-sampled-random-numbers">Uniform Partially-Sampled Random Numbers</h3>
<p>The most trivial example of a PSRN is that of the uniform distribution.</p>
<p>In this document, a <strong>uniform PSRN</strong> is a PSRN that represents a uniform random variate in a given interval. For example, if the PSRN is 3.356…, then it represents a uniformly distributed random variate in the interval [3.356, 3.357]. A uniform PSRN has the property that each additional digit of its fractional part (in this example, <code>.356...</code>) is sampled simply by setting it to an independent uniform random digit, an observation that dates from von Neumann (1951)<sup id="fnref:9" role="doc-noteref"><a href="#fn:9" class="footnote" rel="footnote">9</a></sup> in the binary case.<sup id="fnref:10" role="doc-noteref"><a href="#fn:10" class="footnote" rel="footnote">10</a></sup></p>
<ul>
<li>Flajolet et al. (2010)<sup id="fnref:7:1" role="doc-noteref"><a href="#fn:7" class="footnote" rel="footnote">7</a></sup> use the term <em>geometric bag</em> to refer to a uniform PSRN in the interval [0, 1] that stores binary (base-2) digits, some of which may be unsampled. In this case, the PSRN can consist of just a fractional part, which can be implemented as described earlier.</li>
<li>Karney (2016)<sup id="fnref:1:2" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup> uses the term <em>u-rand</em> to refer to uniform PSRNs that can store a sign, integer part, and a fractional part, where the base of the fractional part’s digits is arbitrary, but Karney’s concept only contemplates sampling digits from left to right without any gaps.</li>
</ul>
<p>A uniform PSRN remains a uniform PSRN even if it was generated using a nonuniform random sampling algorithm (such as Karney’s algorithm for the normal distribution).</p>
<p><a id="Exponential_Partially_Sampled_Random_Numbers"></a></p>
<h3 id="exponential-partially-sampled-random-numbers">Exponential Partially-Sampled Random Numbers</h3>
<p>In this document, an <strong>exponential PSRN</strong> (or <strong><em>e-rand</em></strong>, named similarly to Karney’s “u-rands” for partially-sampled uniform random variates (Karney 2016)<sup id="fnref:1:3" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>) samples each bit that, when combined with the existing bits, results in an exponentially-distributed random variate of the given rate. Also, because <code>-ln(1 - X)</code>, where <code>X</code> is a uniform random variate between 0 and 1, is exponentially distributed, e-rands can also represent the natural logarithm of a partially-sampled uniform random variate in (0, 1]. The difference here is that additional bits are sampled with an equal probability of being 1 or 0 each, but with a probability that gets closer to an equal probability as more bits are sampled. (More specifically, an exponential PSRN generally represents an exponentially-distributed random variate in a given interval.)</p>
<p>Algorithms for sampling e-rands are given in the section “Algorithms for the Beta and Exponential Distributions”.</p>
<p><a id="Other_Distributions"></a></p>
<h3 id="other-distributions">Other Distributions</h3>
<p>PSRNs of other distributions can be implemented via rejection from the uniform distribution. Examples include the following:</p>
<ul>
<li>The beta and continuous Bernoulli distributions, as discussed later in this document.</li>
<li>The standard normal distribution, as shown in (Karney 2016)<sup id="fnref:1:4" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup> by running Karney’s Algorithm N and filling unsampled digits uniformly at random, or as shown in an improved version of that algorithm by Du et al. (2020)<sup id="fnref:11" role="doc-noteref"><a href="#fn:11" class="footnote" rel="footnote">11</a></sup>.</li>
<li>Sampling uniform distributions in [0, <em>n</em>) (not just [0, 1]), is described later in “<a href="#Sampling_Uniform_PSRNs"><strong>Sampling Uniform PSRNs</strong></a>”.)</li>
</ul>
<p>For all these distributions, the PSRN’s unsampled trailing digits converge to the uniform distribution, and this also applies to any continuous distribution with a continuous probability density function (or, more generally, to so-called “absolutely continuous”<sup id="fnref:12" role="doc-noteref"><a href="#fn:12" class="footnote" rel="footnote">12</a></sup> distributions) (Oberhoff 2018)<sup id="fnref:13" role="doc-noteref"><a href="#fn:13" class="footnote" rel="footnote">13</a></sup>, (Hill and Schürger 2005, Corollary 4.4)<sup id="fnref:14" role="doc-noteref"><a href="#fn:14" class="footnote" rel="footnote">14</a></sup>.</p>
<p>PSRNs could also be implemented via rejection from the exponential distribution.</p>
<p><a id="Properties"></a></p>
<h3 id="properties">Properties</h3>
<p>An algorithm that samples from a nondiscrete distribution<sup id="fnref:15" role="doc-noteref"><a href="#fn:15" class="footnote" rel="footnote">15</a></sup> using PSRNs has the following properties:</p>
<ol>
<li>For randomness, the algorithm relies only on a <em>fair coin</em> (something that chooses either 1 or 0 with equal probability, independently of other choices).</li>
<li>The algorithm does not rely on floating-point arithmetic or fixed-precision approximations of irrational numbers or transcendental functions. (The algorithm may calculate approximations that converge to an irrational number, as long as those approximations are rational numbers of arbitrary precision. However, the more implicitly the algorithm works with irrational numbers or transcendental functions, the better.)</li>
<li>The algorithm may use rational arithmetic (such as <code>Fraction</code> in Python or <code>Rational</code> in Ruby), as long as the arithmetic is exact.</li>
<li>If the algorithm gives out a PSRN, the number represented by the sampled digits must follow a distribution that is close to the algorithm’s ideal distribution by a distance of not more than 1/(<em>b</em><sup><em>m</em></sup>), where <em>b</em> is the PSRN’s base, or radix (such as 2 for binary), and <em>m</em> is the position, starting from 1, of the rightmost sampled digit of the PSRN’s fractional part. ((Devroye and Gravel 2020)<sup id="fnref:3:1" role="doc-noteref"><a href="#fn:3" class="footnote" rel="footnote">3</a></sup> suggests Wasserstein distance, or “earth-mover distance”, as the distance to use for this purpose.)</li>
<li>After the algorithm gives out a PSRN, if the algorithm’s caller fills the PSRN’s unsampled fractional digits at random (for example, uniformly at random in the case of a uniform PSRN), so that the PSRN’s fractional part’s first <em>m</em> digits are sampled, the PSRN’s distribution must remain close to the algorithm’s ideal distribution by a distance of not more than 1/(<em>b</em><sup><em>m</em></sup>).</li>
</ol>
<blockquote>
<p><strong>Example</strong>: Suppose an algorithm samples from a normal distribution using base-2 uniform PSRNs. If it gives out a PSRN whose fractional part has three sampled bits (and no unsampled bits before the rightmost sampled bit), the PSRN’s distribution must be within a distance of—</p>
<ul>
<li>1/(2<sup>3</sup>) to the ideal normal distribution at the time the PSRN is output, and</li>
<li>1/(2<sup>6</sup>) to the ideal normal distribution if the caller later inserts three more uniform random bits to the end of the PSRN’s fractional part.</li>
</ul>
<p><strong>Notes:</strong></p>
<ol>
<li>It is not easy to turn a sampler for a nondiscrete distribution into an algorithm that meets these properties. Some reasons for this are given in the section “<a href="#Discussion"><strong>Discussion</strong></a>” later in this document.</li>
<li>The <em>exact rejection sampling</em> algorithm described by Oberhoff (2018)<sup id="fnref:13:1" role="doc-noteref"><a href="#fn:13" class="footnote" rel="footnote">13</a></sup> produces samples that act like PSRNs. However, in general, the algorithm doesn’t have the properties described in this section because some of its operations can introduce numerical error unless care is taken, and these operations include calculating minimums and maximums of probabilities. Moreover, the algorithm’s progression depends on the value of previously sampled bits, not just on the position of those bits as with the uniform and exponential distributions (see also (Thomas and Luk 2008)<sup id="fnref:4:1" role="doc-noteref"><a href="#fn:4" class="footnote" rel="footnote">4</a></sup>). For completeness, Oberhoff’s method appears in the appendix.</li>
</ol>
</blockquote>
<p><a id="Limitations"></a></p>
<h3 id="limitations">Limitations</h3>
<p>Because a PSRN stores a random variate in a certain interval, PSRNs are not well suited for representing numbers in zero-volume sets. Such sets include:</p>
<ul>
<li>Sets of integers or rational numbers.</li>
<li>Sets of individual points.</li>
<li>Curves on two- or higher-dimensional real number space.</li>
<li>Surfaces on three- or higher-dimensional real number space.</li>
</ul>
<p>In the case of curves and surfaces, a PSRN can’t directly store the coordinates, in space, of a random point on that curve or surface (because the exact value of those coordinates may be an irrational number that no computer can store, and no interval can bound those exact coordinates “tightly” enough), but the PSRN <em>can</em> store upper and lower bounds that indirectly give that point’s position on that curve or surface.</p>
<blockquote>
<p><strong>Examples:</strong></p>
<ol>
<li>To represent a point on the edge of a circle, a PSRN can store a random variate in the interval [0, 2*<em>π</em>), via the <strong>RandUniformFromReal</strong> method, given later, for 2*<em>π</em> (for example, it can store an integer part of 2 and a fractional part of [1, 3, 5] and thus represent a number in the interval [2.135, 2.136]), and the number stored this way indicates the distance on the circular arc relative to its starting position. A program that cares about the point’s-, x-, and y-coordinates can then generate enough digits of the PSRN to compute an approximation of cos(<em>P</em>) and sin(<em>P</em>), respectively, to the desired accuracy, where <em>P</em> is the number stored by the PSRN. (However, the direct use of mathematical functions such as <code>cos</code> and <code>sin</code> is outside the scope of this document.)</li>
<li>Example 1 is trivial, because each point on the interval maps evenly to a point on the circle. But this is not true in general: an interval’s or box’s points don’t map evenly to points on a curve or surface in general. For example, take two PSRNs describing the u- and v-coordinates of a 3 dimensional cone’s surface: [1.135, 1.136] for U and [0.288, 0.289] for V, and the cone’s-coordinates are X = U*cos(V), Y = U*sin(V), Z = U. In this example, the PSRNs form a box that’s mapped to a small part of the cone surface. However, the points in the box don’t map to the cone evenly this way, so generating enough digits to calculate X, Y, and Z to the desired accuracy will not sample uniformly from that part of the cone without more work (see Williamson (1987)<sup id="fnref:16" role="doc-noteref"><a href="#fn:16" class="footnote" rel="footnote">16</a></sup> for one solution).</li>
<li>For algorithms that sample uniformly on a general curve with arbitrary accuracy, see Chalkis et al. (2022)<sup id="fnref:17" role="doc-noteref"><a href="#fn:17" class="footnote" rel="footnote">17</a></sup></li>
</ol>
</blockquote>
<p><a id="Relation_to_Constructive_Reals"></a></p>
<h3 id="relation-to-constructive-reals">Relation to Constructive Reals</h3>
<p>Partially-sampled random numbers are related to a body of work dealing with so-called “constructive reals” or “recursive reals”, or operations on real numbers that compute an approximation of the exact result to a user-specified number of digit places. For example, in Hans-J. Boehm’s implementation (Boehm 2020)<sup id="fnref:18" role="doc-noteref"><a href="#fn:18" class="footnote" rel="footnote">18</a></sup>, (Boehm 1987)<sup id="fnref:19" role="doc-noteref"><a href="#fn:19" class="footnote" rel="footnote">19</a></sup>, each operation on “constructive reals” (such as addition, multiplication, <code>exp</code>, <code>ln</code>, and so on) is associated with a function <code>f(n)</code> (where <code>n</code> is usually 0 or greater) that returns an integer <code>m</code> such that <code>abs(m/pow(2, n) - x) < 1/pow(2, n)</code>, where <code>x</code> is the exact result of the operation. In addition, comparisons such as “less than” or “greater than” can operate on “constructive reals”. As suggested in Goubault-Larrecq et al. (2021)<sup id="fnref:20" role="doc-noteref"><a href="#fn:20" class="footnote" rel="footnote">20</a></sup>, there can also be an operation that samples the digits of a uniform random variate between 0 and 1 and gives access to approximations of that variate, sampling random digits as necessary. Similarly, operations of this kind can be defined to access approximations of the value stored in a PSRN (including a uniform or exponential PSRN), sampling digits for the PSRN as necessary.</p>
<p>Details on “constructive real” operations and comparisons are outside the scope of this document.</p>
<p><a id="Sampling_Uniform_and_Exponential_PSRNs"></a></p>
<h2 id="sampling-uniform-and-exponential-psrns">Sampling Uniform and Exponential PSRNs</h2>
<p> </p>
<p><a id="Sampling_Uniform_PSRNs"></a></p>
<h3 id="sampling-uniform-psrns">Sampling Uniform PSRNs</h3>
<p>There are several algorithms for sampling uniform partially-sampled random numbers given another number.</p>
<p>The <strong>RandUniform</strong> algorithm generates a uniformly distributed PSRN (<strong>a</strong>) that is greater than 0 and less than another PSRN (<strong>b</strong>) with probability 1. This algorithm samples digits of <strong>b</strong>’s fractional part as necessary. This algorithm should not be used if <strong>b</strong> is known to be a real number rather than a partially-sampled random number, since this algorithm could overshoot the value <strong>b</strong> had (or appeared to have) at the beginning of the algorithm; instead, the <strong>RandUniformFromReal</strong> algorithm, given later, should be used. (For example, if <strong>b</strong> is 3.425…, one possible result of this algorithm is <strong>a</strong> = 3.42574… and <strong>b</strong> = 3.42575… Note that in this example, 3.425… is not considered an exact number.)</p>
<ol>
<li>Create an empty uniform PSRN <strong>a</strong>. Let <em>β</em> be the base (or radix) of digits stored in <strong>b</strong>’s fractional part (for example, 2 for binary or 10 for decimal). If <strong>b</strong>’s integer part or sign is unsampled, or if <strong>b</strong>’s sign is negative, return an error.</li>
<li>(We now set <strong>a</strong>’s integer part and sign.) Set <strong>a</strong>’s sign to positive and <strong>a</strong>’s integer part to an integer chosen uniformly at random in [0, <em>bi</em>], where <em>bi</em> is <strong>b</strong>’s integer part (note that <em>bi</em> is included). If <strong>a</strong>’s integer part is less than <em>bi</em>, return <strong>a</strong>.</li>
<li>(We now sample <strong>a</strong>’s fractional part.) Set <em>i</em> to 0.</li>
<li>If <strong>b</strong>’s integer part is 0 and <strong>b</strong>’s fractional part begins with a sampled 0-digit, set <em>i</em> to the number of sampled zeros at the beginning of <strong>b</strong>’s fractional part. A nonzero digit or an unsampled digit ends this sequence. Then append <em>i</em> zeros to <strong>a</strong>’s fractional part. (For example, if <strong>b</strong> is 5.000302 or 4.000 or 0.0008, there are three sampled zeros that begin <strong>b</strong>’s fractional part, so <em>i</em> is set to 3 and three zeros are appended to <strong>a</strong>’s fractional part.)</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is unsampled, set the digit at that position to a base-_β_ digit chosen uniformly at random (if <em>β</em> is 2, for example, this is a number that is 1 or 0 with equal probability). (Positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc.)</li>
<li>If the digit at position <em>i</em> of <strong>b</strong>’s fractional part is unsampled, sample the digit at that position according to the kind of PSRN <strong>b</strong> is. (For example, if <strong>b</strong> is a uniform PSRN and <em>β</em> is 2, this can be done by setting the digit at that position to either 1 or 0 with equal probability.)</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is less than the corresponding digit for <strong>b</strong>, return <strong>a</strong>.</li>
<li>If that digit is greater, then discard <strong>a</strong>, then create a new empty uniform PSRN <strong>a</strong>, then go to step 2.</li>
<li>Add 1 to <em>i</em> and go to step 5.</li>
</ol>
<blockquote>
<p><strong>Notes:</strong></p>
<ol>
<li>Karney (2014, end of sec. 4)<sup id="fnref:1:5" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup> discusses how even the integer part can be partially sampled rather than generating the whole integer as in step 2 of the algorithm. However, incorporating this suggestion will add a nontrivial amount of complexity to the algorithm given earlier.</li>
<li>The <strong>RandUniform</strong> algorithm is equivalent to generating the product of a random variate (<strong>b</strong>) and a uniform random variate between 0 and 1.</li>
<li>If <strong>b</strong> is a uniform PSRN with a positive sign, an integer part of 0, and an empty fractional part, the <strong>RandUniform</strong> algorithm is equivalent to generating the product of two uniform random variates between 0 and 1.</li>
</ol>
</blockquote>
<p>The <strong>RandUniformInRangePositive</strong> algorithm generates a uniformly distributed PSRN (<strong>a</strong>) that is greater than one nonnegative real number <strong>bmin</strong> and less than another positive real number <strong>bmax</strong> with probability 1. This algorithm works whether <strong>bmin</strong> or <strong>bmax</strong> is known to be a rational number or not (for example, either number can be the result of an expression such as <code>exp(-2)</code> or <code>ln(20)</code>), but the algorithm notes how it can be more efficiently implemented if <strong>bmin</strong> or <strong>bmax</strong> is known to be a rational number.</p>
<ol>
<li>If <strong>bmin</strong> is greater than or equal to <strong>bmax</strong>, if <strong>bmin</strong> is less than 0, or if <strong>bmax</strong> is 0 or less, return an error.</li>
<li>Create an empty uniform PSRN <strong>a</strong>.</li>
<li>Special case: If <strong>bmax</strong> is 1 and <strong>bmin</strong> is 0, set <strong>a</strong>’s sign to positive, set <strong>a</strong>’s integer part to 0, and return <strong>a</strong>.</li>
<li>Special case: If <strong>bmax</strong> and <strong>bmin</strong> are rational numbers and each of their denominators is a power of <em>β</em>, including 1 (where <em>β</em> is the desired digit base, or radix, of the uniform PSRN, such as 10 for decimal or 2 for binary), then do the following:
<ol>
<li>Let <em>denom</em> be <strong>bmax</strong>’s or <strong>bmin</strong>’s denominator, whichever is greater.</li>
<li>Set <em>c1</em> to floor(<strong>bmax</strong>*<em>denom</em>) and <em>c2</em> to floor((<strong>bmax</strong>−<strong>bmin</strong>)*<em>denom</em>).</li>
<li>If <em>c2</em> is greater than 1, add to <em>c1</em> an integer chosen uniformly at random in [0, <em>c2</em>) (note that <em>c2</em> is excluded).</li>
<li>Let <em>d</em> be the base-_β_ logarithm of <em>denom</em> (this is equivalent to finding the minimum number of base-_β_ digits needed to store <em>denom</em> and subtracting 1). Transfer <em>c1</em>’s least significant digits to <strong>a</strong>’s fractional part; the variable <em>d</em> tells how many digits to transfer to each PSRN this way. Then set <strong>a</strong>’s sign to positive and <strong>a</strong>’s integer part to floor(<em>c1</em>/<em>β</em><sup><em>d</em></sup>). (For example, if <em>β</em> is 10, <em>d</em> is 3, and <em>c1</em> is 7342, set <strong>a</strong>’s fractional part to [3, 4, 2] and <strong>a</strong>’s integer part to 7.) Finally, return <strong>a</strong>.</li>
</ol>
</li>
<li>Calculate floor(<strong>bmax</strong>), and set <em>bmaxi</em> to the result. Likewise, calculate floor(<strong>bmin</strong>) and set <em>bmini</em> to the result.</li>
<li>If <em>bmini</em> is equal to <strong>bmin</strong> and <em>bmaxi</em> is equal to <strong>bmax</strong>, set <strong>a</strong>’s sign to positive and <strong>a</strong>’s integer part to an integer chosen uniformly at random in [<em>bmini</em>, <em>bmaxi</em>) (note that <em>bmaxi</em> is excluded), then return <strong>a</strong>. (It should be noted that determining whether a real number is equal to another is undecidable in general.)</li>
<li>(We now set <strong>a</strong>’s integer part and sign.) Set <strong>a</strong>’s sign to positive and <strong>a</strong>’s integer part to an integer chosen uniformly at random in the interval [<em>bmini</em>, <em>bmaxi</em>] (note that <em>bmaxi</em> is included). If <em>bmaxi</em> is equal to <strong>bmax</strong>, the integer is chosen from the interval [<em>bmini</em>, <em>bmaxi</em>−1] instead. Return <strong>a</strong> if—
<ul>
<li><strong>a</strong>’s integer part is greater than <em>bmini</em> and less than <em>bmaxi</em>, or</li>
<li><em>bmini</em> is equal to <strong>bmin</strong>, and <strong>a</strong>’s integer part is equal to <em>bmini</em> and less than <em>bmaxi</em>.</li>
</ul>
</li>
<li>(We now sample <strong>a</strong>’s fractional part.) Set <em>i</em> to 0 and <em>istart</em> to 0. ( Then, <em>if <strong>bmax</strong> is known rational:</em> set <em>bmaxf</em> to <strong>bmax</strong> minus <em>bmaxi</em>, and <em>if <strong>bmin</strong> is known rational</em>, set <em>bminf</em> to <strong>bmin</strong> minus <em>bmini</em>.)</li>
<li>(This step is not crucial for correctness, but helps improve its efficiency. It sets <strong>a</strong>’s fractional part to the initial digits shared by <strong>bmin</strong> and <strong>bmax</strong>.) If <strong>a</strong>’s integer part is equal to <em>bmini</em> and <em>bmaxi</em>, then do the following in a loop:
1. Calculate the base-_β_ digit at position <em>i</em> of <strong>bmax</strong>’s and <strong>bmin</strong>’s fractional parts, and set <em>dmax</em> and <em>dmin</em> to those digits, respectively. (<em>If <strong>bmax</strong> is known rational:</em> Do this step by setting <em>dmax</em> to floor(<em>bmaxf</em>*<em>β</em>) and <em>dmin</em> to floor(<em>bminf</em>*<em>β</em>).)
2. If <em>dmin</em> equals <em>dmax</em>, append <em>dmin</em> to <strong>a</strong>’s fractional part, then add 1 to <em>i</em> (and, if <strong>bmax</strong>, <strong>bmin</strong>, or both are known to be rational, set <em>bmaxf</em> to <em>bmaxf</em>*<em>β</em>−<em>d</em> and set <em>bminf</em> to <em>bminf</em>*<em>β</em>−<em>d</em>). Otherwise, break from this loop and set <em>istart</em> to <em>i</em>.</li>
<li>(Ensure the fractional part is greater than <strong>bmin</strong>’s.) Set <em>i</em> to <em>istart</em>, then if <strong>a</strong>’s integer part is equal to <em>bmini</em>:
<ol>
<li>Calculate the base-_β_ digit at position <em>i</em> of <strong>bmin</strong>’s fractional part, and set <em>dmin</em> to that digit.</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is unsampled, set the digit at that position to a base-_β_ digit chosen uniformly at random (for example, if <em>β</em> is 2, this is a number that is 1 or 0 with equal probability). (Positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc.)</li>
<li>Let <em>ad</em> be the digit at position <em>i</em> of <strong>a</strong>’s fractional part. If <em>ad</em> is greater than <em>dmin</em>, abort these substeps and go to step 11.</li>
<li>Discard <strong>a</strong>, create a new empty uniform PSRN <strong>a</strong>, and abort these substeps and go to step 7 if <em>ad</em> is less than <em>dmin</em>.</li>
<li>Add 1 to <em>i</em> and go to the first substep.</li>
</ol>
</li>
<li>(Ensure the fractional part is less than <strong>bmax</strong>’s.) Set <em>i</em> to <em>istart</em>, then if <strong>a</strong>’s integer part is equal to <em>bmaxi</em>:
<ol>
<li>If <em>bmaxi</em> is 0 and not equal to <strong>bmax</strong>, and if <strong>a</strong> has no digits in its fractional part, then do the following in a loop:
<ol>
<li>Calculate the base-_β_ digit at position <em>i</em> of <strong>bmax</strong>’s fractional part, and set <em>d</em> to that digit. (<em>If <strong>bmax</strong> is known rational:</em> Do this step by setting <em>d</em> to floor(<em>bmaxf</em>*<em>β</em>).)</li>
<li>If <em>d</em> is 0, append a 0-digit to <strong>a</strong>’s fractional part, then add 1 to <em>i</em> (and, if <strong>bmax</strong> is known to be rational, set <em>bmaxf</em> to <em>bmaxf</em>*<em>β</em>−<em>d</em>). Otherwise, break from this loop.</li>
</ol>
</li>
<li>Calculate the base-_β_ digit at position <em>i</em> of <strong>bmax</strong>’s fractional part, and set <em>dmax</em> to that digit. (<em>If <strong>bmax</strong> is known rational:</em> Do this step by multiplying <em>bmaxf</em> by <em>β</em>, then setting <em>dmax</em> to floor(<em>bmaxf</em>), then subtracting <em>dmax</em> from <em>bmaxf</em>.)</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is unsampled, set the digit at that position to a base-_β_ digit chosen uniformly at random.</li>
<li>Let <em>ad</em> be the digit at position <em>i</em> of <strong>a</strong>’s fractional part. Return <strong>a</strong> if <em>ad</em> is less than <em>dmax</em>.</li>
<li>Discard <strong>a</strong>, create a new empty uniform PSRN <strong>a</strong>, and abort these substeps and go to step 7 if—
<ul>
<li><em><strong>bmax</strong> is not known to be rational</em>, and either <em>ad</em> is greater than <em>dmax</em> or all the digits after the digit at position <em>i</em> of <strong>bmax</strong>’s fractional part are zeros, or</li>
<li><em><strong>bmax</strong> is known to be rational</em>, and either <em>ad</em> is greater than <em>dmax</em> or <em>bmaxf</em> is 0</li>
</ul>
</li>
<li>Add 1 to <em>i</em> and go to the second substep.</li>
</ol>
</li>
<li>Return <strong>a</strong>.</li>
</ol>
<p>The <strong>RandUniformInRange</strong> algorithm generates a uniformly distributed PSRN (<strong>a</strong>) that is greater than one real number <strong>bmin</strong> and less than another real number <strong>bmax</strong> with probability 1. It works for both positive and negative real numbers, but it’s specified separately from <strong>RandUniformInRangePositive</strong> to reduce clutter.</p>
<ol>
<li>If <strong>bmin</strong> is greater than or equal to <strong>bmax</strong>, return an error. If <strong>bmin</strong> and <strong>bmax</strong> are each 0 or greater, return the result of <strong>RandUniformInRangePositive</strong>.</li>
<li>If <strong>bmin</strong> and <strong>bmax</strong> are each 0 or less, call <strong>RandUniformInRangePositive</strong> with <strong>bmin</strong> = abs(<strong>bmax</strong>) and <strong>bmax</strong> = abs(<strong>bmin</strong>), set the result’s fractional part to negative, and return the result.</li>
<li>(At this point, <strong>bmin</strong> is less than 0 and <strong>bmax</strong> is greater than 0.) Set <em>bmaxi</em> to either floor(<strong>bmax</strong>) if <strong>bmax</strong> is 0 or greater, or −ceil(abs(<strong>bmax</strong>)) otherwise, and set <em>bmini</em> to either floor(<strong>bmin</strong>) if <strong>bmin</strong> is 0 or greater, or −ceil(abs(<strong>bmin</strong>)) otherwise. (Described this way to keep implementers from confusing floor with the integer part.)</li>
<li>Set <em>ipart</em> to an integer chosen uniformly at random in the interval [<em>bmini</em>, <em>bmaxi</em>] (note that <em>bmaxi</em> is included). If <em>bmaxi</em> is equal to <strong>bmax</strong>, the integer is chosen from the interval [<em>bmini</em>, <em>bmaxi</em>−1] instead.</li>
<li>If <em>ipart</em> is neither <em>bmini</em> nor <em>bmaxi</em>, create a uniform PSRN <strong>a</strong> with an empty fractional part; then set <strong>a</strong>’s sign to either positive if <em>ipart</em> is 0 or greater, or negative otherwise; then set <strong>a</strong>’s integer part to abs(<em>ipart</em>+1) if <em>ipart</em> is less than 0, or <em>ipart</em> otherwise; then return <strong>a</strong>.</li>
<li>If <em>ipart</em> is <em>bmini</em>, then create a uniform PSRN <strong>a</strong> with a positive sign, an integer part of abs(<em>ipart</em>+1), and an empty fractional part; then run <strong>URandLessThanReal</strong> with <strong>a</strong> = <strong>a</strong> and <strong>b</strong> = abs(<strong>bmin</strong>). If the result is 1, set <strong>a</strong>’s sign to negative and return <strong>a</strong>. Otherwise, go to step 3.</li>
<li>If <em>ipart</em> is <em>bmaxi</em>, then create a uniform PSRN <strong>a</strong> with a positive sign, an integer part of <em>ipart</em>, and an empty fractional part; then run <strong>URandLessThanReal</strong> with <strong>a</strong> = <strong>a</strong> and <strong>b</strong> = <strong>bmax</strong>. If the result is 1, return <strong>a</strong>. Otherwise, go to step 3.</li>
</ol>
<p>The <strong>RandUniformFromReal</strong> algorithm generates a uniformly distributed PSRN (<strong>a</strong>) that is greater than 0 and less than a real number <strong>b</strong> with probability 1. It is equivalent to the <strong>RandUniformInRangePositive</strong> algorithm with <strong>a</strong> = <strong>a</strong>, <strong>bmin</strong> = 0, and <strong>bmax</strong> = <strong>b</strong>.</p>
<p>The <strong>UniformComplement</strong> algorithm generates 1 minus the value of a uniform PSRN (<strong>a</strong>) as follows:</p>
<ol>
<li>If <strong>a</strong>’s sign is negative or its integer part is other than 0, return an error.</li>
<li>For each sampled digit in <strong>a</strong>’s fractional part, set it to <em>base</em>−1−<em>digit</em>, where <em>digit</em> is the digit and <em>base</em> is the base of digits stored by the PSRN, such as 2 for binary.</li>
<li>Return <strong>a</strong>.</li>
</ol>
<p><a id="Sampling_E_rands"></a></p>
<h3 id="sampling-e-rands">Sampling E-rands</h3>
<p><strong>Sampling an e-rand</strong> (a exponential PSRN) makes use of two observations, based on the parameter <em>λ</em> of the exponential distribution (with <em>λ</em> greater than 0):</p>
<ul>
<li>While a coin flip with probability of heads of exp(−<em>λ</em>) is heads, the exponential random variate is increased by 1.</li>
<li>If a coin flip with probability of heads of 1/(1+exp(<em>λ</em>/2<sup><em>prec</em></sup>)) is heads, the exponential random variate is increased by 2<sup>−<em>prec</em></sup>, where <em>prec</em> > 0 is an integer.</li>
</ul>
<p>Devroye and Gravel (2020, sec. 3.8)<sup id="fnref:3:2" role="doc-noteref"><a href="#fn:3" class="footnote" rel="footnote">3</a></sup> already made these observations, but only for <em>λ</em> = 1.</p>
<p>To implement these probabilities using just random bits, the sampler uses the <strong>ExpMinus</strong> and <strong>LogisticExp</strong> algorithms from “<a href="https://peteroupc.github.io/bernoulli.html"><strong>Bernoulli Factory Algorithms</strong></a>”.</p>
<blockquote>
<p><strong>Note:</strong> An exponential PSRN is an exponential random variate built up digit by digit, but an exponential random variate can also be stored in a <em>uniform PSRN</em> and generated by other algorithms.</p>
</blockquote>
<p><a id="Arithmetic_and_Comparisons_with_PSRNs"></a></p>
<h2 id="arithmetic-and-comparisons-with-psrns">Arithmetic and Comparisons with PSRNs</h2>
<p>This section describes addition, subtraction, multiplication, reciprocal, and division involving uniform PSRNs, and discusses other aspects of arithmetic involving PSRNs.</p>
<p><a id="Addition_and_Subtraction"></a></p>
<h3 id="addition-and-subtraction">Addition and Subtraction</h3>
<p>The following algorithm (<strong>UniformAdd</strong>) shows how to add two uniform PSRNs (<strong>a</strong> and <strong>b</strong>) that store digits of the same base (radix) in their fractional parts, and get a uniform PSRN as a result. The input PSRNs may have a positive or negative sign, and it is assumed that their integer parts and signs were sampled. <em>Python code implementing this algorithm is given later in this document.</em></p>
<ol>
<li>If <strong>a</strong> has unsampled digits before the last sampled digit in its fractional part, set each of those unsampled digits to a digit chosen uniformly at random. Do the same for <strong>b</strong>.</li>
<li>If <strong>a</strong> has fewer digits in its fractional part than <strong>b</strong> (or vice versa), sample enough digits (by setting them to uniform random digits, such as numbers that are each 1 or 0 with equal probability if <strong>a</strong> and <strong>b</strong> store binary, or base-2, digits) so that both PSRNs’ fractional parts have the same number of digits. Now, let <em>digitcount</em> be the number of digits in <strong>a</strong>’s fractional part.</li>
<li>Let <em>asign</em> be −1 if <strong>a</strong>’s sign is negative, or 1 otherwise. Let <em>bsign</em> be −1 if <strong>b</strong>’s sign is negative, or 1 otherwise. Let <em>afp</em> be <strong>a</strong>’s integer and fractional parts packed into an integer, as explained in the example, and let <em>bfp</em> be <strong>b</strong>’s integer and fractional parts packed the same way. (For example, if <strong>a</strong> represents the number 83.12344…, <em>afp</em> is 8312344.) Let <em>base</em> be the base of digits stored by <strong>a</strong> and <strong>b</strong>, such as 2 for binary or 10 for decimal.</li>
<li>Calculate the following four numbers:
<ul>
<li><em>afp</em>*<em>asign</em> + <em>bfp</em>*<em>bsign</em>.</li>
<li><em>afp</em>*<em>asign</em> + (<em>bfp</em>+1)*<em>bsign</em>.</li>
<li>(<em>afp</em>+1)*<em>asign</em> + <em>bfp</em>*<em>bsign</em>.</li>
<li>(<em>afp</em>+1)*<em>asign</em> + (<em>bfp</em>+1)*<em>bsign</em>.</li>
</ul>
</li>
<li>Set <em>minv</em> to the minimum and <em>maxv</em> to the maximum of the four numbers just calculated. These are lower and upper bounds to the result of applying interval addition to the PSRNs <strong>a</strong> and <strong>b</strong>. (For example, if <strong>a</strong> is 0.12344… and <strong>b</strong> is 0.38925…, their fractional parts are added to form <strong>c</strong> = 0.51269…., or the interval [0.51269, 0.51271].) However, the resulting PSRN is not uniformly distributed in its interval, and this is what the rest of this algorithm will solve, since in fact, the distribution of numbers in the interval resembles the distribution of the sum of two uniform random variates.</li>
<li>Once the four numbers are sorted from lowest to highest, let <em>midmin</em> be the second number in the sorted order, and let <em>midmax</em> be the third number in that order.</li>
<li>Set <em>x</em> to a uniform random integer in the interval [0, <em>maxv</em>−<em>minv</em>). If <em>x</em> < <em>midmin</em>−<em>minv</em>, set <em>dir</em> to 0 (the left side of the sum density). Otherwise, if <em>x</em> > <em>midmax</em>−<em>minv</em>, set <em>dir</em> to 1 (the right side of the sum density). Otherwise, do the following:
<ol>
<li>Set <em>s</em> to <em>minv</em> + <em>x</em>.</li>
<li>Create a new uniform PSRN, <em>ret</em>. If <em>s</em> is less than 0, set <em>s</em> to abs(1 + <em>s</em>) and set <em>ret</em>’s sign to negative. Otherwise, set <em>ret</em>’s sign to positive.</li>
<li>Transfer the <em>digitcount</em> least significant digits of <em>s</em> to <em>ret</em>’s fractional part. (Note that <em>ret</em>’s fractional part stores digits from most to least significant.) Then set <em>ret</em>’s integer part to floor(<em>s</em>/<em>base</em><sup><em>digitcount</em></sup>), then return <em>ret</em>. (For example, if <em>base</em> is 10, <em>digitcount</em> is 3, and <em>s</em> is 34297, then <em>ret</em>’s fractional part is set to [2, 9, 7], and <em>ret</em>’s integer part is set to 34.)</li>
</ol>
</li>
<li>If <em>dir</em> is 0 (the left side), set <em>pw</em> to <em>x</em> and <em>b</em> to <em>midmin</em>−<em>minv</em>. Otherwise (the right side), set <em>pw</em> to <em>x</em>−(<em>midmax</em>−<em>minv</em>) and <em>b</em> to <em>maxv</em>−<em>midmax</em>.</li>
<li>Set <em>newdigits</em> to 0, then set <em>y</em> to a uniform random integer in the interval [0, <em>b</em>).</li>
<li>If <em>dir</em> is 0, set <em>lower</em> to <em>pw</em>. Otherwise, set <em>lower</em> to <em>b</em>−1−<em>pw</em>.</li>
<li>If <em>y</em> is less than <em>lower</em>, do the following:
<ol>
<li>Set <em>s</em> to <em>minv</em> if <em>dir</em> is 0, or <em>midmax</em> otherwise, then set <em>s</em> to <em>s</em>*(<em>base</em><sup><em>newdigits</em></sup>) + <em>pw</em>.</li>
<li>Create a new uniform PSRN, <em>ret</em>. If <em>s</em> is less than 0, set <em>s</em> to abs(1 + <em>s</em>) and set <em>ret</em>’s sign to negative. Otherwise, set <em>ret</em>’s sign to positive.</li>
<li>Transfer the <em>digitcount</em> + <em>newdigits</em> least significant digits of <em>s</em> to <em>ret</em>’s fractional part, then set <em>ret</em>’s integer part to floor(<em>s</em>/<em>base</em><sup><em>digitcount</em> + <em>newdigits</em></sup>), then return <em>ret</em>.</li>
</ol>
</li>
<li>If <em>y</em> is greater than <em>lower</em> + 1, go to step 7. (This is a rejection event.)</li>
<li>Multiply <em>pw</em>, <em>y</em>, and <em>b</em> each by <em>base</em>, then add a digit chosen uniformly at random to <em>pw</em>, then add a digit chosen uniformly at random to <em>y</em>, then add 1 to <em>newdigits</em>, then go to step 10.</li>
</ol>
<p>The following algorithm (<strong>UniformAddRational</strong>) shows how to add a uniform PSRN (<strong>a</strong>) and a rational number <strong>b</strong>. The input PSRN may have a positive or negative sign, and it is assumed that its integer part and sign were sampled. Similarly, the rational number may be positive, negative, or zero. <em>Python code implementing this algorithm is given later in this document.</em></p>
<ol>
<li>Let <em>ai</em> be <strong>a</strong>’s integer part. Special cases:
<ul>
<li>If <strong>a</strong>’s sign is positive and has no sampled digits in its fractional part, and if <strong>b</strong> is an integer 0 or greater, return a uniform PSRN with a positive sign, an integer part equal to <em>ai</em> + <strong>b</strong>, and an empty fractional part.</li>
<li>If <strong>a</strong>’s sign is negative and has no sampled digits in its fractional part, and if <strong>b</strong> is an integer less than 0, return a uniform PSRN with a negative sign, an integer part equal to <em>ai</em> + abs(<strong>b</strong>), and an empty fractional part.</li>
<li>If <strong>a</strong>’s sign is positive, has an integer part of 0, and has no sampled digits in its fractional part, and if <strong>b</strong> is an integer, return a uniform PSRN with an empty fractional part. If <strong>b</strong> is less than 0, the PSRN’s sign is negative and its integer part is abs(<strong>b</strong>)−1. If <strong>b</strong> is 0 or greater, the PSRN’s sign is positive and its integer part is abs(<strong>b</strong>).</li>
<li>If <strong>b</strong> is 0, return a copy of <strong>a</strong>.</li>
</ul>
</li>
<li>If <strong>a</strong> has unsampled digits before the last sampled digit in its fractional part, set each of those unsampled digits to a digit chosen uniformly at random. Now, let <em>digitcount</em> be the number of digits in <strong>a</strong>’s fractional part.</li>
<li>Let <em>asign</em> be −1 if <strong>a</strong>’s sign is negative or 1 otherwise. Let <em>base</em> be the base of digits stored in <strong>a</strong>’s fractional part (such as 2 for binary or 10 for decimal). Set <em>absfrac</em> to abs(<strong>b</strong>), then set <em>fraction</em> to <em>absfrac</em> − floor(<em>absfrac</em>).</li>
<li>Let <em>afp</em> be <strong>a</strong>’s integer and fractional parts packed into an integer, as explained in the example. (For example, if <strong>a</strong> represents the number 83.12344…, <em>afp</em> is 8312344.) Let <em>asign</em> be −1 if</li>
<li>Set <em>ddc</em> to <em>base</em><sup><em>dcount</em></sup>, then set <em>lower</em> to ((<em>afp</em>*<em>asign</em>)/<em>ddc</em>)+<strong>b</strong> (using rational arithmetic), then set <em>upper</em> to (((<em>afp</em>+1)*<em>asign</em>)/<em>ddc</em>)+<strong>b</strong> (again using rational arithmetic). Set <em>minv</em> to min(<em>lower</em>, <em>upper</em>), and set <em>maxv</em> to min(<em>lower</em>, <em>upper</em>).</li>
<li>Set <em>newdigits</em> to 0, then set <em>b</em> to 1, then set <em>ddc</em> to <em>base</em><sup><em>dcount</em></sup>, then set <em>mind</em> to floor(abs(<em>minv</em>*<em>ddc</em>)), then set <em>maxd</em> to floor(abs(<em>maxv</em>*<em>ddc</em>)). (Outer bounds): Then set <em>rvstart</em> to <em>mind</em>−1 if <em>minv</em> is less than 0, or <em>mind</em> otherwise, then set <em>rvend</em> to <em>maxd</em> if <em>maxv</em> is less than 0, or <em>maxd</em>+1 otherwise.</li>
<li>Set <em>rv</em> to a uniform random integer in the interval [0, <em>rvend</em>−<em>rvstart</em>), then set <em>rvs</em> to <em>rv</em> + <em>rvstart</em>.</li>
<li>(Inner bounds.) Set <em>innerstart</em> to <em>mind</em> if <em>minv</em> is less than 0, or <em>mind</em>+1 otherwise, then set <em>innerend</em> to <em>maxd</em>−1 if <em>maxv</em> is less than 0, or <em>maxd</em> otherwise.</li>
<li>If <em>rvs</em> is greater than <em>innerstart</em> and less than <em>innerend</em>, then the algorithm is almost done, so do the following:
<ol>
<li>Create an empty uniform PSRN, call it <em>ret</em>. If <em>rvs</em> is less than 0, set <em>rvs</em> to abs(1 + <em>rvs</em>) and set <em>ret</em>’s sign to negative. Otherwise, set <em>ret</em>’s sign to positive.</li>
<li>Transfer the <em>digitcount</em> + <em>newdigits</em> least significant digits of <em>rvs</em> to <em>ret</em>’s fractional part, then set <em>ret</em>’s integer part to floor(<em>rvs</em>/<em>base</em><sup><em>digitcount</em> + <em>newdigits</em></sup>), then return <em>ret</em>.</li>
</ol>
</li>
<li>If <em>rvs</em> is equal to or less than <em>innerstart</em> and (<em>rvs</em>+1)/<em>ddc</em> (calculated using rational arithmetic) is less than or equal to <em>minv</em>, go to step 6. (This is a rejection event.)</li>
<li>If <em>rvs</em>/<em>ddc</em> (calculated using rational arithmetic) is greater than or equal to <em>maxv</em>, go to step 6. (This is a rejection event.)</li>
<li>Add 1 to <em>newdigits</em>, then multiply <em>ddc</em>, <em>rvstart</em>, <em>rv</em>, and <em>rvend</em> each by <em>base</em>, then set <em>mind</em> to floor(abs(<em>minv</em>*<em>ddc</em>)), then set <em>maxd</em> to floor(abs(<em>maxv</em>*<em>ddc</em>)), then add a digit chosen uniformly at random to <em>rv</em>, then set <em>rvs</em> to <em>rv</em>+<em>rvstart</em>, then go to step 8.</li>
</ol>
<p><a id="Multiplication"></a></p>
<h3 id="multiplication">Multiplication</h3>
<p>The following algorithm (<strong>UniformMultiplyRational</strong>) shows how to multiply a uniform PSRN (<strong>a</strong>) by a nonzero rational number <strong>b</strong>. The input PSRN may have a positive or negative sign, and it is assumed that its integer part and sign were sampled. <em>Python code implementing this algorithm is given later in this document.</em></p>
<ol>
<li>If <strong>a</strong> has unsampled digits before the last sampled digit in its fractional part, set each of those unsampled digits to a digit chosen uniformly at random. Now, let <em>digitcount</em> be the number of digits in <strong>a</strong>’s fractional part.</li>
<li>Create a uniform PSRN, call it <em>ret</em>. Set <em>ret</em>’s sign to be −1 if <strong>a</strong>’s sign is positive and <strong>b</strong> is less than 0 or if <strong>a</strong>’s sign is negative and <strong>b</strong> is 0 or greater, or 1 otherwise, then set <em>ret</em>’s integer part to 0. Let <em>base</em> be the base of digits stored in <strong>a</strong>’s fractional part (such as 2 for binary or 10 for decimal). Set <em>absfrac</em> to abs(<strong>b</strong>), then set <em>fraction</em> to <em>absfrac</em> − floor(<em>absfrac</em>).</li>
<li>Let <em>afp</em> be <strong>a</strong>’s integer and fractional parts packed into an integer, as explained in the example. (For example, if <strong>a</strong> represents the number 83.12344…, <em>afp</em> is 8312344.)</li>
<li>Set <em>dcount</em> to <em>digitcount</em>, then set <em>ddc</em> to <em>base</em><sup><em>dcount</em></sup>, then set <em>lower</em> to (<em>afp</em>/<em>ddc</em>)*<em>absfrac</em> (using rational arithmetic), then set <em>upper</em> to ((<em>afp</em>+1)/<em>ddc</em>)*<em>absfrac</em> (again using rational arithmetic).</li>
<li>Set <em>rv</em> to a uniform random integer in the interval [floor(<em>lower</em>*<em>ddc</em>), floor(<em>upper</em>*<em>ddc</em>)).</li>
<li>Set <em>rvlower</em> to <em>rv</em>/<em>ddc</em> (as a rational number), then set <em>rvupper</em> to (<em>rv</em>+1)/<em>ddc</em> (as a rational number).</li>
<li>If <em>rvlower</em> is greater than or equal to <em>lower</em> and <em>rvupper</em> is less than <em>upper</em>, then the algorithm is almost done, so do the following: Transfer the <em>dcount</em> least significant digits of <em>rv</em> to <em>ret</em>’s fractional part (note that <em>ret</em>’s fractional part stores digits from most to least significant), then set <em>ret</em>’s integer part to floor(<em>rv</em>/<em>base</em><sup><em>dcount</em></sup>), then return <em>ret</em>. (For example, if <em>base</em> is 10, <em>dcount</em> is 4, and <em>rv</em> is 342978, then <em>ret</em>’s fractional part is set to [2, 9, 7, 8], and <em>ret</em>’s integer part is set to 34.)</li>
<li>If <em>rvlower</em> is greater than <em>upper</em> or if <em>rvupper</em> is less than <em>lower</em>, go to step 4.</li>
<li>Multiply <em>rv</em> and <em>ddc</em> each by <em>base</em>, then add 1 to <em>dcount</em>, then add a digit chosen uniformly at random to <em>rv</em>, then go to step 6.</li>
</ol>
<p>Another algorithm (<strong>UniformMultiply</strong>) shows how to multiply two uniform PSRNs (<strong>a</strong> and <strong>b</strong>) is given in the appendix — the algorithm is complicated and it may be simpler to instead connect PSRNs with “constructive reals” (described earlier) that implement multiplication to arbitrary precision.</p>
<blockquote>
<p><strong>Note:</strong> Let <em>b</em>>0, <em>c</em>≥0, and <em>d</em>>0 be rational numbers where <em>d</em>><em>c</em>. To generate the product of two uniform variates, one in [0, <em>b</em>] and the other in [<em>c</em>, <em>d</em>], the following algorithm can be used.<br />(1) Generate a uniform PSRN using <strong>RandUniformFromReal</strong> with parameter <em>b</em>*(<em>d</em>−<em>c</em>), call it <strong>K</strong>;<br />(2) Get the result of <strong>UniformAddRational</strong> with parameters <strong>K</strong> and <em>b</em>*<em>c</em>, call it <strong>M</strong>;<br />(3) Generate a uniform PSRN using <strong>RandUniform</strong> with parameter <strong>M</strong>; return the PSRN.<br />Broadly speaking: “generate a uniform(0, <em>b</em>*(<em>d</em>−<em>c</em>)) random variate <em>X</em>, then return a uniform(0, <em>X</em>+<em>b</em>*<em>c</em>) random variate”. See the <a href="#Uniform_of_Uniforms_Produces_a_Product_of_Uniforms"><strong>appendix</strong></a> for evidence that this algorithm works, at least when <em>c</em> = 0.</p>
</blockquote>
<p><a id="Reciprocal_and_Division"></a></p>
<h3 id="reciprocal-and-division">Reciprocal and Division</h3>
<p>The following algorithm (<strong>UniformReciprocal</strong>) generates 1/<strong>a</strong>, where <strong>a</strong> is a uniform PSRN, and generates a new uniform PSRN with that reciprocal. The input PSRN may have a positive or negative sign, and it is assumed that its integer part and sign were sampled. All divisions mentioned here should be done using rational arithmetic. <em>Python code implementing this algorithm is given later in this document.</em></p>
<ol>
<li>If <strong>a</strong> has unsampled digits before the last sampled digit in its fractional part, set each of those unsampled digits to a digit chosen uniformly at random. Now, let <em>digitcount</em> be the number of digits in <strong>a</strong>’s fractional part.</li>
<li>Create a uniform PSRN, call it <em>ret</em>. Set <em>ret</em>’s sign to <strong>a</strong>’s sign. Let <em>base</em> be the base of digits stored in <strong>a</strong>’s fractional part (such as 2 for binary or 10 for decimal).</li>
<li>If <strong>a</strong> has no nonzero digit in its fractional part, and has an integer part of 0, then append a digit chosen uniformly at random to <strong>a</strong>’s fractional part. If that digit is 0, repeat this step. (This step is crucial for correctness when both PSRNs’ intervals cover the number 0, since the distribution of their product is different from the usual case.)</li>
<li>Let <em>afp</em> be <strong>a</strong>’s integer and fractional parts packed into an integer, as explained in the example. (For example, if <strong>a</strong> represents the number 83.12344…, <em>afp</em> is 8312344.)</li>
<li>(Calculate lower and upper bounds of 1/<strong>a</strong>, disregarding <strong>a</strong>’s sign.) Set <em>dcount</em> to <em>digitcount</em>, then set <em>ddc</em> to <em>base</em><sup><em>dcount</em></sup>, then set <em>lower</em> to (<em>ddc</em>/(<em>afp</em>+1)), then set <em>upper</em> to (<em>ddc</em>/<em>afp</em>).</li>
<li>Set <em>lowerdc</em> to floor(<em>lower</em>*<em>ddc</em>). If <em>lowerdc</em> is 0, add 1 to <em>dcount</em>, multiply <em>ddc</em> by <em>base</em>, then repeat this step. (This step too is important for correctness.)</li>
<li>(<em>rv</em> represents a tight interval between the lower and upper bounds, and <em>rv2</em> represents a uniform random variate between 0 and 1 to compare with the density function for the reciprocal.) Set <em>rv</em> to a uniform random integer in the interval [<em>lowerdc</em>, floor(<em>upper</em>*<em>ddc</em>)). Set <em>rv2</em> to a uniform random integer in the interval [0, <em>lowerdc</em>).</li>
<li>(Get the bounds of the tight interval <em>rv</em>.) Set <em>rvlower</em> to <em>rv</em>/<em>ddc</em>, then set <em>rvupper</em> to (<em>rv</em>+1)/<em>ddc</em>.</li>
<li>If <em>rvlower</em> is greater than or equal to <em>lower</em> and <em>rvupper</em> is less than <em>upper</em>:
<ol>
<li>Set <em>rvd</em> to <em>lowerdc</em>/<em>ddc</em>, then set <em>rvlower2</em> to <em>rv2</em>/<em>lowerdc</em>, then set <em>rvupper2</em> to (<em>rv2</em>+1)/<em>lowerdc</em>. (<em>rvlower2</em> and <em>rvupper2</em> are bounds of the uniform(0, 1) variate.)</li>
<li>(Compare with upper bounded density: <em>y</em><sup>2</sup>/<em>x</em><sup>2</sup>, where <em>y</em> is the lower bound of 1/<strong>a</strong> and <em>x</em> is between the lower and upper bounds.) If <em>rvupper2</em> is less than (<em>rvd</em>*<em>rvd</em>)/(<em>rvupper</em>*<em>rvupper</em>), then the algorithm is almost done, so do the following: Transfer the <em>dcount</em> least significant digits of <em>rv</em> to <em>ret</em>’s fractional part (note that <em>ret</em>’s fractional part stores digits from most to least significant), then set <em>ret</em>’s integer part to floor(<em>rv</em>/<em>base</em><sup><em>dcount</em></sup>), then return <em>ret</em>. (For example, if <em>base</em> is 10, <em>dcount</em> is 4, and <em>rv</em> is 342978, then <em>ret</em>’s fractional part is set to [2, 9, 7, 8], and <em>ret</em>’s integer part is set to 34.)</li>
<li>(Compare with lower bounded density.) If <em>rvlower2</em> is greater than (<em>rvd</em>*<em>rvd</em>)/(<em>rvlower</em>*<em>rvlower</em>), then abort these substeps and go to step 5. (This is a rejection event.)</li>
</ol>
</li>
<li>If <em>rvlower</em> is greater than <em>upper</em> or if <em>rvupper</em> is less than <em>lower</em>, go to step 5. (This is a rejection event.)</li>
<li>Multiply <em>rv</em>, <em>rv2</em>, <em>lowerdc</em>, and <em>ddc</em> each by <em>base</em>, then add 1 to <em>dcount</em>, then add a digit chosen uniformly at random to <em>rv</em>, then add a digit chosen uniformly at random to <em>rv2</em>, then go to step 8.</li>
</ol>
<p>With this algorithm it’s now trivial to describe an algorithm for dividing one uniform PSRN <strong>a</strong> by another uniform PSRN <strong>b</strong>, here called <strong>UniformDivide</strong>:</p>
<ol>
<li>Run the <strong>UniformReciprocal</strong> algorithm on <strong>b</strong> to create a new uniform PSRN <strong>c</strong>.</li>
<li>If <strong>c</strong>’s fractional part has no digits or all of them are zeros, append uniform random digits to <strong>c</strong> until a nonzero digit is appended this way.</li>
<li>Run the <strong>UniformMultiply</strong> algorithm (given in the appendix) on <strong>a</strong> and <strong>b</strong>, in that order, and return the result of that algorithm.</li>
</ol>
<p>It’s likewise trivial to describe an algorithm for multiplying a uniform PSRN <strong>a</strong> by a nonzero rational number <strong>b</strong>, here called <strong>UniformDivideRational</strong>:</p>
<ol>
<li>If <strong>b</strong> is 0, return an error.</li>
<li>Run the <strong>UniformMultiplyRational</strong> algorithm on <strong>a</strong> and 1/<strong>b</strong>, in that order, and return the result of that algorithm.</li>
</ol>
<p><a id="Using_the_Arithmetic_Algorithms"></a></p>
<h3 id="using-the-arithmetic-algorithms">Using the Arithmetic Algorithms</h3>
<p>The algorithms given earlier for addition and multiplication are useful for scaling and shifting PSRNs. For example, they can transform a normally-distributed PSRN into one with an arbitrary mean and standard deviation (by first multiplying the PSRN by the standard deviation, then adding the mean). Here is a sketch of a procedure that achieves this, given two parameters, <em>location</em> and <em>scale</em>, that are both rational numbers.</p>
<ol>
<li>Generate a uniform PSRN, then transform it into a variate of the desired distribution via an algorithm that employs rejection from the uniform distribution (such as Karney’s algorithm for the standard normal distribution (Karney 2016)<sup id="fnref:1:6" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>)). This procedure won’t work for exponential PSRNs (e-rands).</li>
<li>Run the <strong>UniformMultiplyRational</strong> algorithm to multiply the uniform PSRN by the rational parameter <em>scale</em> to get a new uniform PSRN.</li>
<li>Run the <strong>UniformAddRational</strong> algorithm to add the new uniform PSRN and the rational parameter <em>location</em> to get a third uniform PSRN. Return this third PSRN.</li>
</ol>
<p>See also the section “Discussion” later in this article.</p>
<p><a id="Comparisons"></a></p>
<h3 id="comparisons">Comparisons</h3>
<p>Two PSRNs, each of a different distribution but storing digits of the same base (radix), can be exactly compared to each other using algorithms similar to those in this section.</p>
<p>The <strong>RandLess</strong> algorithm compares two PSRNs, <strong>a</strong> and <strong>b</strong> (and samples additional bits from them as necessary) and returns 1 if <strong>a</strong> turns out to be less than <strong>b</strong> with probability 1, or 0 otherwise (see also (Karney 2016)<sup id="fnref:1:7" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>)).</p>
<ol>
<li>If <strong>a</strong>’s integer part wasn’t sampled yet, sample <strong>a</strong>’s integer part according to the kind of PSRN <strong>a</strong> is. Do the same for <strong>b</strong>.</li>
<li>If <strong>a</strong>’s sign is different from <strong>b</strong>’s sign, return 1 if <strong>a</strong>’s sign is negative and 0 if it’s positive. If <strong>a</strong>’s sign is positive, return 1 if <strong>a</strong>’s integer part is less than <strong>b</strong>’s, or 0 if greater. If <strong>a</strong>’s sign is negative, return 0 if <strong>a</strong>’s integer part is less than <strong>b</strong>’s, or 1 if greater.</li>
<li>Set <em>i</em> to 0.</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is unsampled, set the digit at that position according to the kind of PSRN <strong>a</strong> is. (Positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc.) Do the same for <strong>b</strong>.</li>
<li>Let <em>da</em> be the digit at position <em>i</em> of <strong>a</strong>’s fractional part, and let <em>db</em> be <strong>b</strong>’s corresponding digit.</li>
<li>If <strong>a</strong>’s sign is positive, return 1 if <em>da</em> is less than <em>db</em>, or 0 if <em>da</em> is greater than <em>db</em>.</li>
<li>If <strong>a</strong>’s sign is negative, return 0 if <em>da</em> is less than <em>db</em>, or 1 if <em>da</em> is greater than <em>db</em>.</li>
<li>Add 1 to <em>i</em> and go to step 4.</li>
</ol>
<p><strong>URandLess</strong> is a version of <strong>RandLess</strong> that involves two uniform PSRNs. The algorithm for <strong>URandLess</strong> samples digit <em>i</em> in step 4 by setting the digit at position <em>i</em> to a digit chosen uniformly at random. (For example, if <strong>a</strong> is a uniform PSRN that stores base-2 or binary digits, this can be done by setting the digit at that position to either 1 or 0 with equal probability.)</p>
<blockquote>
<p><strong>Note</strong>: To sample the <strong>maximum</strong> of two uniform random variates between 0 and 1, or the <strong>square root</strong> of a uniform random variate between 0 and 1: (1) Generate two uniform PSRNs <strong>a</strong> and <strong>b</strong> each with a positive sign, an integer part of 0, and an empty fractional part. (2) Run <strong>RandLess</strong> on <strong>a</strong> and <strong>b</strong> in that order. If the call returns 0, return <strong>a</strong>; otherwise, return <strong>b</strong>.</p>
</blockquote>
<p>The <strong>RandLessThanReal</strong> algorithm compares a PSRN <strong>a</strong> with a real number <strong>b</strong> and returns 1 if <strong>a</strong> turns out to be less than <strong>b</strong> with probability 1, or 0 otherwise. This algorithm samples digits of <strong>a</strong>’s fractional part as necessary. This algorithm works whether <strong>b</strong> is known to be a rational number or not (for example, <strong>b</strong> can be the result of an expression such as <code>exp(-2)</code> or <code>ln(20)</code>), but the algorithm notes how it can be more efficiently implemented if <strong>b</strong> is known to be a rational number.</p>
<ol>
<li>If <strong>a</strong>’s integer part or sign is unsampled, return an error.</li>
<li>Set <em>bs</em> to −1 if <strong>b</strong> is less than 0, or 1 otherwise. Calculate floor(abs(<strong>b</strong>)), and set <em>bi</em> to the result. (<em>If <strong>b</strong> is known rational:</em> Then set <em>bf</em> to abs(<strong>b</strong>) minus <em>bi</em>.)</li>
<li>If <strong>a</strong>’s sign is different from <em>bs</em>’s sign, return 1 if <strong>a</strong>’s sign is negative and 0 if it’s positive. If <strong>a</strong>’s sign is positive, return 1 if <strong>a</strong>’s integer part is less than <em>bi</em>, or 0 if greater. (Continue if both are equal.) If <strong>a</strong>’s sign is negative, return 0 if <strong>a</strong>’s integer part is less than <em>bi</em>, or 1 if greater. (Continue if both are equal.)</li>
<li>Set <em>i</em> to 0.</li>
<li>If the digit at position <em>i</em> of <strong>a</strong>’s fractional part is unsampled, set the digit at that position according to the kind of PSRN <strong>a</strong> is. (Positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc.)</li>
<li>Calculate the base-_β_ digit at position <em>i</em> of <strong>b</strong>’s fractional part, and set <em>d</em> to that digit. (<em>If <strong>b</strong> is known rational:</em> Do this step by multiplying <em>bf</em> by <em>β</em>, then setting <em>d</em> to floor(<em>bf</em>), then subtracting <em>d</em> from <em>bf</em>.)</li>
<li>Let <em>ad</em> be the digit at position <em>i</em> of <strong>a</strong>’s fractional part.</li>
<li>Return 1 if—
<ul>
<li><em>ad</em> is less than <em>d</em> and <strong>a</strong>’s sign is positive,</li>
<li><em>ad</em> is greater than <em>d</em> and <strong>a</strong>’s sign is negative, or</li>
<li><em>ad</em> is equal to <em>d</em>, <strong>a</strong>’s sign is negative, and—
<ul>
<li><em><strong>b</strong> is not known to be rational</em> and all the digits after the digit at position <em>i</em> of <strong>b</strong>’s fractional part are zeros (indicating <strong>a</strong> is less than <strong>b</strong> with probability 1), or</li>
<li><em><strong>b</strong> is known to be rational</em> and <em>bf</em> is 0 (indicating <strong>a</strong> is less than <strong>b</strong> with probability 1).</li>
</ul>
</li>
</ul>
</li>
<li>Return 0 if—
<ul>
<li><em>ad</em> is less than <em>d</em> and <strong>a</strong>’s sign is negative,</li>
<li><em>ad</em> is greater than <em>d</em> and <strong>a</strong>’s sign is positive, or</li>
<li><em>ad</em> is equal to <em>d</em>, <strong>a</strong>’s sign is positive, and—
<ul>
<li><em><strong>b</strong> is not known to be rational</em> and all the digits after the digit at position <em>i</em> of <strong>b</strong>’s fractional part are zeros (indicating <strong>a</strong> is greater than <strong>b</strong> with probability 1), or</li>
<li><em><strong>b</strong> is known to be rational</em> and <em>bf</em> is 0 (indicating <strong>a</strong> is greater than <strong>b</strong> with probability 1).</li>
</ul>
</li>
</ul>
</li>
<li>Add 1 to <em>i</em> and go to step 5.</li>
</ol>
<p>An alternative version of steps 6 through 9 in the algorithm above are as follows (see also (Brassard et al. 2019)<sup id="fnref:21" role="doc-noteref"><a href="#fn:21" class="footnote" rel="footnote">21</a></sup>):</p>
<ul>
<li>(6.) Calculate <em>bp</em>, which is an approximation to <strong>b</strong> such that abs(<strong>b</strong> − <em>bp</em>) <= <em>β</em><sup>−<em>i</em> − 1</sup>, and such that <em>bp</em> has the same sign as <strong>b</strong>. Let <em>bk</em> be <em>bp</em>’s digit expansion up to the <em>i</em> + 1 digits after the point (ignoring its sign). For example, if <strong>b</strong> is π or −π, <em>β</em> is 10, and <em>i</em> is 4, one possibility is <em>bp</em> = 3.14159 and <em>bk</em> = 314159.</li>
<li>(7.) Let <em>ak</em> be <strong>a</strong>’s digit expansion up to the <em>i</em> + 1 digits after the point (ignoring its sign).</li>
<li>(8.) If <em>ak</em> <= <em>bk</em> − 2, return either 1 if <strong>a</strong>’s sign is positive or 0 otherwise.</li>
<li>(9.) If <em>ak</em> >= <em>bk</em> + 1, return either 1 if <strong>a</strong>’s sign is negative or 0 otherwise.<sup id="fnref:22" role="doc-noteref"><a href="#fn:22" class="footnote" rel="footnote">22</a></sup></li>
</ul>
<p><strong>URandLessThanReal</strong> is a version of <strong>RandLessThanReal</strong> in which <strong>a</strong> is a uniform PSRN. The algorithm for <strong>URandLessThanReal</strong> samples digit <em>i</em> in step 4 by setting the digit at position <em>i</em> to a digit chosen uniformly at random.</p>
<p>The following shows how to implement <strong>URandLessThanReal</strong> when <strong>b</strong> is a fraction known by its numerator and denominator, <em>num</em>/<em>den</em>.</p>
<ol>
<li>If <strong>a</strong>’s integer part or sign is unsampled, or if <em>den</em> is 0, return an error. Then, if <em>num</em> and <em>den</em> are both less than 0, set them to their absolute values. Then if <strong>a</strong>’s sign is positive, its integer part is 0, and <em>num</em> is 0, return 0. Then if <strong>a</strong>’s sign is positive, its integer part is 0, and <em>num</em>’s sign is different from <em>den</em>’s sign, return 0.</li>
<li>Set <em>bs</em> to −1 if <em>num</em> or <em>den</em>, but not both, is less than 0, or 1 otherwise, then set <em>den</em> to abs(<em>den</em>), then set <em>bi</em> to floor(abs(<em>num</em>)/<em>den</em>), then set <em>num</em> to rem(abs(<em>num</em>), <em>den</em>).</li>
<li>If <strong>a</strong>’s sign is different from <em>bs</em>’s sign, return 1 if <strong>a</strong>’s sign is negative and 0 if it’s positive. If <strong>a</strong>’s sign is positive, return 1 if <strong>a</strong>’s integer part is less than <em>bi</em>, or 0 if greater. (Continue if both are equal.) If <strong>a</strong>’s sign is negative, return 0 if <strong>a</strong>’s integer part is less than <em>bi</em>, or 1 if greater. (Continue if both are equal.) If <em>num</em> is 0 (indicating the fraction is an integer), return 0 if <strong>a</strong>’s sign is positive and 1 otherwise.</li>
<li>Set <em>pt</em> to <em>base</em>, and set <em>i</em> to 0. (<em>base</em> is the base, or radix, of <strong>a</strong>’s digits, such as 2 for binary or 10 for decimal.)</li>
<li>Set <em>d1</em> to the digit at the <em>i</em><sup>th</sup> position (starting from 0) of <strong>a</strong>’s fractional part. If the digit at that position is unsampled, put a digit chosen uniformly at random at that position and set <em>d1</em> to that digit.</li>
<li>Set <em>c</em> to 1 if <em>num</em> * <em>pt</em> >= <em>den</em>, and 0 otherwise.</li>
<li>Set <em>d2</em> to floor(<em>num</em> * <em>pt</em> / <em>den</em>). (In base 2, this is equivalent to setting <em>d2</em> to <em>c</em>.)</li>
<li>If <em>d1</em> is less than <em>d2</em>, return either 1 if <strong>a</strong>’s sign is positive, or 0 otherwise. If <em>d1</em> is greater than <em>d2</em>, return either 0 if <strong>a</strong>’s sign is positive, or 1 otherwise.</li>
<li>If <em>c</em> is 1, set <em>num</em> to <em>num</em> * <em>pt</em> − <em>den</em> * <em>d2</em>, then multiply <em>den</em> by <em>pt</em>.</li>
<li>If <em>num</em> is 0, return either 0 if <strong>a</strong>’s sign is positive, or 1 otherwise.</li>
<li>Multiply <em>pt</em> by <em>base</em>, add 1 to <em>i</em>, and go to step 5.</li>
</ol>
<p><a id="Discussion"></a></p>
<h3 id="discussion">Discussion</h3>
<p>This section discusses issues involving arithmetic with PSRNs.</p>
<p><strong>Uniform PSRN arithmetic produces nonuniform distributions in general.</strong> As can be seen in the arithmetic algorithms earlier in this section (such as <strong>UniformAdd</strong> and <strong>UniformMultiplyRational</strong>), addition, multiplication, and other arithmetic operations with PSRNs (see also (Brassard et al., 2019)<sup id="fnref:21:1" role="doc-noteref"><a href="#fn:21" class="footnote" rel="footnote">21</a></sup>) are not as trivial as adding, multiplying, etc. their integer and fractional parts. A uniform PSRN is ultimately a uniform random variate inside an interval (this is its nature), yet arithmetic on random variates does not produce a uniform distribution in general.</p>
<p>An example illustrates this. Say we have two uniform PSRNs: <em>A</em> = 0.12345… and <em>B</em> = 0.38901…. They represent random variates in the intervals <em>AI</em> = [0.12345, 0.12346] and <em>BI</em> = [0.38901, 0.38902], respectively. Adding two uniform PSRNs is akin to adding their intervals (using interval arithmetic), so that in this example, the result <em>C</em> lies in <em>CI</em> = [0.12345 + 0.38901, 0.12346 + 0.38902] = [0.51246, 0.51248]. However, the resulting random variate is <em>not</em> uniformly distributed in [0.51246, 0.51248], so that simply choosing a uniform random variate in the interval won’t work. (This is true in general for other arithmetic operations besides addition.) This can be demonstrated by generating many pairs of uniform random variates in the intervals <em>AI</em> and <em>BI</em>, summing the numbers in each pair, and building a histogram using the sums (which will all lie in the interval <em>CI</em>). In this case, the histogram will show a triangular distribution that peaks at 0.51247.</p>
<p>The example applies in general to most other math operations besides addition (including multiplication, division, <code>log</code>, <code>sin</code>, and so on): do the math operation on the intervals <em>AI</em> and <em>BI</em>, and build a histogram of random results (products, quotients, etc.) that lie in the resulting interval to find out what distribution forms this way.</p>
<p><strong>Implementing other operations.</strong> In contrast to addition, multiplication, and division, certain other math operations are trivial to carry out in PSRNs. They include negation, as mentioned in (Karney 2016)<sup id="fnref:1:8" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>, and operations affecting the PSRN’s integer part only.</p>
<p>A promising way to connect PSRNs with other math operations (such as multiplication, <code>ln</code>, and <code>exp</code>) is to use “constructive reals” or “recursive reals”. See the section “Relation to Constructive Reals”, earlier.</p>
<p>A sampler can be created that uses the probabilities of getting each digit under the target distribution. But if the distribution is nondiscrete:</p>
<ul>
<li>These probabilities will depend on previous digits except for a very limited class of distributions (including uniform and exponential); see the <a href="#Setting_Digits_by_Digit_Probabilities"><strong>appendix</strong></a> for details.</li>
<li>For distributions outside that limited class, the sampler will be <em>limited-precision</em> (not <em>arbitrary-precision</em>) in practice, since it can hold only so many digit probabilities. For example, the works (Habibizad Navin et al., 2007)<sup id="fnref:23" role="doc-noteref"><a href="#fn:23" class="footnote" rel="footnote">23</a></sup>, (Nezhad et al., 2013)<sup id="fnref:24" role="doc-noteref"><a href="#fn:24" class="footnote" rel="footnote">24</a></sup> point to building a “tree” of such digit probabilities. <sup id="fnref:25" role="doc-noteref"><a href="#fn:25" class="footnote" rel="footnote">25</a></sup></li>
</ul>
<p>Finally, arithmetic with PSRNs may be possible if the result of the arithmetic is distributed with a known probability density function (PDF), allowing for an algorithm that implements rejection from the uniform or exponential distribution. An example of this is found in the <strong>UniformReciprocal</strong> algorithm above. However, that PDF may have an unbounded peak, thus ruling out rejection sampling in practice. For example, if <em>X</em> is a uniform PSRN in the interval [0, 1], then the distribution of <em>X</em><sup>3</sup> has the PDF <code>(1/3) / pow(X, 2/3)</code>, which has an unbounded peak at 0. While this rules out plain rejection samplers for <em>X</em><sup>3</sup> in practice, it’s still possible to sample powers of uniforms using PSRNs, which will be described later in this article.</p>
<p><strong>Reusing PSRNs.</strong> The arithmetic algorithms in this section may give incorrect results if the <em>same PSRN</em> is used more than once in different runs of these algorithms.</p>
<p>This issue happens in general when the original sampler uses the same random variate for different purposes in the algorithm (an example is “<em>W</em>*<em>Y</em>, (1−<em>W</em>)*<em>Y</em>”, where <em>W</em> and <em>Y</em> are independent random variates (Devroye 1986, p. 394)<sup id="fnref:26" role="doc-noteref"><a href="#fn:26" class="footnote" rel="footnote">26</a></sup>). In this case, if one PSRN spawns additional PSRNs (so that they become <em>dependent</em> on the first), those additional PSRNs may become inaccurate once additional digits of the first PSRN are sampled uniformly at random. (This is not always the case, but it’s hard to characterize when the additional PSRNs become inaccurate this way and when not.)</p>
<p>This issue is easy to see for the <strong>UniformAddRational</strong> or <strong>UniformMultiplyRational</strong> algorithm when it’s called more than once with the same PSRN and the same rational number: although the same random variate ought to be returned each time, in reality different variates will be generated this way with probability 1, especially when additional digits are sampled from them afterwards.</p>
<p>It might be believed that the issue just described could be solved by the algorithm below:</p>
<p><em>Assume we want to multiply the same PSRN by different numbers. Let vec be a vector of rational numbers to multiply the same PSRN by, and let vec[i] be the rational number at position i of the vector (positions start at 0).</em></p>
<ol>
<li><em>Set i to 0, set <strong>a</strong> to the input PSRN, set num to vec[i], and set ‘output’ to an empty list.</em></li>
<li><em>Set ret to the result of <strong>UniformMultiplyRational</strong> with the PSRN <strong>a</strong> and the rational number num.</em></li>
<li><em>Add a pointer to ret to the list ‘output’. If vec[i] was the last number in the vector, stop this algorithm.</em></li>
<li><em>Set <strong>a</strong> to point to ret, then add 1 to i, then set num to vec[i]/vec[i−1], then go to step 2.</em></li>
</ol>
<p>However, even this algorithm doesn’t ensure that the output PSRNs will be exactly proportional to the same random variate. An example: Let <strong>a</strong> be the PSRN 0…. (or the interval [0.0, 1.0]), then let <strong>b</strong> be the result of <strong>UniformMultiplyRational</strong>(<strong>a</strong>, 1/2), then let <strong>c</strong> be the result of <strong>UniformMultiplyRational</strong>(<strong>b</strong>, 1/3). One possible result for <strong>b</strong> is 0.41… and for <strong>c</strong> is 0.138…. Now we fill <strong>a</strong>, <strong>b</strong>, and <strong>c</strong> with uniform random bits. Thus, as one possible result, <strong>a</strong> is now 0.13328133…, <strong>b</strong> is now 0.41792367…, and <strong>c</strong> is now 0.13860371…. Here, however, <strong>c</strong> divided by <strong>b</strong> is not exactly 1/3, although it’s close, and <strong>b</strong> divided by <strong>a</strong> is far from 1/2 (especially since <strong>a</strong> was very coarse to begin with). Although this example shows PSRNs with decimal digits, the situation is worse with binary digits.</p>
<p><a id="Building_Blocks"></a></p>
<h2 id="building-blocks">Building Blocks</h2>
<p>This document relies on several building blocks described in this section.</p>
<p>One of them is the “geometric bag” technique by Flajolet and others (2010)<sup id="fnref:7:2" role="doc-noteref"><a href="#fn:7" class="footnote" rel="footnote">7</a></sup>, which generates heads or tails with a probability that is built up digit by digit.</p>
<p><a id="SampleGeometricBag"></a></p>
<h3 id="samplegeometricbag">SampleGeometricBag</h3>
<p>The algorithm <strong>SampleGeometricBag</strong> returns 1 with a probability built up by a uniform PSRN’s fractional part. (Flajolet et al., 2010)<sup id="fnref:7:3" role="doc-noteref"><a href="#fn:7" class="footnote" rel="footnote">7</a></sup> described an algorithm for the base-2 (binary) case, but that algorithm is difficult to apply to other digit bases. Thus the following is a general version of the algorithm for any digit base. For convenience, this algorithm ignores the PSRN’s integer part and sign.</p>
<ol>
<li>Set <em>i</em> to 0, and set <strong>b</strong> to a uniform PSRN with a positive sign and an integer part of 0.</li>
<li>If the item at position <em>i</em> of the input PSRN’s fractional part is unsampled (that is, not set to a digit), set the item at that position to a digit chosen uniformly at random, increasing the fractional part’s capacity as necessary (positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc.), and append the result to that fractional part’s digit expansion. Do the same for <strong>b</strong>.</li>
<li>Let <em>da</em> be the digit at position <em>i</em> of the input PSRN’s fractional part, and let <em>db</em> be the corresponding digit for <strong>b</strong>. Return 0 if <em>da</em> is less than <em>db</em>, or 1 if <em>da</em> is greater than <em>db</em>.</li>
<li>Add 1 to <em>i</em> and go to step 2.</li>
</ol>
<p>For base 2, the following <strong>SampleGeometricBag</strong> algorithm can be used, which is closer to the one given in the Flajolet paper. It likewise ignores the input PSRN’s integer part and sign.</p>
<ol>
<li>Set <em>N</em> to 0.</li>
<li>With probability 1/2, go to the next step. Otherwise, add 1 to <em>N</em> and repeat this step. (When the algorithm moves to the next step, <em>N</em> is what the Flajolet paper calls a <em>geometric</em> random variate (with parameter 1/2), hence the name “geometric bag”, but the terminology “geometric random variate” is avoided in this article since it has several conflicting meanings in academic works.)</li>
<li>If the item at position <em>N</em> in the uniform PSRN’s fractional part (positions start at 0) is not set to a digit (for example, 0 or 1 for base 2), set the item at that position to a digit chosen uniformly at random (for example, either 0 or 1 for base 2), increasing the fractional part’s capacity as necessary. (As a result of this step, there may be “gaps” in the uniform PSRN where no digit was sampled yet.)</li>
<li>Return the item at position <em>N</em>.</li>
</ol>
<p>For more on why these two algorithms are equivalent, see the appendix.</p>
<p><strong>SampleGeometricBagComplement</strong> is the same as the <strong>SampleGeometricBag</strong> algorithm, except the return value is 1 minus the original return value. The result is that if <strong>SampleGeometricBag</strong> outputs 1 with probability <em>U</em>, <strong>SampleGeometricBagComplement</strong> outputs 1 with probability 1 − <em>U</em>.</p>
<p><a id="FillGeometricBag"></a></p>
<h3 id="fillgeometricbag">FillGeometricBag</h3>
<p><strong>FillGeometricBag</strong> takes a uniform PSRN and generates a number whose fractional part has <code>p</code> digits as follows:</p>
<ol>
<li>For each position in [0, <code>p</code>), if the item at that position in the uniform PSRN’s fractional part is unsampled, set the item there to a digit chosen uniformly at random (for example, either 0 or 1 for binary), increasing the fractional part’s capacity as necessary. (Positions start at 0 where 0 is the most significant digit after the point, 1 is the next, etc. See also (Oberhoff 2018, sec. 8)<sup id="fnref:13:2" role="doc-noteref"><a href="#fn:13" class="footnote" rel="footnote">13</a></sup>.)</li>
<li>Let <code>sign</code> be -1 if the PSRN’s sign is negative, or 1 otherwise; let <code>ipart</code> be the PSRN’s integer part; and let <code>bag</code> be the PSRN’s fractional part. Take the first <code>p</code> digits of <code>bag</code> and return <code>sign</code> * (<code>ipart</code> + bag[0] * <em>b</em><sup>−0−1</sup> + bag[1] * <em>b</em><sup>−1−1</sup> + … + bag[<code>p</code>−1] * <em>b</em><sup>−(<code>p</code>−1)−1</sup>), where <em>b</em> is the base, or radix.</li>
</ol>
<p>After step 1, if it somehow happens that digits beyond <code>p</code> in the PSRN’s fractional part were already sampled (that is, they were already set to a digit), then the implementation could choose instead to—</p>
<ul>
<li>fill all unsampled digits between the first and the last set digit,</li>
<li>round the number represented by the PSRN to a number whose fractional part has <code>p</code> digits, with a rounding mode of choice (and without further modifying the PSRN), and</li>
<li>return the rounded number.</li>
</ul>
<p>For example, if <code>p</code> is 4, <em>b</em> is 10, and the PSRN is 0.3437500… or 0.3438500…, the implementation could use a round-to-nearest mode to round the number that the PSRN represents to the number 0.3438 or 0.3439, respectively, and return the rounded number; because this is a PSRN with an “infinite” but unsampled digit expansion, there is no tie-breaking such as “ties to even” applied here.</p>
<p><a id="kthsmallest"></a></p>
<h3 id="kthsmallest">kthsmallest</h3>
<p>The <strong>kthsmallest</strong> method generates the ‘k’th smallest ‘bitcount’-digit uniform random variate in the interval [0, 1] out of ‘n’ of them (also known as the ‘n’th <em>order statistic</em>), is also relied on by this beta sampler. It is used when both <code>a</code> and <code>b</code> are integers, based on the known property that a beta random variate in this case is the <code>a</code>th smallest uniform random variate between 0 and 1 out of <code>a + b - 1</code> of them (Devroye 1986, p. 431)<sup id="fnref:26:1" role="doc-noteref"><a href="#fn:26" class="footnote" rel="footnote">26</a></sup>.</p>
<p><strong>kthsmallest</strong>, however, doesn’t simply generate ‘n’ ‘bitcount’-digit numbers and then sort them. Rather, it builds up their digit expansions digit by digit, via PSRNs. It uses the observation that (in the binary case) each uniform random variate between 0 and 1 is either less than half or greater than half with equal probability; thus, the number of uniform numbers that are less than half vs. greater than half follows a binomial(n, 1/2) distribution (and of the numbers less than half, say, the less-than-one-quarter vs. greater-than-one-quarter numbers follows the same distribution, and so on). Thanks to this observation, the algorithm can generate a sorted sample “on the fly”. A similar observation applies to other bases than base 2 if we use the multinomial distribution instead of the binomial distribution. I am not aware of any other article or paper (besides one by me) that describes the <strong>kthsmallest</strong> algorithm given here.</p>
<p>The algorithm is as follows:</p>
<ol>
<li>Create <code>n</code> uniform PSRNs with positive sign and an integer part of 0.</li>
<li>Set <code>index</code> to 1.</li>
<li>If <code>index <= k</code> and <code>index + n >= k</code>:
<ol>
<li>Generate <strong>v</strong>, a multinomial random vector with <em>b</em> probabilities equal to 1/<em>b</em>, where <em>b</em> is the base, or radix (for the binary case, <em>b</em> = 2, so this is equivalent to generating <code>LC</code>, a binomial random variate with parameters <code>n</code> and 0.5, and setting <strong>v</strong> to {<code>LC</code>, <code>n - LC</code>}).</li>
<li>Starting at <code>index</code>, append the digit 0 to the first <strong>v</strong>[0] PSRNs, a 1 digit to the next <strong>v</strong>[1] PSRNs, and so on to appending a <em>b</em> − 1 digit to the last <strong>v</strong>[<em>b</em> − 1] PSRNs (for the binary case, this means appending a 0 bit to the first <code>LC</code> PSRNs and a 1 bit to the next <code>n - LC</code> PSRNs).</li>
<li>For each integer <em>i</em> in [0, <em>b</em>): If <strong>v</strong>[<em>i</em>] > 1, repeat step 3 and these substeps with <code>index</code> = <code>index</code> + <strong>v</strong>[0] + <strong>v</strong>[1] + … + <strong>v</strong>[<em>i</em> − 1] and <code>n</code> = <strong>v</strong>[<em>i</em>]. (For the binary case, this means: If <code>LC > 1</code>, repeat step 3 and these substeps with the same <code>index</code> and <code>n = LC</code>; then, if <code>n - LC > 1</code>, repeat step 3 and these substeps with <code>index = index + LC</code>, and <code>n = n - LC</code>).</li>
</ol>
</li>
<li>Take the <code>k</code>th PSRN (starting at 1), then optionally fill it with uniform random digits as necessary to give its fractional part <code>bitcount</code> many digits (similarly to <strong>FillGeometricBag</strong> above), then return that number. (Note that the beta sampler described later chooses to fill the PSRN this way via this algorithm.)</li>
</ol>
<p><a id="Power_of_Uniform_Subalgorithm"></a></p>
<h3 id="power-of-uniform-subalgorithm">Power-of-Uniform Subalgorithm</h3>
<p>The power-of-uniform subalgorithm is used for certain cases of the beta sampler below. It returns <em>U</em><sup><em>px</em>/<em>py</em></sup>, where <em>U</em> is a uniform random variate in the interval [0, 1] and <em>px</em>/<em>py</em> is greater than 1, but unlike the naïve algorithm it supports an arbitrary precision, uses only random bits, and avoids floating-point arithmetic. It also uses a <em>complement</em> flag to determine whether to return 1 minus the result.</p>
<p>It makes use of a number of algorithms as follows:</p>
<ul>
<li>It uses an algorithm for <a href="https://peteroupc.github.io/randmisc.html"><strong>sampling unbounded monotone PDFs</strong></a>, which in turn is similar to the inversion-rejection algorithm in (Devroye 1986, ch. 7, sec. 4.4)<sup id="fnref:26:2" role="doc-noteref"><a href="#fn:26" class="footnote" rel="footnote">26</a></sup>. This is needed because when <em>px</em>/<em>py</em> is greater than 1, the distribution of <em>U</em><sup><em>px</em>/<em>py</em></sup> has the PDF <code>(py/px) / pow(U, 1-py/px)</code>, which has an unbounded peak at 0.</li>
<li>It uses a number of Bernoulli factory algorithms, including <strong>SampleGeometricBag</strong> and some algorithms described in “<a href="https://peteroupc.github.io/bernoulli.html"><strong>Bernoulli Factory Algorithms</strong></a>”.</li>
</ul>
<p>However, this algorithm currently only supports generating a PSRN with base-2 (binary) digits in its fractional part.</p>
<p>The power-of-uniform algorithm is as follows:</p>
<ol>
<li>Set <em>i</em> to 1.</li>
<li>Call the <strong>algorithm for (<em>a</em>/<em>b</em>)<sup><em>z</em></sup></strong> described in “<a href="https://peteroupc.github.io/bernoulli.html"><strong>Bernoulli Factory Algorithms</strong></a>”, with parameters <code>a = 1, b = 2, z = py/px</code>. If the call returns 1 and <em>i</em> is less than <em>n</em>, add 1 to <em>i</em> and repeat this step. If the call returns 1 and <em>i</em> is <em>n</em> or greater, return 1 if the <em>complement</em> flag is 1 or 0 otherwise (or return a uniform PSRN with a positive sign, an integer part of 0, and a fractional part filled with exactly <em>n</em> ones or zeros, respectively).</li>
<li>As a result, we will now sample a number in the interval [2<sup>−<em>i</em></sup>, 2<sup>−(<em>i</em> − 1)</sup>). We now have to generate a uniform random variate <em>X</em> in this interval, then accept it with probability (<em>py</em> / (<em>px</em> * 2<sup><em>i</em></sup>)) / <em>X</em><sup>1 − <em>py</em> / <em>px</em></sup>; the 2<sup><em>i</em></sup> in this formula is to help avoid very low probabilities for sampling purposes. The following steps will achieve this without having to use floating-point arithmetic.</li>
<li>Create a positive-sign zero-integer-part uniform PSRN, then create a <em>geobag</em> input coin that returns the result of <strong>SampleGeometricBag</strong> on that PSRN.</li>
<li>Create a <em>powerbag</em> input coin that does the following: “Call the <strong>algorithm for <em>λ</em><sup><em>x</em>/<em>y</em></sup></strong>, described in ‘<a href="https://peteroupc.github.io/bernoulli.html#lambda__x___y"><strong>Bernoulli Factory Algorithms</strong></a>’, using the <em>geobag</em> input coin and with <em>x</em>/<em>y</em> = 1 − <em>py</em> / <em>px</em>, and return the result.”</li>
<li>Append <em>i</em> − 1 zero-digits followed by a single one-digit to the PSRN’s fractional part. This will allow us to sample a uniform random variate limited to the interval mentioned earlier.</li>
<li>Call the <strong>algorithm for ϵ / λ</strong>, described in “<a href="https://peteroupc.github.io/bernoulli.html#x03F5_lambda"><strong>Bernoulli Factory Algorithms</strong></a>”, using the <em>powerbag</em> input coin (which represents <em>b</em>) and with ϵ = <em>py</em>/(<em>px</em> * 2<sup><em>i</em></sup>) (which represents <em>a</em>), thus returning 1 with probability <em>a</em>/<em>b</em>. If the call returns 1, the PSRN was accepted, so do the following:
<ol>
<li>If the <em>complement</em> flag is 1, make each zero-digit in the PSRN’s fractional part a one-digit and vice versa.</li>
<li>Optionally, fill the PSRN with uniform random digits as necessary to give its fractional part <em>n</em> digits (similarly to <strong>FillGeometricBag</strong> above), where <em>n</em> is a precision parameter. Then, return the PSRN.</li>
</ol>
</li>
<li>If the call to the algorithm for ϵ / λ returns 0, remove all but the first <em>i</em> digits from the PSRN’s fractional part, then go to step 7.</li>
</ol>
<p><a id="Algorithms_for_the_Beta_and_Exponential_Distributions"></a></p>
<h2 id="algorithms-for-the-beta-and-exponential-distributions">Algorithms for the Beta and Exponential Distributions</h2>
<p> </p>
<p><a id="Beta_Distribution"></a></p>
<h3 id="beta-distribution">Beta Distribution</h3>
<p>All the building blocks are now in place to describe a <em>new</em> algorithm to sample the beta distribution, described as follows. It takes three parameters: <em>a</em> >= 1 and <em>b</em> >= 1 (or one parameter is 1 and the other is greater than 0 in the binary case) are the parameters to the beta distribution, and <em>p</em> > 0 is a precision parameter.</p>
<ol>
<li>Special cases:
<ul>
<li>If <em>a</em> = 1 and <em>b</em> = 1, return a positive-sign zero-integer-part uniform PSRN.</li>
<li>If <em>a</em> and <em>b</em> are both integers, return the result of <strong>kthsmallest</strong> with <code>n = a - b + 1</code> and <code>k = a</code></li>
<li>In the binary case, if <em>a</em> is 1 and <em>b</em> is less than 1, call the <strong>power-of-uniform subalgorithm</strong> described below, with <em>px</em>/<em>py</em> = 1/<em>b</em>, and the <em>complement</em> flag set to 1, and return the result of that algorithm as is (without filling it as described in substep 7.2 of that algorithm).</li>
<li>In the binary case, if <em>b</em> is 1 and <em>a</em> is less than 1, call the <strong>power-of-uniform subalgorithm</strong> described below, with <em>px</em>/<em>py</em> = 1/<em>a</em>, and the <em>complement</em> flag set to 0, and return the result of that algorithm as is (without filling it as described in substep 7.2 of that algorithm).</li>
</ul>
</li>
<li>If <em>a</em> > 2 and <em>b</em> > 2, do the following steps, which split <em>a</em> and <em>b</em> into two parts that are faster to simulate (and implement the generalized rejection strategy in (Devroye 1986, top of page 47)<sup id="fnref:26:3" role="doc-noteref"><a href="#fn:26" class="footnote" rel="footnote">26</a></sup>):
<ol>
<li>Set <em>aintpart</em> to floor(<em>a</em>) − 1, set <em>bintpart</em> to floor(<em>b</em>) − 1, set <em>arest</em> to <em>a</em> − <em>aintpart</em>, and set <em>brest</em> to <em>b</em> − <em>bintpart</em>.</li>
<li>Do a separate (recursive) run of this algorithm, but with <em>a</em> = <em>aintpart</em> and <em>b</em> = <em>bintpart</em>. Set <em>bag</em> to the PSRN created by the run.</li>
<li>Create an input coin <em>geobag</em> that returns the result of <strong>SampleGeometricBag</strong> using the given PSRN. Create another input coin <em>geobagcomp</em> that returns the result of <strong>SampleGeometricBagComplement</strong> using the given PSRN.</li>
<li>Call the <strong>algorithm for <em>λ</em><sup><em>x</em>/<em>y</em></sup></strong>, described in “<a href="https://peteroupc.github.io/bernoulli.html"><strong>Bernoulli Factory Algorithms</strong></a>”, using the <em>geobag</em> input coin and <em>x</em>/<em>y</em> = <em>arest</em>/1, then call the same algorithm using the <em>geobagcomp</em> input coin and <em>x</em>/<em>y</em> = <em>brest</em>/1. If both calls return 1, return <em>bag</em>. Otherwise, go to substep 2.</li>
</ol>
</li>
<li>Create an positive-sign zero-integer-part uniform PSRN. Create an input coin <em>geobag</em> that returns the result of <strong>SampleGeometricBag</strong> using the given PSRN. Create another input coin <em>geobagcomp</em> that returns the result of <strong>SampleGeometricBagComplement</strong> using the given PSRN.</li>
<li>Remove all digits from the PSRN’s fractional part. This will result in an “empty” uniform random variate between 0 and 1, <em>U</em>, for the following steps, which will accept <em>U</em> with probability <em>U</em><sup>a−1</sup>*(1−<em>U</em>)<sup>b−1</sup> (the proportional probability for the beta distribution), as <em>U</em> is built up.</li>
<li>Call the <strong>algorithm for <em>λ</em><sup><em>x</em>/<em>y</em></sup></strong>, described in “<a href="https://peteroupc.github.io/bernoulli.html"><strong>Bernoulli Factory Algorithms</strong></a>”, using the <em>geobag</em> input coin and <em>x</em>/<em>y</em> = <em>a</em> − 1)/1 (thus returning with probability <em>U</em><sup>a−1</sup>). If the result is 0, go to step 4.</li>
<li>Call the same algorithm using the <em>geobagcomp</em> input coin and <em>x</em>/<em>y</em> = (<em>b</em> − 1)/1 (thus returning 1 with probability (1−<em>U</em>)<sup>b−1</sup>). If the result is 0, go to step 4. (Note that this step and the previous step don’t depend on each other and can be done in either order without affecting correctness, and this is taken advantage of in the Python code below.)</li>
<li><em>U</em> was accepted, so return the result of <strong>FillGeometricBag</strong>.</li>
</ol>
<p>Once a PSRN is accepted by the steps above, optionally fill the unsampled digits of the PSRN’s fractional part with uniform random digits as necessary to give the number a <em>p</em>-digit fractional part (similarly to <strong>FillGeometricBag</strong>), then return the resulting number.</p>
<blockquote>
<p><strong>Notes:</strong></p>
<ul>
<li>A beta random variate with parameters 1/<em>x</em> and 1 is the same as a uniform random variate in [0, 1] raised to the power of <em>x</em>.</li>
<li>For the beta distribution, the bigger <code>alpha</code> or <code>beta</code> is, the smaller the area of acceptance becomes (and the greater the probability that random variates get rejected by steps 5 and 6, raising its run-time). This is because <code>max(u^(alpha-1)*(1-u)^(beta-1))</code>, the peak of the PDF, approaches 0 as the parameters get bigger. To deal with this, step 2 was included, which under certain circumstances breaks the PDF into two parts that are relatively trivial to sample (in terms of bit complexity).</li>
</ul>
</blockquote>
<p><a id="Exponential_Distribution"></a></p>
<h3 id="exponential-distribution">Exponential Distribution</h3>
<p>We also have the necessary building blocks to describe how to sample e-rands. An e-rand consists of four numbers: the first is a multiple of 1/(2<sup><em>k</em></sup>), the second is <em>k</em>, the third is the integer part (initially −1 to indicate the integer part wasn’t sampled yet), and the fourth, <em>λ</em>, is the rate parameter of the exponential distribution (<em>λ</em>>0). (Because exponential random variates are always 0 or greater, the e-rand’s sign is implicitly positive.) In the Python code, e-rands are as described, except <em>λ</em> must be a rational number and its numerator and denominator take up a parameter each.</p>
<p>To sample bit <em>k</em> after the binary point of an exponential random variate with rate <em>λ</em> (where <em>k</em> = 1 means the first digit after the point, <em>k</em> = 2 means the second, etc.), call the <strong>LogisticExp</strong> algorithm (see “<a href="https://peteroupc.github.io/bernoulli.html""><strong>Bernoulli Factory Algorithms</strong></a> with <em>z</em>=<em>λ</em> and <em>prec</em> = <em>k</em>.</p>
<p>The <strong>ExpRandLess</strong> algorithm is a special case of the general <strong>RandLess</strong> algorithm given earlier. It compares two e-rands <strong>a</strong> and <strong>b</strong> (and samples additional bits from them as necessary) and returns 1 if <strong>a</strong> turns out to be less than <strong>b</strong>, or 0 otherwise. (Note that <strong>a</strong> and <strong>b</strong> are allowed to have different <em>λ</em> parameters.)</p>
<ol>
<li>If <strong>a</strong>’s integer part wasn’t sampled yet, call the <strong>ExpMinus</strong> algorithm (see “<a href="https://peteroupc.github.io/bernoulli.html""><strong>Bernoulli Factory Algorithms</strong></a> with <em>z</em>=<em>λ</em>, until the call returns 0, then set the integer part to the number of times 1 was returned this way. Do the same for <strong>b</strong>.</li>
<li>Return 1 if <strong>a</strong>’s integer part is less than <strong>b</strong>’s, or 0 if <strong>a</strong>’s integer part is greater than <strong>b</strong>’s.</li>
<li>Set <em>i</em> to 0.</li>
<li>If <strong>a</strong>’s fractional part has <em>i</em> or fewer bits, call the <strong>LogisticExp</strong> algorithm (see “<a href="https://peteroupc.github.io/bernoulli.html""><strong>Bernoulli Factory Algorithms</strong></a> with <em>z</em>=<em>λ</em> and <em>prec</em> = <em>i</em> + 1, and append the result to that fractional part’s binary expansion. (For example, if the implementation stores the binary expansion as a packed integer and a size, the implementation can shift the packed integer by 1, add the result of the algorithm to that integer, then add 1 to the size.) Do the same for <strong>b</strong>.</li>
<li>Return 1 if <strong>a</strong>’s fractional part is less than <strong>b</strong>’s, or 0 if <strong>a</strong>’s fractional part is greater than <strong>b</strong>’s.</li>
<li>Add 1 to <em>i</em> and go to step 4.</li>
</ol>
<p>The <strong>ExpRandFill</strong> algorithm takes an e-rand and generates a number whose fractional part has <code>p</code> digits as follows:</p>
<ol>
<li>For each position <em>i</em> in [0, <code>p</code>), if the item at that position in the e-rand’s fractional part is unsampled, call the <strong>LogisticExp</strong> algorithm (see “<a href="https://peteroupc.github.io/bernoulli.html""><strong>Bernoulli Factory Algorithms</strong></a> with <em>z</em>=<em>λ</em> and <em>prec</em> = <em>i</em> + 1, and set the item at position <em>i</em> to the result (which will be either 0 or 1), increasing the fractional part’s capacity as necessary. (Bit positions start at 0 where 0 is the most significant bit after the point, 1 is the next, etc. See also (Oberhoff 2018, sec. 8)<sup id="fnref:13:3" role="doc-noteref"><a href="#fn:13" class="footnote" rel="footnote">13</a></sup>.)</li>
<li>Let <code>sign</code> be -1 if the e-rand’s sign is negative, or 1 otherwise; let <code>ipart</code> be the e-rand’s integer part; and let <code>bag</code> be the PSRN’s fractional part. Take the first <code>p</code> digits of <code>bag</code> and return <code>sign</code> * (<code>ipart</code> + bag[0] * 2<sup>−0−1</sup> + bag[1] * 2<sup>−1−1</sup> + … + bag[<code>p</code>−1] * 2<sup>−(<code>p</code>−1)−1</sup>).</li>
</ol>
<p>See the discussion in <strong>FillGeometricBag</strong> for advice on how to handle the case when it somehow happens that bits beyond <code>p</code> in the PSRN’s fractional part were already sampled (that is, they were already set to a digit) after step 2 of this algorithm.</p>
<p>Here is a third algorithm (called <strong>ExpRand</strong>) that generates a <em>uniform PSRN</em>, rather than an e-rand, that follows the exponential distribution. In the algorithm, the rate <em>λ</em> is given as a rational number greater than 0. The method is based on von Neumann’s algorithm (von Neumann 1951)<sup id="fnref:9:1" role="doc-noteref"><a href="#fn:9" class="footnote" rel="footnote">9</a></sup>.</p>
<ol>
<li>Set <em>recip</em> to 1/<em>λ</em>, and set <em>highpart</em> to 0.</li>
<li>Set <em>u</em> to the result of <strong>RandUniformFromReal</strong> with the parameter <em>recip</em>.</li>
<li>Set <em>val</em> to point to the same value as <em>u</em>, and set <em>accept</em> to 1.</li>
<li>Set <em>v</em> to the result of <strong>RandUniformFromReal</strong> with the parameter <em>recip</em>.</li>
<li>Run the <strong>URandLess</strong> algorithm on <em>u</em> and <em>v</em>, in that order. If the call returns 0, set <em>u</em> to <em>v</em>, then set <em>accept</em> to 1 minus <em>accept</em>, then go to step 4.</li>
<li>If <em>accept</em> is 1, add <em>highpart</em> to <em>val</em> via the <strong>UniformAddRational</strong> algorithm given earlier, then return <em>val</em>.</li>
<li>Add <em>recip</em> to <em>highpart</em> and go to step 2.</li>
</ol>
<p>The following alternative version of the previous algorithm (called <strong>ExpRand2</strong>) includes Karney’s improvement to the von Neumann algorithm (Karney 2016)<sup id="fnref:1:9" role="doc-noteref"><a href="#fn:1" class="footnote" rel="footnote">1</a></sup>, namely a so-called “early rejection step”. The algorithm here allows an arbitrary rate parameter (<em>λ</em>), given as a rational number greater than 0, unlike with the von Neumann and Karney algorithms, where <em>λ</em> is 1.</p>
<ol>
<li>Set <em>recip</em> to 1/<em>λ</em>, and set <em>highpart</em> to 0.</li>
<li>Set <em>u</em> to the result of <strong>RandUniformFromReal</strong> with the parameter <em>recip</em>.</li>
<li>Run the <strong>URandLessThanReal</strong> algorithm on <em>u</em> with the parameter <em>recip</em>/2. If the call returns 0, add <em>recip</em>/2 to <em>highpart</em> and go to step 2. (This is Karney’s “early rejection step”, where the parameter is 1/2 when <em>λ</em> is 1. However, Fan et al. (2019)<sup id="fnref:27" role="doc-noteref"><a href="#fn:27" class="footnote" rel="footnote">27</a></sup> point out that the parameter 1/2 in Karney’s “early rejection step” is not optimal.)</li>
<li>Set <em>val</em> to point to the same value as <em>u</em>, and set <em>accept</em> to 1.</li>
<li>Set <em>v</em> to the result of <strong>RandUniformFromReal</strong> with the parameter <em>recip</em>.</li>
<li>Run the <strong>URandLess</strong> algorithm on <em>u</em> and <em>v</em>, in that order. If the call returns 0, set <em>u</em> to <em>v</em>, then set <em>accept</em> to 1 minus <em>accept</em>, then go to step 5.</li>
<li>If <em>accept</em> is 1, add <em>highpart</em> to <em>val</em> via the <strong>UniformAddRational</strong> algorithm given earlier, then return <em>val</em>.</li>
<li>Add <strong><em>recip</em>/2</strong> to <em>highpart</em> and go to step 2.</li>
</ol>
<blockquote>
<p><strong>Note:</strong> A Laplace (double exponential) random variate is then implemented by giving the PSRN returned by <strong>ExpRand</strong> or <strong>ExpRand2</strong> a random sign (with equal probability, the PSRN’s sign is either positive or negative).</p>
</blockquote>
<p><a id="Sampler_Code"></a></p>
<h2 id="sampler-code">Sampler Code</h2>
<p>The following Python code implements the beta sampler described in this document. It relies on two Python modules I wrote:</p>
<ul>
<li>”<a href="https://github.com/peteroupc/peteroupc.github.io/blob/master/bernoulli.py"><strong>bernoulli.py</strong></a>”, which collects a number of Bernoulli factories, some of which are relied on by the following code.</li>
<li>”<a href="https://github.com/peteroupc/peteroupc.github.io/blob/master/randomgen.py"><strong>randomgen.py</strong></a>”, which collects a number of random variate generation methods, including <code>kthsmallest</code>, as well as the <code>RandomGen</code> class.</li>
</ul>
<p>Note that the code uses floating-point arithmetic only to convert the result of the sampler to a convenient form, namely a floating-point number.</p>
<p>This code is far from fast, though, at least in Python.</p>
<p>The Python code below supports only rational-valued <em>λ</em> parameters in the exponential sampler.</p>
<p>```
import math
import random
import bernoulli
from randomgen import RandomGen
from fractions import Fraction</p>
<p>def _toreal(ret, precision):
# NOTE: Although we convert to a floating-point
# number here, this is not strictly necessary and
# is merely for convenience.
return ret*1.0/(1«precision)</p>
<p>def _urand_to_geobag(bag):
return [(bag[0]»(bag[1]-1-i))&1 for i in range(bag[1])]</p>
<p>def _power_of_uniform_greaterthan1(bern, power, complement=False, precision=53):
return bern.fill_geometric_bag(
_power_of_uniform_greaterthan1_geobag(bern, power, complement), precision
)</p>
<p>def _power_of_uniform_greaterthan1_geobag(bern, power, complement=False, precision=53):
if power<1:
raise ValueError(“Not supported”)
if power==1:
return [] # Empty uniform random variate
i=1
powerfrac=Fraction(power)