\\ Copyright (c) 2003 Henri Darmon and Adam Logan
/* Given an element of P^1(Q(sqrt(29)), gives a path from it to [1,0], the
cusp at infinity, by way of adjacent cusps.  Method: the Euclidean algorithm.
We just need to look at nearby lattice points, luckily.  */

/* First a function to calculate a continued fraction for such an element. */

cf(v)= 
	{
	local(x,y,z);
	if (v[2]==0,return([]),quotrem=divide(v);x=quotrem[1];y=cf([v[2],quotrem[2]]);z=concat(x,y);return(z));
	}

addhelp(cf,"cf(v) gives a continued fraction of m/n, where v=[m,n].  Assumes that m and n are quads in Q(sqrt(29)), Q(sqrt(37)), or Q(sqrt(41)); otherwise may not work");
/* Now the problem is the division.  We just have to get something of smaller
norm, so if we have a smaller numerator already, don't bother.  Otherwise
divide and look at nearby lattice points.... This was originally written for
Q(sqrt(29)), but amazingly it should work with only minor modifications for
Q(sqrt(37)) and Q(sqrt(41)).

This function takes a two-component vector and returns another such, namely
[quotient, remainder]. */

divide(v)=
	{
	local(a,b,c,r,s,ccomps,rcomps,closest);
	if(v[2]==0,error("Division by zero in reduce!"),);
	a=v[1];b=v[2];
	s=abs(norm(b));
	if((abs(norm(a))<s),return([0,a]),);
	/* Don't bother if the numerator already has smaller norm. */
	c=a/b;ccomps=[component(c,2),component(c,3)];
	/* ccomps now holds the coefficients of 1 and w of a/b. */
	rcomps=round(ccomps);r=rcomps[1]+rcomps[2]*w;
	/* r is our first guess as to the quotient. */
	closest=r;
	for(i=-2,2,for(j=-2,2,d=r+i+j*w;if(abs(norm(a-d*b))<abs(norm(a-closest*b)),closest=d,)));
	if(abs(norm(a-b*closest))>=s,error("Problem using Euclidean algorithm"),);
	return([closest,a-b*closest])
	}

addhelp(divide,"divide(v): where v = [m,n], returns [q,r] with m=qn+r and |norm(r)| < |norm(n)|");

/* Now a function that produces the nth partial quotient coming from
a list, presumably one representing a continued fraction.  Then one
giving them all.*/

pquotn(v,n)=
	{
	local(s,u);
	u=v[n];
	forstep(s=n-1,1,-1,u=v[s]+1/u);
	return([numerator(u),denominator(u)]);
	}

addhelp(pquotn,"pquotn(v,n): produces the nth convergent from the list of denominators v");

pquot(v)=
	return(vector(length(v),a,pquotn(v,a)))

addhelp(pquot,"pquot(v): produces the list of convergents from the list of denominators v");

/* Finally we can give a path going from an arbitrary cusp to another by
way of adjacent cusps. */

path(a,b)=
	{
	local(pathinftoa,r,pathinftob,pathatoinf);
	pathinftoa=pquot(cf(a));r=length(pathinftoa);
	pathinftob=pquot(cf(b));
	pathatoinf=vector(r,i,pathinftoa[r+1-i]);
	return(concat(concat(pathatoinf,[[1,0]]),pathinftob));
	}

addhelp(path,"path(a,b): gives a path from the cusp a to the cusp b by way of adjacent cusps");

/* A function to find the gcd of quads in the appropriate field. */

qgcd(a,b)=
	{
	local(c);
	if(a==0 || b==0,,if(abs(norm(a))<abs(norm(b)),return(qgcd(b,a)),return(qgcd(b,divide([a,b])[2]))));
	if(a==0,c=b,c=a);
	u=quadunit(nfdisc(component(w,1)));\\i.e., fundamental unit
	u0=quadarch(u)[1];
	if(u0<0,u=conj(u);u0=-1/u0,);
	c0=quadarch(c)[1];
	if(c0<0,c=-c/u;c0=-c0/u0,);
	return(c/u^(round(log(c0)/log(u0))));
	}

addhelp(qgcd,"qgcd(a,b): where a and b are quads, finds their gcd.");