Another modification to RNG generation to be both constant time, secure, and easy to implement as a circuit.

risc0/zkp/src/field/baby_bear.rs

@@ -109,18 +109,35 @@ impl field::Elem for Elem {
 
     /// Generate a random value within the Baby Bear field
     fn random(rng: &mut impl rand_core::RngCore) -> Self {
-        let val: u64 = (rng.next_u32() as u64) << 32 | (rng.next_u32() as u64);
-        // To make sure we are evenly distributed we pull a u64 divide it into
-        // baby-bear sized parts.  If it's not in one of those (chance of less than 6 *
-        // 10^-11) we implode.  In practice, running the proof again will succeed
-        // due to new ZK padding, and even at 1 proof-per-second, mean time to
-        // failure is > 2 years. In a better world, we might propagate this
-        // failure and retry at the proof level
-        const REJECT_CUTOFF: u64 = (u64::MAX / (P as u64)) * (P as u64);
-        if val >= REJECT_CUTOFF {
-            panic!("Random number generator got very unlucky");
+        // Normally, we would use rejection sampling here, but our specialized
+        // verifier circuit really wants an O(1) solution to sampling.  So instead, we
+        // sample [0, 2^192) % P.  This is very close to uniform, as we have 2^192 / P
+        // full copies of P, with only 2^192%P left over elements in the 'partial' copy
+        // (which we would normally reject with rejection sampling).
+        //
+        // Even if we imagined that this failure to reject totally destroys soundess,
+        // the probablity of it occuring even once during proving is vanishingly low
+        // (for the about 50 samples our current verifier pulls and at a probability of
+        // less than2^-161 per sample, this is less than 2^-155).  Even if we target
+        // a soundness of 128 bits, we are millions of times more likely to let an
+        // invalid proof by due to normal low probability events which are part of
+        // soundess analysis than due to imperfect sampling.
+        //
+        // Finally, from an implementation perspective, we can view generating a number
+        // in the [0, 2^192) range as using a linear combination of uniform u32s, r0,
+        // r1, etc and the following formula:
+        // u192 = r0 + 2^32 * r1 + 2^64 * r2 + ... + 2^160 * r5
+        // This is turn can be computed as:
+        // u192 = 2^32*(2^32*(2^32*(2^32*(2^32*(r5) + r4) + r3) + r2) + r1) + r0.
+        // Since we only need the final result modulo P, we can compute the entire
+        // expression above modulo P, and get the following implementation:
+        let mut val: u64 = 0;
+        for _ in 0..6 {
+            val <<= 32;
+            val += rng.next_u32() as u64;
+            val %= P as u64;
         }
-        Elem::from((val % (P as u64)) as u32)
+        Elem::from(val as u32)
     }
 
     fn from_u64(val: u64) -> Self {