function [p,q] = adi_para_sylv(a,b,c,d,n)
% adi_para_sylv() computes shift parameters pj, qj for the general
% two-paramter ADI minimax problem
%
% min_{p,q} max{x in [a,b], y in [c,d]} prod_j|(x-qj)(y-pj)/(x+pj)(y+qj)|
%
% to be used in the Sylvester-ADI method for 
%  AX + XB + FG' = 0 and AXC + MXB + FG' = 0.
% The generated parameters for A=A', M=M' with spec(A,M) in [a,b]
% and B=B', C=C' with spec(B,C) in [c,d].
% Both spectra should be either in the left or right open real axis.
%
% Reference: Algorithm 2.1 in [J. Sabino. Solution of Large-Scale Lyapunov Equations via the 
% Block Modified Smith Method. PhD thesis, Rice University, Houston, Texas, June 2007. 
% Available from: https://hdl.handle.net/1911/102054.]
%       
%
% INPUT:
% a,b,c,d - spectral intervals [a,b] and [c,d] 
% n - number of desired parameters, integer >= 1
%
% OUTPUT:
% p,q - n by 1 vectors of the parameters
%
sgn = sign(a);
if (((b/a) > 1e8) || ((d/c) > 1e8)) %simplification for large cond.number
    p = logspace(log10(a),log10(b),n);
    q = logspace(log10(c),log10(d),n);
else
    mhat = 2*(b-a)*(d-c)/(a+c)/(b+d);
    kp = 1/(1 + mhat + sqrt(mhat*(mhat+2)));
    if (a==c) && (b==d) % essentially Lyapunov-case
        g = 0;
    else
        g = 2*(kp*(b+d)-(a+c))/((a+c)*(b-d)+kp*(b+d)*(c-a));
    end
    f = (2+g*(b-d))/(b+d);
    h = kp*(c-a+2*a*c*g)/(a+c);
    m = 1-kp^2; 
    K = ellipke(m); % complete ellipt. integr. of 1st kind
    u = (2*[1:n]-1)/(2*n) * K; 
    [~,~,w] = ellipj(u,m); % Jacobi elliptic function
    p = (h+w)./(w*g+f);
    q = (h-w)./(w*g-f);
end
p = sgn*p;
q = sgn*q;
