function ar = cot_formula(p)
ar = p*pi/2 - angle(agm(1,(1-i)*sqrt(sin(pi*p))*exp(i*p*pi/2)));
return;
%-----------------------------------------------------------
function m = agm(a,b)
% computes AGM for complex input taking the precautions from
% Borwein/Borwein p. 16.
m = b;
M = a;
while max(abs(M-m)/m) > eps
    m0 = sqrt(m*M);
    M = (m+M)/2;
    if abs(M-m0) > abs(M+m0)
        m = -m0;
    elseif abs(M-m0) == abs(M+m0)
        m = sign(imag(m0/M))*m0;
    else
        m = m0;
    end
end
return