\\ This file includes algorithms concerning complex continued fractions.

\\ The following matrices are generators of the Euclidean Bianchi groups:

V(x) = [x,0;0,1];
T(x) = [1,x;0,1];
U = [0,1;1,0];


\\ The action of the matrix at the boundary of upper half space.

matinfty(A)=
{
	my(x);
	if(A[2,1]==0,x = "infty";,x=A[1,1]/A[2,1];);
	x
}

matzero(A) =
{
	my(x);
	if(A[2,2]==0,x = "infty";,x=A[1,2]/A[2,2];);
	x
}

\\ The action of a matrix to elements of the closure of upper half plane.

matactboundary(M,x)=
{
	(M[1,1]*x + M[1,2])/(M[2,1]*x + M[2,2])
}

\\ The algorithm that we are going to employ is Hurwitz' nearest integer algorithm. 

\\ the following algorithm takes d (to specify the number field) and z (the complex number) as input and outputs the closest element of the lattice.

nearest(d,z) = 
{
	my(x,y,xf,xc,yf,yc,xr,yr,v);
	if(d == 1 
		, xf = floor(real(z)); xc = ceil(real(z)); yf = floor(imag(z)); yc = ceil(imag(z));
		if( ( (real(z)-xf)<= (xc-real(z)) ), x = xf; xr = real(z)-xf; , x = xc; xr = real(z)-xc; );
		if( ( (imag(z)-yf)<= (yc-imag(z)) ), y = yf; yr = imag(z)-yf; , y = yc; yr = imag(z)-yc;);
		v = [x + y*I, xr + yr*I];
		, 
		if(d == 2
		, xf = floor(real(z)); xc = ceil(real(z)); yf = floor(imag(z)/sqrt(2)); yc = ceil(imag(z)/sqrt(2));
		if( ( (real(z)-xf)<= (xc-real(z)) ), x = xf; xr = real(z)-xf; , x = xc; xr = real(z)-xc; );
		if( ( (imag(z)/sqrt(2)-yf)<= (yc-imag(z)/sqrt(2)) ), y = yf; yr = imag(z)-yf; , y = yc; yr = imag(z)-yc;);
		v = [x + y*sqrt(2)*I, xr + yr*sqrt(2)*I];
		,
		print("algorithm is not implemented yet!"); return;
		););
	v
}

CContFrac(d,z) =
{
	my( cf = List([]);, M,rem, a , word, v );
	M = [1,0;0,1];
	word = "";
	a = nearest(d,z);
	listput(cf,a[1]); word = concat(word,"T^("); word = concat(word,a[1]); word = concat(word,")*U*"); M = M*T(a[1])*U;
	if(a[2] == 0
	,  return;
	, rem = a[2]; while( abs(rem)>10^(-100), v = nearest(d,1/rem); rem = v[2]; listput(cf,v[1]); word = concat(word,"T^("); word = concat(word,v[1]); word = concat(word,")*U*"); M = M*T(v[1])*U; ););
	[cf,word,M]
}

CCFracToZ(v) = 
{
	my(z,l);
	l = length(v);z = v[l];
	for(i = 1,l-1, z = v[l-i] + 1/z; );
	print("The corresponding complex number is "z);
	z
}

CCFracToM(v) = 
{
	my(M,l);
	M = [1,0;0,1]; l = length(v);
	for(i = 1,l, M = M*T(v[i])*U;);
	M;
}

MatrixFixedPoint(M) = 
{
	my(D,r1,r2);
	D = (M[1,1] - M[2,2])^2 + 4*M[1,2]*M[2,1];
	r1 = (M[1,1] - M[2,2] + sqrt(D))/(2*M[2,1]);
	r2 = (M[1,1] - M[2,2] - sqrt(D))/(2*M[2,1]);
	print("The roots are " r1 " and "r2);
	[r1,r2]
}