1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.rng.sampling.distribution;
18
19 import java.util.function.Supplier;
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.core.source64.SplitMix64;
22 import org.apache.commons.rng.sampling.RandomAssert;
23 import org.apache.commons.rng.sampling.distribution.StableSampler.Beta0CMSStableSampler;
24 import org.apache.commons.rng.sampling.distribution.StableSampler.Beta0WeronStableSampler;
25 import org.apache.commons.rng.sampling.distribution.StableSampler.CMSStableSampler;
26 import org.apache.commons.rng.sampling.distribution.StableSampler.SpecialMath;
27 import org.apache.commons.rng.sampling.distribution.StableSampler.WeronStableSampler;
28 import org.apache.commons.rng.simple.RandomSource;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.Test;
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52 class StableSamplerTest {
53
54 private static final double PI_2 = Math.PI / 2;
55
56 private static final double PI_4 = Math.PI / 4;
57
58 private static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
59
60
61 private static final double DU = 0x1.0p-53;
62
63 private static final double SMALL_W = 6.564735882096453E-19;
64
65 private static final double TAIL_W = 7.569274694148063;
66
67 private static final double LARGE_W = 4 * TAIL_W;
68
69 private static final double SMALLEST_ALPHA = 1.0 - Math.nextDown(1.0);
70
71 private static final double VALID_ALPHA = 1.23;
72 private static final double VALID_BETA = 0.23;
73 private static final double VALID_GAMMA = 2.34;
74 private static final double VALID_DELTA = 3.45;
75
76 @Test
77 void testAlphaZeroThrows() {
78 assertConstructorThrows(0.0, VALID_BETA, VALID_GAMMA, VALID_DELTA);
79 }
80
81 @Test
82 void testAlphaBelowZeroThrows() {
83 assertConstructorThrows(Math.nextDown(0.0), VALID_BETA, VALID_GAMMA, VALID_DELTA);
84 }
85
86 @Test
87 void testAlphaTooCloseToZeroThrows() {
88
89
90
91 final UniformRandomProvider rng = new SplitMix64(0L);
92 StableSampler s = StableSampler.of(rng, SMALLEST_ALPHA, VALID_BETA, VALID_GAMMA, VALID_DELTA);
93 Assertions.assertNotNull(s);
94
95
96 final double alphaTooSmall = SMALLEST_ALPHA / 2;
97 Assertions.assertNotEquals(0.0, alphaTooSmall, "Expected alpha to be positive");
98 Assertions.assertEquals(1.0, 1 - alphaTooSmall, "Expected rounding to 1");
99
100
101 assertConstructorThrows(alphaTooSmall, VALID_BETA, VALID_GAMMA, VALID_DELTA);
102 }
103
104 @Test
105 void testAlphaAboveTwoThrows() {
106 assertConstructorThrows(Math.nextUp(2.0), VALID_BETA, VALID_GAMMA, VALID_DELTA);
107 }
108
109 @Test
110 void testAlphaNaNThrows() {
111 assertConstructorThrows(Double.NaN, VALID_BETA, VALID_GAMMA, VALID_DELTA);
112 }
113
114 @Test
115 void testBetaBelowMinusOneThrows() {
116 assertConstructorThrows(VALID_ALPHA, Math.nextDown(-1.0), VALID_GAMMA, VALID_DELTA);
117 }
118
119 @Test
120 void testBetaAboveOneThrows() {
121 assertConstructorThrows(VALID_ALPHA, Math.nextUp(1.0), VALID_GAMMA, VALID_DELTA);
122 }
123
124 @Test
125 void testBetaNaNThrows() {
126 assertConstructorThrows(VALID_ALPHA, Double.NaN, VALID_GAMMA, VALID_DELTA);
127 }
128
129 @Test
130 void testGammaNotStrictlyPositiveThrows() {
131 assertConstructorThrows(VALID_ALPHA, VALID_BETA, 0.0, VALID_DELTA);
132 }
133
134 @Test
135 void testGammaInfThrows() {
136 assertConstructorThrows(VALID_ALPHA, VALID_BETA, Double.POSITIVE_INFINITY, VALID_DELTA);
137 }
138
139 @Test
140 void testGammaNaNThrows() {
141 assertConstructorThrows(VALID_ALPHA, VALID_BETA, Double.NaN, VALID_DELTA);
142 }
143
144 @Test
145 void testDeltaInfThrows() {
146 assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.POSITIVE_INFINITY);
147 }
148
149 @Test
150 void testDeltaNegInfThrows() {
151 assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.NEGATIVE_INFINITY);
152 }
153
154 @Test
155 void testDeltaNaNThrows() {
156 assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.NaN);
157 }
158
159
160
161
162
163
164
165
166
167 private static void assertConstructorThrows(double alpha, double beta, double gamma, double delta) {
168 final UniformRandomProvider rng = new SplitMix64(0L);
169 Assertions.assertThrows(IllegalArgumentException.class,
170 () -> StableSampler.of(rng, alpha, beta, gamma, delta));
171 }
172
173
174
175
176
177
178
179
180
181
182 @Test
183 void testTauLimits() {
184
185 final double beta = 1;
186
187
188
189 Assertions.assertEquals(0.0, CMSStableSampler.getTau(2, beta));
190 Assertions.assertEquals(0.0, CMSStableSampler.getTau(0, beta));
191
192
193 for (int i = 0; i <= 512; i++) {
194
195 final double alpha = (double) i / 256;
196 final double tau = CMSStableSampler.getTau(alpha, beta);
197 final double expected = getTauOriginal(alpha, beta);
198 Assertions.assertEquals(expected, tau, 1e-15);
199
200
201 Assertions.assertEquals(tau, CMSStableSampler.getTau(2 - alpha, beta));
202 }
203
204
205 final double limit = beta / PI_2;
206 Assertions.assertEquals(limit, CMSStableSampler.getTau(1, beta));
207 for (double alpha : new double[] {1.01, 1 + 1e-6, 1, 1 - 1e-6, 0.99}) {
208 final double tau = CMSStableSampler.getTau(alpha, beta);
209 final double expected = getTauOriginal(alpha, beta);
210 Assertions.assertEquals(expected, tau, 1e-15);
211
212 Assertions.assertEquals(limit, tau, Math.abs(1 - alpha) + 1e-15);
213 }
214
215
216
217
218 Assertions.assertEquals(0.0, CMSStableSampler.getTau(1.3, 0.0));
219 Assertions.assertEquals(0.0, CMSStableSampler.getTau(1.5, Double.MIN_VALUE));
220 Assertions.assertNotEquals(0.0, CMSStableSampler.getTau(1.0, Double.MIN_VALUE));
221
222
223 Assertions.assertEquals(0.5, CMSStableSampler.getTau(1.5, beta));
224 Assertions.assertEquals(0.5, CMSStableSampler.getTau(0.5, beta));
225 Assertions.assertEquals(-0.5, CMSStableSampler.getTau(1.5, -beta));
226 Assertions.assertEquals(-0.5, CMSStableSampler.getTau(0.5, -beta));
227
228
229 final double tau1 = CMSStableSampler.getTau(Math.nextDown(1.5), 1);
230 final double tau2 = CMSStableSampler.getTau(1.5, 1);
231 final double tau3 = CMSStableSampler.getTau(Math.nextUp(1.5), 1);
232 Assertions.assertTrue(tau1 > tau2);
233 Assertions.assertTrue(tau2 > tau3);
234
235 Assertions.assertEquals(tau1, CMSStableSampler.getTau(2 - Math.nextDown(1.5), 1));
236 Assertions.assertEquals(tau2, CMSStableSampler.getTau(0.5, 1));
237 Assertions.assertEquals(tau3, CMSStableSampler.getTau(2 - Math.nextUp(1.5), 1));
238 }
239
240
241
242
243
244
245
246
247
248
249 private static double getTauOriginal(double alpha, double beta) {
250 final double eps = 1 - alpha;
251
252 double tau;
253
254
255
256
257
258
259 if (eps > -0.99) {
260
261 final double tan2 = eps == 0 ? 1 : Math.tan(eps * PI_2) / (eps * PI_2);
262 tau = beta / (tan2 * PI_2);
263 } else {
264
265 final double meps1 = 1 - eps;
266 final double tan2 = Math.tan(meps1 * PI_2) / (meps1 * PI_2);
267 tau = beta * PI_2 * eps * meps1 * tan2;
268 }
269
270 return tau;
271 }
272
273
274
275
276
277
278 @Test
279 void testA2IsNotZero() {
280
281
282 final double p0 = getU(Long.MIN_VALUE);
283 Assertions.assertEquals(-PI_4, p0);
284
285
286 final double p1 = getU(Long.MIN_VALUE + (1 << 10));
287 final double p2 = getU(Long.MAX_VALUE);
288 Assertions.assertNotEquals(-PI_4, p1);
289 Assertions.assertNotEquals(PI_4, p2);
290 Assertions.assertEquals(-PI_4 + PI_4 * DU, p1);
291 Assertions.assertEquals(PI_4 - PI_4 * DU, p2);
292
293 for (double phiby2 : new double[] {p1, p2}) {
294
295
296 final double a = phiby2 * SpecialMath.tan2(phiby2);
297 Assertions.assertEquals(Math.copySign(Math.nextDown(1.0), phiby2), a);
298 final double da = a * a;
299 final double a2 = 1 - da;
300
301 Assertions.assertNotEquals(0.0, a2);
302
303 Assertions.assertEquals(0x1.0p-52, a2);
304 }
305 }
306
307
308
309
310
311
312
313
314
315
316
317 @Test
318 void testZIsNotAlwaysAboveZero() {
319
320
321 final long x00 = Long.MIN_VALUE;
322 final long x0 = Long.MIN_VALUE + (1 << 10);
323 final long x1 = Long.MAX_VALUE;
324 Assertions.assertEquals(-PI_4, getU(x00));
325 Assertions.assertEquals(-PI_4 + DU * PI_4, getU(x0));
326 Assertions.assertEquals(PI_4 - DU * PI_4, getU(x1));
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 Assertions.assertEquals(0.0, computeNumerator(0.859375, 1, x00));
361
362 Assertions.assertTrue(0.0 > computeNumerator(0.9375, 1, x00));
363 Assertions.assertTrue(0.0 > computeNumerator(1.90625, 1, x00));
364
365
366
367
368 Assertions.assertTrue(0.0 < computeNumerator(0.859375, 1, x0));
369 Assertions.assertTrue(0.0 < computeNumerator(0.9375, 1, x0));
370 Assertions.assertTrue(0.0 < computeNumerator(1.90625, 1, x0));
371
372
373
374
375
376 Assertions.assertTrue(0.0 > computeNumerator(0.828125, 1, x0));
377 Assertions.assertTrue(0.0 > computeNumerator(1.291015625, -1, x1));
378
379
380
381
382
383 Assertions.assertEquals(-1, SpecialMath.tan2(getU(x00)) * getU(x00));
384
385
386 Assertions.assertTrue(-1 < SpecialMath.tan2(getU(x0)) * getU(x0));
387 Assertions.assertTrue(1 > SpecialMath.tan2(getU(x1)) * getU(x1));
388
389 final double beta = 0;
390 Assertions.assertEquals(0.0, computeNumerator(2, beta, x00));
391 Assertions.assertTrue(0.0 < computeNumerator(2, beta, x0));
392 Assertions.assertTrue(0.0 < computeNumerator(2, beta, x1));
393 Assertions.assertTrue(0.0 < computeNumerator(Math.nextDown(2), beta, x0));
394 Assertions.assertTrue(0.0 < computeNumerator(Math.nextDown(2), beta, x1));
395
396
397
398
399
400
401
402
403
404
405
406 final double alpha = 1;
407 Assertions.assertEquals(0.0, computeNumerator(alpha, 1, x00));
408
409 Assertions.assertTrue(0.0 < computeNumerator(alpha, 1, x0));
410 Assertions.assertTrue(0.0 < computeNumerator(alpha, -1, x1));
411
412
413
414 Assertions.assertTrue(0.0 < computeNumerator(alpha, Math.nextUp(-1), x00));
415 }
416
417
418
419
420
421
422
423
424
425 private static double computeNumerator(double alpha, double beta, long x) {
426 final double phiby2 = getU(x);
427 final double eps = 1 - alpha;
428 final double tau = CMSStableSampler.getTau(alpha, beta);
429
430 final double bb = SpecialMath.tan2(eps * phiby2);
431 final double b = eps * phiby2 * bb;
432
433 final double db = b * b;
434 final double b2 = 1 - db;
435
436 return b2 + 2 * phiby2 * bb * tau;
437 }
438
439
440
441
442
443
444
445 @Test
446 void testComputeDWhenZIsFiniteNonZero() {
447 final double[] zs = {Double.MIN_VALUE, Double.MAX_VALUE};
448
449 final double[] alphas = {2, 1.5, 1 + 1e-6, 1, 1 - 1e-6, 0.5, 0.01, 1e-10, SMALLEST_ALPHA};
450 for (final double alpha : alphas) {
451
452 for (final double z : zs) {
453
454 Assertions.assertNotEquals(Double.NaN, computeD(alpha, z));
455 }
456
457
458
459
460
461 final double d0 = computeD(alpha, 0);
462 if (alpha < 1) {
463
464 Assertions.assertEquals(Double.NaN, d0);
465 } else if (alpha == 1) {
466
467 Assertions.assertEquals(Double.NaN, d0);
468 } else {
469
470
471 Assertions.assertEquals(Double.NEGATIVE_INFINITY, d0);
472 }
473
474
475 final double di = computeD(alpha, Double.POSITIVE_INFINITY);
476 if (alpha < 1) {
477
478 Assertions.assertEquals(Double.POSITIVE_INFINITY, di);
479 } else if (alpha == 1) {
480
481 Assertions.assertEquals(Double.NaN, di);
482 } else {
483
484
485 Assertions.assertEquals(Double.NaN, di);
486 }
487 }
488 }
489
490
491
492
493
494
495
496
497 private static double computeD(double alpha, double z) {
498 final double alogz = Math.log(z);
499 final double eps = 1 - alpha;
500 final double meps1 = 1 - eps;
501 return SpecialMath.d2(eps * alogz / meps1) * (alogz / meps1);
502 }
503
504
505
506
507
508
509 @Test
510 void testSinAlphaPhiMinusAtanZeta() {
511
512
513
514
515
516
517 for (double alpha : new double[] {0.25, 0.125}) {
518 for (double phi : new double[] {PI_4, PI_4 / 2}) {
519 double beta = Math.tan(-alpha * phi) / Math.tan(alpha * PI_2);
520 double zeta = -beta * Math.tan(alpha * PI_2);
521 double atanZeta = Math.atan(-zeta);
522 Assertions.assertEquals(0.0, alpha * phi + atanZeta);
523 }
524 }
525 }
526
527
528
529
530
531
532 @Test
533 void testCosPhiMinusAlphaPhiXi() {
534
535 final double cosPi2 = Math.cos(PI_2);
536
537 Assertions.assertEquals(cosPi2, Math.cos(-PI_2));
538
539 Assertions.assertTrue(cosPi2 > 0);
540
541 final UniformRandomProvider rng = RandomAssert.createRNG();
542
543
544 final double[] alphas = {1, Math.nextDown(1), 0.99, 0.5, 0.1, 0.05, 0.01, DU};
545
546
547 final long[] xs = {0, 1 << 10, Long.MIN_VALUE >>> 1, Long.MAX_VALUE};
548 for (final double alpha : alphas) {
549 for (final long x : xs) {
550 assertCosPhiMinusAlphaPhiXi(alpha, x);
551 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
552 }
553 for (int j = 0; j < 1000; j++) {
554 final long x = rng.nextLong();
555 assertCosPhiMinusAlphaPhiXi(alpha, x);
556 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
557 }
558 }
559
560 for (int i = 0; i < 1000; i++) {
561 final double alpha = rng.nextDouble();
562 for (final long x : xs) {
563 assertCosPhiMinusAlphaPhiXi(alpha, x);
564 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
565 }
566 for (int j = 0; j < 1000; j++) {
567 final long x = rng.nextLong();
568 assertCosPhiMinusAlphaPhiXi(alpha, x);
569 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
570 }
571 }
572
573
574 for (int i = 0; i <= 1023; i++) {
575 final double alpha = (double) i / 1023;
576 for (final long x : xs) {
577 assertCosPhiMinusAlphaPhiXi(alpha, x);
578 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
579 }
580 }
581 }
582
583
584
585
586
587
588
589
590
591 private static void assertCosPhiMinusAlphaPhiXi(double alpha, long x) {
592
593 final double eps = 1 - alpha;
594 final double meps1 = 1 - eps;
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618 final double alphaPi2;
619 if (meps1 > 1) {
620
621 alphaPi2 = -(2 - meps1) * PI_2;
622 } else {
623 alphaPi2 = meps1 * PI_2;
624 }
625
626
627
628 double phi = getU(x) * 2;
629 double value = eps * phi + alphaPi2;
630 Assertions.assertTrue(value <= PI_2);
631 Assertions.assertTrue(value >= -PI_2);
632 value = eps * phi - alphaPi2;
633 Assertions.assertTrue(value <= PI_2);
634 Assertions.assertTrue(value >= -PI_2);
635
636
637 phi = -phi;
638 value = eps * phi + alphaPi2;
639 Assertions.assertTrue(value <= PI_2);
640 Assertions.assertTrue(value >= -PI_2);
641 value = eps * phi - alphaPi2;
642 Assertions.assertTrue(value <= PI_2);
643 Assertions.assertTrue(value >= -PI_2);
644 }
645
646
647
648
649
650 @Test
651 void testSinAlphaPhi() {
652
653
654 for (final double phi : new double[] {getU(-1) * 2, getU(1 << 10) * 2}) {
655 final double x = Math.sin(SMALLEST_ALPHA * phi);
656 Assertions.assertNotEquals(0.0, x);
657
658 Assertions.assertEquals(1.9361559566769725E-32, Math.abs(x));
659 }
660 }
661
662
663
664
665
666
667
668 @Test
669 void testExpM1() {
670
671 Assertions.assertEquals(d2(0.5), d2b(0.5));
672
673 Assertions.assertTrue(d2(Math.nextDown(0.5)) <= d2b(0.5));
674 Assertions.assertEquals(d2(-0.5), d2b(-0.5));
675
676 Assertions.assertTrue(d2(-Math.nextDown(0.5)) >= d2b(-0.5));
677
678
679 Assertions.assertFalse(d2(Math.nextDown(0.25)) <= d2b(0.25));
680 }
681
682
683
684
685
686
687
688
689
690
691 void expm1ULPReport() {
692
693 final UniformRandomProvider rng = RandomAssert.createRNG();
694
695 final int size = 1 << 30;
696
697 final long mask = (1L << 52) - 1;
698
699
700
701 Assertions.assertEquals(((double) (1L << 54)) - 1, (double) (1L << 54));
702
703
704
705 for (int signedExp = 6; signedExp >= -4; signedExp--) {
706
707 final long exp = (signedExp + 1023L) << 52;
708
709 Assertions.assertEquals(signedExp, Math.getExponent(Double.longBitsToDouble(exp)));
710 Assertions.assertEquals(signedExp, Math.getExponent(Double.longBitsToDouble((-1 & mask) | exp)));
711
712 long sum1 = 0;
713 long sum2 = 0;
714 long max1 = 0;
715 long max2 = 0;
716 for (int i = size; i-- > 0;) {
717 final long bits = rng.nextLong() & mask;
718 final double x = Double.longBitsToDouble(bits | exp);
719 final double x1 = d2(x);
720 final double x2 = d2b(x);
721 final double x1b = d2(-x);
722 final double x2b = d2b(-x);
723 final long ulp1 = Math.abs(Double.doubleToRawLongBits(x1) - Double.doubleToRawLongBits(x2));
724 final long ulp2 = Math.abs(Double.doubleToRawLongBits(x1b) - Double.doubleToRawLongBits(x2b));
725 sum1 += ulp1;
726 sum2 += ulp2;
727 if (max1 < ulp1) {
728 max1 = ulp1;
729 }
730 if (max2 < ulp2) {
731 max2 = ulp2;
732 }
733 }
734
735 System.out.printf("%-6s %2d %-24s (%d) %-24s (%d)%n",
736 Double.longBitsToDouble(exp), signedExp,
737 (double) sum1 / size, max1, (double) sum2 / size, max2);
738
739 }
740 }
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759 private static double d2(double x) {
760
761 final double d2 = Math.expm1(x) / x;
762 if (Double.isNaN(d2)) {
763
764 if (x == 0) {
765 return 1.0;
766 }
767
768 return x;
769 }
770 return d2;
771 }
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790 private static double d2b(double x) {
791
792 final double d2 = (Math.exp(x) - 1) / x;
793 if (Double.isNaN(d2)) {
794
795 if (x == 0) {
796 return 1.0;
797 }
798
799 return x;
800 }
801 return d2;
802 }
803
804
805
806
807
808 @Test
809 void testD2() {
810 for (final double x : new double[] {Double.MAX_VALUE, Math.log(Double.MAX_VALUE), 10, 5, 1, 0.5, 0.1, 0.05, 0.01}) {
811 Assertions.assertEquals(Math.expm1(x) / x, SpecialMath.d2(x), 1e-15);
812 Assertions.assertEquals(Math.expm1(-x) / -x, SpecialMath.d2(-x), 1e-15);
813 }
814
815
816 Assertions.assertEquals(0.0, Math.expm1(Double.NEGATIVE_INFINITY) / Double.NEGATIVE_INFINITY);
817 Assertions.assertEquals(0.0, SpecialMath.d2(Double.NEGATIVE_INFINITY));
818
819
820 Assertions.assertEquals(Double.NaN, SpecialMath.d2(Double.NaN));
821
822
823 Assertions.assertEquals(Double.NaN, Math.expm1(0) / 0.0);
824 Assertions.assertEquals(Double.NaN, Math.expm1(Double.POSITIVE_INFINITY) / Double.POSITIVE_INFINITY);
825
826 Assertions.assertEquals(1.0, SpecialMath.d2(0.0));
827 Assertions.assertEquals(Double.POSITIVE_INFINITY, SpecialMath.d2(Double.POSITIVE_INFINITY));
828 }
829
830
831
832
833 @Test
834 void testTan2() {
835
836 for (final long x : new long[] {Long.MIN_VALUE + (1 << 10), Long.MAX_VALUE}) {
837 final double phiby2 = getU(x);
838 Assertions.assertEquals(PI_4 - DU * PI_4, Math.abs(phiby2));
839 final double a = phiby2 * SpecialMath.tan2(phiby2);
840
841 Assertions.assertNotEquals(1, Math.abs(a));
842 Assertions.assertTrue(Math.abs(a) < 1.0);
843 }
844
845
846 final double pi = Math.PI;
847 for (final double x : new double[] {pi, pi / 2, pi / 3.99, pi / 4, pi / 4.01, pi / 8, pi / 16}) {
848 final double y = Math.tan(x) / x;
849 Assertions.assertEquals(y, SpecialMath.tan2(x), Math.ulp(y));
850 }
851
852
853
854
855
856 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create(0x1647816481684L);
857 int count = 0;
858 long ulp = 0;
859 long max = 0;
860 long ulp2 = 0;
861 long max2 = 0;
862 for (int i = 0; i < 1000; i++) {
863 final double x = rng.nextDouble() * PI_4;
864 count++;
865 final double tanx = Math.tan(x);
866 final double tan2x = SpecialMath.tan2(x);
867
868 double y = x * tan2x;
869 if (y != tanx) {
870 final long u = Math.abs(Double.doubleToRawLongBits(tanx) - Double.doubleToRawLongBits(y));
871 if (max < u) {
872 max = u;
873 }
874 ulp += u;
875
876
877 Assertions.assertEquals(tanx, y, 4 * Math.ulp(tanx));
878 }
879
880 y = tanx / x;
881 if (y != tan2x) {
882 final long u = Math.abs(Double.doubleToRawLongBits(tan2x) - Double.doubleToRawLongBits(y));
883 if (max2 < u) {
884 max2 = u;
885 }
886 ulp2 += u;
887
888 Assertions.assertEquals(y, tan2x, 3 * Math.ulp(y));
889 }
890 }
891
892
893
894
895
896 Assertions.assertTrue((double) ulp / count < 0.6, "Mean ULP to tan(x) is too high");
897 Assertions.assertTrue((double) ulp2 / count < 0.45, "Mean ULP to tan(x) / x is too high");
898
899 Assertions.assertEquals(1.0, SpecialMath.tan2(0.0), "Must be exact tan(x) / x at x=0");
900 Assertions.assertEquals(4 / Math.PI, SpecialMath.tan2(PI_4), Math.ulp(4 / Math.PI));
901 Assertions.assertEquals(1.0, PI_4 * SpecialMath.tan2(PI_4), Math.ulp(1.0));
902
903 Assertions.assertTrue(1.0 >= PI_4 * SpecialMath.tan2(PI_4));
904 Assertions.assertTrue(1.0 >= PI_4 * SpecialMath.tan2(Math.nextDown(PI_4)));
905
906 Assertions.assertTrue(SpecialMath.tan2(Math.nextUp(PI_4)) >= SpecialMath.tan2(PI_4));
907 }
908
909
910
911
912
913
914 @Test
915 void testSamplesWithAlphaNot1() {
916
917 final double[] alphas = {0.3, 0.9, 1.1, 1.5};
918 final double[] betas = {-1, -0.5, -0.3, 0};
919 final double[] ws = {0.1, 1, 3};
920 final double[] us = {0.1, 0.25, 0.5, 0.8};
921
922 final double relative = 1e-5;
923 final double absolute = 1e-10;
924 for (final double alpha : alphas) {
925 for (final double beta : betas) {
926 for (final double w : ws) {
927 for (final double u : us) {
928 final double x = sampleCMS(alpha, beta, w, u);
929 final double y = sampleWeronAlphaNot1(alpha, beta, w, u);
930 Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative));
931
932 final double z = sampleCMS(alpha, -beta, w, 1 - u);
933 Assertions.assertEquals(x, -z, 0.0);
934 }
935 }
936 }
937 }
938 }
939
940
941
942
943
944
945 @Test
946 void testSamplesWithAlpha1() {
947
948 final double[] betas = {-1, -0.5, -0.3, 0};
949 final double[] ws = {0.1, 1, 3};
950 final double[] us = {0.1, 0.25, 0.5, 0.8};
951
952 final double relative = 1e-5;
953 final double absolute = 1e-10;
954 final double alpha = 1;
955 for (final double beta : betas) {
956 for (final double w : ws) {
957 for (final double u : us) {
958 final double x = sampleCMS(alpha, beta, w, u);
959 final double y = sampleWeronAlpha1(beta, w, u);
960 Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative));
961
962 final double z = sampleCMS(alpha, -beta, w, 1 - u);
963 Assertions.assertEquals(x, -z, 0.0);
964 }
965 }
966 }
967 }
968
969
970
971
972
973
974 @Test
975 void testConvergenceWithAlphaCloseTo1() {
976 final double[] betas = {-1, -0.5, 0, 0.3, 1};
977 final double[] ws = {0.1, 1, 10};
978 final double[] us = {0.1, 0.25, 0.5, 0.8};
979 final int steps = 30;
980
981
982
983
984
985 for (double deltaStart : new double[] {-0.0625, 0.0625}) {
986
987
988 int cmsCount = 0;
989 int weronCount = 0;
990
991 for (final double beta : betas) {
992 for (final double w : ws) {
993 for (final double u : us) {
994
995 double x0 = sampleCMS(1, beta, w, u);
996 Assertions.assertTrue(Double.isFinite(x0), "Target must be finite");
997
998
999 double delta = deltaStart;
1000 double dx = Math.abs(x0 - sampleCMS(1 + delta, beta, w, u));
1001 for (int i = 0; i < steps; i++) {
1002 delta /= 2;
1003 final double dx2 = Math.abs(x0 - sampleCMS(1 + delta, beta, w, u));
1004 if (dx2 > dx) {
1005 cmsCount++;
1006 }
1007 dx = dx2;
1008 }
1009
1010
1011 x0 = sampleWeronAlpha1(beta, w, u);
1012 Assertions.assertTrue(Double.isFinite(x0), "Target must be finite");
1013
1014
1015 delta = deltaStart;
1016 dx = Math.abs(x0 - sampleWeronAlphaNot1(1 + delta, beta, w, u));
1017 for (int i = 0; i < steps; i++) {
1018 delta /= 2;
1019 final double dx2 = Math.abs(x0 - sampleWeronAlphaNot1(1 + delta, beta, w, u));
1020 if (dx2 > dx) {
1021 weronCount++;
1022 }
1023 dx = dx2;
1024 }
1025 }
1026 }
1027 }
1028
1029
1030 Assertions.assertEquals(0, cmsCount);
1031
1032
1033 Assertions.assertTrue(weronCount > 200);
1034 }
1035 }
1036
1037
1038
1039
1040
1041
1042
1043 @Test
1044 void testExtremeInputsToSample() {
1045
1046 Assertions.assertEquals(Double.NaN, sampleCMS(1.3, 0.7, 0, 0.25));
1047 Assertions.assertTrue(Double.isFinite(sampleCMS(1.3, 0.7, SMALL_W, 0.25)));
1048
1049
1050 Assertions.assertEquals(Double.NaN, sampleCMS(1.1, 1.0, 0.1, 0));
1051 Assertions.assertTrue(Double.isFinite(sampleCMS(1.1, 1.0, 0.1, DU)));
1052
1053
1054
1055
1056 Assertions.assertEquals(Double.NaN, sampleCMS(0.01, 0.7, SMALL_W, 0.5));
1057
1058
1059
1060 Assertions.assertEquals(Double.NaN, sampleCMS(1e-5, 0.7, 1.0, 1e-4));
1061 Assertions.assertEquals(Double.NaN, sampleCMS(1e-5, -0.7, 1.0, 1 - 1e-4));
1062
1063 final double[] alphas = {Math.nextDown(2), 1.3, 1.1, Math.nextUp(1), 1, Math.nextDown(1), 0.7, 0.1, 0.05, 0.01, 0x1.0p-16};
1064 final double[] betas = {1, 0.9, 0.001, 0};
1065
1066
1067
1068 final double[] ws = {0, SMALL_W, 0.001, 1, 10, LARGE_W};
1069
1070
1071
1072 final double[] us = {DU, 2 * DU, 0.0001, 0.5 - DU, 0.5};
1073
1074 int nan1 = 0;
1075
1076 for (final double alpha : alphas) {
1077 for (final double beta : betas) {
1078 if (alpha == 1 && beta == 0) {
1079
1080 continue;
1081 }
1082
1083
1084 final double[] support = getSupport(alpha, beta);
1085 final double lower = support[0];
1086 final double upper = support[1];
1087 for (final double w : ws) {
1088 for (final double u : us) {
1089 final double x1 = sampleCMS(alpha, beta, w, u);
1090 final double x2 = sampleWeron(alpha, beta, w, u);
1091
1092 if (Double.isNaN(x1)) {
1093 nan1++;
1094 }
1095
1096 Assertions.assertNotEquals(Double.NaN, x2);
1097
1098
1099
1100 Assertions.assertEquals(x1, 0.0 - sampleCMS(alpha, -beta, w, 1 - u), 0.0);
1101 Assertions.assertEquals(x2, 0.0 - sampleWeron(alpha, -beta, w, 1 - u), 0.0);
1102
1103 if (Double.isInfinite(x1) && x1 != x2) {
1104
1105
1106
1107 Assertions.assertTrue(lower <= x2 && x2 <= upper);
1108 }
1109 }
1110 }
1111 }
1112 }
1113
1114
1115 Assertions.assertNotEquals(0, nan1);
1116 }
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130 private static double sampleCMS(double alpha, double beta, double w, double u) {
1131 final double phiby2 = PI_2 * (u - 0.5);
1132 final double eps = 1 - alpha;
1133
1134
1135 final double meps1 = 1 - eps;
1136
1137
1138 final double tau = CMSStableSampler.getTau(alpha, beta);
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151 final double a = phiby2 * SpecialMath.tan2(phiby2);
1152
1153 final double bb = SpecialMath.tan2(eps * phiby2);
1154
1155 final double b = eps * phiby2 * bb;
1156
1157 final double da = a * a;
1158 final double db = b * b;
1159
1160 final double a2 = 1 - da;
1161
1162 final double a2p = 1 + da;
1163
1164 final double b2 = 1 - db;
1165
1166 final double b2p = 1 + db;
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181 final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
1182
1183 final double alogz = Math.log(z);
1184 final double d = SpecialMath.d2(eps * alogz / meps1) * (alogz / meps1);
1185
1186
1187 return (1 + eps * d) *
1188 (2 * ((a - b) * (1 + a * b) - phiby2 * tau * bb * (b * a2 - 2 * a))) /
1189 (a2 * b2p) + tau * d;
1190 }
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206 private static double sampleWeron(double alpha, double beta, double w, double u) {
1207 return alpha == 1 ? sampleWeronAlpha1(beta, w, u) : sampleWeronAlphaNot1(alpha, beta, w, u);
1208 }
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232 private static double sampleWeronAlphaNot1(double alpha, double beta, double w, double u) {
1233
1234 final double eps = 1 - alpha;
1235 final double meps1 = 1 - eps;
1236
1237 double zeta;
1238 if (meps1 > 1) {
1239 zeta = beta * Math.tan((2 - meps1) * PI_2);
1240 } else {
1241 zeta = -beta * Math.tan(meps1 * PI_2);
1242 }
1243
1244 final double scale = Math.pow(1 + zeta * zeta, 0.5 / meps1);
1245 final double invAlpha = 1.0 / meps1;
1246 final double invAlphaM1 = invAlpha - 1;
1247
1248 final double phi = Math.PI * (u - 0.5);
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273 final double atanZeta = Math.atan(-zeta);
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290 double t1 = Math.sin(meps1 * phi + atanZeta);
1291 if (t1 == 0) {
1292 return zeta;
1293 }
1294
1295
1296 t1 /= Math.pow(Math.cos(phi), invAlpha);
1297
1298
1299
1300
1301
1302 final double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, invAlphaM1);
1303 if (t2 == 0) {
1304 return zeta;
1305 }
1306
1307 return t1 * t2 * scale + zeta;
1308 }
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325 private static double sampleWeronAlpha1(double beta, double w, double u) {
1326
1327 final double phi = Math.PI * (u - 0.5);
1328
1329
1330 final double betaPhi = PI_2 + beta * phi;
1331 return (betaPhi * Math.tan(phi) -
1332 beta * Math.log(PI_2 * w * Math.cos(phi) / betaPhi)) / PI_2;
1333 }
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343 @Test
1344 void testSamplesWithZBelow0() {
1345
1346
1347
1348
1349 final long[] longs = {Long.MAX_VALUE, -6261465550279131136L};
1350
1351 final double phiby2 = PI_4 - PI_4 * DU;
1352 final double w = 5.0;
1353 assertUWSequence(new double[] {
1354 phiby2, w,
1355 }, longs);
1356
1357
1358
1359 final double alpha = 1.291015625;
1360 final double beta = -1;
1361 Assertions.assertTrue(0.0 > computeNumerator(alpha, beta, Long.MAX_VALUE));
1362
1363
1364
1365 final double eps = 1 - alpha;
1366 final double tau = CMSStableSampler.getTau(alpha, beta);
1367 final double a = phiby2 * SpecialMath.tan2(phiby2);
1368 final double bb = SpecialMath.tan2(eps * phiby2);
1369 final double b = eps * phiby2 * bb;
1370 final double da = a * a;
1371 final double db = b * b;
1372 final double a2 = 1 - da;
1373 final double a2p = 1 + da;
1374 final double b2 = 1 - db;
1375 final double b2p = 1 + db;
1376 final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
1377 Assertions.assertTrue(0.0 > z);
1378
1379 final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
1380
1381 Assertions.assertTrue(Double.isFinite(sampler.sample()), "Sampler did not recover");
1382 }
1383
1384
1385
1386
1387
1388 @Test
1389 void testSamplesWithZInfinite() {
1390
1391 final long[] longs = {Long.MIN_VALUE >>> 1, 0};
1392
1393 assertUWSequence(new double[] {
1394 PI_4 / 2, 0,
1395 }, longs);
1396
1397 for (final double alpha : new double[] {0.789, 1, 1.23}) {
1398
1399 for (final double beta : new double[] {-0.56, 0, 0.56}) {
1400
1401 if (alpha == 1 && beta == 0) {
1402 continue;
1403 }
1404 final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
1405 final double x = sampler.sample();
1406
1407 Assertions.assertFalse(Double.isNaN(x), "Sampler did not recover");
1408 if (beta != 0) {
1409
1410 if (alpha < 0) {
1411
1412 Assertions.assertEquals(Math.copySign(Double.POSITIVE_INFINITY, beta), x);
1413 } else if (alpha > 1) {
1414
1415 final double[] support = getSupport(alpha, beta);
1416 final double mu = support[2];
1417 Assertions.assertEquals(mu, x);
1418 }
1419 }
1420 }
1421 }
1422 }
1423
1424
1425
1426
1427
1428 @Test
1429 void testSamplesWithDInfinite() {
1430
1431
1432 testSamplesWithDInfinite(0.01);
1433 testSamplesWithDInfinite(-0.01);
1434 }
1435
1436
1437
1438
1439
1440 @Test
1441 void testBeta0SamplesWithDInfinite() {
1442 testSamplesWithDInfinite(0.0);
1443 }
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453 private static void testSamplesWithDInfinite(double beta) {
1454
1455
1456 final long xuLo = Long.MIN_VALUE + (1024 << 10);
1457 final long xuHi = Long.MAX_VALUE - (1023 << 10);
1458
1459
1460 final long x = 3L;
1461 final long[] longs = {xuLo, x, xuHi, x, 0, x};
1462
1463 assertUWSequence(new double[] {
1464 -PI_4 + 1024 * DU * PI_4, SMALL_W,
1465 PI_4 - 1024 * DU * PI_4, SMALL_W,
1466 0.0, SMALL_W
1467 }, longs);
1468
1469
1470
1471
1472
1473
1474 final double alpha = 0.03;
1475 final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
1476 final double x1 = sampler.sample();
1477 final double x2 = sampler.sample();
1478 final double x3 = sampler.sample();
1479
1480 final double max = Double.POSITIVE_INFINITY;
1481 Assertions.assertEquals(-max, x1);
1482 Assertions.assertEquals(max, x2);
1483
1484 Assertions.assertNotEquals(Double.NaN, x3);
1485
1486 if (beta == 0) {
1487
1488 Assertions.assertEquals(0.0, x3);
1489 } else {
1490
1491 Assertions.assertEquals(Math.copySign(max, beta), x3);
1492 }
1493 }
1494
1495
1496
1497
1498
1499 @Test
1500 void testAlpha1SamplesWithExtremePhi() {
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510 final long[] longs1 = {Long.MIN_VALUE + (1 << 10), 2703662416942444033L};
1511 assertUWSequence(new double[] {
1512 -PI_4 + PI_4 * DU, 1.0,
1513 }, longs1);
1514 final StableSampler sampler1 = StableSampler.of(createRngWithSequence(longs1), 1.0, 1.0);
1515 final double x1 = sampler1.sample();
1516 Assertions.assertTrue(Double.isFinite(x1), "Sampler did not recover");
1517
1518
1519 final long[] longs2 = {Long.MAX_VALUE, 2703662416942444033L};
1520 assertUWSequence(new double[] {
1521 PI_4 - PI_4 * DU, 1.0,
1522 }, longs2);
1523 final StableSampler sampler2 = StableSampler.of(createRngWithSequence(longs2), 1.0, -1.0);
1524 final double x2 = sampler2.sample();
1525 Assertions.assertTrue(Double.isFinite(x2), "Sampler did not recover");
1526
1527
1528 Assertions.assertEquals(x1, -x2);
1529 }
1530
1531
1532
1533
1534
1535
1536 @Test
1537 void testSupport() {
1538 testSupport(1.0, 0.0);
1539 }
1540
1541
1542
1543
1544
1545
1546 @Test
1547 void testSupportWithTransformation() {
1548
1549 for (final double gamma : new double[] {0.78, 1.23, Double.MAX_VALUE, Double.MIN_VALUE}) {
1550 for (final double delta : new double[] {0.43, 12.34, Double.MAX_VALUE}) {
1551 testSupport(gamma, delta);
1552 testSupport(gamma, -delta);
1553 }
1554 }
1555 }
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566 private static void testSupport(double gamma, double delta) {
1567
1568 final double[] alphas = {2.0, 1.5, 1.0, Math.nextDown(1), 0.99, 0.75, 0.5, 0.25, 0.1, 0.01};
1569 for (final double alpha : alphas) {
1570 testSupport(alpha, 1, gamma, delta);
1571 testSupport(alpha, -1, gamma, delta);
1572 }
1573 }
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583 private static void testSupport(double alpha, double beta, double gamma, double delta) {
1584
1585 final double[] support = getSupport(alpha, beta);
1586
1587 double lower;
1588 if (support[0] == -Double.MAX_VALUE) {
1589 lower = Double.NEGATIVE_INFINITY;
1590 } else {
1591 lower = support[0] * gamma + delta;
1592 }
1593 double upper;
1594 if (support[1] == Double.MAX_VALUE) {
1595 upper = Double.POSITIVE_INFINITY;
1596 } else {
1597 upper = support[1] * gamma + delta;
1598 }
1599
1600
1601
1602 final long[] longs = new long[] {
1603
1604
1605 Long.MIN_VALUE,
1606
1607
1608 Long.MIN_VALUE + (1 << 10), 0,
1609
1610 Long.MIN_VALUE + (1 << 10), -1, -1, -1, -1, -1, -1, -1, -1, 0,
1611
1612 Long.MAX_VALUE, 0,
1613
1614 Long.MAX_VALUE, -1, -1, -1, -1, -1, -1, -1, -1, 0,
1615
1616 0, 0,
1617
1618 0, -1, -1, -1, -1, -1, -1, -1, -1, 0,
1619
1620
1621
1622 Long.MIN_VALUE + (1 << 10), 2703662416942444033L,
1623
1624 Long.MAX_VALUE, 2703662416942444033L,
1625
1626 0, 2703662416942444033L,
1627
1628
1629
1630 Long.MIN_VALUE >> 1, 0,
1631
1632 Long.MIN_VALUE >> 1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
1633
1634 Long.MIN_VALUE >>> 1, 0,
1635
1636 Long.MIN_VALUE >>> 1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
1637 };
1638
1639
1640 final double phiby2low = -PI_4 + PI_4 * DU;
1641 final double phiby2high = PI_4 - PI_4 * DU;
1642 assertUWSequence(new double[] {
1643 phiby2low, 0,
1644 phiby2low, LARGE_W,
1645 phiby2high, 0,
1646 phiby2high, LARGE_W,
1647 0, 0,
1648 0, LARGE_W,
1649 phiby2low, 1.0,
1650 phiby2high, 1.0,
1651 0, 1.0,
1652 -PI_4 / 2, 0,
1653 -PI_4 / 2, LARGE_W,
1654 PI_4 / 2, 0,
1655 PI_4 / 2, LARGE_W,
1656 }, longs);
1657
1658 final StableSampler sampler = StableSampler.of(
1659 createRngWithSequence(longs), alpha, beta, gamma, delta);
1660 for (int i = 0; i < 100; i++) {
1661 final double x = sampler.sample();
1662 if (!(lower <= x && x <= upper)) {
1663 Assertions.fail(String.format("Invalid sample. alpha=%s, beta=%s, gamma=%s, delta=%s [%s, %s] x=%s",
1664 alpha, beta, gamma, delta, lower, upper, x));
1665 }
1666 }
1667 }
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684 private static double[] getSupport(double alpha, double beta) {
1685
1686 double eps = 1 - alpha;
1687 double meps1 = 1 - eps;
1688
1689
1690
1691 double mu;
1692 if (alpha > 1) {
1693 mu = beta * Math.tan((2 - meps1) * PI_2);
1694 } else {
1695
1696
1697 if (alpha == 0.5) {
1698 mu = -beta;
1699 } else {
1700 mu = -beta * Math.tan(meps1 * PI_2);
1701 }
1702 }
1703
1704
1705 double lower = -Double.MAX_VALUE;
1706 double upper = Double.MAX_VALUE;
1707 if (meps1 < 1) {
1708 if (beta == 1) {
1709
1710 lower = mu;
1711 } else if (beta == -1) {
1712
1713 upper = mu;
1714 }
1715 }
1716 return new double[] {lower, upper, mu};
1717 }
1718
1719
1720
1721
1722
1723 @Test
1724 void testRandomDeviatesUandW() {
1725
1726 final double d = DU * PI_4;
1727
1728 Assertions.assertNotEquals(-PI_4, getU(createRngWithSequence(Long.MIN_VALUE)));
1729 Assertions.assertEquals(-PI_4 + d, getU(createRngWithSequence(Long.MIN_VALUE + (1 << 10))));
1730 Assertions.assertEquals(-PI_4 / 2, getU(createRngWithSequence(Long.MIN_VALUE >> 1)));
1731 Assertions.assertEquals(-d, getU(createRngWithSequence(-1)));
1732 Assertions.assertEquals(0.0, getU(createRngWithSequence(0)));
1733 Assertions.assertEquals(d, getU(createRngWithSequence(1 << 10)));
1734 Assertions.assertEquals(PI_4 / 2, getU(createRngWithSequence(Long.MIN_VALUE >>> 1)));
1735 Assertions.assertEquals(PI_4 - d, getU(createRngWithSequence(Long.MAX_VALUE)));
1736
1737
1738 Assertions.assertEquals(0, ZigguratSampler.Exponential.of(
1739 createRngWithSequence(0L)).sample());
1740 Assertions.assertEquals(SMALL_W, ZigguratSampler.Exponential.of(
1741 createRngWithSequence(3)).sample());
1742 Assertions.assertEquals(0.5, ZigguratSampler.Exponential.of(
1743 createRngWithSequence(1446480648965178882L)).sample());
1744 Assertions.assertEquals(1.0, ZigguratSampler.Exponential.of(
1745 createRngWithSequence(2703662416942444033L)).sample());
1746 Assertions.assertEquals(2.5, ZigguratSampler.Exponential.of(
1747 createRngWithSequence(6092639261715210240L)).sample());
1748 Assertions.assertEquals(5.0, ZigguratSampler.Exponential.of(
1749 createRngWithSequence(-6261465550279131136L)).sample());
1750 Assertions.assertEquals(TAIL_W, ZigguratSampler.Exponential.of(
1751 createRngWithSequence(-1, -1, 0)).sample());
1752 Assertions.assertEquals(3 * TAIL_W, ZigguratSampler.Exponential.of(
1753 createRngWithSequence(-1, -1, -1, -1, -1, -1, 0)).sample(), 1e-14);
1754 }
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765 private static double getU(UniformRandomProvider rng) {
1766 final double x = getU(rng.nextLong());
1767 if (x == -PI_4) {
1768 return getU(rng);
1769 }
1770 return x;
1771 }
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796 private static double getU(long x) {
1797 return (x >> 10) * PI_4_SCALED;
1798 }
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845 private static UniformRandomProvider createRngWithSequence(final long... longs) {
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873 final UniformRandomProvider rng = RandomAssert.seededRNG();
1874
1875
1876 return new SplitMix64(0L) {
1877 private int l;
1878
1879 @Override
1880 public long nextLong() {
1881 if (l == longs.length) {
1882 return rng.nextLong();
1883 }
1884 return longs[l++];
1885 }
1886 };
1887 }
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900 private static void assertUWSequence(double[] expected, long[] longs) {
1901 final UniformRandomProvider rng = createRngWithSequence(longs);
1902
1903
1904 final SharedStateContinuousSampler exp = ZigguratSampler.Exponential.of(rng);
1905 for (int i = 0; i < expected.length; i += 2) {
1906 final int j = i / 2;
1907 Assertions.assertEquals(expected[i], getU(rng), () -> j + ": Incorrect u");
1908 if (i + 1 < expected.length) {
1909 Assertions.assertEquals(expected[i + 1], exp.sample(), () -> j + ": Incorrect w");
1910 }
1911 }
1912 }
1913
1914
1915
1916
1917
1918
1919 @Test
1920 void testSamplerOutputIsContinuousFunction() {
1921
1922 for (final double beta : new double[] {0.5, 0.2, 0.1, 0.001}) {
1923 testSamplerOutputIsContinuousFunction(1 + 8096 * DU, beta, 1.0, beta, 1 - 8096 * DU, beta, 0);
1924 testSamplerOutputIsContinuousFunction(1 + 1024 * DU, beta, 1.0, beta, 1 - 1024 * DU, beta, 1);
1925
1926 testSamplerOutputIsContinuousFunction(1 + 128 * DU, beta, 1.0, beta, 1 - 128 * DU, beta, 1);
1927 testSamplerOutputIsContinuousFunction(1 + 16 * DU, beta, 1.0, beta, 1 - 16 * DU, beta, 4);
1928
1929
1930 testSamplerOutputIsContinuousFunction(1 + DU, beta, 1.0, beta, 1 - DU, beta, 0);
1931 }
1932
1933 for (final double alpha : new double[] {1.5, 1.2, 1.1, 1.001}) {
1934 testSamplerOutputIsContinuousFunction(alpha, 8096 * DU, alpha, 0, alpha, -8096 * DU, 0);
1935 testSamplerOutputIsContinuousFunction(alpha, 1024 * DU, alpha, 0, alpha, -1024 * DU, 0);
1936 testSamplerOutputIsContinuousFunction(alpha, 128 * DU, alpha, 0, alpha, -128 * DU, 1);
1937
1938 testSamplerOutputIsContinuousFunction(alpha, 16 * DU, alpha, 0, alpha, -16 * DU, 64);
1939 testSamplerOutputIsContinuousFunction(alpha, DU, alpha, 0, alpha, -DU, 4);
1940 }
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950 }
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964 private static void testSamplerOutputIsContinuousFunction(double alpha1, double beta1,
1965 double alpha2, double beta2,
1966 double alpha3, double beta3,
1967 int ulp) {
1968 final long seed = 0x62738468L;
1969 final UniformRandomProvider rng1 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
1970 final UniformRandomProvider rng2 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
1971 final UniformRandomProvider rng3 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
1972 final StableSampler sampler1 = StableSampler.of(rng1, alpha1, beta1);
1973 final StableSampler sampler2 = StableSampler.of(rng2, alpha2, beta2);
1974 final StableSampler sampler3 = StableSampler.of(rng3, alpha3, beta3);
1975 final Supplier<String> msg = () -> String.format("alpha=%s, beta=%s", alpha2, beta2);
1976 for (int i = 0; i < 1000; i++) {
1977 final double x1 = sampler1.sample();
1978 final double x2 = sampler2.sample();
1979 final double x3 = sampler3.sample();
1980
1981 if (x3 > x1) {
1982 if (x2 > x3) {
1983
1984 Assertions.assertEquals(x3, x2, ulp * Math.ulp(x3), msg);
1985 } else if (x2 < x1) {
1986 Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1), msg);
1987 }
1988 } else if (x3 < x1) {
1989 if (x2 < x3) {
1990
1991 Assertions.assertEquals(x3, x2, ulp * Math.ulp(x3), msg);
1992 } else if (x2 > x1) {
1993 Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1), msg);
1994 }
1995 }
1996 }
1997 }
1998
1999
2000
2001
2002 @Test
2003 void testSharedStateSampler() {
2004
2005 testSharedStateSampler(2.0, 0.0);
2006
2007 testSharedStateSampler(1.0, 0.0);
2008
2009 testSharedStateSampler(0.5, 1.0);
2010
2011 testSharedStateSampler(0.5, 0.1);
2012
2013 testSharedStateSampler(1.3, 0.0);
2014
2015 testSharedStateSampler(1.0, 0.23);
2016
2017 testSharedStateSampler(Math.nextUp(1.0), 0.23);
2018
2019 testSharedStateSampler(1.3, 0.1);
2020
2021 testSharedStateSampler(1e-5, 0.1);
2022 testSharedStateSampler(1e-5, 0.0);
2023
2024
2025 testSharedStateSampler(1.99, 0.1);
2026 }
2027
2028
2029
2030
2031
2032
2033
2034
2035 private static void testSharedStateSampler(double alpha, double beta) {
2036 final UniformRandomProvider rng1 = RandomAssert.seededRNG();
2037 final UniformRandomProvider rng2 = RandomAssert.seededRNG();
2038 StableSampler sampler1 = StableSampler.of(rng1, alpha, beta);
2039 StableSampler sampler2 = sampler1.withUniformRandomProvider(rng2);
2040 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
2041
2042 sampler1 = StableSampler.of(rng1, alpha, beta, 1.3, 13.2);
2043 sampler2 = sampler1.withUniformRandomProvider(rng2);
2044 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
2045 }
2046
2047
2048
2049
2050 @Test
2051 void testTransformedSampler() {
2052
2053
2054
2055
2056 testTransformedSampler(2.0, 0.0, 1);
2057
2058 testTransformedSampler(1.0, 0.0);
2059
2060 testTransformedSampler(0.5, 1.0);
2061
2062 testTransformedSampler(1.3, 0.0);
2063
2064 testTransformedSampler(1.0, 0.23);
2065
2066 testTransformedSampler(Math.nextUp(1.0), 0.23);
2067
2068 testTransformedSampler(1.3, 0.1);
2069
2070 testTransformedSampler(1e-5, 0.1);
2071
2072
2073 testTransformedSampler(1.99, 0.1);
2074 }
2075
2076
2077
2078
2079
2080
2081
2082
2083 private static void testTransformedSampler(double alpha, double beta) {
2084 testTransformedSampler(alpha, beta, 0);
2085 }
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095 private static void testTransformedSampler(double alpha, double beta, int ulp) {
2096 final UniformRandomProvider[] rngs = RandomAssert.createRNG(2);
2097 final UniformRandomProvider rng1 = rngs[0];
2098 final UniformRandomProvider rng2 = rngs[1];
2099 final double gamma = 3.4;
2100 final double delta = -17.3;
2101 final StableSampler sampler1 = StableSampler.of(rng1, alpha, beta, gamma, delta);
2102 final ContinuousSampler sampler2 = createTransformedSampler(rng2, alpha, beta, gamma, delta);
2103 if (ulp == 0) {
2104 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
2105 } else {
2106 for (int i = 0; i < 10; i++) {
2107 final double x1 = sampler1.sample();
2108 final double x2 = sampler2.sample();
2109 Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1));
2110 }
2111 }
2112 }
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125 private static ContinuousSampler createTransformedSampler(UniformRandomProvider rng,
2126 double alpha, double beta,
2127 final double gamma, final double delta) {
2128 final StableSampler delegate = StableSampler.of(rng, alpha, beta);
2129 return new ContinuousSampler() {
2130 @Override
2131 public double sample() {
2132 return gamma * delegate.sample() + delta;
2133 }
2134 };
2135 }
2136
2137
2138
2139
2140 @Test
2141 void testSymmetry() {
2142 final byte[] seed = RandomSource.KISS.createSeed();
2143 for (final double alpha : new double[] {1e-4, 0.78, 1, 1.23}) {
2144 for (final double beta : new double[] {-0.43, 0.23}) {
2145 for (final double gamma : new double[] {0.78, 1, 1.23}) {
2146 for (final double delta : new double[] {-0.43, 0, 0.23}) {
2147
2148
2149
2150
2151
2152
2153
2154
2155 final UniformRandomProvider rng1 = RandomSource.KISS.create(seed);
2156 final UniformRandomProvider rng2 = RandomSource.KISS.create(seed);
2157
2158
2159 final UniformRandomProvider forward = new SplitMix64(0) {
2160 private int i;
2161 @Override
2162 public long nextLong() {
2163
2164 if ((i++ & 0x1) == 0) {
2165
2166
2167 final long x = rng1.nextLong() >>> 10 | 1L;
2168
2169 return x << 10;
2170 }
2171
2172 long x;
2173 do {
2174 x = rng1.nextLong();
2175 } while ((x & 0xff) >= 252);
2176 return x;
2177 }
2178 };
2179
2180
2181 final UniformRandomProvider reverse = new SplitMix64(0) {
2182 private final long upper = 1L << 54;
2183 private int i;
2184 @Override
2185 public long nextLong() {
2186
2187 if ((i++ & 0x1) == 0) {
2188
2189
2190 final long x = rng2.nextLong() >>> 10 | 1L;
2191
2192 return (upper - x) << 10;
2193 }
2194
2195 long x;
2196 do {
2197 x = rng2.nextLong();
2198 } while ((x & 0xff) >= 252);
2199 return x;
2200 }
2201 };
2202
2203 final StableSampler s1 = StableSampler.of(forward, alpha, beta, gamma, delta);
2204
2205 final StableSampler s2 = StableSampler.of(reverse, alpha, -beta, gamma, -delta);
2206 for (int i = 0; i < 100; i++) {
2207 Assertions.assertEquals(s1.sample(), -s2.sample());
2208 }
2209 }
2210 }
2211 }
2212 }
2213 }
2214
2215
2216
2217
2218 @Test
2219 void testSymmetryLevy() {
2220 final double alpha = 0.5;
2221 final double beta = 1.0;
2222 final byte[] seed = RandomSource.KISS.createSeed();
2223 final UniformRandomProvider rng1 = RandomSource.KISS.create(seed);
2224 final UniformRandomProvider rng2 = RandomSource.KISS.create(seed);
2225 for (final double gamma : new double[] {0.78, 1, 1.23}) {
2226 for (final double delta : new double[] {-0.43, 0, 0.23}) {
2227 final StableSampler s1 = StableSampler.of(rng1, alpha, beta, gamma, delta);
2228
2229 final StableSampler s2 = StableSampler.of(rng2, alpha, -beta, gamma, -delta);
2230 for (int i = 0; i < 100; i++) {
2231 Assertions.assertEquals(s1.sample(), -s2.sample());
2232 }
2233 }
2234 }
2235 }
2236
2237
2238
2239
2240
2241
2242 @Test
2243 void testToString() {
2244 final UniformRandomProvider rng = RandomAssert.seededRNG();
2245 for (final double[] p : new double[][] {
2246 {1.3, 0.1},
2247 {2.0, 0.0},
2248 {1.0, 0.0},
2249 {0.5, 1.0},
2250 {1e-5, 0},
2251 {1e-5, 0.1},
2252 {0.7, 0.1, 3.0, 4.5},
2253 }) {
2254 StableSampler sampler;
2255 if (p.length == 2) {
2256 sampler = StableSampler.of(rng, p[0], p[1]);
2257 } else {
2258 sampler = StableSampler.of(rng, p[0], p[1], p[2], p[3]);
2259 }
2260 final String s = sampler.toString().toLowerCase();
2261 Assertions.assertTrue(s.contains("stable"));
2262 }
2263 }
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275 @Test
2276 void testImplementationsMatch() {
2277
2278
2279
2280
2281 final long unsetHighBit = ~(1L << 54);
2282 final long setLowBit = 1L << 53;
2283 final double hi = getU(Long.MAX_VALUE & unsetHighBit);
2284 final double lo = getU(Long.MIN_VALUE | setLowBit);
2285
2286 Assertions.assertEquals(PI_4, hi, 2e-3);
2287 Assertions.assertEquals(-PI_4, lo, 2e-3);
2288 Assertions.assertEquals(0.0, lo + hi, 1e-3);
2289
2290
2291 final UniformRandomProvider rng = createRngWithSequence(setLowBit);
2292 final double w = ZigguratSampler.Exponential.of(rng).sample();
2293 Assertions.assertNotEquals(0.0, w);
2294
2295 Assertions.assertEquals(0.0036959349092519837, w);
2296
2297 final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
2298 final long seed = 0x83762b3daf1c43L;
2299 final UniformRandomProvider rng1 = new SplitMix64(0L) {
2300 private UniformRandomProvider delegate = source.create(seed);
2301 @Override
2302 public long next() {
2303 final long x = delegate.nextLong();
2304 return (x & unsetHighBit) | setLowBit;
2305 }
2306 };
2307 final UniformRandomProvider rng2 = new SplitMix64(0L) {
2308 private UniformRandomProvider delegate = source.create(seed);
2309 @Override
2310 public long next() {
2311 final long x = delegate.nextLong();
2312 return (x & unsetHighBit) | setLowBit;
2313 }
2314 };
2315
2316
2317 final double[] alphas = {0.3, 0.5, 1.2, 1.5};
2318 final double[] betas = {-0.5, -0.3, -0.1, 0};
2319
2320 final double relative = 1e-5;
2321 final double absolute = 1e-10;
2322
2323 for (final double alpha : alphas) {
2324 for (final double beta : betas) {
2325 final Supplier<String> msg = () -> String.format("alpha=%s, beta=%s", alpha, beta);
2326
2327
2328
2329 StableSampler s1;
2330 StableSampler s2;
2331 if (beta == 0) {
2332 s1 = new Beta0CMSStableSampler(rng1, alpha);
2333 s2 = new Beta0WeronStableSampler(rng2, alpha);
2334 } else {
2335 s1 = new CMSStableSampler(rng1, alpha, beta);
2336 s2 = new WeronStableSampler(rng2, alpha, beta);
2337 }
2338 for (int i = 0; i < 1000; i++) {
2339 final double x = s1.sample();
2340 final double y = s2.sample();
2341 Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative), msg);
2342 }
2343 }
2344 }
2345 }
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355 @Test
2356 void testSpecializedBeta0CMSImplementation() {
2357 final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
2358
2359 final byte[] seed = source.createSeed();
2360 final UniformRandomProvider rng1 = source.create(seed);
2361 final UniformRandomProvider rng2 = source.create(seed);
2362
2363 final double[] alphas = {0.3, 0.5, 1.2, 1.5};
2364 for (final double alpha : alphas) {
2365
2366
2367
2368 final StableSampler sampler1 = new CMSStableSampler(rng1, alpha, 0.0);
2369 final StableSampler sampler2 = new Beta0CMSStableSampler(rng2, alpha);
2370 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
2371 }
2372 }
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382 @Test
2383 void testSpecializedBeta0WeronImplementation() {
2384 final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
2385
2386 final byte[] seed = source.createSeed();
2387 final UniformRandomProvider rng1 = source.create(seed);
2388 final UniformRandomProvider rng2 = source.create(seed);
2389
2390 final double[] alphas = {0.3, 0.5, 1.2, 1.5};
2391 for (final double alpha : alphas) {
2392
2393
2394
2395 final StableSampler sampler1 = new WeronStableSampler(rng1, alpha, 0.0);
2396 final StableSampler sampler2 = new Beta0WeronStableSampler(rng2, alpha);
2397 RandomAssert.assertProduceSameSequence(sampler1, sampler2);
2398 }
2399 }
2400
2401
2402
2403
2404
2405
2406
2407 @Test
2408 void testWeronImplementationEdgeCase() {
2409 double alpha = 0.25;
2410
2411 double beta = -0.48021693505171;
2412
2413
2414 final long x = Long.MIN_VALUE >>> 1;
2415 final long[] longs = new long[] {
2416
2417 x, 0,
2418
2419 x, -1, -1, -1, -1, -1, -1, -1, -1, 0,
2420
2421 x, 2703662416942444033L,
2422 };
2423
2424
2425 assertUWSequence(new double[] {
2426 PI_4 / 2, 0,
2427 PI_4 / 2, LARGE_W,
2428 PI_4 / 2, 1.0,
2429 }, longs);
2430
2431 final double zeta = -beta * Math.tan(alpha * PI_2);
2432 Assertions.assertEquals(0.0, alpha * PI_4 + Math.atan(-zeta));
2433
2434 final UniformRandomProvider rng = createRngWithSequence(longs);
2435 final StableSampler sampler = new WeronStableSampler(rng, alpha, beta);
2436
2437
2438 Assertions.assertEquals(zeta, sampler.sample());
2439 Assertions.assertEquals(zeta, sampler.sample());
2440 Assertions.assertEquals(zeta, sampler.sample());
2441 }
2442 }