/************************************************************/
/* Authors : Gustavo Rama, Gonzalo Tornaria
 * Date    : June 2020
 *
 * Computation of the L-function associated to a paramodular form
 * of weight 3 and square-free level N
 *
 * See files "form_61a.gp" and "form_167a.gp" for examples
 *
 * 2023-04-19:
 *
 *     Added support for non-squarefree levels
 *     See "form_76.gp" and "form_96.gp" for examples
 */

lfunparamodular(data)=
/* data = [N, class, dim, root_numbers, ap_values, ap2_values] */
{
  N = data[1];
  dim = data[3];
  ep_s = Map(Mat(data[4]~));
  e = -vecprod(Mat(data[4]~)[,2]);
  ap_s = Map(Mat(data[5]~));
  pmax = vecmax(Mat(data[5]~)[,1]);
  ap2_s = Map(Mat(data[6]~));

  if(dim != 1, error("only ok for dim = 1"));

  my(euler(p)=
      if(!mapisdefined(ap_s, p),
          /* no information at p */
          return(1 + O(X)));
      ap=mapget(ap_s,p);
      vp=valuation(N, p);
      if(vp==0,
          if(mapisdefined(ap2_s, p),
              ap2 = mapget(ap2_s, p);
              /* Euler factor at a good prime */
              return(1 - ap * X + p*(ap2+1+p^2) * X^2 - ap*p^3*X^3 + p^6*X^4)));
      if(vp==1,
          if(mapisdefined(ep_s, p),
              ep = mapget(ep_s, p);
              if(mapisdefined(ap2_s, p),
                  ap2 = mapget(ap2_s, p);
                  if(ap*ep + ap2 + p != 0,
                    error("check: p=",p," expected mu=",-p-ap*ep)));
              /* Euler factor at a prime exactly dividing N */
              return((1 + ep*p*X) * (1 - (ap+ep*p)*X + p^3*X^2))));
      if(vp>=2,
          if(mapisdefined(ap2_s, p),
              ap2 = mapget(ap2_s, p);
              /* Euler factor at a prime of higher valuation */
              return(1 - ap*X + p*(ap2+p^2)*X^2)));
      /* partial information at p */
      return(1 - ap*X + O(X^2)));

  dirichlet = direuler(p=2,nextprime(pmax+1)-1,1/euler(p));

  return( lfuncreate([dirichlet, 0, [-1,0,0,1], 4, N, e]) );
}