function a = find_a(ld, Q, prec) % Q = -2*P
arguments
    ld
    Q
    prec = 32;
end
    digits(prec)
    C = Q * ld^(2*ld-2);
    syms x
    a = vpasolve(C * (1-x)^ld == x, x, [0 C]);
end

% C = Q * ld^(2*ld-2);
% eqn = @(x) ld * hypergeom([-ld+1 1], 2, C * (1 - C * x)) .* (1 - C * x) - x;
% a1 = fzero(eqn, 0.0);
% a = C * (1 - C * a1);