| 61 | } |
| 62 | |
| 63 | BlockNumber |
| 64 | BlockSampler_Next(BlockSampler bs) |
| 65 | { |
| 66 | BlockNumber K = bs->N - bs->t; /* remaining blocks */ |
| 67 | int k = bs->n - bs->m; /* blocks still to sample */ |
| 68 | double p; /* probability to skip block */ |
| 69 | double V; /* random */ |
| 70 | |
| 71 | Assert(BlockSampler_HasMore(bs)); /* hence K > 0 and k > 0 */ |
| 72 | |
| 73 | if ((BlockNumber) k >= K) |
| 74 | { |
| 75 | /* need all the rest */ |
| 76 | bs->m++; |
| 77 | return bs->t++; |
| 78 | } |
| 79 | |
| 80 | /*---------- |
| 81 | * It is not obvious that this code matches Knuth's Algorithm S. |
| 82 | * Knuth says to skip the current block with probability 1 - k/K. |
| 83 | * If we are to skip, we should advance t (hence decrease K), and |
| 84 | * repeat the same probabilistic test for the next block. The naive |
| 85 | * implementation thus requires a sampler_random_fract() call for each |
| 86 | * block number. But we can reduce this to one sampler_random_fract() |
| 87 | * call per selected block, by noting that each time the while-test |
| 88 | * succeeds, we can reinterpret V as a uniform random number in the range |
| 89 | * 0 to p. Therefore, instead of choosing a new V, we just adjust p to be |
| 90 | * the appropriate fraction of its former value, and our next loop |
| 91 | * makes the appropriate probabilistic test. |
| 92 | * |
| 93 | * We have initially K > k > 0. If the loop reduces K to equal k, |
| 94 | * the next while-test must fail since p will become exactly zero |
| 95 | * (we assume there will not be roundoff error in the division). |
| 96 | * (Note: Knuth suggests a "<=" loop condition, but we use "<" just |
| 97 | * to be doubly sure about roundoff error.) Therefore K cannot become |
| 98 | * less than k, which means that we cannot fail to select enough blocks. |
| 99 | *---------- |
| 100 | */ |
| 101 | V = sampler_random_fract(bs->randstate); |
| 102 | p = 1.0 - (double) k / (double) K; |
| 103 | while (V < p) |
| 104 | { |
| 105 | /* skip */ |
| 106 | bs->t++; |
| 107 | K--; /* keep K == N - t */ |
| 108 | |
| 109 | /* adjust p to be new cutoff point in reduced range */ |
| 110 | p *= 1.0 - (double) k / (double) K; |
| 111 | } |
| 112 | |
| 113 | /* select */ |
| 114 | bs->m++; |
| 115 | return bs->t++; |
| 116 | } |
| 117 | |
| 118 | /* |
| 119 | * This is used for sampling AO/CO row numbers, in the flattened |
no test coverage detected