Probabilistic Sieve of Eratosthenes

algorithm, math

Solution

Let's define a series of random variables recursively:

Let Xk,r denote the indicator variable, taking value `1` iff `X[k] == true` by the end of the iteration in which the variable `i` took value `r`.

In order to have fewer symbols and since it makes more intuitive sense with the code, we'll just write Xk,i which is valid although would have been confusing in the definition since `i` taking value `i` is confusing when the first refers to the variable in the loop and the latter to the value of the variable.

Now we note that:

P(Xk,i ~ 0) = P(Xk,i-1 ~ 0) + P(Xk,i-1 ~ 1) * P(Xk-1,i-1 ~ 1) * 1/i

(~ is used in place of = just to make it understandable, since = would otherwise take two separate meanings and looks confusing).

This equality holds by virtue of the fact that either `X[k]` was `false` at the end of the `i` iteration either because it was false at the end of the `i-1`, or it was `true` at that point, but in that last iteration `X[k-1]` was `true` and so we entered the loop and changed `X[k]` with probability of 1/i. The events are mutually exclusive, so there is no intersection.

The base of the recursion is simply the fact that P(Xk,1 ~ 1) = 1 and P(X2,i ~ 1) = 1.

Lastly, we note simply that P(`X[k] == true`) = P(Xk,k-1 ~ 1).

This can be programmed rather easily. Here's a javascript implementation that employs memoisation (you can benchmark if using nested indices is better than string concatenation for the dictionary index, you could also redesign the calculation to maintain the same runtime complexity but not run out of stack size by building bottom-up and not top-down). Naturally the implementation will have a runtime complexity of `O(k^2)` so it's not practical for arbitrarily large numbers:

function P(k) {
   if (k<2 || k!==Math.round(k)) return -1;
   var _ = {};
   function _P(n,i) {
      if(n===2) return 1;
      if(i===1) return 1;
      var $ = n+'_'+i;
      if($ in _) return _[$];
      return _[$] =  1-(1-_P(n,i-1) + _P(n,i-1)*_P(n-1,i-1)*1/i);
   }
   return _P(k,k-1);
}

P(1000); // 0.12274162882390949

More interesting would be how the 1/i probability changes things. I.e. whether or not the probability converges to 0 or to some other value, and if so, how changing the 1/i affects that.

Of course if you ask on mathSE you might get a better answer - this answer is pretty simplistic, I'm sure there is a way to manipulate it to acquire a direct formula.

Problem

Consider the following algorithm. ``` function Rand(): return a uniformly random real between 0.0 and 1.0 function Sieve(n): assert(n >= 2) for i = 2 to n X[i] = true for i = 2 to n if (X[i]) for j = i+1 to n if (Rand() < 1/i) X[j] = false return X[n] ``` What is the probability that Sieve(k) returns `true` as a function of k ?

Original source