your table look like this:
a b | sab
c d | scd
----------
sac sbd| S
and you have 4 unknowns with 4 constraints (forget about the maximum constraints on a,b,c and d for a moment):
a+b=sab
a+c=sac
c+d=scd
b+d=sbd
the 4 constraints are not independent (otherwise you would have only one possible solution!), and a
bit of algebra shows that the matrix of this linear system has rank 3. so you have one degree of
freedom to play with. pick a for example, and then vary a from 0 to its maximum value. For each value
of a then compute b, c and d using the row and column sum constraints, and check that they satisy the
positivity and maximum constraints.
The R code for your example is as follows:
sab <- 59
scd <- 64
sac <- 31
sbd <- sab + scd - sac ### this is always true
amax <- 20
bmax <- 40
cmax <- 12
dmax <- 70
### let us vary a, our only degree of freedom
for (a in 0:amax){
### let us compute b, c and d by satisfying row and column sum constraints
b <- sab - a
c <- sac - a
d <- sbd - b
### let us check inequality constraints
if (b <= bmax && b>= 0 && c <= cmax && c >= 0 && d <= dmax && d >= 0){
cat("\nSolution:\n")
print(m <- rbind(c(a,b),c(c,d)))
cat("\nrowSums:", rowSums(m))
cat("\ncolsums:", colSums(m))
cat("\n---------------\n")
if (! identical(rowSums(m), c(sab,scd)))
stop("\nrow sum is not right!\n")
if (! identical(colSums(m), c(sac,sbd)))
stop("\ncolumns sum is not right!\n")
}
}