\\ Copyright (c) 2003 Henri Darmon and Adam Logan
\r ellpn.gp
\r helpers.gp
/* A file with the functions that will calculate the integrals we need. */

/* General function for printing debugging messages. */

silentprint(strvec,silent,priority,newline=1)=
	if(silent>priority,,for(i=1,length(strvec),print1(strvec[i]));if(newline,print("")));

addhelp(silentprint,"silentprint(strings,i,j,n = 1): print the (vector of) strings if i > j, followed by a newline if n is not 0.  Used to print optional messages")

/* Let's start with an easy (if not short) one.  Given an elliptic curve over
a quadratic field (the field in nfinit form, as usual, then the 5-component
vector) and a positive integer, let's find all the a_i for i an element with
that norm and pair them with the totally positive elements having that norm. */

ans(bnf,v,n)=
	{
	local(fact,apsvec,d,i,pdecs,apivec,pivec,c0,c1,c2,j,indexvec,normnvec,anvec);
	if(n==1, return([[1], [1]]),);
	fact=factor(n)~;
	apsvec=pdecs=vector(length(fact));
	/* Now we have a vector which we will fill with the a_p for the
	various primes dividing n. */
	d=bnf.disc;
	for(i=1,length(fact),
		if(kronecker(d,fact[1,i])==-1,if(Mod(fact[2,i],2),return([]),apsvec[i]=fact[1,i]^2+1-find(coefsvec,fact[1,i],Pi);if(apsvec[i]==fact[1,i]^2+1-Pi,error("Coefficient not found!"),));pdecs[i]=[fact[1,i]],fdec=idealprimedec(bnf,fact[1,i]);apsvec[i]=fpcoefs(bnf,v,fdec)[1];pdecs[i]=vector(length(fdec),j,positivize(bnf,nfbasistoalg(bnf,bnfisprincipal(bnf,fdec[j])[2])))));
\\	print("apsvec has ", apsvec);
\\	print("pdecs has ", pdecs);
	/* The return ([]) is for the case where n has p^odd, where p is a
	prime inert in the field.  Such a thing can't be a norm.
	Anyway, now we have a vector with the a_p for the primes dividing
	n.  Next goal is to produce the a_{p^i}. */
	apivec=vector(length(fact));
	for(i=1,length(fact),
		if(kronecker(d,fact[1,i])==-1,
			c0=1;c2=c1=apsvec[i];for(j=1,fact[2,i]/2-1,c2=apsvec[i]*c1-fact[1,i]^2*c0;c0=c1;c1=c2);apivec[i]=[c2],);
		if(kronecker(d,fact[1,i])==0,
			c0=1;c2=c1=apsvec[i][1];for(j=1,fact[2,i]-1,c2=apsvec[i][1]*c1-fact[1,i]*c0;c0=c1;c1=c2);apivec[i]=[c2],);
		/* If inert or split, we just have to work out the
		linear recurrence.  Otherwise, a bit more work... */
		if(kronecker(d,fact[1,i])==1,
			tvec=vector(2,x,vector(fact[2,i]+1));
			for(y=1,2,for(j=1,fact[2,i]+1,if(j==1,tvec[y][j]=1,if(j==2,tvec[y][j]=apsvec[i][y],tvec[y][j]=tvec[y][j-1]*tvec[y][2]-fact[1,i]*tvec[y][j-2]))));
			apivec[i]=vector(fact[2,i]+1,j,tvec[1][j]*tvec[2][fact[2,i]+2-j])));
\\	print("apivec has ", apivec);
	/* apivec is now a vector with 1 component for each prime dividing
	the order, giving the a_i for elements of norm p^n.  Time now for
	the task of producing the generators for the powers of primes. */
	pivec=vector(length(fact));
	for(i=1,length(fact),
		if(kronecker(d,fact[1,i])==1,pivec[i]=vector(fact[2,i]+1,j,pdecs[i][1]^(j-1)*pdecs[i][2]^(fact[2,i]+1-j)),if(kronecker(d,fact[1,i])==0,pivec[i]=[pdecs[i][1]^fact[2,i]],pivec[i]=[pdecs[i][1]^(fact[2,i]/2)])));
\\	print("pivec has ", pivec);
	/* Now it's just two forvecs to run through them all at once. */
	indexvec = vector(length(apivec),i,[1,length(apivec[i])]);
	anvec=[];
	forvec(ind=indexvec,anvec=concat(anvec,prod(i=1,length(ind),apivec[i][ind[i]])));
	normnvec=[];
	forvec(ind=indexvec,normnvec=concat(normnvec,prod(i=1,length(ind),pivec[i][ind[i]])));
	return([anvec,normnvec]);
	}

addhelp(ans,"ans(bnf, ev, n): where bnf is the bnfinit of a field (restricted to be quadratic, at least for now), ev a 5-component vector defining an elliptic curve, and n a positive integers, this gives the vector of [nu, a_nu] for a complete set of totally positive nu of norm n up to squares");

/* Now let's do the inner loop.  Given a positive integer, we'll evaluate
the sum given in the display below (8.1) in Darmon over all elements with
that norm.  We'll need the answer from ans(n) above, a requested precision
in number of digits, and the big number field (convenient way to include
the discriminant and the fundamental unit).  Also, of course, we'll need
the limits of integration: xi and yi are expected to be 2-component vectors.
Finally, we have to keep track of which place is which, so send along a
vector of roots.

Basically, we evaluate the sum over a class of representatives, stopping when
the terms become too tiny to matter. */

innern(bnf,ansvec,xi,yi,dif,pr)=
	{
	local(fusq,di,total,n,ni,n2,j,tempterm,firstterm,inter);
	if(length(ansvec)==0,return(0),);
	if(bnf.clgp[1] != 1, error("Field doesn't have class number 1!"),);
\\	print(ansvec);
	n=norm(ansvec[2][1]);
\\	print("working with norm ", n);
	fusq=arch(bnf.fu[1]^2,bnf.roots);
	di=arch(dif,bnf.roots);
\\	print(fusq,di);
	total=0;
	for(i=1,length(ansvec[1]),
		if(ansvec[1][i],
			n2=ansvec[2][i];
			ni=arch(n2,bnf.roots);
			inter=sum(j=-5,5,term(xi,yi,di,ni,fusq,j));
			j=6;tempterm=term(xi,yi,di,ni,fusq,5);
			while(1,
				lastterm=tempterm;
				tempterm=term(xi,yi,di,ni,fusq,j);
				inter=inter+tempterm;
				if((abs(tempterm)<10^-(pr+log(pr)) && abs(tempterm)< abs(lastterm)/100)||abs(tempterm)<10^(-2*pr),break,j=j+1));
			j=-6;tempterm=term(xi,yi,di,ni,fusq,-5);
			while(1,
				lastterm=tempterm;
				tempterm=term(xi,yi,di,ni,fusq,j);
				inter=inter+tempterm;
				if((abs(tempterm)<10^-(pr+log(pr)) && abs(tempterm)< abs(lastterm)/100)||tempterm==0,break,j=j-1));
			inter=inter*bnf.disc*ansvec[1][i]/n;
			total=total+inter,));
	return(total);
	}

addhelp(innern,"not really intended to be called directly");

/* Isolating the computation of an individual term. */

term(x,y,d,n,fu,k)=(e2(n[1]*y[1]*fu[1]^k/d[1])-e2(n[1]*x[1]*fu[1]^k/d[1]))*(e2(n[2]*y[2]*fu[2]^k/d[2])-e2(n[2]*x[2]*fu[2]^k/d[2]));

addhelp(term,"calculates a term in the sum that will be evaluated in an integral.  Not really intended for users");

/* Attempting an improved version that knows how to use the invariance to
move the limits around.  */

integrate(bnf,v,xi,yi,pr,silent=0)=
	{
	local(dif,val,an,tmp,upto,newlimits,integ);
	silentprint(["integrate called with parameters ",v," ",xi," ",yi," ",pr],silent,1);
	val=0;integ=0;
	newlimits=adjustlimits(bnf,xi,yi);
	silentprint(["found new limits"],silent,0);
	dif=positivize(bnf,nfbasistoalg(bnf,bnfisprincipal(bnf,bnf.diff)[2]));
	for(j=1,length(newlimits),
	silentprint(["entering main loop"],silent,1);
	for(i=1,500,
		an=ans(bnf,v,i);
		if(Mod(i,100),,silentprint([" ... ",i],silent,0,0));
		val=val+innern(bnf,an,newlimits[j][1],newlimits[j][2],dif,pr));
\\	print(val);
	upto=500;
	while(1,
		tmp=0;
		for(i=upto+1,upto+500,
			if(Mod(i,100),,silentprint([" ... ",i],silent,0,0));
			an=ans(bnf,v,i);
			tmp=tmp+innern(bnf,an,newlimits[j][1],newlimits[j][2],dif,pr));
		silentprint([" ",tmp," ",val],silent,0);
		if(abs(val)<10^-pr,silentprint(["\ndone one integral, returning 0"],silent,1);break(),);
		if(abs(tmp/val)<10^-pr,silentprint(["\ndone one integral, returning ", tmp+val],silent,1);integ=tmp+val;break,val=tmp+val;upto=upto+500)));
	silentprint(["total value: ", integ],silent,1);
	return(integ);
	}

addhelp(integrate,"integrate(bnf,v,xi,yi,pr,silent=0): attempts to work out a nonnormalized integral of a Hilbert modular form over bnf corresponding to the elliptic curve whose defining equation is given in v, with the limits xi, yi, to the (slightly optimistic) precision pr, and silently if silent = 0.  Most users would be better off to call normalint or normalintminus instead of this");
	

/* normalint.  Call syntax: bnf is the bnfinit of the number field; v is
the 5-component vector for the elliptic curve, with x being a root of the
polynomial defining bnf; xi and yi are limits; pr is the requested precision,
in digits. */

normalint(bnf,v,xi,yi,pr,silent=0)=
	{
	local(ellvec,cst,gu);
	varch=arch(v,bnf.roots);
	for(i=1,length(varch),varch[i]=realize(varch[i]));
	ellvec=vector(length(bnf.roots),i,ellinit(varch[i]));
	cst=sqrt(1/bnf.disc)/prod(i=2,length(ellvec),ellvec[i].omega[1]/2);
	gu=goodunit(bnf);
\\	print("good unit: ", gu);
\\	print("multiplying by constant: ", cst);
	return(cst*(integrate(bnf,v,xi,yi,pr,silent)+integrate(bnf,v,adjust(xi,gu),adjust(yi,gu),pr,silent)));
	}

addhelp(normalint,"normalint(bnf,v,xi,yi,pr,silent=0): bnf is the bnfinit of a number field, v a 5-component vector for an elliptic curve, symbolically in terms of a generator of v, xi and yi are vectors of limits, pr is the precision (supposedly), and values of silent greater than 0 make the program quieter.  This computes the integral of Darmon's normalized differential form omega_+ for the modular form believed to come from v");

/* normalintminus.  Call syntax: as for normalint. */

normalintminus(bnf,v,xi,yi,pr,silent=0)=
	{
	local(ellvec,cst,gu);
	varch=arch(v,bnf.roots);
	for(i=1,length(varch),varch[i]=realize(varch[i]));
	ellvec=vector(length(bnf.roots),i,ellinit(varch[i]));
	cst=sqrt(1/bnf.disc)/(I*imag(ellvec[1].omega[2]));
	gu=goodunit(bnf);
\\	print("good unit: ", gu);
\\	print("multiplying by constant: ", cst);
	return(cst*(integrate(bnf,v,xi,yi,pr,silent)-integrate(bnf,v,adjust(xi,gu),adjust(yi,gu),pr,silent)));
	}

addhelp(normalintminus,"normalintminus(bnf,v,xi,yi,pr,silent=0): with notation as for normalint, this computes the integral of Darmon's normalized differential form omega_-");

normalintminus2(bnf,v,xi,yi,pr,silent=0)=
	{
	local(ellvec,cst,gu);
	varch=arch(v,bnf.roots);
	for(i=1,length(varch),varch[i]=realize(varch[i]));
	ellvec=vector(length(bnf.roots),i,ellinit(varch[i]));
	cst=sqrt(1/bnf.disc)/prod(i=2,length(ellvec),ellvec[i].omega[1]/2);
	gu=goodunit(bnf);
\\	print("good unit: ", gu);
\\	print("multiplying by constant: ", cst);
	return(cst*(integrate(bnf,v,xi,yi,pr,silent)-integrate(bnf,v,adjust(xi,gu),adjust(yi,gu),pr,silent)));
	}

normalintminus3(bnf,v,xi,yi,pr,silent=0)=
	{
	local(ellvec,cst,gu);
	varch=arch(v,bnf.roots);
	for(i=1,length(varch),varch[i]=realize(varch[i]));
	ellvec=vector(length(bnf.roots),i,ellinit(varch[i]));
	cst=sqrt(1/bnf.disc)/(I*imag(ellvec[1].omega[2]));
	gu=-goodunit(bnf);
\\	print("good unit: ", gu);
\\	print("multiplying by constant: ", cst);
	return(cst*(integrate(bnf,v,xi,yi,pr,silent)-integrate(bnf,v,adjust2(xi,gu),adjust2(yi,gu),pr,silent)));
	}

/* A little one to do integrals from 0 to infinity.  Now improved to manage
the limits appropriately.  Note that multiplying and dividing by units is
completely useless, so that is now eliminated. */

normalint0toinf(bnf,v,x0,y0,pr,silent=0)=
	{
	local(g1,g2);
	g1 = abs(1/imag(1/sqrt(x0*y0)))*I/2; g2 = -.1/g1;
	silentprint(["using ",[g1, g2]," to fiddle the limits"],silent,1);
	return(normalint(a,ev,[x0,10^pr*I],[g1,I],pr,silent)+normalint(a,ev,[-1/x0,I],[-1/g1,10^pr*I],pr,silent)+normalint(a,ev,[g1,10^pr*I],[y0,g2],pr,silent)+normalint(a,ev,[-1/g1,-1/g2],[-1/y0,10^pr*I],pr,silent));
	}

addhelp(normalint0toinf,"normalint0toinf(bnf,v,x0,y0,pr,silent=0): calculates the integral of Darmon's normalized differential form for the elliptic curve given by v over the bnf, from x0 to y0 and from 0 to infinity, to precision pr, and quietly if silent > 0");

/* Analogy test: normalint0toinfminus is to normalint0toinf as normalintminus
is to ...? */

normalint0toinfminus(bnf,v,x0,y0,pr,silent=0)=
	{
	local(g1,g2);
	g1 = abs(1/imag(1/sqrt(x0*y0)))*I/2; g2 = -.1/g1;
	silentprint(["using ",[g1, g2]," to fiddle the limits"],silent,1);
	return(normalintminus(a,ev,[x0,10^pr*I],[g1,I],pr,silent)+normalintminus(a,ev,[-1/x0,I],[-1/g1,10^pr*I],pr,silent)+normalintminus(a,ev,[g1,10^pr*I],[y0,g2],pr,silent)+normalintminus(a,ev,[-1/g1,-1/g2],[-1/y0,10^pr*I],pr,silent));
	}

addhelp(normalint0toinfminus,"normalint0toinfminus(bnf,v,x0,y0,pr,silent=0): calculates the integral of Darmon's normalized differential form for the elliptic curve given by v over the bnf, from x0 to y0 and from 0 to infinity, to precision pr, and quietly if silent > 0");

/* invtest tests the SL_2-invariance of the integral.  Seems to work...

Notation as above, except that m is a matrix in GL_2(Q(x)), representing
a matrix in SL_2(O_F) by putting x equal to the root of bnf.pol. */

invtest(bnf,v,xi,yi,pr,m,silent=0)=
	{
	local(ui,vi,a1,a2,mi);
	mi=[subst(m,x,bnf.roots[1]),subst(m,x,bnf.roots[2])];
	a1=normalint(bnf,v,xi,yi,pr,silent);
	ui=vector(2,i,(slact(mi[i],xi[i])));
	vi=vector(2,i,(slact(mi[i],yi[i])));
	a2=normalint(bnf,v,ui,vi,pr,silent);
	print("1st answer is ", a1, "\n2nd answer is ", a2);
	}

addhelp(invtest, "invtest(bnf,v,xi,yi,pr,m,silent=0): reassures suspicious users that the integral of the form corresponding to the elliptic curve v over bnf with limits xi, yi is invariant with respect to the matrix m, whose entries should be polynomials in the generator x of the bnf.  The output gives the integral and the integral of the form in the limits obtained by transforming the given ones by m");

cocycle(bnf,v,tau,X,g0,g1,pr,silent=0)=
	{
	local(gamma0,gamma1);
	gamma0=[subst(g0,x,bnf.roots[1]),subst(g0,x,bnf.roots[2])];
	gamma1=[subst(g1,x,bnf.roots[1]),subst(g1,x,bnf.roots[2])];
	return(normalint(bnf,v,[tau,slact(gamma0[2],X)],[slact(gamma0[1],tau),slact(gamma0[2],slact(gamma1[2],X))],pr,silent));
	}

addhelp(cocycle,"cocycle: evaluates Darmon's cocycle (see chapter 8 of the CBMS notes for details) at a particular gamma_0, gamma_1.  Not maintained, so use is now deprecated");


/* How do we find the nice unit in F(sqrt(whatever)) anyway?  This function
does the trick.  Input is a polynomial in some variable (not x!) indicating
the thing to square root on the integral basis of F; one can also input a
new variable name if you want the square root of a rational.  At present
the program will laugh at you, though.  Also, the field a (a global variable)
has to be quadratic, otherwise the unit computation won't work.  Output is
the continued fraction we'll use to make the integral. 

Finally, the method for choosing the best continued fraction is somewhat
crude and could presumably be improved. */

makeunitcf(quad,newvar=tttt,silent=0)=
	{
	local(var,btemp,bbtemp,units,relnorms,umodsq,changes,k,goodunits,bestsofar,currunit,fixmat,cftry,optimizer,relgen,scl);
	if(type(quad)=="t_POL",atemp=subst(a,x,variable(quad)),atemp=subst(a,x,newvar));
	var=variable(atemp);
	if(quadarch(subst(quad,var,w))[1]>0,quad=subst(quad,var,1-var);silentprint(["(using conjugate quad ", quad, ")"],silent,2),);
	tautmp=sqrt(conj(subst(quad,var,w)));if(real(tautmp),error("bug in makeunitcf"),);silentprint(["(using tau = ", tautmp,")"],silent,1);
	btemp=rnfinit(atemp,x^2-Mod(quad,atemp.pol));
	bbtemp=bnfinit(btemp[11][1]);
	if(bbtemp.sign[2]==1,,error("not an ATR extension"));
	if((scl=bnrclass(bbtemp,[1,vector(poldegree(bbtemp.pol)-2,x,1)])[1])>1,print("*** Warning: strict class number is ", scl, " ***"),);
	optimizer=optimize(atemp,btemp,rnfbasis(atemp,btemp[1])); \\find an optimal embedding;
	opt0=matrix(2,2,i,j,quadarch(optimizer[i,j])[1]);opt1=matrix(2,2,i,j,quadarch(optimizer[i,j])[2]);
	tau2=slact(opt0,tautmp);\\or should it be opt1?
	silentprint(["new tautmp is ", tau2],silent,1);
\\	if(imag(tau2)<0,print("oops, wrong sign, conjugating by [0, 1; 1, 0], tautmp now ", 1/tau2);tau2=1/tau2;optimizer=[0,1;1,0]*optimizer,);
	units = Mod(bbtemp.fu,bbtemp.pol);
	silentprint(["units of extension: ", units],silent,2);
	relnorms=vector(length(units),i,relnorm(atemp,btemp,lift(units[i])));
	silentprint(["found relative norms: ", relnorms],silent,1);
\\We've found the units; now let's locate their norms in the base field mod squares.
\\(first component is the unit, second is +- 1)
	umodsq=[1,-1,atemp.fu[1],-atemp.fu[1]];
	changes=vector(length(bb),i,
		k=0;
		silentprint(["analyzing unit ", i],silent,0);
		for(l=1,length(umodsq),
			utmp=relnorms[i]*umodsq[l];
			if(nfissquare(atemp,utmp),k=l,));
		if(k==0,error("Error dealing with norms of units!"),k));
	goodunits=vector(length(units),i,units[i]^2);
	if(changes[1]==changes[2],goodunits[1]=units[1]*units[2]);
	for(i=1,2,if(changes[i]==1,goodunits[i]=units[i]));
\\This would be more complicated over fields of higher degree.
\\Now we will find the (ignorable) unit that is already in the base field.
\\Or actually, we don't need to....it seems that GP's units will always
\\include such, or give it as product or quotient.
	silentprint(["Generators of units of square norm: ", goodunits],silent,0);
	relgen=1;
	if(isinbase(goodunits[1]),relgen=goodunits[2],);
	for(i=-1,1,if(isinbase(goodunits[1]^i*goodunits[2]),relgen=goodunits[1],));
	if(relgen==1,error("bug finding relative generator"),silentprint(["found relative generator ", relgen, " for units"],silent,1));
\\That should be fixed...
	bestsofar=[0,1000,10^-100,[]];
	for(i=-3,3,  \\arbitrary choice
		if(i==0, next,);
		currunit=lift(relgen^i);
		silentprint(["now trying ", currunit, " with exponent ", i],silent,0);
		fixmat=subst([quad*polcoeff(currunit,2)+polcoeff(currunit,0),quad*(quad*polcoeff(currunit,3)+polcoeff(currunit,1));(quad*polcoeff(currunit,3)+polcoeff(currunit,1)),quad*polcoeff(currunit,2)+polcoeff(currunit,0)],variable(quad),w);\\should it be optimizer or conjugate thereof?
\\Kludgy fix.  We really should choose the right one from the beginning.
		silentprint(["starting from matrix ", fixmat],silent,0);
		fm=vector(2,zz,matrix(2,2,z1,z2,quadarch(fixmat[z1,z2])[zz]));if(abs(slact(fm[1],tautmp)-tautmp)<10^-10,fixmat=conj(optimizer*fixmat*optimizer^-1),if(abs(slact(fm[2],tautmp)-tautmp)<10^-10,fixmat=(optimizer*fixmat*optimizer^-1),print(fm,tautmp);error("bug in makeunitcf")));
		silentprint(["using matrix ", fixmat],silent,0);
		if(fixmat[2,1]==0,silentprint(["rejecting unit--in base field"],silent,0);next,);
		cftry=conj(cf(fixmat[,1]~));
		silentprint(["found continued fraction ", cftry],silent,0);
		limvec = unpack(makepath(cftry,tau2,silent));recvec=vector(length(limvec),i,-1/limvec[i]);closest=min(vecmin(imag(limvec)),vecmin(imag(recvec)));
		silentprint(["comparing to previous best: distance ", closest, " from real axis"],silent,0);
		if((closest^length(cftry)>=bestsofar[3]^bestsofar[2]),bestsofar=[i,length(cftry),closest,cftry],));
	silentprint(["Best unit found: ", relgen^bestsofar[1], " giving rise to continued fraction ", bestsofar[4], " with minimum distance ", bestsofar[3], " from the real axis"],silent,1);
	return([bestsofar[4],tau2]);
	}

addhelp(makeunitcf, "makeunitcf(quad,newvariable,silent): given a quad, thought of as an element of Q(sqrt(n)), finds a good unit in the quadratic extension obtained by adjoining a square root of the quad, in the sense that it gives a path that stays far from the boundary of the upper half-plane so that integrals will be manageable.  See Darmon's CBMS lecture notes for some insights on this");

/* The optimize function produces the optimal embedding of a number field
(whence the name...).  The input is the bnf of the quadratic field and
a relative extension, assumed casually to be a quadratic extension of a
field with class number 1 defined by a quadratic polynomial of the form
x^2-f(y).  Output is a matrix M such that M^-1 phi M is an optimal embedding,
where phi is the standard embedding taking f(y) to [0,f(y);1,0].  Uses
various properties, etc.  */

optimize(q,qr,relok)=
	{
	local(okgen,matroot,matgen,a2,c1,c2,cg,r,s);
	okgen=relok[,2];
	matroot=[0,-subst(lift(polcoeff(qr[1],0)),variable(polcoeff(qr[1],0)),w);1,0];
	matgen=matroot*(okgen[2][1]+okgen[2][2]*w) + matid(2)*(okgen[1][1]+okgen[1][2]*w);
\\matgen now gives a relative generator of O_K in the standard embedding.
	silentprint(["starting with ", matgen],silent,0);
	if(isintegral(trace(matgen)) && isintegral(matdet(matgen)),,error("bug1 in optimize"));
\\now must choose r such that a - rc is integral.
	a2=denominator(matgen[1,1]);if(a2!=2 && a2!=1,error("bug2 in optimize"),);
	c1=numerator(matgen[2,1]);c2=denominator(matgen[2,1]);if(type(matgen[2,1])=="t_QUAD",cg=c2/qgcd(c1,c2),cg=c2);
\\	if(isintegral(cg/a2),print("no bug 3"),error("bug3 in optimize"));
\\We'll catch that later--it has to be tolerated for now, in some cases.
	r=cg/a2;cp=matgen[2,1]*r;
\\r is now such that cp is half-integral.  This is ugly, but...
	s=Pi;\\sentinel
	for(i=-1,0,for(j=-1,0,if(isintegral(matgen[1,1]+(i+j*w)*cp),s=(i+j*w),)));
	if(s != Pi,,error("bug4 in optimize"));
	matgen2=[1,-s;0,1]*matgen*[1,s;0,1];
\\	print("first conjugating by ", [1,s;0,1]," to ", matgen2);
	matgen3=[1,0;0,cg]*matgen2*[1,0;0,1/cg];
\\	print("...then by ", [1,0;0,1/cg]," to ", matgen3);
	for(i=1,2,for(j=1,2,if(isintegral(matgen3[i,j]),,error("bug5 in optimize"))));
	silentprint(["returning ", [1,0;0,cg]*[1,-s;0,1], " from optimize"],silent,0);
	return([1,0;0,cg]*[1,-s;0,1]);
	}

addhelp(optimize,"optimize: converts the obvious embedding of F(sqrt(a)) into an optimal one.  Not of much interest to casual users");

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");

/* A function to find points on the relevant elliptic curve over the
relevant number field.  Input is the 5-component vector of the elliptic
curve; output is, with luck, the z-value of the point we're looking for.
Assumes a mildly special format for ev. */

findpoint(ev,quad,llim=1,ulim=100,silent=0,multiple=0)=
	{
	local(evsymb,ellpol,elldisc,qg,qnorm,found,try,dtry,xarch,qa,xgood,ygood,atemp,btemp,bbtemp,foundpt,e,foundz,fracp,foundx,foundy);
	default(realprecision,100);
	qg=subst(quad,variable(quad),w);
	qnorm=norm(qg);
	qa=quadarch(qg);
	\\if(qa[1]>0,qg=conj(qg);silentprint(["(using conjugate of suggested quad)"],silent,2),);
	evsymb=subst(ev,x,w);
	ellpol=y^2+x*y*evsymb[1]+y*evsymb[3]-(x^3+x^2*evsymb[2]+x*evsymb[4]+evsymb[5]);
	elldisc=poldisc(ellpol,y);
	silentprint(["using discriminant ", elldisc],silent,1);
\\This is the discriminant that must be a square in the extension for the 
\\point to be in the right field.
\\We'll search in larger and larger radii until we get lucky.
	found=0;
	for(i=llim,ulim,
		for(a1=-2*i^2,2*i^2,for(a2=-2*i^2,2*i^2,for(a3=1,i,
			if((abs(a1)<=2*(i-1)^2 && abs(a2)<=2*(i-1)^2 && abs(a3) <=(i-1))||gcd(gcd(a1,a2),a3)>1,next,);
			try=(a1+a2*w)/a3;dtry=subst(elldisc,x,try);
			if(dtry,,silentprint(["2-torsion point with x = ", try],silent,1);next); \\gets rid of 2-torsion
			if(!issquare(qnorm*norm(dtry)),next,);
			silentprint(["Norm is right when trying x-coord of ", try],silent,1);
			if(quadissquare(dtry*qg)!=[],silentprint(["found point with x-coordinate ",try],silent,2);found=1;if(torsord=istorsion(ev,try),silentprint(["point is torsion of order ", torsord],silent,2),break(3+(multiple==0))),)))));
	if(!found, print("point not found; consider longer search?");return([0]),);
\\OK, now we have the x-coordinate, but what we really need is the z-value.
\\We also need to be careful to get the right archimedean embedding, and
\\we want the y-coordinate as well.
	foundx=toquartic(try,qg);
	ellpy=subst(ellpol,x,try);
	if(quadissquare(poldisc(ellpy)/qg)!=[],sqd=subst(lift(quadissquare(poldisc(ellpy)/qg)[1]),y,w),error("bug in findpoint"));
	foundy=(toquartic(-polcoeff(ellpy,1),qg)+toquartic(sqd,qg)*x)/2;
	if(ellisoncurve(ellinit(subst(ev,x,Mod(x,algdep(w,2)))),[foundx,foundy]),silentprint(["symbolically, point is ", [foundx,foundy]],silent,2),error("bug in findpoint"));
	xarch=quadarch(try);
	e=ellinit(subst(ev,x,a.roots[1]));
	if(xarch[1]<0,xgood=xarch[1],xgood=xarch[2]);
	ygood=ellordinate(e,xgood)[1];
	if(type(quad)=="t_POL",atemp=subst(a,x,variable(quad)),atemp=subst(a,x,newvar));
	btemp=rnfinit(atemp,x^2-Mod(quad,atemp.pol));
	bbtemp=bnfinit(btemp[11][1]);
	if(nfdisc(algdep(ygood,4))!=bbtemp.disc,print("*** warning: possible bug or insufficient precision in findpoint ***");,);
\\	print(" "); \\uncomment this if the programs crash near the end....
	foundpt=[xgood,ygood];print("numerically, point is ",foundpt);
	\\silentprint(["trying to divide point..."],silent,0);
	if(ellisoncurve(e,foundpt),foundz=ellpointtoz(e,foundpt),print("*** warning: probable bug in findpoint ***"));
	\\foundpt=dividept(e,foundpt,quad,silent);
	realpart=real(foundz)/e.omega[1];fracp=algdep(realpart,1);if(max(polcoeff(fracp,0),polcoeff(fracp,1))>20,print("*** warning: point doesn't seem to have a small multiple in minus part.  ", realpart, " is supposed to be a rational number. ***"),);
	silentprint(["found z value ", foundz, ", returning imaginary part only"],silent,1);
	return(imag(foundz)*I);
	}

addhelp(findpoint,"findpoint(ev,quad,lower,upper,silent,multiple): tries to find a point on an elliptic curve (given by ev) over a quadratic extension of a quadratic extension of Q (given by the square root of quad).  This function is rather unsophisticated and could presumably be improved.  Note that the function quits after finding a single point if multiple is 0, and otherwise continues through the range");

/* istorsion: determines whether a number is the x-coordinate of a
torsion point. */

istorsion(esymb, xcoord)=
	{
	local(wreal,ereal,xreal,yreal,xyz,v);
	wreal=quadarch(w)[1];
	ereal=ellinit(subst(esymb,variable(esymb),wreal));
	xreal=quadarch(xcoord)[1];
	yreal=ellordinate(ereal,xreal)[1];
	xyz=ellpointtoz(ereal,[xreal,yreal]);
	v=lindep([ereal.omega[1],ereal.omega[2],xyz]);
	if(abs(v[3])<100,return(abs(v[3])),return(0));
	}

addhelp(istorsion,"istorsion(e,x): where e is a 5-component vector defining an elliptic curve over a quadratic field and x a number, returns the order of a torsion point with x as its x-coordinate, or 0 if the point is not torsion");

/* dividept: tells whether a point can be divided on the elliptic curve. */

dividept(ell,pt,quartic,silent)=
	{
	local();
	divpoint=pt;
	ptz=ellpointtoz(ell,pt);
	nftry=nfinit(algdep(sqrt(quadarch(subst(quartic,variable(quartic),w))[1]),4));
	trydisc=nftry.disc;
	forprime(i=2,13,silentprint(["entering main loop at ", i],silent,1);
		for(r1=0,i-1,for(r2=0,i-1,
			trydiv=ellztopoint(ell,(ptz+r1*ell.omega[1]+r2*ell.omega[2])/i);
			tryxalg=algdep(trydiv[1],4);
			if(poldegree(tryxalg) == 3 ||
	(poldegree(tryxalg)==2 && nfdisc(tryxalg)!=a.disc) ||
	(poldegree(tryxalg)==4 && !issquare(poldisc(tryxalg)/trydisc)),next,);
			if(poldegree(tryxalg) == 4 && nfdisc(tryxalg) != trydisc, next,);
			if(!nfisincl(tryxalg,nftry.pol),next,);
			tryyalg=algdep(trydiv[2],4);
			if(poldegree(tryyalg) == 3 ||
	(poldegree(tryyalg)==2 && nfdisc(tryyalg)!=a.disc) ||
	(poldegree(tryyalg)==4 && !issquare(poldisc(tryyalg)/trydisc)),next,);
			if(poldegree(tryyalg) == 4 && nfdisc(tryyalg) != trydisc,next,);
			if(!nfisincl(tryyalg,nftry.pol),next,);
			divpoint=trydiv;div=i;next(3);)));
	if(divpoint==pt,return(divpoint),silentprint(["managed to divide point by ", i,"; trying again"],silent,1);return(dividept(ell,divpoint,quartic,silent)));
	/* can't hurt to make it recursive... */
	}

addhelp(dividept,"dividept tries to divide a point on an elliptic curve over a quartic field, but may not work");

makepath(cfvec,ta,silent)=
	{
	local(l,t,u);
	l=length(cfvec);
	t=vector(l);
	silentprint(["entering makepath with parameters ", cfvec, " ", ta],silent,0);
	if(cfvec[1],
		t[1]=ta-quadarch(conj(cfvec[1]))[2];  \\make sure we get the right conjugate!
		for(i=2,l,t[i]=-1/t[i-1]+(-1)^i*quadarch(conj(cfvec[i]))[2]),
		t[2]=-1/ta+quadarch(conj(cfvec[2]))[2];
		for(i=3,l,t[i]=-1/t[i-1]+(-1)^i*quadarch(conj(cfvec[i]))[2]);
		t[1]=t[2];
		t[2]=ta);
	silentprint(["built vector of limits: ", t],silent,0);
\\At this point, the vector t is such that the integral we're looking
\\for is the sum of the semi-indefinite integrals from infinity to 0.
\\Now, semiindef(a)+semiindef(b)=int(a to -1/b) int(0 to inf), so we'll
\\combine them in adjacent pairs and run the first one by itself if there's
\\an odd number.  Alternatively, we can always use b to -1/a, whichever
\\seems better.
	u=vector(ceil(l/2));
	if(Mod(l,2),u[1]=[t[1]-1,1-1/t[1]],); \\reverse, because it's from infinity to 0.  But should it be?
	for(i=1+lift(Mod(l,2)),length(u),
		if(imag(-1/t[2*i-1-lift(Mod(l,2))])*imag(t[2*i-lift(Mod(l,2))])<imag(-1/t[2*i-lift(Mod(l,2))])*imag(t[2*i-1-lift(Mod(l,2))]),u[i]=[t[2*i-1-lift(Mod(l,2))],-1/t[2*i-lift(Mod(l,2))]],u[i]=[t[2*i-lift(Mod(l,2))],-1/t[2*i-1-lift(Mod(l,2))]]));
	return(u);
	}

addhelp(makepath,"makepath(cf,tau): makes the path corresponding to a continued fraction cf for tau, the generator of an ATR extension.  See Darmon's CBMS lecture notes for details");

/* Finally we have it all and can realize our grandiose plans.
quad is the thing we're taking a square root of to generate our field.
test produces a point on the curve, the period, and a point found by
integration, and looks for a relation with manageable coefficients. */

test(quad,prec=14,silent=0)=
	{
	local(cfr,tau,taupath,fz,kz,period,temp);
	default(realprecision,2*prec);
	a=bnfinit(a);
	e0=ellinit(subst(ev,x,a.roots[1]));
	period=if(real(e0.omega[2]),2*imag(e0.omega[2])*I,imag(e0.omega[2])*I);
	tau=sqrt(quadarch(subst(quad,variable(quad),w)));
	if(imag(tau[1])==0,tau=tau[2],tau=tau[1]);
	temp=makeunitcf(quad,,silent);
	cfr=temp[1];tau=temp[2];
	silentprint(["continued fraction: ", cfr],silent,1);
	taupath=makepath(cfr,tau,silent);
	silentprint(["found integrals to evaluate: ", taupath],silent,0);
	fz=sum(i=1,length(taupath),normalint0toinf(a,ev,taupath[i][1],taupath[i][2],prec,silent));
	realpart=real(fz)/e0.omega[1];fracp=algdep(realpart,1);if(max(polcoeff(fracp,0),polcoeff(fracp,1))>20 && abs(realpart)>10^(3-prec),print("*** warning: point doesn't seem to have a small multiple in minus part.  ", realpart, " is supposed to be a rational number. ***"),);
	silentprint(["found z value ", fz, " by integration, using imaginary part only"], silent,2);
	fz=imag(fz)*I;
	kz=findpoint(ev,quad,silent);
	silentprint(["found point over field giving rise to z-value ", kz], silent,2);
	default(realprecision,prec-2);
	print("integration produced z-value ", fz, "\nsearching produced z-value ", kz, "\nimaginary period is ", period);
	print("suggested linear dependence: ", lindep([fz,kz,period]));
	}

addhelp(test,"test(quad,precision,silent=0) tests Darmon's Stark-Heegner point conjecture over the field generated by the square root of quad to the given precision, silently if requested (default 0: chatty, 1: quieter, 2: very quiet).  One should read a data file giving the field and curve first");

testminus(quad,prec=14,silent=0)=
	{
	local(cfr,tau,taupath,fz,kz,period,temp,e1);
	default(realprecision,2*prec);
	a=bnfinit(a);
	e1=ellinit(subst(ev,x,a.roots[2]));
	tau=sqrt(quadarch(subst(quad,variable(quad),w)));
	if(imag(tau[1])==0,tau=tau[2],tau=tau[1]);
	temp=makeunitcf(quad,,silent);
	cfr=temp[1];tau=temp[2];
	silentprint(["continued fraction: ", cfr],silent,1);
	taupath=makepath(cfr,tau,silent);
	silentprint(["found integrals to evaluate: ", taupath],silent,0);
	fz=sum(i=1,length(taupath),normalint0toinfminus(a,ev,taupath[i][1],taupath[i][2],prec,silent));
	print("found z value ", fz, " by integration");
	default(realprecision,prec-2);
	print("periods: ", e1.omega);
	print("suggested linear dependence on real and imaginary periods: ", lindep([fz,e1.omega[1],e1.omega[2]]));
	}

addhelp(testminus, "testminus(quad,precision,silent=0) tests Darmon's prediction that omega^- integrates to 0 on the same path that produces points by integrating omega^+, at least in situations where ordinary and narrow class number are equal.  Usage notes as for test");

testminus2(quad,prec=14,silent=0)=
	{
	local(cfr,tau,taupath,fz,kz,period,temp);
	default(realprecision,2*prec);
	a=bnfinit(a);
	e0=ellinit(subst(ev,x,a.roots[1]));
	tau=sqrt(quadarch(subst(quad,variable(quad),w)));
	if(imag(tau[1])==0,tau=tau[2],tau=tau[1]);
	temp=makeunitcf(quad,,silent);
	cfr=temp[1];tau=temp[2];
	silentprint(["continued fraction: ", cfr],silent,1);
	taupath=makepath(cfr,tau,silent);
	silentprint(["found integrals to evaluate: ", taupath],silent,0);
	fz=sum(i=1,length(taupath),normalint0toinfminus(a,ev,taupath[i][1],taupath[i][2],prec,silent));
	print("found z value ", fz, " by integration");
	default(realprecision,prec-2);
	print("periods: ", e0.omega);
	print("suggested linear dependence on real and imaginary periods: ", lindep([fz,e0.omega[1],e0.omega[2]]));
	}

testminus3(quad,prec=14,silent=0)=
	{
	local(cfr,tau,taupath,fz,kz,period,temp,e1);
	default(realprecision,2*prec);
	a=bnfinit(a);
	e1=ellinit(subst(ev,x,a.roots[2]));
	tau=sqrt(quadarch(subst(quad,variable(quad),w)));
	if(imag(tau[1])==0,tau=tau[2],tau=tau[1]);
	temp=makeunitcf(quad,,silent);
	cfr=temp[1];tau=temp[2];
	silentprint(["continued fraction: ", cfr],silent,1);
	taupath=makepath(cfr,tau,silent);
	print(taupath);
silentprint(["found integrals to evaluate: ", taupath],silent,0);
	fz=sum(i=1,length(taupath),normalint0toinfminus3(a,ev,taupath[i][1],taupath[i][2],prec,silent));
	print("found z value ", fz, " by integration");
	default(realprecision,prec-2);
	print("periods: ", e1.omega);
	print("suggested linear dependence on real and imaginary periods: ", lindep([fz,e1.omega[1],e1.omega[2]]));
	}

normalint0toinfminus3(bnf,v,x0,y0,pr,silent=0)=
	{
	local(g1,g2);
	g1 = abs(1/imag(1/sqrt(x0*y0)))*I/2; g2 = -.1/g1;
	silentprint(["using ",[g1, g2]," to fiddle the limits"],silent,1);
	return(normalintminus3(a,ev,[x0,10^pr*I],[g1,I],pr,silent)+normalintminus3(a,ev,[-1/x0,I],[-1/g1,10^pr*I],pr,silent)+normalintminus3(a,ev,[g1,10^pr*I],[y0,g2],pr,silent)+normalintminus3(a,ev,[-1/g1,-1/g2],[-1/y0,10^pr*I],pr,silent));
	}