# performs multiton subsampling on a count vector
# written by John Alroy on 21 March 2017

multiton<-function(n,target)	{
	S <- length(n)
	singletons <- length(which(n == 1))
	N <- sum(n)
	if (N / (singletons + 1) * (S - singletons) / S < target)
		return(NA)
	lastS <- 1
	for (i in 1:N)	{
		not <- 0
		thisSingletons <- 0
		all <- lchoose(N , i)
		for (j in 1:length(n))	{
			not <- not + exp(lchoose(N - n[j] , i) - all)
			thisSingletons <- thisSingletons + exp(log(n[j]) + lchoose(N - n[j] , i - 1) - all)
		}
		thisS <- S - not
		if (i / (thisSingletons + 1) * (thisS - thisSingletons) / thisS > target)
			break
		lastS <- thisS
	}
	return(as.numeric((lastS + thisS) / 2))
}