function b = find_b(ld, Q, prec) % Q = -2*P
arguments
    ld
    Q
    prec = 32;
end
    digits(prec)
    C = ld^2 * Q^(1/(ld-1));
    syms x
    b1 = vpasolve(x^2 + C == x^(2*ld), x, [1 1+C/(2*ld-2)]);
    b = vpa(Q^(1/(2*ld-2))/b1);
end

% C = ld^2 * Q^(1/(ld-1));
% eqn = @(x) powm1(-C*x, 2*ld-2) .* x .* (1 + C*x).^2 - 1;
% b1 = fzero(eqn, 1/(2*ld-2));
% b = Q^(1/(2*ld-2)) * 1/(1 + C * b1);