Haskell, 67 bytes

a#b|h<-a/2,c<-sqrt$(1+h)^2-b,c>0=h+c:[h-c|h>c]|l<-(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 < e^a\$, which follows directly from \$\sum r_i=a\$ and \$\prod\left(r_i+1\right)=b\$ along with \$e^{r_i}>r_i+1\$. Since \$e^a\$ is the increasing limit of \$(1+a/N)^N\$ as \$N \to \infty\$, the bound \$b < e^a\$ implies that \$b < (1+a/N)^N\$, or equivalently that \$b^{1/N} < 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'} < 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}$$

Haskell, ^102… 95 bytes

• -6 bytes thanks to xnor (also, go upvote his amazing answer!).

a#b|m<-until(\m->(1+a/m)**m>=b)(+2)2=filter(>0)[a/m+sqrt((1+a/m)^2-b**(2/m))*(-1)**i|i<-[1..m]]

Try it online!

The relevant function is (#) which takes as input a and b as described in the statement.

Charcoal, 52 bytes

ＮθＮη≔¹ζＷ›ηＸ⊕∕θ⊗ζ⊗ζ≦⊕ζＦζＩΦＥ²⁺∕θ⊗ζ×∨λ±¹₂⁻Ｘ⊕∕θ⊗ζ²Ｘη∕¹ζκ

Try it online! Link is to verbose version of code. Explanation: Basically another port of @xnor's answer.

ＮθＮη≔¹ζ

Input \$ a \$ and \$ b \$ and start counting the number of repetitions of the output that we need.

Ｗ›ηＸ⊕∕θ⊗ζ⊗ζ≦⊕ζ

Keep incrementing (doubling also works of course, but it's not strictly necessary) the number of repetitions while \$ b > (1 + \frac a {2n}) ^ {2n} \$. (This explains why \$ b < e^a \$.)

Ｆζ

Repeat the output the appropriate number of times...

ＩΦＥ²⁺∕θ⊗ζ×∨λ±¹₂⁻Ｘ⊕∕θ⊗ζ²Ｘη∕¹ζκ

... output up to two non-zero solutions to the quadratic.

Jelly, 32 31 bytes

ḷ‘ɼ_’;N};1Ær½ßH}¥/}x®$ÆiṪṪƊ?,ḟ0

Try it online!

A full program taking two arguments, b and a in that order. Implements @xnor’s algorithm, so be sure to upvote that one too! Uses recursion to find real roots to a quadratic equation of the form \$x^2 - w x + (y - w - 1) = 0\$ where \$w = \frac{a}{2 ^ n}\$ and \$y = b ^ \frac{1}{2n} \$ for some non-negative integer \$n\$.

Explanation

ḷ‘ɼ                              | Increment the register, then revert to the original left argument (i.e. b)
   _                             | b - a
    ’                            | Subtract 1
     ;N}                         | Concatenate to -a
        ;1                       | Concatenate to 1
          Ær                     | Roots of polynomial with terms (b - a - 1, -a, 1)
                      ÆiṪṪƊ?     | Does the last root havecan imaginary component?
                ¥/}         ,    | - If yes, do the following using the original arguments to this link:
            ½ßH}                 |   - Call the current link recursively with the square root of the left argument and half the right argument
                   x®$           | - If not, repeat each term the number of times indicated by the register
                             ḟ0  | Filter out zeros

JavaScript (ES10), 93 bytes

This is really just a port of @xnor's answer.

Expects (a)(b).

a=>g=(b,r=1,x=(q=a/=2)+(d=(++q*q-b)**.5))=>x?Array(r).fill(a-d?[x,a-d]:x).flat():g(b**.5,r*2)

Try it online!