/* Copyright (C) 2014 The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */

#include "pari.h"

/* Return 1 if the point Q is a Weierstrass (2-torsion) point of the
 * curve E, return 0 otherwise */
static long
ellisweierstrasspoint(GEN E, GEN Q)
{ return ell_is_inf(Q) || gequal0(ec_dmFdy_evalQ(E, Q)); }


/* Given an elliptic curve E = [a1, a2, a3, a4, a6] and t,w in the ring of
 * definition of E, return the curve
 *  E' = [a1, a2, a3, a4 - 5t, a6 - (E.b2 t + 7w)] */
static GEN
make_velu_curve(GEN E, GEN t, GEN w)
{
  GEN A4, A6, a1 = ell_get_a1(E), a2 = ell_get_a2(E), a3 = ell_get_a3(E);
  A4 = gsub(ell_get_a4(E), gmulsg(5LL, t));
  A6 = gsub(ell_get_a6(E), gadd(gmul(ell_get_b2(E), t), gmulsg(7LL, w)));
  return mkvec5(a1,a2,a3,A4,A6);
}

/* If phi = (f(x)/h(x)^2, g(x,y)/h(x)^3) is an isogeny, return the
 * variables x and y in a vecsmall */
INLINE GEN
get_isog_vars(GEN phi)
{
  long vx = varn(gel(phi, 1));
  long vy = varn(gel(phi, 2));
  if (vy == vx) vy = varn(leading_term(gel(phi, 2)));
  return mkvecsmall2(vx, vy);
}

INLINE GEN
substvec(GEN target, GEN vars, GEN in)
{
  long nvars = lg(vars) - 1, i;
  GEN res = target;
  for (i = 1; i <= nvars; ++i)
    res = gsubst(res, vars[i], gel(in, i));
  return res;
}
/* Given isogenies F:E' -> E and G:E'' -> E', return the composite
 * isogeny F o G:E'' -> E */
static GEN
ellcompisog(GEN F, GEN G)
{
  pari_sp av = avma;
  GEN Gh, Gh2, Gh3, f, g, h, z, numz, denz, denz2, in, comp;

  checkellisog(F);
  checkellisog(G);
  Gh = gel(G,3); Gh2 = gsqr(Gh); Gh3 = gmul(Gh, Gh2);
  in = mkvec2(gdiv(gel(G,1), Gh2), gdiv(gel(G,2), Gh3));
  comp = substvec(F, get_isog_vars(F), in);
  f = gel(comp,1);
  g = gel(comp,2);
  z = gel(comp,3); numz = numer(z); denz = denom(z);
  if (!issquareall(denom(f), &h))
    pari_err_BUG("ellcompisog (denominator of composite abscissa not square)");
  h  = RgX_mul(h, numz);
  h = RgX_Rg_div(h, leading_term(h));

  denz2 = gsqr(denz);
  f = RgX_mul(numer(f), denz2);
  g = RgX_mul(numer(g), gmul(denz,denz2));
  return gerepilecopy(av, mkvec3(f,g,h));
}

/* Given an isogeny phi from ellisogeny() and a point P in the domain of phi,
 * return phi(P) */
GEN
ellisogenyapply(GEN phi, GEN P)
{
  pari_sp ltop = avma;
  GEN f, g, h, img_f, img_g, img_h, img_h2, img_h3, img, vars, tmp;
  long vx, vy;
  if (lg(P) == 4) return ellcompisog(phi,P);
  checkellisog(phi);
  checkellpt(P);
  if (ell_is_inf(P)) return ellinf();
  f = gel(phi, 1);
  g = gel(phi, 2);
  h = gel(phi, 3);
  vars = get_isog_vars(phi);
  vx = vars[1];
  vy = vars[2];
  img_h = poleval(h, gel(P, 1));
  if (gequal0(img_h)) { avma = ltop; return ellinf(); }

  img_h2 = gsqr(img_h);
  img_h3 = gmul(img_h, img_h2);
  img_f = poleval(f, gel(P, 1));
  /* FIXME: This calculation of g is perhaps not as efficient as it could be */
  tmp = gsubst(g, vx, gel(P, 1));
  img_g = gsubst(tmp, vy, gel(P, 2));
  img = cgetg(3, t_VEC);
  gel(img, 1) = gdiv(img_f, img_h2);
  gel(img, 2) = gdiv(img_g, img_h3);
  return gerepileupto(ltop, img);
}

/* isog = [f, g, h] = [x, y, 1] = identity */
static GEN
isog_identity(long vx, long vy)
{ return mkvec3(pol_x(vx), pol_x(vy), pol_1(vx)); }

/* Returns an updated value for isog based (primarily) on tQ and uQ. Used to
 * iteratively compute the isogeny corresponding to a subgroup generated by a
 * given point. Ref: Equation 8 in Velu's paper.
 * isog = NULL codes the identity */
static GEN
update_isogeny_polys(GEN isog, GEN E, GEN Q, GEN tQ, GEN uQ, long vx, long vy)
{
  pari_sp ltop = avma, av;
  GEN xQ = gel(Q, 1), yQ = gel(Q, 2);
  GEN rt = deg1pol_shallow(gen_1, gneg(xQ), vx);
  GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);

  GEN gQx = ec_dFdx_evalQ(E, Q);
  GEN gQy = ec_dFdy_evalQ(E, Q);
  GEN tmp1, tmp2, tmp3, tmp4, f, g, h, rt_sqr, res;

  /* g -= uQ * (2 * y + E.a1 * x + E.a3)
   *   + tQ * rt * (E.a1 * rt + y - yQ)
   *   + rt * (E.a1 * uQ - gQx * gQy) */
  av = avma;
  tmp1 = gmul(uQ, gadd(deg1pol_shallow(gen_2, gen_0, vy),
                       deg1pol_shallow(a1, a3, vx)));
  tmp1 = gerepileupto(av, tmp1);
  av = avma;
  tmp2 = gmul(tQ, gadd(gmul(a1, rt),
                       deg1pol_shallow(gen_1, gneg(yQ), vy)));
  tmp2 = gerepileupto(av, tmp2);
  av = avma;
  tmp3 = gsub(gmul(a1, uQ), gmul(gQx, gQy));
  tmp3 = gerepileupto(av, tmp3);

  if (!isog) isog = isog_identity(vx,vy);
  f = gel(isog, 1);
  g = gel(isog, 2);
  h = gel(isog, 3);
  rt_sqr = gsqr(rt);
  res = cgetg(4, t_VEC);
  av = avma;
  tmp4 = gdiv(gadd(gmul(tQ, rt), uQ), rt_sqr);
  gel(res, 1) = gerepileupto(av, gadd(f, tmp4));
  av = avma;
  tmp4 = gadd(tmp1, gmul(rt, gadd(tmp2, tmp3)));
  gel(res, 2) = gerepileupto(av, gsub(g, gdiv(tmp4, gmul(rt, rt_sqr))));
  av = avma;
  gel(res, 3) = gerepileupto(av, gmul(h, rt));
  return gerepileupto(ltop, res);
}

/* Given a point P on E, return the curve E/<P> and, if only_image is zero,
 * the isogeny pi: E -> E/<P>. The variables vx and vy are used to describe
 * the isogeny (ignored if only_image is zero) */
static GEN
isogeny_from_kernel_point(GEN E, GEN P, int only_image, long vx, long vy)
{
  pari_sp av = avma;
  GEN isog, EE, f, g, h, h2, h3;
  GEN Q = P, t = gen_0, w = gen_0;
  long c;
  if (!oncurve(E,P))
    pari_err_DOMAIN("isogeny_from_kernel_point", "point", "not on", E, P);
  if (ell_is_inf(P))
  {
    if (only_image) return E;
    return mkvec2(E, isog_identity(vx,vy));
  }

  isog = NULL; c = 1;
  for (;;)
  {
    GEN tQ, xQ = gel(Q,1), uQ = ec_2divpol_evalx(E, xQ);
    int stop = 0;
    if (ellisweierstrasspoint(E,Q))
    { /* ord(P)=2c; take Q=[c]P into consideration and stop */
      tQ = ec_dFdx_evalQ(E, Q);
      stop = 1;
    }
    else
      tQ = ec_half_deriv_2divpol_evalx(E, xQ);
    t = gadd(t, tQ);
    w = gadd(w, gadd(uQ, gmul(tQ, xQ)));
    if (!only_image) isog = update_isogeny_polys(isog, E, Q,tQ,uQ, vx,vy);
    if (stop) break;

    Q = elladd(E, P, Q);
    ++c;
    /* IF x([c]P) = x([c-1]P) THEN [c]P = -[c-1]P and [2c-1]P = 0 */
    if (gequal(gel(Q,1), xQ)) break;
    if (gc_needed(av,1))
    {
      if(DEBUGMEM>1) pari_warn(warnmem,"isogeny_from_kernel_point");
      gerepileall(av, isog? 4: 3, &Q, &t, &w, &isog);
    }
  }

  EE = make_velu_curve(E, t, w);
  if (only_image) return EE;

  if (!isog) isog = isog_identity(vx,vy);
  f = gel(isog, 1);
  g = gel(isog, 2);
  if ( ! (typ(f) == t_RFRAC && typ(g) == t_RFRAC))
    pari_err_BUG("isogeny_from_kernel_point (f or g has wrong type)");

  /* Clean up the isogeny polynomials f, g and h so that the isogeny
   * is given by (x,y) -> (f(x)/h(x)^2, g(x,y)/h(x)^3) */
  h = gel(isog, 3);
  h2 = gsqr(h);
  h3 = gmul(h, h2);
  f = gmul(f, h2);
  g = gmul(g, h3);
  if (typ(f) != t_POL || typ(g) != t_POL)
    pari_err_BUG("isogeny_from_kernel_point (wrong denominator)");
  return mkvec2(EE, mkvec3(f,g, gel(isog,3)));
}

/* Given a t_POL x^n - s1 x^{n-1} + s2 x^{n-2} - s3 x^{n-3} + ...
 * return the first three power sums (Newton's identities):
 *   p1 = s1
 *   p2 = s1^2 - 2 s2
 *   p3 = (s1^2 - 3 s2) s1 + 3 s3 */
static void
first_three_power_sums(GEN pol, GEN *p1, GEN *p2, GEN *p3)
{
  long d = degpol(pol);
  GEN s1, s2, ms3;

  *p1 = s1 = gneg(RgX_coeff(pol, d-1));

  s2 = RgX_coeff(pol, d-2);
  *p2 = gsub(gsqr(s1), gmulsg(2LL, s2));

  ms3 = RgX_coeff(pol, d-3);
  *p3 = gadd(gmul(s1, gsub(*p2, s2)), gmulsg(-3LL, ms3));
}


/* Let E and a t_POL h of degree 1 or 3 whose roots are 2-torsion points on E.
 * - if only_image != 0, return [t, w] used to compute the equation of the
 *   quotient by the given 2-torsion points
 * - else return [t,w, f,g,h], along with the contributions f, g and
 *   h to the isogeny giving the quotient by h. Variables vx and vy are used
 *   to create f, g and h, or ignored if only_image is zero */

/* deg h = 1; 2-torsion contribution from Weierstrass point */
static GEN
contrib_weierstrass_pt(GEN E, GEN h, long only_image, long vx, long vy)
{
  GEN p = ellbasechar(E);
  GEN a1 = ell_get_a1(E);
  GEN a3 = ell_get_a3(E);
  GEN x0 = gneg(constant_term(h)); /* h = x - x0 */
  GEN b = gadd(gmul(a1,x0), a3);
  GEN y0, Q, t, w, t1, t2, f, g;

  if (!equalis(p, 2LL)) /* char(k) != 2 ==> y0 = -b/2 */
    y0 = gmul2n(gneg(b), -1);
  else
  { /* char(k) = 2 ==> y0 = sqrt(f(x0)) where E is y^2 + h(x) = f(x). */
    if (!gequal0(b)) pari_err_BUG("two_torsion_contrib (a1*x0+a3 != 0)");
    y0 = gsqrt(ec_f_evalx(E, x0), 0);
  }
  Q = mkvec2(x0, y0);
  t = ec_dFdx_evalQ(E, Q);
  w = gmul(x0, t);
  if (only_image) return mkvec2(t,w);

  /* Compute isogeny, f = (x - x0) * t */
  f = deg1pol_shallow(t, gmul(t, gneg(x0)), vx);

  /* g = (x - x0) * t * (a1 * (x - x0) + (y - y0)) */
  t1 = deg1pol_shallow(a1, gmul(a1, gneg(x0)), vx);
  t2 = deg1pol_shallow(gen_1, gneg(y0), vy);
  g = gmul(f, gadd(t1, t2));
  return mkvec5(t, w, f, g, h);
}
/* deg h =3; full 2-torsion contribution. NB: assume h is monic; base field
 * characteristic is odd or zero (cannot happen in char 2).*/
static GEN
contrib_full_tors(GEN E, GEN h, long only_image, long vx, long vy)
{
  GEN p1, p2, p3, half_b2, half_b4, t, w, f, g;
  first_three_power_sums(h, &p1,&p2,&p3);
  half_b2 = gmul2n(ell_get_b2(E), -1);
  half_b4 = gmul2n(ell_get_b4(E), -1);

  /* t = 3*(p2 + b4/2) + p1 * b2/2 */
  t = gadd(gmulsg(3LL, gadd(p2, half_b4)), gmul(p1, half_b2));

  /* w = 3 * p3 + p2 * b2/2 + p1 * b4/2 */
  w = gadd(gmulsg(3LL, p3), gadd(gmul(p2, half_b2),
                                gmul(p1, half_b4)));
  if (only_image) return mkvec2(t,w);

  /* Compute isogeny */
  {
    GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), t1, t2;
    GEN s1 = gneg(RgX_coeff(h, 2));
    GEN dh = RgX_deriv(h);
    GEN psi2xy = gadd(deg1pol_shallow(a1, a3, vx),
                      deg1pol_shallow(gen_2, gen_0, vy));

    /* f = -3 (3 x + b2/2 + s1) h + (3 x^2 + (b2/2) x + (b4/2)) h'*/
    t1 = RgX_mul(h, gmulsg(-3, deg1pol(stoi(3), gadd(half_b2, s1), vx)));
    t2 = mkpoln(3, stoi(3), half_b2, half_b4);
    setvarn(t2, vx);
    t2 = RgX_mul(dh, t2);
    f = RgX_add(t1, t2);

    /* 2g = psi2xy * (f'*h - f*h') - (a1*f + a3*h) * h; */
    t1 = RgX_sub(RgX_mul(RgX_deriv(f), h), RgX_mul(f, dh));
    t2 = RgX_mul(h, RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(h, a3)));
    g = RgX_divs(gsub(gmul(psi2xy, t1), t2), 2LL);

    f = RgX_mul(f, h);
    g = RgX_mul(g, h);
  }
  return mkvec5(t, w, f, g, h);
}

/* Given E and a t_POL T whose roots define a subgroup G of E, return the factor
 * of T that corresponds to the 2-torsion points E[2] \cap G in G */
INLINE GEN
two_torsion_part(GEN E, GEN T)
{ return RgX_gcd(T, elldivpol(E, 2, varn(T))); }

/* Return the jth Hasse derivative of the polynomial f = \sum_{i=0}^n a_i x^i,
 * i.e. \sum_{i=j}^n a_i \binom{i}{j} x^{i-j}. It is a derivation even when the
 * coefficient ring has positive characteristic */
static GEN
derivhasse(GEN f, ulong j)
{
  ulong i, d = degpol(f);
  GEN df;
  if (gequal0(f) || d == 0) return pol_0(varn(f));
  if (j == 0) return gcopy(f);
  df = cgetg(2 + (d-j+1), t_POL);
  df[1] = f[1];
  for (i = j; i <= d; ++i) gel(df, i-j+2) = gmul(binomialuu(i,j), gel(f, i+2));
  return normalizepol(df);
}

static GEN
non_two_torsion_abscissa(GEN E, GEN h0, GEN x)
{
  GEN mp1, dh0, ddh0, t, u, t1, t2, t3;
  long m = degpol(h0);
  mp1 = gel(h0, m + 1); /* negative of first power sum */
  dh0 = RgX_deriv(h0);
  ddh0 = RgX_deriv(dh0);
  t = ec_2divpol_evalx(E, x);
  u = ec_half_deriv_2divpol_evalx(E, x);
  t1 = RgX_sub(RgX_sqr(dh0), RgX_mul(ddh0, h0));
  t2 = RgX_mul(u, RgX_mul(h0, dh0));
  t3 = RgX_mul(RgX_sqr(h0),
               deg1pol_shallow(stoi(2*m), gmulsg(2LL, mp1), varn(x)));
  /* t * (dh0^2 - ddh0*h0) - u*dh0*h0 + (2*m*x - 2*s1) * h0^2); */
  return RgX_add(RgX_sub(RgX_mul(t, t1), t2), t3);
}

static GEN
isog_abscissa(GEN E, GEN kerp, GEN h0, GEN x, GEN two_tors)
{
  GEN f0, f2, h2, t1, t2, t3;
  f0 = (degpol(h0) > 0)? non_two_torsion_abscissa(E, h0, x): pol_0(varn(x));
  f2 = gel(two_tors, 3);
  h2 = gel(two_tors, 5);

  /* Combine f0 and f2 into the final abcissa of the isogeny. */
  t1 = RgX_mul(x, RgX_sqr(kerp));
  t2 = RgX_mul(f2, RgX_sqr(h0));
  t3 = RgX_mul(f0, RgX_sqr(h2));
  /* x * kerp^2 + f2 * h0^2 + f0 * h2^2 */
  return RgX_add(t1, RgX_add(t2, t3));
}

static GEN
non_two_torsion_ordinate_char_not2(GEN E, GEN f, GEN h, GEN psi2)
{
  GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
  GEN df = RgX_deriv(f), dh = RgX_deriv(h);
  /* g = df * h * psi2/2 - f * dh * psi2
   *   - (E.a1 * f + E.a3 * h^2) * h/2 */
  GEN t1 = RgX_mul(df, RgX_mul(h, RgX_divs(psi2, 2LL)));
  GEN t2 = RgX_mul(f, RgX_mul(dh, psi2));
  GEN t3 = RgX_mul(RgX_divs(h, 2LL),
                   RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(RgX_sqr(h), a3)));
  return RgX_sub(RgX_sub(t1, t2), t3);
}

/* h = kerq */
static GEN
non_two_torsion_ordinate_char2(GEN E, GEN h, GEN x, GEN y)
{
  GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), a4 = ell_get_a4(E);
  GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
  GEN h2, dh, dh2, ddh, D2h, D2dh, H, psi2, u, t, alpha;
  GEN p1, t1, t2, t3, t4;
  long m, vx = varn(x);

  h2 = RgX_sqr(h);
  dh = RgX_deriv(h);
  dh2 = RgX_sqr(dh);
  ddh = RgX_deriv(dh);
  H = RgX_sub(dh2, RgX_mul(h, ddh));
  D2h = derivhasse(h, 2);
  D2dh = derivhasse(dh, 2);
  psi2 = deg1pol_shallow(a1, a3, vx);
  u = mkpoln(3, b2, gen_0, b6);
  setvarn(u, vx);
  t = deg1pol_shallow(b2, b4, vx);
  alpha = mkpoln(4, a1, a3, gmul(a1, a4), gmul(a3, a4));
  setvarn(alpha, vx);
  m = degpol(h);
  p1 = RgX_coeff(h, m-1); /* first power sum */

  t1 = gmul(gadd(gmul(a1, p1), gmulgs(a3, m)), RgX_mul(h,h2));

  t2 = gmul(a1, gadd(gmul(a1, gadd(y, psi2)), RgX_add(RgX_Rg_add(RgX_sqr(x), a4), t)));
  t2 = gmul(t2, gmul(dh, h2));

  t3 = gadd(gmul(y, t), RgX_add(alpha, RgX_Rg_mul(u, a1)));
  t3 = gmul(t3, RgX_mul(h, H));

  t4 = gmul(u, psi2);
  t4 = gmul(t4, RgX_sub(RgX_sub(RgX_mul(h2, D2dh), RgX_mul(dh, H)),
                        RgX_mul(h, RgX_mul(dh, D2h))));

  return gadd(t1, gadd(t2, gadd(t3, t4)));
}

static GEN
isog_ordinate(GEN E, GEN kerp, GEN kerq, GEN x, GEN y, GEN two_tors, GEN f)
{
  GEN g;
  if (! equalis(ellbasechar(E), 2LL)) {
    /* FIXME: We don't use (hence don't need to calculate)
     * g2 = gel(two_tors, 4) when char(k) != 2. */
    GEN psi2 = ec_dmFdy_evalQ(E, mkvec2(x, y));
    g = non_two_torsion_ordinate_char_not2(E, f, kerp, psi2);
  } else {
    GEN h2 = gel(two_tors, 5);
    GEN g2 = gmul(gel(two_tors, 4), RgX_mul(kerq, RgX_sqr(kerq)));
    GEN g0 = non_two_torsion_ordinate_char2(E, kerq, x, y);
    g0 = gmul(g0, RgX_mul(h2, RgX_sqr(h2)));
    g = gsub(gmul(y, RgX_mul(kerp, RgX_sqr(kerp))), gadd(g2, g0));
  }
  return g;
}

/* Given an elliptic curve E and a polynomial kerp whose roots give the
 * x-coordinates of a subgroup G of E, return the curve E/G and,
 * if only_image is zero, the isogeny pi:E -> E/G. Variables vx and vy are
 * used to describe the isogeny (and are ignored if only_image is zero). */
static GEN
isogeny_from_kernel_poly(GEN E, GEN kerp, long only_image, long vx, long vy)
{
  long m;
  GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
  GEN p1, p2, p3, x, y, f, g, two_tors, EE, t, w;
  GEN kerh = two_torsion_part(E, kerp);
  GEN kerq = RgX_divrem(kerp, kerh, ONLY_DIVIDES);
  if (!kerq) pari_err_BUG("isogeny_from_kernel_poly");
  /* isogeny degree: 2*degpol(kerp)+1-degpol(kerh) */
  m = degpol(kerq);

  kerp = RgX_Rg_div(kerp, leading_term(kerp));
  kerq = RgX_Rg_div(kerq, leading_term(kerq));
  kerh = RgX_Rg_div(kerh, leading_term(kerh));
  switch(degpol(kerh))
  {
  case 0:
    two_tors = mkvec5(gen_0, gen_0, pol_0(vx), pol_0(vx), pol_1(vx));
    break;
  case 1:
    two_tors = contrib_weierstrass_pt(E, kerh, only_image,vx,vy);
    break;
  case 3:
    two_tors = contrib_full_tors(E, kerh, only_image,vx,vy);
    break;
  default:
    two_tors = NULL;
    pari_err_DOMAIN("isogeny_from_kernel_poly", "kernel polynomial",
                    "does not define a subgroup of", E, kerp);
  }
  first_three_power_sums(kerq,&p1,&p2,&p3);
  x = pol_x(vx);
  y = pol_x(vy);

  /* t = 6 * p2 + b2 * p1 + m * b4, */
  t = gadd(gmulsg(6LL, p2), gadd(gmul(b2, p1), gmulsg(m, b4)));

  /* w = 10 * p3 + 2 * b2 * p2 + 3 * b4 * p1 + m * b6, */
  w = gadd(gmulsg(10LL, p3),
           gadd(gmul(gmulsg(2LL, b2), p2),
                gadd(gmul(gmulsg(3LL, b4), p1), gmulsg(m, b6))));

  EE = make_velu_curve(E, gadd(t, gel(two_tors, 1)),
                          gadd(w, gel(two_tors, 2)));
  if (only_image) return EE;

  f = isog_abscissa(E, kerp, kerq, x, two_tors);
  g = isog_ordinate(E, kerp, kerq, x, y, two_tors, f);
  return mkvec2(EE, mkvec3(f,g,kerp));
}

/* Given an elliptic curve E and a subgroup G of E, return the curve
 * E/G and, if only_image is zero, the isogeny corresponding
 * to the canonical surjection pi:E -> E/G. The variables vx and
 * vy are used to describe the isogeny (and are ignored if
 * only_image is zero). The subgroup G may be given either as
 * a generating point P on E or as a polynomial kerp whose roots are
 * the x-coordinates of the points in G */
GEN
ellisogeny(GEN E, GEN G, long only_image, long vx, long vy)
{
  pari_sp av = avma;
  GEN j, z;
  checkell(E);j = ell_get_j(E);
  if (vx < 0) vx = 0;
  if (vy < 0) vy = 1;
  if (varncmp(vx, vy) >= 0) pari_err_PRIORITY("ellisogeny", pol_x(vx), "<=", vy);
  if (varncmp(vy, gvar(j)) >= 0) pari_err_PRIORITY("ellisogeny", j, ">=", vy);
  switch(typ(G))
  {
  case t_VEC:
    checkellpt(G);
    if (!ell_is_inf(G))
    {
      GEN x =  gel(G,1), y = gel(G,2);
      if (varncmp(vy, gvar(x)) >= 0) pari_err_PRIORITY("ellisogeny", x, ">=", vy);
      if (varncmp(vy, gvar(y)) >= 0) pari_err_PRIORITY("ellisogeny", y, ">=", vy);
    }
    z = isogeny_from_kernel_point(E, G, only_image, vx, vy);
    break;
  case t_POL:
    if (varncmp(vy, gvar(constant_term(G))) >= 0)
      pari_err_PRIORITY("ellisogeny", constant_term(G), ">=", vy);
    z = isogeny_from_kernel_poly(E, G, only_image, vx, vy);
    break;
  default:
    z = NULL;
    pari_err_TYPE("ellisogeny", G);
  }
  return gerepilecopy(av, z);
}