/************************************************************/
/* Author : Gonzalo Tornaria <mailto:tornaria@math.utexas.edu>
 * Date   : December 2003
 */

/*
 * Explicit routines to compute the central value of the L-series of
 * quadratic twists using the data in <data_imag.gp> and <data_real.gp>.
 * NOT practical to compute lots of twists of the same elliptic curve:
 *     all the coefficients of each theta series can and should be computed
 *     at the same time, in a optimized routine.
 *     In pari, qfrep, but only available from version 2.2.7
 * (anyway, this is LOTS faster than using elllseries on the twist!!)
 *
 * Take this as an example of how to use the data, as well as the precise
 * formulas to compute the central value of the twisted L-series.
 * I've checked this at least for all the available data, and for all
 * twists by discriminants less than 1000.
 *
 * Imaginary case based on Gross, 1987 (Heights and the Special Values of L-series)
 * Real case based on joint work with Ariel Pacetti, in preparation.
 *
 */

/************************************************************/

/* function to compute the coefficient of a theta series; 
 * qfminim is not flexible enough to do it all at once,
 * and qfrep is not available until 2.2.7 (alpha)
 */
q3_coeff(q,i)=
{
  local(A);
  A=[2*q[1],q[6],q[5];q[6],2*q[2],q[4];q[5],q[4],2*q[3]];
  qfminim(A,2*i,0)[1]-qfminim(A,2*i-1,0)[1];
}

/* wrappers on the elliptic curve data */
tw_ell(tw)=tw[1][4];
tw_deg(tw)=tw[1][3];
tw_level(tw)=tw[1][2];
tw_name(tw)=tw[1][1];

/* compute a coefficient of the given linear combination of theta series */
tw_coeff(tw,n) = sum(i=1,length(tw[2]),tw[2][i]*q3_coeff(tw[3][i],n));
/* The height of the linear combination, for the formula */
tw_height(tw) = sum(i=1,length(tw[2]),tw[4][i]*tw[2][i]^2);

/************************************************************/

/* Central value of the L-series of the twist -- using theta series */
ctwist(tw,D,epsilon=1E-20)=
{
  if(D>0,
     ctwist_const_r(tw),
     ctwist_const_i(tw)
  ) * ctwist_D(tw,D);
}

/* The variable part */
/* Optimize here !!! */
ctwist_D(tw,D) =
{
  local(locsign,loc);

  if(D%4>1,error("Not a discriminant in ctwist"));

  /* This assumes that N is prime and D is fundamental */
  locsign=if(D>0,ellap(ellinit(tw_ell(tw)),tw_level(tw)),-1);
  loc=(1+locsign*kronecker(D,tw_level(tw)))/2;

  /* loc=0 implies that tw_coeff==0 and the rank of the twist is odd */
  /* Thus, it should never return "N/A" for prime conductor*/
  if(loc==0, return(if(tw_coeff(tw,abs(D))==0, 0, "N/A")));

  tw_coeff(tw,abs(D))^2 / sqrt(abs(D)) / loc;
}

/* The constant part for real twists */
ctwist_const_r(tw, epsilon=1E-20) =
{
  local(e,L);
  e=ellinit(ellquadtwist(tw_ell(tw),-tw_level(tw)));
  L=elllseries(e,1);
  if(abs(L)<epsilon, error("L-series = 0"));
  tw_deg(tw) * (1+tw_level(tw)) /* factor for the degree of the twist */
             * e.area / tw_height(tw) / L
             * sqrt(tw_level(tw));
}

/* The constant part for imaginary twists */
ctwist_const_i(tw, epsilon=1E-20) =
{
  local(e,L);
  e=ellinit(tw_ell(tw));
  L=elllseries(e,1);
  if(abs(L)<epsilon, error("L-series = 0"));
  tw_deg(tw) * e.area / tw_height(tw) / L
}

/************************************************************/

/* equation of the quadratic twist of an elliptic curve */
ellquadtwist(E,D)=
	{if(D%4!=0&&D%4!=1,error("not a discriminant"));
	[E[1]*D,D*(E[2]-E[1]^2/4*(D-1)),E[3]*D^2,D^2*(E[4]-E[1]*E[3]/2*(D-1)),
	D^3*(E[5]-E[3]^2/4*(D-1))]}

/* Central value of the L-series of the twist -- using elllseries */
etwist(tw,D)=elllseries(ellinit(ellquadtwist(tw_ell(tw),D)),1);


/************************************************************/

/* TESTING (see <verify_imag.gp> and <verify_real.gp> */

CRTAB = "
	";
CRTABTAB = "
		";

/* This verifies that ctwist == etwist
 * (for imaginary fundamental discriminants up to bound)
 */
verify_twist_i(tw,bound=100,epsilon=1E-20) =
{
  local(const,c,e,ok,bug,maxdiff);
  const=ctwist_const_i(tw,epsilon);
  for(i=1,bound,D=-i;if(isfundamental(D),
    c=const*ctwist_D(tw,D);
    e=etwist(tw,D);
    ok=abs(e-c)<epsilon;
    if(!ok,bug++;print(CRTAB,[D,e,c]));
    maxdiff=max(maxdiff,abs(e-c));
    print1(CRTAB,D," ");
  ));
  if(!bug,print(CRTABTAB,"ok (max diff = ",maxdiff,")"));
  if(bug,print(CRTABTAB,bug," errors (max diff = ",maxdiff,")"));
}

/* This verifies that ctwist == etwist
 * (for real fundamental discriminants up to bound)
 */
verify_twist_r(tw,bound=100,epsilon=1E-20) =
{
  local(const,c,e,ok,bug,maxdiff);
  const=ctwist_const_r(tw,epsilon);
  for(i=1,bound,D=i;if(isfundamental(D),
    c=const*ctwist_D(tw,D);
    e=etwist(tw,D);
    ok=abs(e-c)<epsilon;
    if(!ok,bug++;print(CRTAB,[D,e,c]));
    maxdiff=max(maxdiff,abs(e-c));
    print1(CRTAB,D," ");
  ));
  if(!bug,print(CRTABTAB,"ok (max diff = ",maxdiff,")"));
  if(bug,print(CRTABTAB,bug," errors (max diff = ",maxdiff,")"));
}

/************************************************************/