efficiently approximate real solution for polynomial function

I want to efficiently solve a degree-7 polynomial in k.

For example, with the following set of 7 unconditional probabilities,

p <- c(0.0496772, 0.04584501, 0.04210299, 0.04026439, 0.03844668, 0.03487194, 0.03137491)

the overall event probability is approximately 25% :

> 1 - prod(1 - p)
[1] 0.2506676

And if I want to approximate a constant k to proportionally change all elements of p so that the overall event probability is now approximately 30%, I can do this using an equation solver (such as Wolfram Alpha), which may use Newton's method or bisection to approximate k in:

here, k is approximately 1.23:

> 1 - prod(1 - 1.23*p)
[1] 0.3000173

But what if I want to solve this for many different overall event probabilities, how can I efficiently do this in R?

I've looked at the function SMfzero in the package NLRoot, but it's still not clear to me how I can achieve it.

EDIT I've benchmarked the solutions so far. On the toy data p above:

Unit: nanoseconds
              expr     min      lq      mean  median      uq      max neval
 approximation_fun     800    1700    3306.7    3100    4400    39500  1000
       polynom_fun 1583800 1748600 2067028.6 1846300 2036300 16332600  1000
      polyroot_fun  596800  658300  863454.2  716250  792100 44709000  1000
         bsoln_fun   48800   59800   87029.6   85100  102350   613300  1000
        find_k_fun   48500   60700   86657.4   85250  103050   262600  1000

NB, I'm not sure if its fair to compare the approximation_fun with the others but I did ask for an approximate solution so it does meet the brief.

The real problem is a degree-52 polynomial in k. Benchmarking on the real data:

Unit: microseconds
              expr     min       lq       mean   median       uq     max neval
 approximation_fun     1.9     3.20     7.8745     5.50    14.50    55.5  1000
       polynom_fun 10177.2 10965.20 12542.4195 11268.45 12149.95 80230.9  1000
         bsoln_fun    52.3    60.95    91.4209    71.80   117.75   295.6  1000
        find_k_fun    55.0    62.80    90.1710    73.10   118.40   358.2  1000

Another option would be to just search for a root on a segment without specifically calculating polynomial coefficients. This can be done e.g. with the uniroot function.

Only one not-so-trivial thing we need to do here is to specify the segment. k is obviously >=0 - so that would be the left point. Then we know that all the k*p values should be probabilities, hence <=1. Therefore k <= 1/max(p) - that's the right point.

And so the code is:

find_k <- function(p, taget_op) {
  f <- function(x) 1 - prod(1 - x*p) - target_op
  max_k <- 1/max(p)
  res <- uniroot(f, c(0, max_k))
  res$root
}

p <- runif(1000, 0, 1)
target_op <- 0.3
k <- find_k(p, target_op)
k
#> [1] 0.000710281

1 - prod(1 - k*p)
#> [1] 0.2985806

^{Created on 2021-04-29 by the reprex package (v1.0.0)}

This works pretty fast even for 1000 probabilities.