\\ Copyright (c) 2003 Henri Darmon and Adam Logan
/* Some code that supports finding a_p for elliptic curves defined over
quadratic fields, and some general nonsense. */

/* First, a little snippet to find a quadratic nonresidue.  Just tries things
until we find one.  Input: p and f, an allegedly irreducible polynomial mod p.
Output: a polymod giving a nonsquare element of the field Z[x]/(p,f). */

nonres(p,f)=
	{
	local(d,v);
	if(p==2,error("everything's a square in that field"),);
	v= (mattranspose(factormod(f,p)));
	if (length(v)>1 || v[2,1]>1,error("reducible polynomial in nonres"),);
	d=poldegree(f);
	while(1,try=Mod(sum(i=0,d-1,variable(f)^i*Mod(random(p),p)),f);if(try^((p^d-1)/2)==p-1,break,));
	return(try)
	}

addhelp(nonres,"nonres(p,pol): assumes that pol is irreducible mod p, and returns a nonsquare in the field Z[x]/(p,f)");

/* Taking square roots in finite fields, given a quadratic nonresidue (see
above).  Is this called Coppersmith's algorithm?  Parameters are polymods
giving the element of the field and a nonresidue, and the order of the field
(just for convenience); return value is the square root, if there is one. */

sqrff(sq, nr, q)=
	{
	local(n,r,s,guess,i);
	if (Mod(q,2),,return(sq^(q/2))); /* char 2 is easy (if not useful)*/
	if (sq == 0, return(0),);
	if (sq^((q-1)/2)-1,error("not a square in sqrff",));
	r=valuation(q-1,2);s=(q-1)/2^r;
	n= nr^s;
	guess = sq^((s+1)/2);
	/* Now we have n, an element of order r, the highest power of
	2 that divides the order of the group of units.  We also have an
	initial guess for what the square root might be; it's off by a
	power of n.  The method now is to adjust our guess at each step
	to bring down the order, if necessary.  We know that a^s has order
	dividing 2^(r-1); if it is 2^(r-1), we need to adjust the guess
	by multiplying it by n.  Then we're off by something of order
	dividing 2^(r-2), etc. */
	forstep(i=r-1,1,-1,
		if((guess^2/sq)^(2^(i-1))-1,guess=guess*n^(2^(r-i-1)),));
	if(guess^2 != sq, error("Bug in sqrff!"),return(guess));
	}

addhelp(sqrff,"sqrff(square,nonres,q): given square, a polymod in a finite field which should be a square, nonres, a nonsquare in the same field (use the nonres() function to find if necessary), and q, the order of a the field, finds a square root of the square");

/* The baby-step-giant-step method for finding the order of a point
on an elliptic curve.  Input is the curve, the point, and the order of
the field.  Output, with luck, is the order. */

bsgs(ec,pt,q)=
	{
	local(wub,steps,babysteps,newstep,giantstep,giant,v);
	if(length(pt)==1,return(1),);
	if(pt[2]==0,return(2),);
	/* Now we're done with the trivial cases.  First priority is to
	compute the Weil bound, then figure out how many steps we need. */
	wub = q + 1 + floor(2*sqrt(q));
	steps = ceil(sqrt(wub));
	babysteps=vector(steps);
	newstep=pt;
	babysteps[1]=vectorize(pt,1);
	for(i=2,steps,newstep=elladd(ec,newstep,pt);if(newstep == [0],return(i),);babysteps[i]=vectorize(newstep,i));
	/* Now the vector babysteps contains the first so-many multiples
	of the point, and if we've found the order already we can go home.
	Time to sort it.  Note that newstep retains steps * initial point. */
	babysteps = vecsort(babysteps,,2);
	/* And now for the giant steps. */
	giantstep=ellpow(ec,newstep,-1);giant=[0];
	for(i=1,steps,
		giant=elladd(ec,giant,giantstep);
		if(giant==[0],return(i*steps),);
		v=vectorize(giant)[1];
		if(found=search(v,babysteps),return(i*steps+found),));
	/* Something's wrong if we haven't yet found it. */
	error("Bug in bsgs routines.");
	}

addhelp(bsgs,"bsgs(ec,point,q): finds the order of a point on an elliptic curve over a finite field (where q is the order of the field)");

/* Sample timing: for a point on a curve over F_1019^2, the naive algorithm
to find the order took 18.4 seconds, this took 0.5. */

/* Now the real program, at last.  Takes a 5-component vector giving the
coefficients (of course, this breaks if a1 or a3 is not 0--shh!), a
nonresidue to twist by, and, for convenience again, the polynomial defining
the field and the residue characteristic.  Output should, with luck, be the
order of the group of points on the curve. 

Uses the Shanks-Mestre algorithm with baby-step-giant-step, and tacitly
assumes that p is not 2. 

Here I try a trick that should speed things up a lot in practice.  The
reality is that a randomly chosen point probably has order N_E/k for some
k a lot smaller than the order of the field, so...well, read the code. 
*/

ellaq(e,n,f,p,slow=0)=
	{
	local(d,curves,orders,tryx,tryy,curvept,a,wlb);
	d=poldegree(f);
	if(p^d<100,return(trivialellaq(e,f,p,d)),);
	curves=[ellinit(e),ellinit(elltwist(e,n))];
	orders=[1,1];
	wlb=p^d+1-2*sqrt(p^d);
	while(1,
		tryx=Mod(sum(i=0,d-1,Mod(random(p),p)*x^i),f);
		tryy=tryx^3+e[2]*tryx^2+e[4]*tryx+e[5];
		/* We've chosen (x,f(x)) now.  If f(x) = 0, we get a 2-torsion
		point on both the curves; deal appropriately.  Otherwise, we
		get a curve-and-point, either the original or the twisted
		one.  Find its order and adjust. */
		if(tryy,if(tryy^((p^d-1)/2)+1,curvept=[1,[tryx,sqrff(tryy,n,p^d)]],curvept=[2,[tryx*n,sqrff(tryy*n^3,n,p^d)]]);orders[curvept[1]]=lcm(orders[curvept[1]],if(slow,bsgs(curves[curvept[1]],curvept[2],p^d),fastorder(curves[curvept[1]],curvept[2],p^d))),orders=lcm(orders,2));
		a=lcm(orders[1],orders[2]);
		if(a>4*sqrt(p^d),cong=chinese(Mod(0,orders[1]),2*Mod(p^d+1,orders[2]));t=floor(wlb/a)*a+lift(cong);if(t<wlb,t=t+a,);return(t),));
		/* Now the order of whichever one it is has been adjusted. 
		The above checks whether we're done.*/
	}
	
/* Estimated timings: when q = 10007^2 this takes about 220 seconds for one
given curve.  The naive method without precomputation of squares would
take about 39000 seconds.  I couldn't even allocate a p by p matrix for
the precomputation. */

addhelp(ellaq, "ellaq(e,n,f,p,slow=0): uses the baby-step-giant-step method to count points on an elliptic curve over a finite field.  Here e gives the curve, n a nonresidue, the field is Z[x]/(f,p), and slow not being 0 causes the use of a much slower method, not recommended.");

/* Let's now try a totally different method for finding the order of a
point on an elliptic curve.  First, check whether its order is small.
If not, it's enough to find a multiple of the order within the Weil bounds:
that then has to be the number of points on the curve! */

mediumorder(ec,pt,q)=
	{
	local(wlb,wub,steps,babysteps,newstep,giant,v,i);
	if(length(pt)==1,return(1),);
	if(pt[2]==0,return(2),);
	/* Now we're done with the trivial cases.  First priority is to
	compute the Weil bounds, then figure out how many steps we need. */
	wub = q + 1 + floor(2*sqrt(q));
	wlb = 2*q + 2 - wub;
	steps = wub - wlb + 1;
	babysteps=vector(steps);
	newstep=pt;
	babysteps[1]=vectorize(pt,1);
	for(i=2,steps,newstep=elladd(ec,newstep,pt);if(newstep == [0],return(i),);babysteps[i]=vectorize(newstep,i));
	babysteps = vecsort(babysteps,,2);
	giant = ellpow(ec,pt,wub);
	if(giant==[0],return(wub),);
	v=vectorize(giant)[1];
	if(found=search(v,babysteps),return(wub-found),error("Bug in mediumorder."));
	}

addhelp(mediumorder, "mediumorder: obsolete function, included for completeness");
/* A hybrid method that should be faster than the above.  We use the
baby-step-giant-step idea to check whether the order is at most the
difference between lower and upper bounds.  If so, that's all;
if not, there is a unique multiple of the order between lower and
upper bounds.  Find it! */

fastorder(ec,pt,q)=
	{
	local();
	if(length(pt)==1,return(1),);
	if(pt[2]==0,return(2),);
	/* Now we're done with the trivial cases.  First priority is to
	compute the Weil bounds, then figure out how many steps we need. */
	wub = q + 1 + floor(2*sqrt(q));
	wlb = 2*q + 2 - wub;
	range = wub-wlb+1;
	steps = floor(sqrt(range))+1;
\\	print("Dividing range ", [wlb, wub], " into steps of size ", steps);
	babysteps = vector(steps);
	newstep = pt;
	babysteps[1]=vectorize(pt,1);
	multiple = -1;
	monster=ellpow(ec,pt,wub);
\\	print("finished setting up");
	if (monster == [0], multiple = wub,);
	if (monster == pt, multiple = wub-1);
	for(i=2,steps,newstep=elladd(ec,newstep,pt);if(newstep == [0],return(i),);if(monster == newstep, multiple = wub - i);babysteps[i]=vectorize(newstep,i));
	/* Now the vector babysteps contains the first so-many multiples
	of the point, and if we've found the order already we can go home.
	Time to sort it.  Note that newstep retains steps * initial point. */
	babysteps = vecsort(babysteps,,2);
\\	print("finished sorting");
	/* And now for the giant steps. */
	giantstep=ellpow(ec,newstep,-1);giant=[0];
	/* We'll work down from 'monster', thus discovering the order if
	it isn't too small.  In this loop, we count up by giant steps,
	checking whether the point has small order.  In case it doesn't,
	we're also counting down from the top. */
	for(i=1,steps,
\\		if(Mod(i,20),,print("20 medium steps"));
		giant=elladd(ec,giant,giantstep);
		if(giant==[0],return(i*steps),);
		if(giant==ellpow(ec,pt,-i*steps),,error("bug1 in fastorder"));
		v=vectorize(giant)[1];
		if(found=search(v,babysteps),return(i*steps+found),);
		monster = elladd(ec,monster,giantstep);
		if(monster == ellpow(ec,pt,wub-i*steps),,error("bug2 in fastorder"));
		if(monster==[0],multiple = wub - i*steps,);
		v=vectorize(monster)[1];
		if(found=search(v,babysteps),multiple=wub-i*steps-found,));
	if(multiple==-1,error("bug3 in fastorder"),/*print("point has large order: only multiple is ", multiple)*/;return(multiple);)	
	}		

addhelp(fastorder,"fastorder(ec,point,q): fast function for finding the order of a point on an elliptic curve");

/* It might be even better to replace the binary search with a hashing
algorithm.  As long as we're dealing with only q^1/4 points, we should
be able to build a hash table with very few collisions.... */

/* trivialellaq: what we do when it's too small for b-s-g-s.  Just go
through and check...*/

trivialellaq(e,f,p,d)=
	{
	local(i,j,tx,tcx,total);
	total=1;
	forvec(i=vector(d,a,[0,p-1]),tx=Mod(sum(j=1,d,i[j]*variable(f)^(j-1)),f);tcx=tx^3+tx^2*e[2]+tx*e[4]+e[5];if(tcx,if(tcx^((p^d-1)/2)+1,total=total+2,),total=total+1));
	return(total);
	}

addhelp(trivialellaq,"trivialellaq(ec,poly,p,d): counts the points on the elliptic curve ec over Z[x]/(p,poly), where poly is of degree d; useful when p^d is small or p = 2, else use apvec");

/* Now a useful function in dealing with number fields.  Takes an nf,
a vector of polynomials representing elements of a number field, and the
idealprimedec of a rational prime, and returns all the reductions of the
vector at primes above the given one.  Passes along the polynomials to reduce
by, for convenience.*/

reductions(nf,v,pdec)=
	{
	local(nfredp,nfmodppols,i,a,j,w,p);
	p=pdec[1][1];
	nfredp=vector(length(pdec),i,nfmodprinit(nf,pdec[i]));
	/* Time to extract the polynomials to mod out by from the
	nfmodprinit information. */
	nfmodppols=vector(length(pdec),i,a=nfeltpowmodpr(nf,nf.zk[2],pdec[i][4],nfredp[i]);x^pdec[i][4]-sum(j=1,pdec[i][4],x^(j-1)*lift(a[j])));
	w=vector(length(pdec),i,[Mod(modpol(v,p),nfmodppols[i]),nfmodppols[i]]);
	for(i=1,length(w),if(poldegree(nfmodppols[i])==1,w[i]=lift(w[i]),));
	return(w);
	}

addhelp(reductions, "reductions(nf,v,pdec): where pdec is the prime decomposition of p in nf (as found by idealprimedec), reduces v at all primes of nf above p");
/* Now, at last, a function that takes a 5-component vector of coefficients,
the polynomial defining the field, and a prime, and gives as output a
vector of [n_p, f_p], one for each prime of the field over p.  For simplicity
we will avoid bad reduction, or even primes with the same residue
characteristic as such a prime. */

apvec(v,f,p,slow=0)=
	{
	local(e,a,w,i,etest,y,yy,pd,z);
	e=ellinit(v);
	if(p==2||valuation(norm(Mod(e.disc,f)),p),return([]),);
	a=nfinit(f);
	pd=idealprimedec(a,p);
	w=reductions(a,v,pd);
	z=vector(length(w));
	for(i=1,length(w),z[i]=[ellaq(Mod(modpol(v,p),w[i][2]),nonres(p,w[i][2]),w[i][2],p,slow),poldegree(w[i][2])]);
	z
	}

addhelp(apvec,"apvec(v,pol,p,slow=0): gives a vector of [n_p, f_p], one for each prime of Q[x]/(pol) over p, for the elliptic curve defined by v.  If slow is not 0, uses a much slower algorithm.");

/* A wrapper function for calling ellaq.  Given, as above, a 5-component
vector of coefficients, the polynomial, and upper and lower bounds for
primes, we give [n_p, p^d] for primes inert in the field.  This will be
fairly useless unless the polynomial is actually quadratic, but that could
easily be amended. */

fp2coefs(v,f,pmin,pmax,slow=0)=
	{
	local(p,w,nff);
	w=ellconvert(v);
	nff=nfinit(f);
	forprime(p=pmin,pmax,
		if(length(idealprimedec(nff,p))==1 && idealprimedec(nff,p)[1][3]==1,
/*			if(p==2,print([trivialellaq(Mod(modpol(v,2),f),f,2,poldegree(f)),2^poldegree(f)]),nr=nonres(p,f);print([ellaq(Mod(modpol(w,p),f),nr,f,p),p^poldegree(f)]))));  NB that trivialellaq doesn't work for p = 2 yet */
			if(p==2,,nr=nonres(p,f);print([ellaq(Mod(modpol(w,p),f),nr,f,p,slow),p]))));
	}

\\if(p == 2, print([apvec(v,f,2),2^poldegree(f)]), print([apvec(w,f,p)[1][1],p^poldegree(f)])),));

/* Finally, what do we do at primes that split or are inert in the field? 
GP is so fast at these that it probably isn't worth storing them; just
compute on the fly as needed. */

fpcoefs(nf,v,pd)=
	{
	local(i,w,w2);
	w=simplify(reductions(nf,v,pd));
	w2=[];
	for(i=1,length(w),if(poldegree(w[i][2])==1,w2=concat(w2,[w[i][1]]),));
	w2=lift(w2);
	/* Now we've built the list of things to count. */
	return([vector(length(w2),i,ellap(ellinit(w2[i]),pd[1][1])),pd[1][1]]);
	}

addhelp(fpcoefs,"fpcoefs(nf,v,pd): determines a_p for the curve given by v for primes of degree 1 above p in nf");