function  [gamma,alpha,n1vec] = s4vecp(a,b,c,d);
%   Subroutine for SKEWED-MU power method.

% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: 1.1.6.2 $ $Date: 2004/08/17 21:56:05 $

    na = norm(a);
    nb = norm(b);
    nc = norm(c);
    nd = norm(d);
    beta4 = nd^2;
    beta3 = 2 *real(c'*d);
    beta2 = nc^2 + nb^2;
    beta1 = 2*real(a'*b);
    beta0 = na^2 - 1;
    
    if na<1
       %[beta4 beta3 beta2 beta1 beta0]
        tmp = LOCALroots([beta4 beta3 beta2 beta1 beta0],5);
        %tmp = roots([beta4 beta3 beta2 beta1 beta0]);
        [ind] = find( imag(tmp)==0 & tmp>0 );
        if isempty(ind)
            SW= warning('off','all');
           %disp('In Bisection')
             glow = 0;
             ghigh = min([1e6 (na+1)/nb 1/nc (nc+sqrt(nc^2+4*nd^2))/2/nd])^2;
             cnt = 0;
             while (ghigh-glow>1e-6)
                gtry = 0.5*(ghigh+glow);
                sg = sqrt(gtry);
                if norm([a+sg*b;sg*c+gtry*d])<1
                   glow = gtry;
                else
                   ghigh = gtry;
                end
                cnt = cnt + 1;
             end
             beta = sg;
             n1vec = [a+sg*b;sg*c+gtry*d];
             alpha = gtry;
             gamma = 1;
            warning(SW);
        else
            beta = tmp(ind(1));
            n1vec = [a+beta*b;beta*(c+beta*d)];
            gamma = 1;
            alpha = beta^2;
            % else
        %   beta = 100000;
        %    n1vec = [a+beta*b;beta*(c+beta*d)];
        %   gamma = norm(n1vec);
        %    alpha = beta^2;
        % end
        end
    else
        gamma = 1/na;
        alpha = 0;
        n1vec = [a/na;0*d];
    end

function r = LOCALroots(c,n)

%leadz = find(abs(c)<1e-10*max(abs(c)));
leadz = find(abs(c)==0);

if isempty(leadz)
   a = diag(ones(1,n-2),-1);
   a(1,:) = -c(2:n) ./ c(1);
   r = eig(a);
else
   len = length(leadz);
   strp = find((1:len)==leadz);
   c(strp) = [];
   n = length(c);
   a = diag(ones(1,n-2),-1);
   a(1,:) = -c(2:n) ./ c(1);
   r = eig(a);
end