This problem arose in a different context at work, but I have translated it to pizza.

Suppose you have a circular pizza of radius $R$. Upon this disc, $n$ pepperoni will be distributed completely randomly. All pepperoni have the same radius $r$.

A pepperoni is "free" if it does not overlap any other pepperoni.

You are free to choose $n$.

Suppose you choose a small $n$. The chance that any given pepperoni is free are very large. But $n$ is small so the total number of free pepperoni is small. Suppose you choose a large $n$. The chance that any given pepperoni is free are small. But there are a lot of them.

Clearly, for a given $R$ and $r$, there is some optimal $n$ that maximizes the expected number of free pepperoni. How to find this optimum?

**Edit: picking the answer**

So it looks like leonbloy's answer given the best approximation in the cases I've looked at:

```
r/R n* by simulation n_free (sim) (R/2r)^2
0.1581 12 4.5 10
0.1 29 10.4 25
0.01 2550 929.7 2500
```

(There's only a few hundred trials in the r=0.01 sim, so 2550 might not be super accurate.) So I'm going to pick it for the answer. I'd like to thank everyone for their contributions, this has been a great learning experience.

Here are a few pictures of a simulation for r/R = 0.1581, n=12:

**Edit after three answers posted:**

I wrote a little simulation. I'll paste the code below so it can be checked (edit: it's been fixed to correctly pick points randomly on a unit disc). I've looked at ~~two~~ three cases so far. First case, r = 0.1581, R = 1, which is roughly p = 0.1 by mzp's notation. At these parameters I got n* = 12 (free pepperoni = 4.52). Arthur's expression did not appear to be maximized here. leonbloy's answer would give 10. I also did r = 0.1, R = 1. I got n* = 29 (free pepperoni = 10.38) in this case. Arthur's expression was not maximized here and leonbloy's answer would give 25. Finally for r = 0.01 I get roughly n*=2400 as shown here:

Here's my (ugly) code, now edited to properly pick random points on a disc:

```
from __future__ import division
import numpy as np
# the radius of the pizza is fixed at 1
r = 0.1 # the radius of the pepperoni
n_to_try = [1,5,10,20,25,27,28,29,30,31,32,33,35] # the number of pepperoni
trials = 10000# the number of trials (each trial randomly places n pepperoni)
def one_trial():
# place the pepperoni
pepperoni_coords = []
for i in range(n):
theta = np.random.rand()*np.pi*2 # a number between 0 and 2*pi
a = np.random.rand() # a number between 0 and 1
coord_x = np.sqrt(a) * np.cos(theta) # see http://mathworld.wolfram.com/DiskPointPicking.html
coord_y = np.sqrt(a) * np.sin(theta)
pepperoni_coords.append((coord_x, coord_y))
# how many pepperoni are free?
num_free_pepperoni = 0
for i in range(n): # for each pepperoni
pepperoni_coords_copy = pepperoni_coords[:] # copy the list so the orig is not changed
this_pepperoni = pepperoni_coords_copy.pop(i)
coord_x_1 = this_pepperoni[0]
coord_y_1 = this_pepperoni[1]
this_pepperoni_free = True
for pep in pepperoni_coords_copy: # check it against every other pepperoni
coord_x_2 = pep[0]
coord_y_2 = pep[1]
distance = np.sqrt((coord_x_1 - coord_x_2)**2 + (coord_y_1 - coord_y_2)**2)
if distance < 2*r:
this_pepperoni_free = False
break
if this_pepperoni_free:
num_free_pepperoni += 1
return num_free_pepperoni
for n in n_to_try:
results = []
for i in range(trials):
results.append(one_trial())
x = np.average(results)
print "For pizza radius 1, pepperoni radius", r, ", and number of pepperoni", n, ":"
print "Over", trials, "trials, the average number of free pepperoni was", x
print "Arthur's quantity:", x* ((((1-r)/1)**(x-1) - (r/1)) / ((1-r) / 1))
```