\\ %TODO: form to matrix for general forms(the version for reduced forms is already present!) The obstacle here is: you do not know when you are at the main cycle. To understand this, one has to figure out a way to classify those forms on the main cycle. This last point is a work in progress.
\\ %TODO: halka acilmadan once bir suru kontrol koymak gerekli: mesela MtoF fonksiyonuna verilen matris beklentileri saglamayabilir!
\\ %TODO: IsAmbiguous konulabilinir.
\\ %TODO: IsReciprocal konulabilinir.
\\ %TODO: Have to write procedures to connect this theory to continued fractions.
\\ %TODO: IsEquivalent would be a nice function.
\\ %TODO: Have to put a control to the function CarkToForm: If the vector that this function eats is periodic, one has to be careful! The trace and discriminant fails, as well as the cark itself!
\\ %TODO: Reduce does not check for the indefiniteness. Put a control!
\\ %TODO: ...



\\ This code is designed neither for optimal speed, nor for optimal memory management. It just works on my desktop!!!

\\Let's define our fundamental matrices:
R = [1,-1;1,0];
S = [0,1;-1,0];
L = [1,-1;1,0];
LS= L*S;
L2S = L^2*S;
U = [0,1;1,0];
A = [-1,0;0,1];
B = [0,1;1,0];
C = [-1,1;0,1];
R2 = [1,0;2,-1];

\\ this is used to compute gcd of 3 integers
{
gcd3(x,y,z) = 
my(delta);
delta = gcd(y,z);
gcd(x,delta)
}
\\ this turns 4 integers to a form
MtoF(a,b,c,d) = 
{
	my(delta,S);
	delta = gcd3(c,d-a,-b);
	S = sign(a+d);
	print("( " S/delta*c " , " S/delta*(d-a) " , -" S/delta*b " ) is created");
	Qfb(S/delta*c,S/delta*(d-a),-S/delta*b)
}

\\ to ease our computations, here is conjugation

conjugate(M,N) = 
{
	N^(-1)*M*N
}

\\ this turns a matrix into a form.
MtoF(M) = 
{
	my(delta,S);
	delta = gcd3(M[2,1],M[2,2]-M[1,1],-M[1,2]);
	S = sign(M[1,1]+M[2,2]);
	print("( " S/delta*M[2,1] " , " S/delta*(M[2,2]-M[1,1]) " , -" S/delta*M[1,2] " ) is created");
	Qfb(S/delta*M[2,1],S/delta*(M[2,2]-M[1,1]),-S/delta*M[1,2])
}

\\ this is classical discriminant of a form. Can read out coefficients of a form and also vectors by component(f,i)!
{
disc(f) = component(f,2)^2 - 4*component(f,1)*component(f,3)
}

\\ usual trace and determinant functions are trace and matdet!!!
\\ this function is to be used in the reduction process

t(f)=
{
	my(res);
	if (abs(component(f,3))>= sqrt(disc(f)),res = sign(component(f,3))*floor((abs(component(f,3))+component(f,2))/(2*abs(component(f,3)))),res = sign(component(f,3))*floor((sqrt(disc(f))+component(f,2))/(2*abs(component(f,3)))));
	res
}

\\ this is obvious!!
IsReduced(f) =
{
	my(res);
	res = 0;
	if (abs(sqrt(disc(f))-2*abs(component(f,1)))<component(f,2) && component(f,2)<sqrt(disc(f)),res = 1);
	res
}

\\ this is obvious, too. I have to admit that I do not know any use of this!
IsZagierReduced(f) = 
{
	my(res);
	if (sqrt(disc(f))< component(f,2) && component(f,2)<sqrt(disc(f)) + 2*component(f,1) && disc(f)<component(f,2) && component(f,2)<sqrt(disc(f))+2*component(f,3),res = 1,res = 0);
	res
}
\\this is just a subroutine to be used in constructing the reduced form
OneStepReduction(f) = 
{
	my(i,a,b,c,M,N,w1,w,g);
	a = component(f,1);
	b = component(f,2);
	c = component(f,3);
	M = [a,b/2;b/2,c];
	word = "";
	N = S*(L*S)^(t(f));
	M = (N~)*M*N;
	w1 = Str(t(f)); 
	w = concat(" S(LS)^",w1); 
	print("Conjugation is done by: "w); 
	print("The new form is: " M); 
	g = Qfb(M[1,1],2*M[1,2],M[2,2]);
	[g,w]
}
OneStepReductionOhneErklarung(f) = 
{
	my(i,a,b,c,M,N,w1,w,g);
	a = component(f,1);
	b = component(f,2);
	c = component(f,3);
	M = [a,b/2;b/2,c];
	w = "";
	N = S*(L*S)^(t(f));
	M = (N~)*M*N;
	w1 = Str(t(f)); 
	w = concat(" S(LS)^",w1); 
	g = Qfb(M[1,1],2*M[1,2],M[2,2]);
	[g,w]
}


\\the result is a reduced form
Reduce(f) =
{
	my(g);
	g=f;
	until (IsReduced(g)==1,g = OneStepReduction(g)[1];);
	g
}

ReduceOhneErklarung(f) =
{
	my(g);
	g=f;
	until (IsReduced(g)==1,g = OneStepReductionOhneErklarung(g)[1];);
	g
}

\\ a routine to compute an automorphism of a given form f:
\\ important flaw of the above code is the fact that PARI/GP cannot always solve the pell equation, which is crucial in computing the automorphism. However, what you can do is, if the form is reduced, just look at the cark. The corresponding matric is nothing but the generator of the automorphism group!

aut (f) =
{
	my(g,S=vector(2));
	g = Qfb(1,0,-disc(f));
	S = abs(qfbsolve(g,4));
	S = [(S[1] - component(f,2)*S[2])/2, - component(f,3)*S[2];component(f,1)*S[2],(S[1] + component(f,2)*S[2])/2];
	S
}


\\ a sub-routine that gives the vector of left and right turns on the Farey graph!

VectorofReduction(f) =
{
	my(g,n,V,i);
	if (IsReduced(f) == 0, print("Given form is not reduced. Reducing it!"); 	g = ReduceOhneErklarung(f);, g=f;);
	V = [];
	until(g == f, V = concat(V,t(g)); g = OneStepReductionOhneErklarung(g)[1];);
	V
}

\\ feed the following with a reduced form and get the corresponding matrix!
\\ this version however works only for reduced forms!

FtoM(f) = 
{
	if (IsReduced(f) == 0, print("Given form is not reduced. Wait for a new version!"); break;);
	my(V,M);
	V = VectorofReduction(f);
	M = [1,0;0,1];
	for (i=1,#V, M = M * S * (L*S)^(V[i]););
	M
}

\\ this computes the action of the modular group on forms.

act(f,N) =
{
	my(M,K,a,b,c);
	a = component(f,1);
	b = component(f,2);
	c = component(f,3);
	M = [a,b/2;b/2,c];
	K = (N~)*M*N;
	Qfb(K[1,1],2*K[1,2],K[2,2])
}
\\ this is to write the cycle of the form f, using our theory.

Trace(f) = 
{
	my(M);
	M = FtoM(f);
	trace(M)
}
cark(f) =
{
	my(g,h,word,M,V);
	if(IsReduced(f)==0, print("The given form is not reduced. Reducing it:"); g=ReduceOhneErklarung(f);,g=f;);
	print("A reduced element in the class represented by " f " is " g);
	word = "";
	h = g;
	n = 1;
	until(g == h, word = concat(word, OneStepReductionOhneErklarung(h)[2]); h = OneStepReductionOhneErklarung(h)[1]; print("Now reducing " h););
	M = FtoM(g);
	V = VectorofReduction(g);
	print("The corresponding matrix is: " word " = " M);
	print("The trace of this class is: " trace(M));
	print("Discriminant of this class is: "disc(g));
	[h,M,V]
}


briefcark(f) =
{
	if(IsReduced(f)==0, g=ReduceOhneErklarung(f);,g=f;);
	word = "";
	h = g;
	n = 1;
	until(g == h, word = concat(word, OneStepReductionOhneErklarung(h)[2]); h = OneStepReductionOhneErklarung(h)[1];);
	M = FtoM(g);
	print("The corresponding matrix is: " word " = " M);
	[h,M]
}

\\ this is to provide information on the class group. Here, given an element we compute its powers.

PowersInClassGroup(f) = 
{
	my(n,i,j,V,W);
	n = qfbclassno(disc(f));
	if (n == 1, print("The class group is trivial; quitting!");break;);
	print("Working with discriminant: " disc(f));
	print("The class group has "n" elements!");
	V = vector(n,i, qfbpowraw(f,i););
	for(j=1, n, print("Cark of f^" j " is as follows:"); cark(V[j]););
	V;
}

\\ the following eats a vector and spits out a form. But, the vector is an integral 1x2n vector whose entries represent the Farey components. The inner boundary components are to be denoted by "-". 

CarkToForm(V)=
{
	my(i,j,f,l,M);
	l = #V;
	M = [1,0;0,1];
	if (l%2 == 0,,print("The given cark does not have even number of entries; not possible. Quitting!");return;);
	for (j=1,l-1, if(V[j]*V[j+1]<0,,print("Consecutive entries should have opposite sign(to be able to differentiate between the boundary components)"); return;));
	if (V[1]<0, for(i=1,l/2, M = M * (S*L^2)^(-V[2*i-1])*(S*L)^(V[2*i]););,for(i=1,l/2, M = M * (S*L)^(V[2*i-1])*(S*L^2)^(-V[2*i]);););
	if(V[1]<0,	f = MtoF(M); f = act(f,S);,f = MtoF(M););
	print("The form "f" is created!");
	print("The class has trace " trace(M) " and discriminant "disc(f)".");
	f
}

briefCarkToForm(V)=
{
	my(i,j,f,l,M);
	l = #V;
	M = [1,0;0,1];
	if (l%2 == 0,,print("The given cark does not have even number of entries; not possible. Quitting!");return;);
	for (j=1,l-1, if(V[j]*V[j+1]<0,,print("Consecutive entries should have opposite sign(to be able to differentiate between the boundary components)"); return;));
	if (V[1]<0, for(i=1,l/2, M = M * (S*L^2)^(-V[2*i-1])*(S*L)^(V[2*i]););,for(i=1,l/2, M = M * (S*L)^(V[2*i-1])*(S*L^2)^(-V[2*i]);););
	f = MtoF(M);
	f
}

\\ The code down below has a gap. Yes carks have to have even number of components, but they might add up to an odd number, which is in fact the case "half" of the time. Need to correct this.

List2Carks(n)=
{
	my(V,W,U,i,j,f);
	if (n%2 == 0,,print("The carks have to have even number of Farey components; not possible. Quitting!");return;);
	V = vector(n/2);
	W = matrix(n-1,2);
	for(i = 1,n/2, V[i] = partitions(n,n-1)[i];);
	for(j=1,n/2-1, W[2*j-1,] = [V[j][1],- V[j][2]]; W[2*j,] = [V[j][2],-V[j][1]];);
	W[n-1,] = [n/2,-n/2];
	for(i=1,n-1, print("Cark under consideration: " W[i,] "."); f = CarkToForm(W[i,]); print("The class group has: " qfbclassno(disc(f)) " element(s)."); cark(f); U = concat(U,f););
	U
}

\\ We need our own continued fraction algoritm in the end... Of course, there are things to do with this algoritm: it has to return a vector without zeroes, etc.

\\ the following gives the continued fraction expansion of the two fixed points of the hyperbolic element that determines the binary quadratic form given. Remember the continued fraction expansions are modified so that they fit nicely into carks. However the case by case treated continued fractions are somewhat unnatural!

cf(x,n) =
{ 
	my( a = vector(n), t );
	if (x>0, for (i = 1, n, a[i] = floor(x); t = x - a[i]; if (!t || i == n, break); x = 1 / t;); , for (i = 1, n, a[i] = ceil(x); t = x - a[i]; if (!t || i == n, break); x = 1 / t;); );
	a;
}

\\ The following is the minus continued fraction of Zagier, and will hopefully be a better suited option.

cfminus(x,n) =
{ 
	my( a = vector(n), t );
	for (i = 1, n, a[i] = ceil(x); t = a[i] - x; if (!t || i == n, break); x = 1 / t;);
	a;
}

\\ the following is nearest continued fraction expansion.

cfnear(x,n) =
{ 
	my( a = vector(n), t, d,e );
	if (x>0, for (i = 1, n, d= x-floor(x); e=ceil(x)-x; if(d<=e,a[i] = floor(x); t = x - a[i]; if (!t || i == n, break); x = 1 / t;, a[i]=ceil(x);t=x-a[i]; if (!t || i == n, break); x = 1 / t; ););, for (i = 1, n, a[i] = ceil(x); t = x - a[i]; if (!t || i == n, break); x = 1 / t;); );
	a;
}


\\ The following gives the "minus" continued fraction associated to a cark.

CtoCFM(V) =
{
	my(res);
	res = V[#V];
	for(i = 1, #V-1, res = V[#V-i] - 1/(res); print("Result is now " res););
	res;
}

\\ %TODO: It is a good idea also to display the roots not numerically but algebraically!

\\ For some reason that I do not understand the following does not work!

FixedPointAll(f) = 
{
	my(x01,x00,x1,x2);
	x1 = (-component(f,2) - sqrt(disc(f)))/(2*component(f,1));
	x2 = (-component(f,2) + sqrt(disc(f)))/(2*component(f,1));
	if (x1>x2, x01 = x1; x00 = x2, x01 = x2; x00=x1;);
	print("The larger root is "x01" having the following continued fraction expansion : " cf(x01,40));
	print("The larger root is "x01" having the following MINUS continued fraction expansion : " cfminus(x01,40));
	print("The smaller root is "x00" having the following continued fraction expansion : " cf(x00,40));
	print("The smaller root is "x00" having the following MINUS continued fraction expansion : " cfminus(x00,40));
	[cf(x01,20),cf(x00,20)];
}

FixedPoint(f,n) = 
{
	my(x01,x00,x1,x2);
	x1 = (-component(f,2) - sqrt(disc(f)))/(2*component(f,1));
	x2 = (-component(f,2) + sqrt(disc(f)))/(2*component(f,1));
	if (x1>x2, x01 = x1; x00 = x2, x01 = x2; x00=x1;);
	print("The larger root is "x01" having the following PLUS continued fraction expansion : " cf(x01,n));
	print("The smaller root is "x00" having the following PLUS continued fraction expansion : " cf(x00,n));
	[cf(x01,n),cf(x00,n)]
}

FixedPointMinus(f,n) = 
{
	my(x01,x00,x1,x2);
	x1 = (-component(f,2) - sqrt(disc(f)))/(2*component(f,1));
	x2 = (-component(f,2) + sqrt(disc(f)))/(2*component(f,1));
	if (x1>x2, x01 = x1; x00 = x2, x01 = x2; x00=x1;);
	print("The larger root is "x01" having the following MINUS continued fraction expansion : " cfminus(x01,n));
	print("The smaller root is "x00" having the following MINUS continued fraction expansion : " cfminus(x00,n));
	[cf(x01,n),cf(x00,n)]
}

\\ The following two routines compute the exact value of the given continued fraction.

DecodeCFPLUS(V) =
{
	for (i = 2, #V, if (V[i] <= 0,print("There cannot be any non-positive entry(ies). Quitting!");return;););
	my(res);
	res = V[#V];
	for (i = 1, #V-1, res = V[#V-i] + 1/res);
	print(V " is equal to " res);
	res
}

DecodeCFMINUS(V) =
{
	for (i = 2, #V, if (V[i] <= 0,print("There cannot be any non-positive entries. Quitting!");return;););
	my(res);
	res = V[#V];
	for (i = 1, #V-1, res = V[#V-i] - 1/res);
	print(V " is equal to " res);
	res
}

\\ Now I start reduction theory. For this, I need the value of a form at the given lattice points.

Eval(f,x,y)=
{
	component(f,1)*x^2+ component(f,2)*x*y + component(f,3)*y^2
}

\\ investigate the lengths of carks of non-squarefree discriminants. The input is a cark and the type of the matrix which determines those forms that lie above the corresponding cark. A is of very special type, given in Buell's book.

A(p) = [1,0;0,p]

Findprimes(C,A,n) = 
{
	my(m,l);
	l = length(C);
	if(1 == Mod(l,2),print("Cark has to be of even length. Quitting!");return;,);
	W = [1,0;0,1];
	if (C[1]<0, for(i=1,l/2, W = W * (L2S)^(-C[2*i-1])*(LS)^(C[2*i]););,for(i=1,l/2, W = W * (LS)^(C[2*i-1])*(L2S)^(-C[2*i]););print(W););
	M = A^(-1)*W^n*A;
	m = M[2,1]*p;
	m = subst(m,p,0); \\ in the above line m is of type t_POL, this way we make it into an integer again! Otherwise factorization does not work.
	factor(m)
}

Nonsquarefree(C,F,A,n) = 
{
	my(W,M,l,r,g);
	l = length(C);
	f = briefCarkToForm(C);
	W = [1,0;0,1];
	if (C[1]<0, for(i=1,l/2, W = W * (L2S)^(-C[2*i-1])*(LS)^(C[2*i]););,for(i=1,l/2, W = W * (LS)^(C[2*i-1])*(L2S)^(-C[2*i]););print(W););
	M = A^(-1)*W^n*A;
	r = matsize(F)[1];
	for(i=1,r, print("The divisors are" F); print("For the prime divisor " F[i,1] " the cark lying above " f " is "); g = MtoF(subst(M,p,F[i,1])); briefcark(g);)
}

Listlifts(C,A,n) = 
{
	my(F);
	for(i = 2,n, print("Now computing primes at " i "th power: "); F = Findprimes(C,A,i); Nonsquarefree(C,F,A,i);)
}

Listprimelifts(C,A,n) = 
{
	my(F,P);
	P = primes(n);
	for(i = 1,length(P), print("Now computing primes at " P[i] "th power: "); F = Findprimes(C,A,P[i]); Nonsquarefree(C,F,A,P[i]);)
}

\\ We now relate the cark to the Picard group. For this, we will need some routines. 
\\ The first routine will admit a form as an argument and produce all spinal forms in the corresponding form.

ListSpinalForms(f) = 
{
	my(fR,Spine = [];);
	if(component(f,1)*component(f,3)>0, print("f is not a spinal form. Reducing it!"); fR =ReduceOhneErklarung(f);, fR = f; );
	until(fR == f, fR = act(fR,LS); if(component(fR,1)*component(fR,3)>0, fR = act(fR,S*LS););  Spine = concat(Spine,fR);  );
	for(i = 1,length(Spine), Spine = concat(Spine, act(Spine[i],S)););
	Spine
}

\\ This is a subroutine that that takes a list of forms as input and returns a vector which contains forms which are acted by LS and L2S

VectorAct(V) = 
{
	my(W);
	W = List();
	for(i = 1,#V,W = concat(W, List([act(V[i],LS),act(V[i],L2S)])));
	W
}

\\ The following takes an integer and a form as input, which is assumed to be the root of a Farey branch. The output is a list of forms having distance n+1 to the spine.

FormsOnFareyBranch(root,n) = 
{
	my(V,W);
	V = List(root); 
	for(i = 1,n, W = List(); for(j= 2^(i-1), #V, W = concat(W,List(V[j])););  V = concat(V,VectorAct(W);););
	V
}

\\ the following routine lists forms up to distance n. However, the ones with S are not inlcluded, in order not to complicate the picture.

ListLevelForms(f,n) = 
{
	my(fR,Spine = List();,root);
	if(component(f,1)*component(f,3)>0, print("f is not a spinal form. Reducing it!"); fR =ReduceOhneErklarung(f);, fR = f; );
	until(fR == f, fR = act(fR,LS); if(component(fR,1)*component(fR,3)>0, root = fR; fR = act(fR,S*LS);, root = act(fR,S*LS));  Spine = concat(Spine,List(fR)); Spine = concat(Spine,FormsOnFareyBranch(root,n))  );
	Spine
}

ListFacialFormswithoutS(f,n) =
{
	my(Spine=[f]);
	for(i = 1,n, Spine = concat(Spine,act(f,LS^i)); Spine = concat(Spine,act(f,(S*L^2)^(i))); );
	Spine
}

\\The following routine admits a form as its argument and produces a list of spinal points in the Picard group. Notice however that the points are represented in complex coordinates!
SpinalPointsinPic(f) = 
{
	my(points = [];,forms,a,b,c,d);
	forms = ListSpinalForms(f);
	d=disc(f);
	for(i = 1,length(forms), a = component(forms[i],1);b = component(forms[i],2);c = component(forms[i],3); points = concat( points, ( abs(b - sqrt(d)))/( 2 * sqrt(abs(a*c))) + I*( abs(b + sqrt(d)))/(2 * sqrt(abs(a*c)))););
	points
}

\\ The following routine takes a list of forms and computes their images in the picard group. The forms are assumed to be equivalent.

PointsinPic(forms) = 
{
	my(points = [];,a,b,c,d);
	d=disc(forms[1]);
	for(i = 1,length(forms), a = component(forms[i],1);b = component(forms[i],2);c = component(forms[i],3); points = concat( points, ( abs(b - sqrt(d)))/( 2 * sqrt(abs(a*c))) + I*( abs(b + sqrt(d)))/(2 * sqrt(abs(a*c)))););
	points
}


\\

WriteSpinalPointsinPic(f) = 
{
	write(picard,ListSpinalForms(f));
	write(picard,SpinalPointsinPic(f));
}


\\ The following routine admits a list of forms as input and calculates the corresponding points in the Picard group.

WritePointsinPic(forms) = 
{
	my(points = [];,a,b,c,d);
	d = disc(forms[1]);
	for(i = 1,length(forms), a = component(forms[i],1);b = component(forms[i],2);c = component(forms[i],3); points = concat( points, ( abs(b - sqrt(d)))/( 2 * sqrt(abs(a*c))) + I*( abs(b + sqrt(d)))/(2 * sqrt(abs(a*c)))););
	write(picard,forms);
	write(picard,points);
}

\\ The following is an auxiliary function that helps to compute the forms on a face when one feeds it with LS*L2S^n and L2S^LS^n. IOne may regard the root form (a,b,c) to be a variable, as well.

f(M) = 
{
	my(p,var);
	var = M*[x;y]; 
	a*var[1,1]^2 + b*var[1,1]*var[2,1]+c*var[2,1]^2
}

g(M) = 
{
	my(p,var);
	var = M*[x;y]; 
	(a+b+c)*var[1,1]^2 + (2*a+3*b+4*c)*var[1,1]*var[2,1]+(a+2*b+4*c)*var[2,1]^2
}

\\ What follows is related to Hecke continued fraction theory
\\The followings are standard matrices to use in computations.

Lq(q) = 
{
	my(K);
	K = [2*cos(Pi/q),-1;1,0];
	K
}

Tq(q) = 
{
	my(K);
	K = [1,2*cos(Pi/q);0,1];
	K
}

l(q) = 
{
	my(k);
	k = 2*cos(Pi/q);
	k
}

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
}

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


\\ The followings are related to the Hecke reduction algorithm of ours, not Rosen's. There are certain simplifications of the algorithm. This way, however, it is easier to read. At the moment, the algorithm works only for positive real numbers. Negatives can be obtained by a very slight modification, which will be dealt with later....

\\ TODO: do the negative integer case, too...

\p 200

heckecf(x,q) =
{
	my(M, diff = List([]),int = List([]), v = List([]),w = List([]),y,n,j,	word, N);
	for(i=1,q-1,listput(int,matzero(Lq(q)^i*S)));
	word = ""; y=x; N = [1,0;0,1];
	while(y > 10^(-30), print("y is " y);
			if( y>l(q), 
					n=0;
					while(((y-n*l(q))*(y-(n+1)*l(q))>0), n=n+1;); 
					listput(w,(Lq(q)*S)^n); 
					y = matactboundary((S*Lq(q)^(q-1))^n,y);
					for(i = 1,length(w),N = N*w[i]); print("the convergent at this step is " matzero(N)); N = [1,0;0,1]; 
					word = concat(word,"(Lq*S)^"); word = concat(word,n);
					print("the word is " word);
					,
					j=2; diff=List([]);
					for(i=1,q-1,listput(diff,y-int[i]);); 
					while((diff[j]*diff[j-1]>10^(-50)),j=j+1;);
					listput(w,(Lq(q)^(j)*S)); 
					for(i = 1,length(w),N = N*w[i]); print("the convergent at this step is " matzero(N)); N = [1,0;0,1];
					y = matactboundary(S*Lq(q)^(q-j),y);
					word = concat(word,"Lq^"); word = concat(word,j); word = concat(word,"S");
					print("the word is "word);
					);
		); 
	M = [1,0;0,1];
	for(i = 1, length(w), M = M*w[i];);
\\	convergents(w);
	[M,w,word]
}

\\ the following is the fast hecke cf algorithm. It contains also a control on the accuracy of the algorithm, the 10^-80 part...

fastheckecf(x,q) =
{
	\\ put a control to see you've got a positive number. We'll deal with negative integers later...
	my(int , y , A, word, M, mat, count);
	int=List([]);word = ""; M = [1,0;0,1]; mat = List([]); count = 1; y=x;
	listput(int,0);for(i=2,q-1,listput(int,matzero(Lq(q)^(q-i)*S));); \\ This is the list of intervals. Note that the entries in this list are increasing.
	\\ for(i = 1,q-1,if(y-int[i]<10^(-100),print("This is trivial!");print("The factorization is Lq^"q-i"*S");return;);); \\ this line checks whether you've hit the jackpot! 
	print("intervals are " int);
	while( ((abs(y)) > 10^(-105)) && ((abs(y))<10^(105)), 
				A = onestepheckecf(y,q); y = A[1]; word = concat(word,A[2]); M = M*(Lq(q)^(A[3])*S)^(A[4]); listput(mat,Lq(q)^(A[3])*S); print("Step " count " is done!"); count +=1; print("A is "A);print("convergents are \n"matzero(M) "\nand\n " matinfty(M));
			);
	print("M is " M " , word is " word);
	[M,word,mat]
}

onestepheckecf(y,q) = 
{
	my(int,x, p, diff, word,n);
	x=y;p = 1;word = "";n=1;int = List([]); diff = List([]);
	listput(int,0);
	for(i=2,q-1,listput(int,matzero(Lq(q)^(q-i)*S));); \\ we know that the entries in this list is increasing.
	for(i = 1,q-1,listput(diff,y-int[i]);); \\ create a new list called diff containing the differences. As a result of the previous two lines, this list is decreasing. 
	while(diff[p]*diff[p+1]>0, p+=1; if(p==q-1,break;)); \\ finds the interval in which y falls. This means we know which element to act on y. 
	if(p==1, 
		n = numberofrightturns(q,x); 
		word = concat(word,"*(Lq^"); word = concat(word,q-1); word = concat(word,"*S)^"); word = concat(word,n); 
		x = matactboundary((S*Lq(q))^n,x);
		, 
		if(p == q-1 ,
			n = numberofleftturns(q,x); 
			word = concat(word,"*(Lq*S)^"); word = concat(word,n); 
			x = matactboundary((S*Lq(q)^(q-1))^n,x);
			, 
			word = concat(word,"*Lq^"); word = concat(word,q-p); word = concat(word,"*S");
			x = matactboundary(S*Lq(q)^p,x);
			);
		); 
	[x,word,q-p,n]
}

numberofleftturns(q,y) = 
{
	my(n);
	n=0;
	while(((y-n*l(q))*(y-(n+1)*l(q))>0), n=n+1;); 
	n
}

numberofrightturns(q,y) = 
{
	my(n);
	n=0;
	while(((1/y-n*l(q))*(1/y-(n+1)*l(q))>0), n=n+1;); 
	n
}

\\This checks, to high presicision (meaning 200 digits), whether your matrix works for the given number with increased precision.

numcheck(M,r) = 
{
	my(x,y);
	x = (matzero(M) - r); y = (matinfty(M) - r);
	if( (abs(x)>10^(-200)) && (abs(y)>10^(-200) ), print("your approximation is bad!");, 
		if( (abs(x)<=10^(-200)), print("You have a good approximation at 0");
			,print("You have a good approximation at infinity");
		);
	);
}

\\ %TODO: an algebraic check is also possible... You must be careful about the arguments though.


\\ The following routine eats the sequence of matrices and spits out the convergents, numerically of course. 

convergents(A) = 
{
	my(M);
	M = [1,0;0,1];
	for(i = 1,length(A),M = M*A[i]; print("the "i"th convergent is " matzero(M)))
}

\\ The following gives the algebraic version of end points.

algLq(a) = [a,-1;1,0];

algboundary(q) =
{
	my(M,vec);
	vec = listcreate();
	for(i = 1,q-1, listput(vec,matzero(algLq(a)^i*S)););
	vec
}
numericalboundary(q) =
{
	my(M,vec);
	vec = listcreate();
	for(i = 1,q-1, listput(vec,matzero(algLq(l(q))^i*S)););
	vec
}

\\ the following two are functions that give the maximum and minimum of lists. 

maximum(list)=
{
	my(m,index);
	m=list[1];index = 1;
	for(i=2,length(list),if(m<=list[i],m=list[i]);index = i;);
	print("The longest interval is " index "th.");
	m
}

minimum(list)=
{
	my(m,index);
	m=list[1];index = 1;
	for(i=2,length(list),if(m>=list[i],m=list[i]);index = i;);
	print("The shortest interval is " index "th.");
	m
}

convergence(boundary) = 
{
	my(diff);
	diff = List([]);
	for(x = 1,length(boundary)-1,listput(diff,boundary[x]-boundary[x+1]););
	print("the maximum length is " maximum(diff));
	print("the minimum length is " minimum(diff));
}

\\ This matrix gives the numerical boundary of intervals.

Nq(q,l) = [sin(Pi*(l+1)/q),-sin(Pi*l/q);sin(Pi*l/q),-sin(Pi*(l-1)/q)]/sin(Pi/q)


\\  writetex(filename,"heckecf(2,9)")
\\  writetex(filename,heckecf(2,9))
\\  are the two commands that we use to write results to a file. It is useful whenever the results are very long...


\\ these are calculations corresponding to the latest conjecture about the trace of the element and the fundamental unit.

integerpartoffundamentalunit(D) = 
{
	my(s,x);
	if(isfundamental(D)==0,s = 0; return(0););
	x = quadunit(D);
	if(norm(x)==1, if (Mod(D,4)==1, s =  real(x) + imag(x)/2; , s = real(x););, if( Mod(D,4)==1,  s = real(x^2) + imag(x^2)/2; , s = real(x^2); ); );
	s
}
quadraticpartoffundamentalunit(D) = 
{
	my(x,s);
	if(isfundamental(D)==0,return(0););
	x = quadunit(D);
	if(norm(x)==1, if (Mod(D,4)==1,s =  imag(x)/2;, s = imag(x);), if(Mod(D,4)==1,  s = imag(x^2)/2;, s = imag(x^2)); );
	s
}

testconjecture(D) = 
{
	my(f,t,a);
	if(isfundamental(D) ==0, return(0););
	f = decidetheform(D);
	t = Trace(f);
	a = integerpartoffundamentalunit(D);
	if(abs(t) - 2*abs(a)==0,,print("For D " D " we have an exception!!!"));
\\	print("Trace is " t " integer part of fundamental unit is " a);
\\	print("Difference is " t^2 - 4 - 4*(a^2 - 1 ) );
}

decidetheform(D)=
{
	my(f);
	if(isfundamental(D) ==0, return(0););
	if(qfbclassno(D)==1, if(Mod(D,4) == 1, f = ReduceOhneErklarung(Qfb(1,1,-(D-1)/4)); , f = ReduceOhneErklarung(Qfb(1,0,-D/4));) , f=ReduceOhneErklarung(quadclassunit(D)[3][1]););
	f
}

decidethefactor(D) =
{
	my(n,x);
	n=1;
	if(isfundamental(D) ==0, return(0););
	until(issquare(n^2*D+4),n = n+1;);
	x = n^2*D+4;
	[n,SQRT(x)]
}

SQRT(n) = 
{
	my(k,s,M);
	s = 1;
	if(n==1,return(1);break;, M = factor(n););
	if(issquare(n)==0, s = sqrt(n); , k = matsize(M)[1]; for(i = 1,k, s = s*M[i,1]^(M[i,2]/2););  );
	s
}

funddisc(n)=
{
	my(V);
	V = Vec(5);
	for(i = 6,n, if(isfundamental(i)==1,V = concat(V,i);););
	V
}

testtable(D)=
{
	my(V,W,f,t);
	V= Vec(D);
	V = concat(V,decidethefactor(D));
	W = Vec([integerpartoffundamentalunit(D),quadraticpartoffundamentalunit(D)]);
	V = concat(V,W);
	if(Mod(D,4)==1,W = Vec([W[1]^2+D*W[2]^2,2*W[1]*W[2]]);, W = Vec([W[1]^2+(D/4)*W[2]^2,2*W[1]*W[2]]););
	V = concat(V,W);
	f = decidetheform(D); t = Trace(f);
	if(Mod(D,4)==1, V = concat(V,[(t^2-2)/2,abs(t)*SQRT((t^2-4)/D)/2]);, V = concat(V,[(t^2-2)/2,abs(t)*SQRT((t^2-4)/D)]););
	V
}