/* 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"
#include "paripriv.h"

/* Not so fast arithmetic with points over elliptic curves over Fl */

/***********************************************************************/
/**                                                                   **/
/**                              Flj                                  **/
/**                                                                   **/
/***********************************************************************/

/* Arithmetic is implemented using Jacobian coordinates, representing
 * a projective point (x : y : z) on E by [z*x , z^2*y , z].  This is
 * probably not the fastest representation available for the given
 * problem, but they're easy to implement and up to 60% faster than
 * the school-book method used in Fle_mulu().
 */

/*
 * Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
 * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl
 */

INLINE void
Flj_dbl_indir_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
{
  ulong X1, Y1, Z1;
  ulong XX, YY, YYYY, ZZ, S, M, T;

  X1 = P[1]; Y1 = P[2]; Z1 = P[3];

  if (Z1 == 0)
  {
    Q[1] = X1; Q[2] = Y1; Q[3] = Z1;
    return;
  }

  XX = Fl_sqr_pre(X1, p, pi);
  YY = Fl_sqr_pre(Y1, p, pi);
  YYYY = Fl_sqr_pre(YY, p, pi);
  ZZ = Fl_sqr_pre(Z1, p, pi);
  S = Fl_double(Fl_sub(Fl_sqr_pre(Fl_add(X1, YY, p), p, pi),
                       Fl_add(XX, YYYY, p), p), p);
  M = Fl_add(Fl_triple(XX, p),
             Fl_mul_pre(a4, Fl_sqr_pre(ZZ, p, pi), p, pi), p);
  T = Fl_sub(Fl_sqr_pre(M, p, pi), Fl_double(S, p), p);
  Q[1] = T;
  Q[2] = Fl_sub(Fl_mul_pre(M, Fl_sub(S, T, p), p, pi),
                Fl_double(Fl_double(Fl_double(YYYY, p), p), p), p);
  Q[3] = Fl_sub(Fl_sqr_pre(Fl_add(Y1, Z1, p), p, pi),
                Fl_add(YY, ZZ, p), p);
}

INLINE void
Flj_dbl_pre_inplace(GEN P, ulong a4, ulong p, ulong pi)
{
  Flj_dbl_indir_pre(P, P, a4, p, pi);
}

GEN
Flj_dbl_pre(GEN P, ulong a4, ulong p, ulong pi)
{
  GEN Q = cgetg(4, t_VECSMALL);
  Flj_dbl_indir_pre(P, Q, a4, p, pi);
  return Q;
}

/*
 * Cost: 11M + 5S + 9add + 4*2.
 * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl
 */

INLINE void
Flj_add_indir_pre(GEN P, GEN Q, GEN R, ulong a4, ulong p, ulong pi)
{
  ulong X1, Y1, Z1, X2, Y2, Z2;
  ulong Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W;
  X1 = P[1]; Y1 = P[2]; Z1 = P[3];
  X2 = Q[1]; Y2 = Q[2]; Z2 = Q[3];

  if (Z2 == 0) {
    R[1] = X1; R[2] = Y1; R[3] = Z1;
    return;
  }

  if (Z1 == 0) {
    R[1] = X2; R[2] = Y2; R[3] = Z2;
    return;
  }

  Z1Z1 = Fl_sqr_pre(Z1, p, pi);
  Z2Z2 = Fl_sqr_pre(Z2, p, pi);
  U1 = Fl_mul_pre(X1, Z2Z2, p, pi);
  U2 = Fl_mul_pre(X2, Z1Z1, p, pi);
  S1 = Fl_mul_pre(Y1, Fl_mul_pre(Z2, Z2Z2, p, pi), p, pi);
  S2 = Fl_mul_pre(Y2, Fl_mul_pre(Z1, Z1Z1, p, pi), p, pi);
  H = Fl_sub(U2, U1, p);
  r = Fl_double(Fl_sub(S2, S1, p), p);

  /* If points are equal we must double. */
  if (H == 0) {
    if (r == 0) {
      /* Points are equal so double. */
      Flj_dbl_indir_pre(P, R, a4, p, pi);
    } else {
      /* Points are opposite so return zero. */
      R[1] = R[2] = 1; R[3] = 0;
    }
    return;
  }
  I = Fl_sqr_pre(Fl_double(H, p), p, pi);
  J = Fl_mul_pre(H, I, p, pi);
  V = Fl_mul_pre(U1, I, p, pi);
  W = Fl_sub(Fl_sqr_pre(r, p, pi), Fl_add(J, Fl_double(V, p), p), p);
  R[1] = W;
  R[2] = Fl_sub(Fl_mul_pre(r, Fl_sub(V, W, p), p, pi),
                Fl_double(Fl_mul_pre(S1, J, p, pi), p), p);
  R[3] = Fl_mul_pre(Fl_sub(Fl_sqr_pre(Fl_add(Z1, Z2, p), p, pi),
                           Fl_add(Z1Z1, Z2Z2, p), p), H, p, pi);
}

INLINE void
Flj_add_pre_inplace(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
{
  Flj_add_indir_pre(P, Q, P, a4, p, pi);
}

GEN
Flj_add_pre(GEN P, GEN Q, ulong a4, ulong p, ulong pi)
{
  GEN R = cgetg(4, t_VECSMALL);
  Flj_add_indir_pre(P, Q, R, a4, p, pi);
  return R;
}

GEN
Flj_neg(GEN Q, ulong p)
{
  return mkvecsmall3(Q[1], Fl_neg(Q[2], p), Q[3]);
}

typedef struct {
  ulong n; /* The number being represented */
  ulong pbits, nbits;  /* Positive bits and negative bits */
  ulong lnzb; /* Leading non-zero bit */
} naf_t;

/*
 * Return the signed binary representation (i.e. the Non-Adjacent Form
 * in base 2) of 0 <= a < 2^63.
 */
static void
naf_repr(naf_t *x, ulong a)
{
  long t, i;
  ulong pbits, nbits;
  ulong c0 = 0, c1, a0;

  x->n = a;
  pbits = nbits = 0;
  for (i = 0; a; a >>= 1, ++i) {
    a0 = a & 1;
    c1 = (c0 + a0 + ((a & 2) >> 1)) >> 1;
    t = c0 + a0 - (c1 << 1);
    if (t < 0)
      nbits |= (1ULL << i);
    else if (t > 0)
      pbits |= (1ULL << i);
    c0 = c1;
  }
  c1 = c0 >> 1;
  t = c0 - (c1 << 1);
  /* Note that we don't need to check whether t < 0, since a >= 0 implies
   * that this most significant signed bit must be non-negative. */
#if 0
  if (t < 0)
    nbits |= (1ULL << i);
  else
#endif
  if (t > 0)
    pbits |= (1ULL << i);

  x->pbits = pbits;
  x->nbits = nbits;
  /* Note that expu returns the least nonzero bit in the argument,
   * like the bit-scan-rev instruction on Intel architectures. */
  /* Using pbits here is justified by the fact that a >= 0, so the
   * most significant bit must be positive. */
  x->lnzb = expu(pbits) - 2;
}

/*
 * Standard left-to-right signed double-and-add to compute [n]P.
 */

static GEN
Flj_mulu_pre_naf(GEN P, ulong n, ulong a4, ulong p, ulong pi, const naf_t *x)
{
  GEN R, Pinv;
  ulong pbits, nbits, lnzb;
  ulong m;

  if (n == 0)
    return mkvecsmall3(1, 1, 0);
  if (n == 1)
    return Flv_copy(P);

  R = Flj_dbl_pre(P, a4, p, pi);
  if (n == 2)
    return R;

  pbits = x->pbits;
  nbits = x->nbits;
  lnzb = x->lnzb;

  Pinv = Flj_neg(P, p);
  m = (1ULL << lnzb);
  for ( ; m; m >>= 1) {
    Flj_dbl_pre_inplace(R, a4, p, pi);
    if (m & pbits)
      Flj_add_pre_inplace(R, P, a4, p, pi);
    else if (m & nbits)
      Flj_add_pre_inplace(R, Pinv, a4, p, pi);
  }
  avma = (pari_sp)R; return R;
}

GEN
Flj_mulu_pre(GEN P, ulong n, ulong a4, ulong p, ulong pi)
{
  naf_t x;
  naf_repr(&x, n);
  return Flj_mulu_pre_naf(P, n, a4, p, pi, &x);
}

ulong
Flj_order_ufact(GEN P, ulong n, GEN F, ulong a4, ulong p, ulong pi)
{
  pari_sp av = avma;
  ulong res = 1;
  long i, nfactors;
  GEN primes, exps;

  primes = gel(F, 1);
  nfactors = lg(primes);
  exps = gel(F, 2);

  for (i = 1; i < nfactors; ++i) {
    ulong q, pp = primes[i];
    long ei = exps[i];
    long j, k;
    naf_t x;
    GEN b;

    for (q = pp, j = 1; j < ei; ++j)
      q *= pp;
    b = Flj_mulu_pre(P, n / q, a4, p, pi);

    naf_repr(&x, pp);
    for (j = 0; j < ei && b[3] != 0; ++j)
      b = Flj_mulu_pre_naf(b, pp, a4, p, pi, &x);
    if (b[3] != 0)
      return 0;
    for (k = 0; k < j; ++k)
      res *= pp;
    avma = av;
  }
  return res;
}

GEN
Fle_to_Flj(GEN P)
{ return ell_is_inf(P) ? mkvecsmall3(1ULL, 1ULL, 0ULL):
                         mkvecsmall3(P[1], P[2], 1ULL);
}

GEN
Flj_to_Fle_pre(GEN P, ulong p, ulong pi)
{
  if (P[3] == 0) return ellinf();
  else
  {
    ulong Z = Fl_inv(P[3], p);
    ulong Z2 = Fl_sqr_pre(Z, p, pi);
    ulong X3 = Fl_mul_pre(P[1], Z2, p, pi);
    ulong Y3 = Fl_mul_pre(P[2], Fl_mul_pre(Z, Z2, p, pi), p, pi);
    return mkvecsmall2(X3, Y3);
  }
}

INLINE void
random_Fle_pre_indir(ulong a4, ulong a6, ulong p, ulong pi,
                     ulong *pt_x, ulong *pt_y)
{
  ulong x, x2, y, rhs;
  do
  {
    x   = random_Fl(p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
    x2  = Fl_sqr_pre(x, p, pi);
    rhs = Fl_add(Fl_mul_pre(x, Fl_add(x2, a4, p), p, pi), a6, p);
  } while ((!rhs && !Fl_add(Fl_triple(x2,p),a4,p)) || krouu(rhs, p) < 0);
  y = Fl_sqrt_pre(rhs, p, pi);
  *pt_x = x; *pt_y = y;
}

GEN
random_Flj_pre(ulong a4, ulong a6, ulong p, ulong pi)
{
  ulong x, y;
  random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
  return mkvecsmall3(x, y, 1);
}

/***********************************************************************/
/**                                                                   **/
/**                              Fle                                  **/
/**                                                                   **/
/***********************************************************************/
GEN
Fle_changepoint(GEN x, GEN ch, ulong p)
{
  ulong p1,u,r,s,t,v,v2,v3;
  GEN z;
  if (ell_is_inf(x)) return x;
  u = ch[1]; r = ch[2];
  s = ch[3]; t = ch[4];
  v = Fl_inv(u, p); v2 = Fl_sqr(v,p); v3 = Fl_mul(v,v2,p);
  p1 = Fl_sub(x[1],r,p);
  z = cgetg(3,t_VECSMALL);
  z[1] = Fl_mul(v2, p1, p);
  z[2] = Fl_mul(v3, Fl_sub(x[2], Fl_add(Fl_mul(s,p1, p),t, p),p),p);
  return z;
}

GEN
Fle_changepointinv(GEN x, GEN ch, ulong p)
{
  ulong u, r, s, t, X, Y, u2, u3, u2X;
  GEN z;
  if (ell_is_inf(x)) return x;
  X = x[1]; Y = x[2];
  u = ch[1]; r = ch[2];
  s = ch[3]; t = ch[4];
  u2 = Fl_sqr(u, p); u3 = Fl_mul(u,u2,p);
  u2X = Fl_mul(u2,X, p);
  z = cgetg(3, t_VECSMALL);
  z[1] = Fl_add(u2X,r,p);
  z[2] = Fl_add(Fl_mul(u3,Y,p), Fl_add(Fl_mul(s,u2X,p), t, p), p);
  return z;
}
static GEN
Fle_dbl_slope(GEN P, ulong a4, ulong p, ulong *slope)
{
  ulong x, y, Qx, Qy;
  if (ell_is_inf(P) || !P[2]) return ellinf();
  x = P[1]; y = P[2];
  *slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
                  Fl_double(y, p), p);
  Qx = Fl_sub(Fl_sqr(*slope, p), Fl_double(x, p), p);
  Qy = Fl_sub(Fl_mul(*slope, Fl_sub(x, Qx, p), p), y, p);
  return mkvecsmall2(Qx, Qy);
}

GEN
Fle_dbl(GEN P, ulong a4, ulong p)
{
  ulong slope;
  return Fle_dbl_slope(P,a4,p,&slope);
}

static GEN
Fle_add_slope(GEN P, GEN Q, ulong a4, ulong p, ulong *slope)
{
  ulong Px, Py, Qx, Qy, Rx, Ry;
  if (ell_is_inf(P)) return Q;
  if (ell_is_inf(Q)) return P;
  Px = P[1]; Py = P[2];
  Qx = Q[1]; Qy = Q[2];
  if (Px==Qx)
    return Py==Qy ? Fle_dbl_slope(P, a4, p, slope): ellinf();
  *slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
  Rx = Fl_sub(Fl_sub(Fl_sqr(*slope, p), Px, p), Qx, p);
  Ry = Fl_sub(Fl_mul(*slope, Fl_sub(Px, Rx, p), p), Py, p);
  return mkvecsmall2(Rx, Ry);
}

GEN
Fle_add(GEN P, GEN Q, ulong a4, ulong p)
{
  ulong slope;
  return Fle_add_slope(P,Q,a4,p,&slope);
}

static GEN
Fle_neg(GEN P, ulong p)
{
  if (ell_is_inf(P)) return P;
  return mkvecsmall2(P[1], Fl_neg(P[2], p));
}

GEN
Fle_sub(GEN P, GEN Q, ulong a4, ulong p)
{
  pari_sp av = avma;
  ulong slope;
  return gerepileupto(av, Fle_add_slope(P, Fle_neg(Q, p), a4, p, &slope));
}

struct _Fle
{
  ulong a4,a6;
  ulong p;
};

static GEN
_Fle_dbl(void *E, GEN P)
{
  struct _Fle *ell = (struct _Fle *) E;
  return Fle_dbl(P, ell->a4, ell->p);
}

static GEN
_Fle_add(void *E, GEN P, GEN Q)
{
  struct _Fle *ell=(struct _Fle *) E;
  return Fle_add(P, Q, ell->a4, ell->p);
}

GEN
Fle_mulu(GEN P, ulong n, ulong a4, ulong p)
{
  ulong pi;
  if (!n || ell_is_inf(P)) return ellinf();
  if (n==1) return zv_copy(P);
  if (n==2) return Fle_dbl(P, a4, p);
  pi = get_Fl_red(p);
  return Flj_to_Fle_pre(Flj_mulu_pre(Fle_to_Flj(P), n, a4, p, pi), p, pi);
}

static GEN
_Fle_mul(void *E, GEN P, GEN n)
{
  pari_sp av = avma;
  struct _Fle *e=(struct _Fle *) E;
  long s = signe(n);
  GEN Q;
  if (!s || ell_is_inf(P)) return ellinf();
  if (s < 0) P = Fle_neg(P, e->p);
  if (is_pm1(n)) return s > 0? zv_copy(P): P;
  Q = (lgefint(n)==3) ? Fle_mulu(P, uel(n,2), e->a4, e->p):
                        gen_pow(P, n, (void*)e, &_Fle_dbl, &_Fle_add);
  return s > 0? Q: gerepileuptoleaf(av, Q);
}

GEN
Fle_mul(GEN P, GEN n, ulong a4, ulong p)
{
  struct _Fle E;
  E.a4 = a4; E.p = p;
  return _Fle_mul(&E, P, n);
}

/* Finds a random non-singular point on E */

GEN
random_Fle_pre(ulong a4, ulong a6, ulong p, ulong pi)
{
  ulong x, y;
  random_Fle_pre_indir(a4, a6, p, pi, &x, &y);
  return mkvecsmall2(x, y);
}

GEN
random_Fle(ulong a4, ulong a6, ulong p)
{ return random_Fle_pre(a4, a6, p, get_Fl_red(p)); }

static GEN
_Fle_rand(void *E)
{
  struct _Fle *e=(struct _Fle *) E;
  return random_Fle(e->a4, e->a6, e->p);
}

static const struct bb_group Fle_group={_Fle_add,_Fle_mul,_Fle_rand,hash_GEN,zv_equal,ell_is_inf,NULL};

GEN
Fle_order(GEN z, GEN o, ulong a4, ulong p)
{
  pari_sp av = avma;
  struct _Fle e;
  e.a4=a4;
  e.p=p;
  return gerepileuptoint(av, gen_order(z, o, (void*)&e, &Fle_group));
}

GEN
Fle_log(GEN a, GEN b, GEN o, ulong a4, ulong p)
{
  pari_sp av = avma;
  struct _Fle e;
  e.a4=a4;
  e.p=p;
  return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &Fle_group));
}

ulong
Fl_ellj(ulong a4, ulong a6, ulong p)
{
  if (SMALL_ULONG(p))
  {
    /* a43 = 4 a4^3 */
    ulong a43 = Fl_double(Fl_double(Fl_mul(a4, Fl_sqr(a4, p), p), p), p);
    /* a62 = 27 a6^2 */
    ulong a62 = Fl_mul(Fl_sqr(a6, p), 27 % p, p);
    ulong z1 = Fl_mul(a43, 1728 % p, p);
    ulong z2 = Fl_add(a43, a62, p);
    return Fl_div(z1, z2, p);
  } else
    return Fl_ellj_pre(a4, a6, p, get_Fl_red(p));
}

void
Fl_ellj_to_a4a6(ulong j, ulong p, ulong *pt_a4, ulong *pt_a6)
{
  ulong zagier = 1728 % p;
  if (j == 0)           { *pt_a4 = 0; *pt_a6 =1; }
  else if (j == zagier) { *pt_a4 = 1; *pt_a6 =0; }
  else
  {
    ulong k = Fl_sub(zagier, j, p);
    ulong kj = Fl_mul(k, j, p);
    ulong k2j = Fl_mul(kj, k, p);
    *pt_a4 = Fl_triple(kj, p);
    *pt_a6 = Fl_double(k2j, p);
  }
}

ulong
Fl_elldisc_pre(ulong a4, ulong a6, ulong p, ulong pi)
{
  /* D = -(4A^3 + 27B^2) */
  ulong t1, t2;
  t1 = Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi);
  t1 = Fl_double(Fl_double(t1, p), p);
  t2 = Fl_mul_pre(27 % p, Fl_sqr_pre(a6, p, pi), p, pi);
  return Fl_neg(Fl_add(t1, t2, p), p);
}

ulong
Fl_elldisc(ulong a4, ulong a6, ulong p)
{
  if (SMALL_ULONG(p))
  {
    /* D = -(4A^3 + 27B^2) */
    ulong t1, t2;
    t1 = Fl_mul(a4, Fl_sqr(a4, p), p);
    t1 = Fl_double(Fl_double(t1, p), p);
    t2 = Fl_mul(27 % p, Fl_sqr(a6, p), p);
    return Fl_neg(Fl_add(t1, t2, p), p);
  } else
    return Fl_elldisc_pre(a4, a6, p, get_Fl_red(p));
}

void
Fl_elltwist(ulong a4, ulong a6, ulong D, ulong p, ulong *pt_a4, ulong *pt_a6)
{
  ulong D2 = Fl_sqr(D, p);
  *pt_a4 = Fl_mul(a4, D2, p);
  *pt_a6 = Fl_mul(a6, Fl_mul(D, D2, p), p);
}