\\ Copyright (c) 2003 Henri Darmon and Adam Logan
/* find: another binary search function.  This one works in our tables of
coefficients, which are vectors of the form [n_p, p].  The third argument
is a sentinel value to be returned in case of failure. */

find(v, n, sentinel)=
	{
	local(lower,upper,current);
	lower = 1; upper = length(v);
	while (1,
		current = floor((lower+upper)/2);
		if(v[current][2] == n, return(v[current][1]),);
		if(v[current][2] < n, lower = current+1, upper = current-1);
		if(lower > upper, return(sentinel),));
	}

addhelp(find,"binary search function, not really intended to be called by users")
/* positive takes a bnf and a vector of polynomials in the variable and
changes the signs if necessary. */

positive(bnf,v)=
	{
	for(i=1,length(v),
		if(subst(v[i],variable(v[i]),polroots(bnf.pol)[1]) < 0, v[i]=-v[i],);
		v[i]=vector(poldegree(bnf.pol),j,subst(v[i],variable(v[i]),polroots(bnf.pol)[j])));
	return(v);
	}

addhelp(positive,"positive(bnf,v): given v, a vector of polynomials, replaces each with the vector of its values at the roots of the polynomial defining bnf, with signs changed so that the first component is positive in each case")

/* positivize takes a bnf and an element in polymod form and returns a
totally positive generator of the same ideal.  You'd have to do something
a little cleverer if the field had higher degree. */

positivize(bnf,c)=
	{
	local(cv);
	cv=avals(c);
	if(sig(cv)==[1,1],return(c),);
	if(sig(cv)==[-1,-1],return(-c),);
\\	print("now trying to work with ", c*bnf.fu[1], " with vector of signs ", sig(avals(c*bnf.fu[1])));
	if(sig(avals(c*bnf.fu[1]))==[1,1],return(c*bnf.fu[1]),);
	if(sig(avals(-c*bnf.fu[1]))==[1,1],return(-c*bnf.fu[1]),);
	error("Field apparently doesn't have strict class number 1 in positivize");
	}

addhelp(positivize,"positivize(bnf,c): given a quadratic bnf and c, a polymod viewed as element of c, returns a totally positive element d of the bnf such that (c) = (d)");

sig(v)=
	return(vector(length(v),i,sign(v[i])));

addhelp(sig, "sig(v): the vector of signs of the elements of v")

avals(pm)=
	{
	local(base,roots);
	base=component(pm,1);roots=polroots(base);
	roots=realize(roots);
	return(vector(length(roots),i,subst(lift(pm),variable(lift(pm)),roots[i])));
	}

addhelp(avals,"avals(polymod): given a polymod Mod(f(t),g(t)), gives the vector of values of f at the roots of g");

e2(x)=return(if(real(2*Pi*I*x)<-500,0,if(real(2*Pi*I*x)>5000,error("too big!"),exp(2*Pi*I*x))));

addhelp(e2,"e2(x): the normalized exponential e^(2 Pi I x)");

arch(c,v)=
	return(vector(length(v),i,subst(lift(c),x,v[i])));

addhelp(arch,"arch(value,vector): substitutes value into the elements of vector.  Typically used when value is the real or complex value of an element of a number field");

/* Helper for normalint. */

adjust(lims,unit)=
	return([lims[1]*unit[1],conj(lims[2])*unit[2]]);

addhelp(adjust,"minor adjustment of limits to help in normalization. No reason to call this directly that I can see");

/* Helper for normalint. */

adjust2(lims,unit)=
	return([conj(lims[1])*unit[1],lims[2]*unit[2]]);

addhelp(adjust,"minor adjustment of limits to help in normalization. No reason to call this directly that I can see");

realize(v)= 
	{
	for(i=1,length(v),if(abs(imag(v[i]))<10^-15,v[i]=real(v[i]),));return(v);
	}

addhelp(realize,"realize(v): replaces elements of v that are very close to real numbers with those real numbers");
/* We also want to compute the normalized version of the integral, display
8.2 in Darmon.  To do this we need a unit negative at the first place
but positive at the second: first, this helper. */

goodunit(bnf)=
	{
	local(c,d);
	c=arch(bnf.fu,bnf.roots);c=[c[1][1],c[2][1]];
\\	print("identified fundamental unit; converting ", c, " to reals");
	d=realize(c);
	if(sig(d)==[1,-1],return(c),return(-c));
	}

addhelp(goodunit,"goodunit(bnf): where bnf is quadratic, returns the archimedean version of a unit of bnf that is positive at the first place and negative at the second");

semiindef(bnf,v,t,pr,silent=0)=normalint0toinf(bnf,v,1-1/t,t-1,pr,silent);

addhelp(semiindef,"semiindef: computes Darmon's semi-indefinite integral with the limit t of the form for the elliptic curve given by v over bnf, to pr places of accuracy or so, silently if silent > 0");

slact(m,a)=(m[1,1]*a+m[1,2])/(m[2,1]*a+m[2,2]);

addhelp(slact, "slact(m,a) is ma for the usual action of SL_2");

/* Helper, converting quadgen into number. */

numb(q,n)=if(type(q)=="t_QUAD",return(component(q,2)+n*component(q,3)),return(q));

/* Helper. */

nfissquare(a,v)=
	{
	local(nff);
	if(v==0,return([0]),nff=mattranspose(nffactor(a,x^2-v));if(length(nff)==1,return([]),return([polcoeff(nff[1,1],0),-polcoeff(nff[1,1],0)])));
	}

addhelp(nfissquare,"nfissquare(a,v): determines whether v is a square in the number field a, and returns the square roots if so.");
/* Another. */

quadissquare(q)=
	if(type(q)=="t_QUAD",return(nfissquare(nfinit(subst(component(q,1),x,y)),component(q,3)*y+component(q,2))),if(issquare(q),return([sqrt(q),-sqrt(q)]),return([])));

addhelp(quadissquare,"quadissquare(q): determines whether the quad q is a square, returning the numerical square roots if so.");

toquartic(q,r)=
	{
	local(pmod,ld);
	pmod=algdep(r,2);
	ld=lindep([q,1,r]);
	return(Mod(-ld[2]/ld[1]-ld[3]/ld[1]*x^2,subst(pmod,x,x^2)));
	}

addhelp(toquartic,"toquartic(q,r): embeds the quad q into the field generated by a square root of r as a polymod");

relnorm(nf,rnf,v)=relchar=rnfcharpoly(nf,rnf[1],v);return(polcoeff(relchar,0));

addhelp(relnorm,"relnorm(nf,rnf,v): returns the relative norm of v from rnf to nf.");

unpack(v) = local(vv=[]);for(i=1,length(v),vv=concat(vv,v[i]));return(vv);

addhelp(unpack,"unpack(v): expands a vector of vectors into an ordinary vector.");

isinbase(u)= if(polcoeff(lift(u),1)==0 && polcoeff(lift(u),3)==0,return(1),return(0));

addhelp(isinbase,"isinbase(u): given a polymod (expected to be modulo a polynomial of degree 4, otherwise this is useless) returns 1 if the coefficients of x and x^3 are 0--that is, if it is in the quadratic subfield");

isintegral(a)= if(type(a)=="t_INT",return(1),);if(type(a)=="t_FRAC",return(0),);if(type(a)=="t_QUAD",if(type(component(a,2))=="t_INT" && type(component(a,3))=="t_INT",return(1),return(0)),error("unknown type in isintegral"));

addhelp(isintegral,"isintegral(a): returns 1 if a is integral, else 0, at least within its range of types");

/* Helper finding arch. embeddings of a quadgen. */

quadarch(q)=if(type(q)=="t_QUAD",realize(arch(component(q,2)+x*component(q,3),polroots(component(q,1)))),[q,q]);

addhelp(quadarch,"quadarch(q): returns the archimedean embeddings of the quad q");

adjustlimits(bnf,xi,yi)=[[xi,yi]];

addhelp(adjustlimits,"adjustlimits: it would be nice if this broke up integrals into shorter but easier-to-evaluate integrals.  Unfortunately, it doesn't");

/* Now a little thing to convert a point on an elliptic curve to a vector,
so that we can track it easily.  When we record the baby steps, we need
to keep track of which ones they are, or that will be lost in sorting. */

vectorize(pt,n=0)=
	{
	local(r,s);
	if (length(pt)==1,return([-1,-1,-1,-1,n]),);
	r=lift(lift(pt[1]));s=lift(lift(pt[2]));
	return([[polcoeff(r,1),polcoeff(r,0),polcoeff(s,1),polcoeff(s,0)],n]);
	}

addhelp(vectorize, "vectorize: puts points on curves over finite fields into predictable form.  Not useful to casual users");

/* Here is the binary searching routine we need for the b-s-g-s above.
The idea is simple enough: just compare until found... Return value is
0 if not found, else the index.  To make sure we don't have to copy the
vector too much, the implementation is not recursive. */

search(v,w)=
	{
	local(left, right, middle);
	left=1;right=length(w);
	while ((left != right),
		middle = floor((left+right)/2);
		if(v==w[middle][1],return(w[middle][2]),);
		if(lex(v,w[middle][1])==1,left = max(left+1,middle), right = middle));
	if(v==w[left][1],return(w[left][2]),return(0))
	}

addhelp(search, "search: binary search function, not recommended for direct use");

/* Produce a twist of an elliptic curve, given a nonresidue to twist it by. */

elltwist(e,n)=
	return([0,n*e[2],0,n^2*e[4],n^3*e[5]]);

addhelp(elltwist,"elltwist(e,n): twists the elliptic curve e, in the form y^2 = x^3 + a[2]x^2 + a[4]x + a[5], by the nonresidue n");

/* modpol takes a vector of polynomials and reduces them mod p. */

modpol(v,p)=
	vector(length(v),i,if(poldegree(v[i])==0,Mod(v[i],p),sum(j=0,poldegree(v[i]),Mod(polcoeff(v[i],j),p)*variable(v[i])^j)));

addhelp(modpol, "modpol(v,p): reduce polynomials in v mod p");

/* A snippet that takes [a,b,c,d,e], data for an elliptic curve, and converts
it to [0,b',0,d',e']. */

ellconvert(v)=
	return(vecextract(ellchangecurve(ellinit(v),[1,0,-v[1]/2,-v[3]/2]),"1..5"));

addhelp(ellconvert,"ellconvert(e): clears a_1 and a_3 from e");