✨ Add exponentialWad to LibPRNG (ethereum#836)

Author: Vectorized

.gas-snapshot

@@ -649,16 +649,17 @@ LibMapTest:testUint8MapSetAndGet() (gas: 883454)
 LibMapTest:testUint8MapSetAndGet(uint256) (runs: 256, μ: 59152, ~: 59103)
 LibMapTest:testUint8MapSetAndGet2(uint256) (runs: 256, μ: 67845, ~: 69101)
 LibMapTest:test__codesize() (gas: 13187)
-LibPRNGTest:testLCGGas() (gas: 20802)
+LibPRNGTest:testExponentialWad() (gas: 4268529)
+LibPRNGTest:testLCGGas() (gas: 20736)
 LibPRNGTest:testPRNGGas() (gas: 25645)
-LibPRNGTest:testPRNGNext() (gas: 16162)
-LibPRNGTest:testPRNGShuffle() (gas: 504485)
+LibPRNGTest:testPRNGNext() (gas: 16184)
+LibPRNGTest:testPRNGShuffle() (gas: 504507)
 LibPRNGTest:testPRNGShuffleBytes() (gas: 222095)
 LibPRNGTest:testPRNGShuffleBytesGas() (gas: 1322452)
 LibPRNGTest:testPRNGShuffleGas() (gas: 1610949)
 LibPRNGTest:testPRNGUniform() (gas: 559341)
-LibPRNGTest:testStandardNormalWad() (gas: 4254566)
-LibPRNGTest:test__codesize() (gas: 5879)
+LibPRNGTest:testStandardNormalWad() (gas: 4260605)
+LibPRNGTest:test__codesize() (gas: 6473)
 LibRLPTest:testComputeAddressDifferential(address,uint256) (runs: 256, μ: 1913, ~: 1633)
 LibRLPTest:testComputeAddressForLargeNonces() (gas: 1771)
 LibRLPTest:testComputeAddressForSmallNonces() (gas: 967)
@@ -858,29 +859,29 @@ MetadataReaderLibTest:testReadStringTruncated(uint256) (runs: 256, μ: 799541, ~
 MetadataReaderLibTest:testReadUint() (gas: 911216)
 MetadataReaderLibTest:testReadUint(uint256) (runs: 256, μ: 46033, ~: 47093)
 MetadataReaderLibTest:test__codesize() (gas: 9186)
-MinHeapLibTest:testHeapEnqueue(uint256) (runs: 256, μ: 188312, ~: 192351)
-MinHeapLibTest:testHeapEnqueue2(uint256) (runs: 256, μ: 628359, ~: 566381)
+MinHeapLibTest:testHeapEnqueue(uint256) (runs: 256, μ: 185713, ~: 187563)
+MinHeapLibTest:testHeapEnqueue2(uint256) (runs: 256, μ: 648386, ~: 573268)
 MinHeapLibTest:testHeapEnqueueGas() (gas: 854801)
-MinHeapLibTest:testHeapEnqueueZeroMaxLengthReverts(uint256) (runs: 256, μ: 47693, ~: 28615)
-MinHeapLibTest:testHeapPSiftTrick(uint256,uint256,uint256) (runs: 256, μ: 716, ~: 877)
-MinHeapLibTest:testHeapPushAndPop(uint256) (runs: 256, μ: 107245, ~: 99564)
-MinHeapLibTest:testHeapPushPop(uint256) (runs: 256, μ: 245886, ~: 244732)
-MinHeapLibTest:testHeapReplace(uint256) (runs: 256, μ: 293870, ~: 302883)
-MinHeapLibTest:testHeapReplaceOrPopEmptyHeapReverts(uint256) (runs: 256, μ: 47585, ~: 50916)
-MinHeapLibTest:testHeapRoot(uint256) (runs: 256, μ: 48520, ~: 47811)
-MinHeapLibTest:testHeapSmallest(uint256) (runs: 256, μ: 1703864, ~: 1271598)
+MinHeapLibTest:testHeapEnqueueZeroMaxLengthReverts(uint256) (runs: 256, μ: 48157, ~: 48176)
+MinHeapLibTest:testHeapPSiftTrick(uint256,uint256,uint256) (runs: 256, μ: 711, ~: 708)
+MinHeapLibTest:testHeapPushAndPop(uint256) (runs: 256, μ: 103034, ~: 99564)
+MinHeapLibTest:testHeapPushPop(uint256) (runs: 256, μ: 246900, ~: 245280)
+MinHeapLibTest:testHeapReplace(uint256) (runs: 256, μ: 298263, ~: 315174)
+MinHeapLibTest:testHeapReplaceOrPopEmptyHeapReverts(uint256) (runs: 256, μ: 46295, ~: 50906)
+MinHeapLibTest:testHeapRoot(uint256) (runs: 256, μ: 49603, ~: 67523)
+MinHeapLibTest:testHeapSmallest(uint256) (runs: 256, μ: 1666715, ~: 1271598)
 MinHeapLibTest:testHeapSmallestGas() (gas: 50027269)
-MinHeapLibTest:testMemHeapEnqueue(uint256) (runs: 256, μ: 857848, ~: 823530)
-MinHeapLibTest:testMemHeapEnqueue2(uint256) (runs: 256, μ: 793762, ~: 741115)
+MinHeapLibTest:testMemHeapEnqueue(uint256) (runs: 256, μ: 897229, ~: 824239)
+MinHeapLibTest:testMemHeapEnqueue2(uint256) (runs: 256, μ: 803087, ~: 735967)
 MinHeapLibTest:testMemHeapEnqueueGas() (gas: 233981)
-MinHeapLibTest:testMemHeapEnqueueZeroMaxLengthReverts(uint256) (runs: 256, μ: 25191, ~: 24159)
-MinHeapLibTest:testMemHeapPushAndPop(uint256) (runs: 256, μ: 824109, ~: 748166)
+MinHeapLibTest:testMemHeapEnqueueZeroMaxLengthReverts(uint256) (runs: 256, μ: 25211, ~: 24170)
+MinHeapLibTest:testMemHeapPushAndPop(uint256) (runs: 256, μ: 795803, ~: 743874)
 MinHeapLibTest:testMemHeapPushGas() (gas: 31982)
-MinHeapLibTest:testMemHeapPushPop(uint256) (runs: 256, μ: 1028744, ~: 775102)
-MinHeapLibTest:testMemHeapReplace(uint256) (runs: 256, μ: 879282, ~: 824296)
-MinHeapLibTest:testMemHeapReplaceOrPopEmptyHeapReverts(uint256) (runs: 256, μ: 25909, ~: 25517)
-MinHeapLibTest:testMemHeapRoot(uint256) (runs: 256, μ: 694847, ~: 644571)
-MinHeapLibTest:testMemHeapSmallest(uint256) (runs: 256, μ: 835597, ~: 744766)
+MinHeapLibTest:testMemHeapPushPop(uint256) (runs: 256, μ: 1061534, ~: 775073)
+MinHeapLibTest:testMemHeapReplace(uint256) (runs: 256, μ: 895879, ~: 845419)
+MinHeapLibTest:testMemHeapReplaceOrPopEmptyHeapReverts(uint256) (runs: 256, μ: 25898, ~: 25371)
+MinHeapLibTest:testMemHeapRoot(uint256) (runs: 256, μ: 720351, ~: 644424)
+MinHeapLibTest:testMemHeapSmallest(uint256) (runs: 256, μ: 807207, ~: 744317)
 MinHeapLibTest:testMemHeapSmallestGas() (gas: 2871518)
 MinHeapLibTest:test__codesize() (gas: 13708)
 MulticallableTest:testMulticallableBenchmark() (gas: 29588)
@@ -948,8 +949,8 @@ ReceiverTest:testOnERC721Received() (gas: 64127)
 ReceiverTest:test__codesize() (gas: 3310)
 RedBlackTreeLibTest:testRedBlackTreeBenchUint160() (gas: 3429423)
 RedBlackTreeLibTest:testRedBlackTreeBenchUint256() (gas: 5843130)
-RedBlackTreeLibTest:testRedBlackTreeInsertAndRemove(uint256) (runs: 256, μ: 627514, ~: 487609)
-RedBlackTreeLibTest:testRedBlackTreeInsertAndRemove2(uint256) (runs: 256, μ: 420920, ~: 390160)
+RedBlackTreeLibTest:testRedBlackTreeInsertAndRemove(uint256) (runs: 256, μ: 624653, ~: 478269)
+RedBlackTreeLibTest:testRedBlackTreeInsertAndRemove2(uint256) (runs: 256, μ: 407187, ~: 348916)
 RedBlackTreeLibTest:testRedBlackTreeInsertAndRemove3() (gas: 21546528)
 RedBlackTreeLibTest:testRedBlackTreeInsertBenchStep() (gas: 3703676)
 RedBlackTreeLibTest:testRedBlackTreeInsertBenchUint160() (gas: 3469294)
@@ -958,16 +959,16 @@ RedBlackTreeLibTest:testRedBlackTreeInsertOneGas() (gas: 45564)
 RedBlackTreeLibTest:testRedBlackTreeInsertTenGas() (gas: 282751)
 RedBlackTreeLibTest:testRedBlackTreeInsertThreeGas() (gas: 96193)
 RedBlackTreeLibTest:testRedBlackTreeInsertTwoGas() (gas: 69797)
-RedBlackTreeLibTest:testRedBlackTreeNearest(uint256) (runs: 256, μ: 228165, ~: 224911)
-RedBlackTreeLibTest:testRedBlackTreeNearestAfter(uint256) (runs: 256, μ: 251011, ~: 248814)
-RedBlackTreeLibTest:testRedBlackTreeNearestBefore(uint256) (runs: 256, μ: 232084, ~: 195909)
+RedBlackTreeLibTest:testRedBlackTreeNearest(uint256) (runs: 256, μ: 234562, ~: 224911)
+RedBlackTreeLibTest:testRedBlackTreeNearestAfter(uint256) (runs: 256, μ: 239070, ~: 225794)
+RedBlackTreeLibTest:testRedBlackTreeNearestBefore(uint256) (runs: 256, μ: 224829, ~: 195447)
 RedBlackTreeLibTest:testRedBlackTreePointers() (gas: 91897)
 RedBlackTreeLibTest:testRedBlackTreeRejectsEmptyValue() (gas: 3238)
 RedBlackTreeLibTest:testRedBlackTreeRemoveViaPointer() (gas: 58172)
 RedBlackTreeLibTest:testRedBlackTreeTreeFullReverts() (gas: 50335)
 RedBlackTreeLibTest:testRedBlackTreeTryInsertAndRemove() (gas: 56139)
 RedBlackTreeLibTest:testRedBlackTreeValues() (gas: 191599)
-RedBlackTreeLibTest:testRedBlackTreeValues(uint256) (runs: 256, μ: 383999, ~: 267070)
+RedBlackTreeLibTest:testRedBlackTreeValues(uint256) (runs: 256, μ: 392283, ~: 267070)
 RedBlackTreeLibTest:test__codesize() (gas: 14087)
 ReentrancyGuardTest:testRecursiveDirectUnguardedCall() (gas: 34254)
 ReentrancyGuardTest:testRecursiveIndirectUnguardedCall() (gas: 47771)

src/utils/LibPRNG.sol

@@ -175,4 +175,37 @@ library LibPRNG {
                 shr(192, mul(s, and(m, mulmod(r3, a, n))))))), 7745966692414833770)
         }
     }
+
+    /// @dev Returns a sample from the unit exponential distribution denominated in `WAD`.
+    function exponentialWad(PRNG memory prng) internal pure returns (uint256 result) {
+        /// @solidity memory-safe-assembly
+        assembly {
+            let n := 0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff43 // Prime.
+            let a := 0x100000000000000000000000000000051 // Prime and a primitive root of `n`.
+            let m := 0x7fffffffffffffffffffffffffffffff
+            let r := keccak256(prng, 0x20)
+            mstore(prng, r)
+            let p := and(m, r)
+            let w := and(m, shr(127, r))
+            if iszero(gt(w, p)) {
+                // Passes the Kolmogorov-Smirnov test for 200k samples.
+                // Gas usage varies, starting from about 199+ gas.
+                for {} 1 {} {
+                    r := mulmod(r, a, n)
+                    let v := and(m, r)
+                    if iszero(lt(v, w)) {
+                        r := mulmod(r, a, n)
+                        result := add(result, 1000000000000000000)
+                        w := and(m, shr(127, r))
+                        p := and(m, r)
+                        if iszero(lt(w, p)) { break }
+                        continue
+                    }
+                    w := and(m, shr(127, r))
+                    if iszero(lt(w, v)) { break }
+                }
+            }
+            result := add(result, div(p, 170141183460469231732))
+        }
+    }
 }

test/LibPRNG.t.sol

@@ -201,4 +201,36 @@ contract LibPRNGTest is SoladyTest {
             assertLt(FixedPointMathLib.abs(sd - wad), uint256(wad / 8));
         }
     }
+
+    function testExponentialWad() public {
+        LibPRNG.PRNG memory prng;
+        unchecked {
+            uint256 n = 1000;
+            int256 oldM;
+            int256 newM;
+            int256 oldS;
+            int256 newS;
+            for (uint256 i; i != n;) {
+                uint256 gasBefore = gasleft();
+                int256 x = int256(prng.exponentialWad());
+                uint256 gasUsed = gasBefore - gasleft();
+                if (++i == 1) {
+                    oldM = (newM = x);
+                } else {
+                    newM = oldM + (x - oldM) / int256(i);
+                    newS = oldS + (x - oldM) * (x - newM);
+                    oldM = newM;
+                    oldS = newS;
+                }
+                emit LogInt("exponentialWad", x);
+                emit LogUint("gasUsed", gasUsed);
+            }
+            int256 wad = int256(FixedPointMathLib.WAD);
+            emit LogInt("mean", newM);
+            int256 sd = int256(FixedPointMathLib.sqrt(uint256(newS / int256(n - 1))));
+            assertLt(FixedPointMathLib.abs(newM - wad), uint256(wad / 8));
+            emit LogInt("standard deviation", sd);
+            assertLt(FixedPointMathLib.abs(sd - wad), uint256(wad / 8));
+        }
+    }
 }