Write $\rho$ for the odds ratio, $\beta=\Pr(E)$, $\gamma=\Pr(O)$. Four independent equations are
$$\cases{a+b=\beta \\ a+c=\gamma \\ a+b+c+d=1 \\ ad = \rho bc.}$$
Adding the first two shows
$$b+c = \beta + \gamma - 2a. \tag{1}$$
Multiplying the first two equations and using $(1)$ yields
$$bc = \beta\gamma - (\beta+\gamma)a + a^2.\tag{2}$$
Multiplying the third equation by $a$, using the fourth to re-express $ad$ and plugging in $(1)$ and $(2)$ gives
$$a^2 + a(\beta+\gamma-2a) + \rho(\beta\gamma - (\beta+\gamma)a + a^2) = a.$$
In a more standard form, this a zero of the quadratic
$$(\rho-1)a^2 + [(\beta+\gamma)(1-\rho)-1]a + \rho\beta\gamma.\tag{3}$$
Provided $\beta, \gamma,\rho$ are consistent with some $2\times 2$ table, there will be at least one real zero of $(3)$, easily found using any quadratic formula. For either zero, solutions for the remaining entries are readily found from the first three of the original equations as
$$\cases{b=\beta-a\\c=\gamma-a\\d = 1+a-\beta-\gamma.}$$
There will be at most one valid solution for $a$, determined by the non-negativity of the coefficients. Here is an R implementation of the solution as the function f along with a test using randomly generated tables. The test outputs the random table, its reconstructed value from f, and a measure of the difference between them. By wrapping the test in replicate, I have run it 10,000 times. The final output gives the largest difference found: up to floating point error it equals zero, demonstrating the correctness of this approach.
f <- function(beta, gamma, rho, eps=1e-15) {
a <- rho-1
b <- (beta+gamma)*(1-rho)-1
c_ <- rho*beta*gamma
if (abs(a) < eps) {
z <- -c_ / b
} else {
d <- b^2 - 4*a*c_
if (d < eps*eps) s <- 0 else s <- c(-1,1)
z <- (-b + s*sqrt(max(0, d))) / (2*a)
}
y <- vapply(z, function(a) zapsmall(matrix(c(a, gamma-a, beta-a, 1+a-beta-gamma), 2, 2)),
matrix(0.0, 2, 2))
i <- apply(y, 3, function(u) all(u >= 0))
return(y[,,i])
}
set.seed(17)
i<-0
sim <- replicate(1e4, {
while(TRUE) {
x <- matrix(round(rexp(4), 2), 2, 2)
if(all(rowSums(x) > 0) && all(colSums(x) > 0) && x[1,2]*x[2,1] > 0) break
}
x <- x / sum(x)
beta <- rowSums(x)[1]
gamma <- colSums(x)[1]
rho <- x[1,1]*x[2,2] / (x[1,2]*x[2,1])
y<- f(beta, gamma, rho)
delta <- try(zapsmall(c(1, sqrt(crossprod(as.vector(x-y)))))[2])
if ("try-error" %in% class(delta)) cat("Error processing ", x, "\n")
delta
})
max(sim)
