Set the sum and product Given two positive reals \$a\$ and \$b\$, output some positive reals \$r_i\$, such that \$\sum r_i=a\$ and \$\prod\left(r_i+1\right)=b\$. You can assume that it's possible. You can also assume that your float type have infinite precision.
Test cases:
2,3 =&gt; 2
2,4 =&gt; 1,1 or 1/2,(sqrt(57)+9)/12,(9-sqrt(57))/12 or etc.
2,5 =&gt; 1/3,1/2,2/3,1/2 or etc.
2,8 =&gt; (undefined behavior)
2,2 =&gt; (undefined behavior)
e,2e =&gt; 1,e-1 or etc. (e is natural base, though any number&gt;1 work)
3,e^3 =&gt; (undefined behavior) (It's possible to get as close as you want to e^3, but not reaching)

Shortest code wins.
Notes

Given the assumption with infinite precision, you can't solve arbitrary equation if your language doesn't have such functions(you can still use functions like Solve in Mathematica to solve such). In this case, some forms may have solution mathematically but you can't work it out, e.g. \$p,p,p,...,p,q\$ where \$p\$ and \$q\$ are reals. (At least software don't provide exact solution for x*x*x*y=5,3*x+y=7)
Haskell, 67 bytes

a#b|h&lt;-a/2,c&lt;-sqrt$(1+h)^2-b,c&gt;0=h+c:[h-c|h&gt;c]|l&lt;-(a/2)#sqrt b=l++l

Try it online!
Based on this solution of Delfad0r, who also saved 1 byte off this
We first try to find a two-element solution with \$r_1+r_2=a\$ and \$(r_1+1)(r_2+1)=b\$, which we do by solving a quadratic equation. If we get real solutions, we're done. (If one solution is zero, we remove it and output a singleton.)
If this quadratic isn't solvable in the reals, we replace \$a \to a/2\$ and \$b \to \sqrt{b}\$ and solve from there. If we find a solution to this, then duplicating the list (that is, repeating each element a second time) gives a solution to the original problem, since this doubles the sum and squares the product of numbers-plus-one. We recursively make these substitutions until we simplify \$(a,b)\$ to ones that give a two-element solution.
We show that this process always reaches a solution. For the problem to be solvable, we must have \$b &lt; e^a\$, which follows directly from \$\sum r_i=a\$ and \$\prod\left(r_i+1\right)=b\$ along with \$e^{r_i}&gt;r_i+1\$. Since \$e^a\$ is the increasing limit of \$(1+a/N)^N\$ as \$N \to \infty\$, the bound \$b &lt; e^a\$ implies that \$b &lt; (1+a/N)^N\$, or equivalently that \$b^{1/N} &lt; 1+a/N\$, for sufficiently large \$N\$. By choosing \$N=2^n\$ as a large power of \$2\$, doing \$n-1\$ steps of the substitution \$a \to a/2\$ and \$b \to \sqrt{b}\$ gives us \$a' = 2a/N\$ and \$b' = b^{2/N}\$, which we've shown satisfy \$\sqrt{b'} &lt; 1 + a'/2\$.
This is exactly what we need to have real solutions for \$r_1+r_2=a\$ and \$(r_1+1)(r_2+1)=b\$, which gives a quadratic with discriminant \$(1+a/2)^2-b\$. It solutions are $$r  = a/2 \pm \sqrt{(1+a/2)^2-b}$$