function[aa,bb,gg]=double3(a,b,g)
%function[aa,bb,gg]=double3(a,b,g)
%   Raw doubling algorithm for repeatedly squaring a symplectic matrix
%   that is represented in the form
%
%          |     a(1)^(-1)          [a(1)^(-1)]*b(1)         |
%       Z =|                                                 |
%          |   g(1)*a(1)^(-1)    a(1)'+g(1)*[(1)^(-1)]b(1)   |
%
%   The inputs are a, b, and g.  The algorithm iterates to convergence
%   of g(j) in the following recursions:
%
%       a(j+1) = a(j)*[(I + b(j)*g(j))^(-1)]*g(j)
%
%       g(j+1) = g(j)+ a(j)'*g(j)*[(I + b(j)*g(j))^(-1)]*a(j)
%
%       b(j+1) = b(j) + a(j)*[(I + b(j)*g(j))^(-1)]*b(j)*a(j)'
%
%   The "limiting values" of a, b, and g are returned in aa, bb, and gg.
%   Note that even though gg converges, aa wil not.
a0=a;
b0=b;
g0=g;
tol=1e-15;
dd=1;
ss=max(size(a));
v=eye(ss);
while dd>tol
a1=a0 *((v+b0*g0)\a0);
b1=b0+a0*((v+b0*g0)\(b0*a0'));
g1=g0+a0'*g0*((v+b0*g0)\a0);
dd=max(max(abs(g1-g0)));
a0=a1;
b0=b1;
g0=g1;
end
aa=a1;bb=b1;gg=g1;
