# 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= 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)