Brain teaser: How to generate 7 integers with equal probability using a biased coin that has a pr(head) = p?

Question

This is a question I found on Glassdoor: How does one generate 7 integers with equal probability using a coin that has a $\mathbb{Pr}(\text{Head}) = p\in(0,1)$?

Basically, you have a coin that may or may not be fair, and this is the only random-number generating process you have, so come up with random number generator that outputs integers from 1 to 7 where the probability of getting each of these integers is 1/7.

Efficiency of the data-generates process matters.

Answer 1

Flip the coin twice. If it lands HH or TT, ignore it and flip it twice again.

Now, the coin has equal probability of coming up HT or TH. If it comes up HT, call this H1. If it comes up TH, call this T1.

Keep obtaining H1 or T1 until you have three in a row. These three results give you a number based on the table below:

H1 H1 H1 -> 1
H1 H1 T1 -> 2
H1 T1 H1 -> 3
H1 T1 T1 -> 4
T1 H1 H1 -> 5
T1 H1 T1 -> 6
T1 T1 H1 -> 7
T1 T1 T1 -> [Throw out all results so far and repeat]

I argue that this would work perfectly fine, although you would have a lot of wasted throws in the process!

Answer 2

Assume that $p \in (0,1)$.

Step 1:. Toss the coin 5 times.

If the outcome is

$(H, H, H, T, T)$, return $1$ and stop.

$(H, H, T, T, H)$, return $2$ and stop.

$(H, T, T, H, H)$, return $3$ and stop.

$(T, T, H, H, H)$, return $4$ and stop.

$(T, H, H, H, T)$, return $5$ and stop.

$(H, H, T, H, T)$, return $6$ and stop.

$(H, T, H, T, H)$, return $7$ and stop.

Step 2:. If the outcome is none of the above, repeat Step 1.

Note that regardless of the value of $p \in (0,1)$, each of the seven outcomes listed above has probability $q = p^3(1-p)^2$, and the expected number of coin tosses is $\displaystyle \frac{5}{7q}$. The tosser doesn't need to know the value of $p$ (except that $p\neq 0$ and $p\neq 1$); it is guaranteed that the seven integers are equally likely to be returned by the experiment when it terminates (and it is guaranteed to end with probability $1$).

Answer 3

Generalizing the case described by Dilip Sarwate

Some of the methods described in the other answers use a scheme in which you throw a sequence of $n$ coins in a 'turn' and depending on the result you choose a number between 1 or 7 or discard the turn and throw again.

The trick is to find in the expansion of possibilities a multiple of 7 outcomes with the same probability $p^k(1-p)^{n-k}$ and match those against each other.

Because the total number of outcomes is not a multiple of 7, we have a few outcomes that we can not assign to a number, and have some probability that we need to discard the outcomes and start over.

---

The case of using 7 coin flips per turn

Intuitively we could say that rolling the dice seven times would be very interesting. Since we only need to throw out 2 out of the $2^7$ possibilities. Namely, the 7 times heads and 0 times heads.

For all other $2^7-2$ possibilities there is always a multiple of 7 cases with the same number of heads. Namely 7 cases with 1 heads, 21 cases with 2 heads, 35 cases with 3 heads, 35 cases with 4 heads, 21 cases with 5 heads, and 7 cases with 6 heads.

So if you compute the number (discarding 0 heads and 7 heads) $$X = \sum_{k=1}^{7} (k-1) \cdot C_k $$

with $C_k$ Bernoulli distributed variables (value 0 or 1), then X modulo 7 is a uniform variable with seven possible results.

---

Comparing different number of coin flips per turn

The question remains what the optimal number of rolls per turn would be. Rolling more dices per turn cost you more, but you reduce the probability to have to roll again.

The image below shows a manual computations for the first few numbers of coin flips per turn. (possibly there might be an analytical solution, but I believe it is safe to say that a system with 7 coin flips provides the best method regarding the expectation value for the necessary number of coin flips)

# plot an empty canvas
plot(-100,-100,
     xlab="flips per turn",
     ylab="E(total flips)",
     ylim=c(7,400),xlim=c(0,20),log="y")
title("expectation value for total number of coin flips
(number of turns times flips per turn)")

# loop 1
# different values p from fair to very unfair 
# since this is symmetric only from 0 to 0.5 is necessary 

# loop 2
# different values for number of flips per turn
# we can only use a multiple of 7 to assign 
#   so the modulus will have to be discarded
#   from this we can calculate the probability that the turn succeeds
#   the expected number of flips is 
#       the flips per turn 
#             divided by 
#       the probability for the turn to succeed 

for (p in c(0.5,0.2,0.1,0.05)) {
  Ecoins <- rep(0,16)
  for (dr in (5:20)){
    Pdiscards = 0
    for (i in c(0:dr)) { 
      Pdiscards = Pdiscards + p^(i)*(1-p)^(dr-i) * (choose(dr,i) %% 7)
    }
    Ecoins[dr-4] = dr/(1-Pdiscards)
  }
  lines(5:20, Ecoins)
  points(5:20, Ecoins, pch=21, col="black", bg="white", cex=0.5)
  text(5, Ecoins[1], paste0("p = ",p), pos=2)
}

---

Using an early stopping rule

note: the calculations below, for the expectation value of number of flips, are for a fair coin $p=0.5$, it would become a mess to do this for different $p$, but the principle remains the same (although different book-keeping of the cases is needed)

We should be able to choose the cases (instead of the formula for $X$) such that we might be able to stop earlier.

• With 5 coin flips we have for the six possible different unordered sets of heads and tails:

  1+5+10+10+5+1 ordered sets

  And we can use the groups with ten cases (that is the group with 2 heads or the group with 2 tails) to choose (with equal probability) a number. This occurs in 14 out of 2^5=32 cases. This leaves us with:

  1+5+3+3+5+1 ordered sets

• With an extra (6-th) coin flip we have for the seven possible different unordered sets of heads and tails:

  1+6+8+6+8+6+1 ordered sets

  And we can use the groups with eight cases (that is the group with 3 heads or the group with 3 tails) to choose (with equal probability) a number. This occurs in 14 out of 2*(2^5-14)=36 cases. This leaves us with:

  1+6+1+6+1+6+1 ordered sets

• With another (7-th) extra coin flip we have for the eight possible different unordered sets of heads and tails:

  1+7+7+7+7+7+7+1 ordered sets

  And we can use the groups with seven cases (all except the all tails and all heads cases) to choose (with equal probability) a number. This occurs in 42 out of 44 cases. This leaves us with:

  1+0+0+0+0+0+0+1 ordered sets

  (we could continue this but only in the 49-th step does this give us an advantage)

So the probability to select a number

• at 5 flips is $\frac{14}{32} = \frac{7}{16}$
• at 6 flips is $\frac{9}{16}\frac{14}{36} = \frac{7}{32}$
• at 7 flips is $\frac{11}{32}\frac{42}{44} = \frac{231}{704}$
• not in 7 flips is $1-\frac{7}{16}-\frac{7}{32}-\frac{231}{704} = \frac{2}{2^7}$

This makes the expectation value for the number of flips in one turn, conditional that there is success and p=0.5:

$$5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704} = 5.796875 $$

The expectation value for the total number of flips (until there is a success), conditional that p=0.5, becomes:

$$\left(5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704}\right) \frac{2^7}{2^7-2} = \frac{53}{9} = 5.88889 $$

---

The answer by NcAdams uses a variation of this stopping-rule strategy (each time come up with two new coin flips) but is not optimally selecting out all the flips.

The answer by Clid might be similar as well although there might be an uneven selection rule that each two coin flips a number might be chosen but not necessarily with equal probability (a discrepancy which is being repaired during later coin flips)

---

Comparison with other methods

Other methods using a similar principle are the one by NcAdams and AdamO.

The principle is: A decision for a number between 1 and 7 is made after a certain number of heads and tails. After an $x$ number of flips, for each decision that leads to a number $i$ there is a similar, equally probable, decision that leads to a number $j$ (the same number of heads and tails but just in a different order). Some series of heads and tails can lead to a decision to start over.

For such type of methods the one placed here is the most efficient because it makes the decisions as early as possible (as soon as there is a possibility for 7 equal probability sequences of heads and tails, after the $x$-th flip, we can use them to make a decision on a number and we do not need to flip further if we encounter one of those cases).

This is demonstrated by the image and simulation below:

#### mathematical part #####
set.seed(1)


#plotting this method
p <- seq(0.001,0.999,0.001)
tot <- (5*7*(p^2*(1-p)^3+p^3*(1-p)^2)+
       6*7*(p^2*(1-p)^4+p^4*(1-p)^2)+
       7*7*(p^1*(1-p)^6+p^2*(1-p)^5+p^3*(1-p)^4+p^4*(1-p)^3+p^5*(1-p)^2+p^6*(1-p)^1)+
        7*1*(0+p^7+(1-p)^7) )/
             (1-p^7-(1-p)^7)
plot(p,tot,type="l",log="y",
     xlab="p",
     ylab="expactation value number of flips"
     )

#plotting method by AdamO
tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
          (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
          (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
          35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
          1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
           *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
           5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
           15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot,col=2,lty=2)

#plotting method by NcAdam
lines(p,3*8/7/(p*(1-p)),col=3,lty=2)

legend(0.2,500,
       c("this method calculation","AdamO","NcAdams","this method simulation"),
       lty=c(1,2,2,0),pch=c(NA,NA,NA,1),col=c(1,2,3,1))


##### simulation part ######

#creating decision table
mat<-matrix(as.numeric(intToBits(c(0:(2^5-1)))),2^5,byrow=1)[,c(1:12)]
colnames(mat) <- c("b1","b2","b3","b4","b5","b6","b7","sum5","sum6","sum7","decision","exit")

# first 5 rolls
mat[,8] <- sapply(c(1:2^5), FUN = function(x) {sum(mat[x,1:5])})

mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 3 heads
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 6th roll
mat <- rbind(mat,mat)
mat[c(33:64),6] <- rep(1,32)
mat[,9] <- sapply(c(1:2^6), FUN = function(x) {sum(mat[x,1:6])})

mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 4 heads
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 7th roll
mat <- rbind(mat,mat)
mat[c(65:128),7] <- rep(1,64)
mat[,10] <- sapply(c(1:2^7), FUN = function(x) {sum(mat[x,1:7])})

for (i in 1:6) {
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],12] = rep(7,7) # we can stop for 7 cases with i heads
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],11] = c(1:7)   
}


mat[1,12] = 7           # when we did not have succes we still need to count the 7 coin tosses
mat[2^7,12] = 7


draws = rep(0,100)
num = rep(0,100)
# plotting simulation
for (p in seq(0.05,0.95,0.05)) {
  n <- rep(0,1000)
  for (i in 1:1000) {
    coinflips <- rbinom(7,1,p)  # draw seven numbers
    I <- mat[,1:7]-matrix(rep(coinflips,2^7),2^7,byrow=1) == rep(0,7)                      # compare with the table
    Imatch = I[,1]*I[,2]*I[,3]*I[,4]*I[,5]*I[,6]*I[,7]        # compare with the table 
      draws[i] <- mat[which(Imatch==1),11]                 # result which number
      num[i]   <- mat[which(Imatch==1),12]                 # result how long it took
  }
  Nturn <- mean(num)                   #how many flips we made
  Sturn <- (1000-sum(draws==0))/1000   #how many numbers we got (relatively)
  points(p,Nturn/Sturn)
}

another image which is scaled by $p*(1-p)$ for better comparison:

zoom in comparing methods described in this post and comments

the 'conditional skipping of the 7-th step' is a slight improvement which can be made on the early stopping rule. In this case you select not groups with equal probabilities after the 6-th flips. You have 6 groups with equal probabilities, and 1 groups with a slightly different probability (for this last group you need to flip one more extra time when you have 6 heads or tails and because you discard the 7 heads or 7 tails, you will end up with the same probability after all)

Answer 4

Divide a box into seven equal-area regions, each labeled with an integer. Throw the coin into the box in such a way that it has equal probability of landing in each region.

This is only half in jest -- it's essentially the same procedure as estimating $\pi$ using a physical Monte Carlo procedure, like dropping rice grains onto paper with a circle drawn on it.

This is one of the only answer that works for the case of $p = 1$ or $p=0$.

Answer 5

EDIT: based on others' feedback.

Here's an interesting thought:

set the list of {1,2,3,4,5,6,7}. Throw the coin for each element in the list sequentially. If it lands head side up for a particular element, remove the number from the list. If all the numbers from a particular iteration of the list are removed, repeat the sampling. Do so until only one number remains.

drop.one <- function(x, p) {
  drop <- runif(length(x)) < p
  if (all(drop))
    return(x)
  return(x[!drop])
}

sample.recur <- function(x, p) {
  if (length(x) > 1)
    return(sample.recur(drop.one(x, p), p))
  return(x)
}

# x <- c(1:7,7:1)
x <- 1:7
p <- 0.01

out <- replicate(1e5, sample.recur(x, p))

round(prop.table(table(out)), 2)

gives me an approximately uniform distribution

> round(prop.table(table(out)), 2)
out
   1    2    3    4    5    6    7 
0.14 0.14 0.15 0.14 0.14 0.14 0.14 

It's interesting to note (if I haven't made a dire mistake) that this produces a different result than generating $N$ binomial outcomes as the sum of 13 tosses of the coin (counting 0 heads as an outcome) and mapping the {0,1,2,...,12,13} index onto the earlier list of {1,2,3,...,3,2,1}. I don't quite know how to prove that my method works.

---

Evaluation of expectation value for number of coin throws

The expectation value for the number coin throws can be calculated using the transition matrix below (answering the question when we start with $x$ non-eliminated numbers then what is the probability to get to $y$ non-eliminated numbers)

$$M = \begin{bmatrix} q^7 & 0 & 0 & 0 &0 & 0 & 1 & 1\\ 7p^1q^6 & q^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 21p^2q^5 & 6p^1q^5 & q^5 & 0 & 0 & 0 & 0 & 0\\ 35 p^3q^4 & 15 p^2q^4 & 5q^4 & q^4 & 0 & 0 & 0 & 0\\ 35 p^4q^3 & 20 p^3q^3 & 10 p^2q^3 & 4 p^1q^3 & q^3 & 0 & 0 & 0\\ 21p^5q^2 & 15 p^4q^2 & 10 p^3q^2 & 6 p^2q^2 & 3 p^1q^2 & q^2 & 0 & 0\\ 7p^6q^1 & 6 p^5q^1 & 5 p^4q^1 & 4 p^3q^1 & 3p^2q^1 & 2p^1q^1 & 0 & 0\\ p^7 & p^6 & p^5 & p^4 & p^3 & p^2 & 0 & 0 \end{bmatrix}$$

The eigenvector associated with the eigenvalue 1 (which can be found by solving $(M-I)v=0$) depicts how much time is relatively spend in what state. Then 7th state is how often you will be able to draw a number from 1 to 7. The other states tell how many coin flips it costs.

Below is the image which compares with the answer from NcAdams which has expectation value for coin throws being $E(n) = \frac{24}{7}p(1-p)$

Remarkable is that the method performs better roughly for $p>2/3$. But also the performance is non-symmetric. A symmetric and better overall performance could be made when a probabilistic switching rule would be made which changes the decision rule from tails to heads when heads happens to be improbable.

Solution found with wxMaxima

M: matrix(
 [(1-p)^7,        0,          0,0,0,0,1,1], 
 [7* p*(1-p)^6,   (1-p)^6,        0,0,0,0,0,0], 
 [21*p^2*(1-p)^5, 6*p*(1-p)^5,    (1-p)^5,0,0,0,0,0], 
 [35*p^3*(1-p)^4, 15*p^2*(1-p)^4, 5*p*(1-p)^4,(1-p)^4,0,0,0,0], 
 [35*p^4*(1-p)^3, 20*p^3*(1-p)^3, 10*p^2*(1-p)^3,4*p*(1-p)^3,(1-p)^3,0,0,0], 
 [21*p^5*(1-p)^2, 15*p^4*(1-p)^2, 10*p^3*(1-p)^2,6*p^2*(1-p)^2,3*p*(1-p)^2,(1-p)^2,0,0], 
 [7* p^6*(1-p)^1, 6*p^5*(1-p),    5*p^4*(1-p),4*p^3*(1-p),3*p^2*(1-p),2*(1-p)*p,0,0], 
 [p^7,        p^6,        p^5,p^4,p^3,p^2,0,0]
);
z: nullspace(M-diagmatrix(8,1));
x : apply (addcol, args (z));
t : [7,6,5,4,3,2,0,0];
plot2d(t.x/x[7],[p,0,1],logy);

Calculations in R

# plotting empty canvas
plot(-100,-100,
     xlab="p",
     ylab="E(total flips)",
     ylim=c(10,1000),xlim=c(0,1),log="y")

# plotting simulation
for (p in seq(0.1,0.9,0.05)) {

  n <- rep(0,10000)
  for (i in 1:10000) {
    success  = 0
    tests = c(1,1,1,1,1,1,1)     # start with seven numbers in the set
    count = 0
    while(success==0) {
      for (j in 1:7)  {
        if (tests[j]==1) {
          count = count + 1
          if  (rbinom(1,1,p) == 1) {
            tests[j] <- 0        # elliminate number when we draw heads
          }
        }
      }
      if (sum(tests)==1) {
        n[i] = count
        success = 1              # end     when 1 is left over
      }
      if (sum(tests)==0) {
        tests = c(1,1,1,1,1,1,1) # restart when 0 are left over
      }
    }
  }
  points(p,mean(n))
}

# plotting formula
p <- seq(0.001,0.999,0.001)

tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
    (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
    (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
    35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
    1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
   *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
 5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
  15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot)

#plotting comparison with alternative method
lines(p,3*8/7/(p*(1-p)),lty=2)

legend(0.2,500,
       c("simulation","calculation","comparison"),
       lty=c(0,1,2),pch=c(1,NA,NA))

Answer 6

The question is a bit ambiguous, is it asking "generate a random integer equal or less than 7 with equal probability", or is it asking "generate 7 random integers with equal probability?" - but what is the space of integers?!?

I'll assume it's the former, but the same logic I'm applying can be extended to the latter case too, once that problem is cleared up.

With a biased coin, you can produce a fair coin by following the following procedure: https://en.wikipedia.org/wiki/Fair_coin#Fair_results_from_a_biased_coin

A number 7 or less can be written in binary as three {0,1} digits. So all one needs to do is follow the above procedure three times, and convert the binary number produced back to decimal.

Answer 7

As mentioned in earlier comments, this puzzle relates to John von Neumann's 1951 paper "Various Techniques Used in Connection With Random Digits" published in the research journal of the National Bureau of Standards:

There is a wider literature about such problems that goes under the name of Bernoulli factory problems, that is, given a coin with tail probability $p$, how to simulate a coin with tail probability $f(p)$. If feasible, since some functions $f$ cannot be used as for instance $f(p)=\min\{1,2p\}$. Nacu and Peres (2005) study fast algorithms for solving [solvable] Bernoulli factory problems where fast means exponential decay of the tail of the distribution of the number $N$ of trials.

Answer 8

A solution that never wastes flips, which helps a lot for very-biased coins.

The disadvantage of this algorithm (as written, at least) is that it's using arbitrary-precision arithmetic. Practically, you probably want to use this until integer overflow, and only then throw it away and start over.

Also, you need to know what the bias is ... which you might not, say, if it is temperature-dependent like most physical phenomena.

---

Assuming the chance of heads is, say, 30%.

• Start with the range [1, 8).
• Flip your coin. If heads, use the left 30%, so your new range is [1, 3.1). Else, use the right 70%, so your new range is [3.1, 8).
• Repeat until the entire range has the same integer part.

---

Full code:

#!/usr/bin/env python3
from fractions import Fraction
from collections import Counter
from random import randrange


BIAS = Fraction(3, 10)
STAT_COUNT = 100000


calls = 0
def biased_rand():
    global calls
    calls += 1
    return randrange(BIAS.denominator) < BIAS.numerator


def can_generate_multiple(start, stop):
    if stop.denominator == 1:
        # half-open range
        stop = stop.numerator - 1
    else:
        stop = int(stop)
    start = int(start)
    return start != stop


def unbiased_rand(start, stop):
    if start < 0:
        # negative numbers round wrong
        return start + unbiased_rand(0, stop - start)
    assert isinstance(start, int) and start >= 0
    assert isinstance(stop, int) and stop >= start
    start = Fraction(start)
    stop = Fraction(stop)
    while can_generate_multiple(start, stop):
        if biased_rand():
            old_diff = stop - start
            diff = old_diff * BIAS
            stop = start + diff
        else:
            old_diff = stop - start
            diff = old_diff * (1 - BIAS)
            start = stop - diff
    return int(start)


def stats(f, *args, **kwargs):
    c = Counter()
    for _ in range(STAT_COUNT):
        c[f(*args, **kwargs)] += 1

    print('stats for %s:' % f.__qualname__)
    for k, v in sorted(c.items()):
        percent = v * 100 / STAT_COUNT
        print('  %s: %f%%' % (k, percent))


def main():
    #stats(biased_rand)
    stats(unbiased_rand, 1, 7+1)
    print('used %f calls at bias %s' % (calls/STAT_COUNT, BIAS))


if __name__ == '__main__':
    main()

Answer 9

This also only works for $p \neq 1$ and $p \neq 0$.

We first turn the (possibly) unfair coin into a fair coin using the process from NcAdams answer:

Flip the coin twice. If it lands HH or TT, ignore it and flip it twice again.

Now, the coin has equal probability of coming up HT or TH. If it comes up HT, call this H1. If it comes up TH, call this T1.

Now we use the fair coin to generate a real number between $0$ and $1$ in binary. Let H1$= 1$ and T1 $= 0$. Start with the string 0., flip the coin and append the resulting digit to the at the end of the string. Repeat with the new string. For example, the sequence H1 H1 T1 would give you the number $0.110$.

$1/7$ is a repeating decimal, and with the right-hand side being in base 2 we have that:

$1/7 = 0.001 001 001 ...$

$2/7 = 0.010 010 010 ...$

$3/7 = 0.011 011 011 ...$

$4/7 = 0.100 100 100 ...$

$5/7 = 0.101 101 101 ...$

$6/7 = 0.110 110 110 ...$

Keep flipping the fair coin to generate the decimal digits until the digits of your sequence does not match one of the above sequences, then chose the number $n$ such that your generated number is less than $n/7$ and greater than $(n-1)/7$. Since each generated number is equally likely we have chosen a number between $1$ and $7$ with equal probability.

Answer 10

Inspired by AdamO's answer, here is a Python solution that avoids bias:

def roll(p, n):
    remaining = range(1,n+1)
    flips = 0
    while len(remaining) > 1:
        round_winners = [c for c in remaining if random.choices(['H','T'], [p, 1.0-p]) == ['H']]
        flips += len(remaining)
        if len(round_winners) > 0:
            remaining = round_winners
        p = 1.0 - p
    return remaining[0], flips

There are two main changes here: The main one is that if all the number are discarded in a round, repeat the round. Also I flip the choice of whether heads or tails means discard every time. This reduces the number of flips needed in cases where p is close to 0 or 1 by ~70% when p=0.999

Answer 11

It appears we are allowed to change the mapping of the outcome of each flip, every time we flip. So, using for convenience the first seven positive integers, we give the following orders:

1st Flip, map $H \to 1$
2nd Flip, map $H \to 2$
...
7th flip, map $H \to 7$
8th flip, map $H \to 1$

etc

Repeat, always in batches of 7 flips. Map the $T$ outcomes to nothing.

SOME REMARKS ON EFFICIENCY
Our RNG, index it by $AP$, will generate zero useful outcomes in one 7-flip batch if we get $T$ in all 7 flips. So

$$P_{AP}(\text{no integers generated}) = (1-p)^7$$

As we run $N_b$ 7-flip batches, the total number of useless flips will tend to

$$\text{Count}_{AP}(\text{useless flips}) \to 7\cdot N_b(1-p)^7$$

Consider now the RNG of @DilipSarwate. There, we use a binomial $B(p,n=5)$ and 5-flip batches. The seven outcomes that generate an integer each has probability of occuring $p^3(1-p)^2$, so, in a 5-flip batch

$$P_{DS}(\text{no integers generated}) = 1-7\cdot p^3(1-p)^2$$

The count of useless flips will here tend to

$$\text{Count}_{DS}(\text{useless flips}) \to 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

For the $AP$ RNG to tend to produce less useless flips, it must be the case that

$$ \text{Count}_{AP}(\text{useless flips}) < \text{Count}_{DS}(\text{useless flips})$$

$$\implies 7\cdot N_b(1-p)^7 < 5\cdot N_b\cdot [1-7\cdot p^3(1-p)^2]$$

$$\implies 7\cdot (1-p)^7 < 5\cdot [1-7\cdot p^3(1-p)^2]$$

Numerical examination shows that if $p>0.0467$, then the $AP$ RNG generates less useless flips.

We also find that the number of useless flips is monotonically decreasing in $p$ for the $AP$ RNG, while for the $DS$ RNG it has a minimum at around $p\approx 0.5967$ and then increases again, while in general it stays high. The ratio

$$\frac{\text{Count}_{AP}(\text{useless flips})}{\text{Count}_{DS}(\text{useless flips})}$$

declines pretty quickly. For example it is equal to $0.67$ for$p=0.1$, equal to $0.3$ for $p=0.2$, equal to $0.127$ for $p=0.4$.