# http://recursed.blogspot.com/2010/05/how-to-test-for-syphilis.html

# values of calc(n,1) are particularly intriguing.

# cache
d = {}
n_max = 0

# any is 0 if we don't know there are cases in this group,
# else 1.
# This function is just a memoizing wrapper.
def calc(n,any):
  # return (#groups, expected #tests)
  try: return d[n,any]
  except:
    d[n,any] = do_calc(n,any)
    return d[n,any]

def do_calc(n,any):
  global n_max
  if n > n_max+15:
    # crawl upward, filling the cache, so as not to blow
    # the stack.
    for m in range(n_max+10,n,10):
      calc(m,any)
  n_max = n
  if n==0: return (0,0) # nothing to do
  if any:
    if n==1: return (0,0) # we already know who the case is
    best = (0,n+1)
    for k in range(1,n): # no point testing whole group, of course
      # first test a group of size k, then continue
      q = prob(k,n) # probability of positive result for the group
      # expected cost:
      #   1 to test the group
      #   with prob q: test the group (knowing it has a positive)
      #                and test the remainder (which now need no longer be positive)
      #   with prob 1-q: test the remainder (still known positive)
      cost = 1 + q*(calc(k,1)[1] + calc(n-k,0)[1]) + (1-q)*calc(n-k,1)[1]
      if cost<best[1]: best = (k,cost)
    return best
  else:
    if n==1: return (1,1) # just test our single subject
    best = (0,n+1)
    for k in range(1,n+1):
      # first test a group of size k, then continue
      q = 1-(1-p)**k # probability of positive result for the group
      cost = 1 + q*calc(k,1)[1] + calc(n-k,0)[1]
      if cost < best[1]: best = (k,cost)
    return best

# probability of >= 1 case in a group of size k,
# given that there's at least one in a containing group of size n
def prob(k,n):
  # = Pr(>=1 in k | >=1 in n)
  # = Pr(>=1 in k) / Pr(>=1 in n)
  return (1-(1-p)**k) / (1-(1-p)**n)