There is surely a more elegant way to do this, but a brute-force approach would be to generate all possible 2x2 tables with the given row and column marginals and then filter the ones that satisfy the joint marginal conditions.
library(tidyverse)
## Constraints
r1_marg <- 20
r2_marg <- 20
c1_marg <- 29
c2_marg <- 11
## Range of values
r1c1 <- c(0:max(c(r1_marg, c1_marg)))
r1c2 <- c(0:max(c(r1_marg, c2_marg)))
r2c1 <- c(0:max(c(c1_marg, r2_marg)))
r2c2 <- c(0:max(c(c2_marg, r2_marg)))
## This is super-inefficient and I am
## sure there is a clever way to do
## this directly.
marg_combos <- expand_grid(r1c1, r1c2,
r2c1, r2c2) %>%
filter(r1c1 + r1c2 == r1_marg &
r2c1 + r2c2 == r2_marg &
r1c1 + r2c1 == c1_marg &
r1c2 + r2c2 == c2_marg) %>%
tibble::rowid_to_column(var = "tab_id") %>%
pivot_longer(r1c1:r2c2, names_to = "pos") %>%
group_by(tab_id) %>%
summarize(mat = list(matrix(value,
nrow = 2, ncol = 2, byrow = TRUE))) %>%
group_by(tab_id) %>%
mutate(fisher = map(mat, fisher.test)) # probably want to read the fisher docs
marg_combos
#> # A tibble: 12 x 3
#> # Groups: tab_id [12]
#> tab_id mat fisher
#> <int> <list> <list>
#> 1 1 <int[,2] [2 × 2]> <htest>
#> 2 2 <int[,2] [2 × 2]> <htest>
#> 3 3 <int[,2] [2 × 2]> <htest>
#> 4 4 <int[,2] [2 × 2]> <htest>
#> 5 5 <int[,2] [2 × 2]> <htest>
#> 6 6 <int[,2] [2 × 2]> <htest>
#> 7 7 <int[,2] [2 × 2]> <htest>
#> 8 8 <int[,2] [2 × 2]> <htest>
#> 9 9 <int[,2] [2 × 2]> <htest>
#> 10 10 <int[,2] [2 × 2]> <htest>
#> 11 11 <int[,2] [2 × 2]> <htest>
#> 12 12 <int[,2] [2 × 2]> <htest>
## E.g.
marg_combos$mat[[1]]
#> [,1] [,2]
#> [1,] 9 11
#> [2,] 20 0
marg_combos$fisher[[1]]
#>
#> Fisher's Exact Test for Count Data
#>
#> data: .x[[i]]
#> p-value = 0.0001453
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#> 0.0000000 0.2397675
#> sample estimates:
#> odds ratio
#> 0
marg_combos$mat[[7]]
#> [,1] [,2]
#> [1,] 15 5
#> [2,] 14 6