/*

  The tate collection

*/

#include "pari.h"

GEN cpnonsquare(GEN);
GEN cpsqrt(GEN);
GEN cprewrite(GEN);
GEN ellchangeofvars(GEN,GEN,long);
GEN elltate(GEN,GEN);

GEN cpnonsquare(GEN Gp) {
  /* this should be modified sometime to handle arbitrary sized primes, and to use gcarreparfait */
  long ltop = avma;
  long i,j,*L,p = itos(Gp);
  if (p == 2) err(talker,"cpgen does work at p=2");
  else if ((p % 4) == 3) return stoi(-1);
  else if ((p % 8) == 5) return gdeux;
  else {
    L = malloc(p*sizeof(long));
    for (i = 0; i < p; i++) L[i] = 0;
    for (i = 0; i < p; i++) L[(i*i) % p] = 1;
    for (i = 1; i < p; i++) {
      if (L[i] = 0) { j = i; break; }
      else if (L[p-i] == 0) { j = -i; break; }
    }
    free(L);
    return(stoi(j));
  }
}

GEN cpsqrt(GEN x) {
  long ltop = avma;
  GEN ns,b,c,bd,col1,col2,outmat;

  if (typ(x) != t_PADIC) err(talker,"input should be p-adic");
  if ((valp(x) % 2) != 0) err(talker,"not an element of C_p");
  ns = cpnonsquare((GEN) x[2]);
  if (kronecker((GEN) x[4],(GEN) x[2]) == -1) {
    c = gsqrt(gdiv(x,ns),0);
    b = gmul(c,ns);
    outmat = cgetg(3,t_MAT);
    col1 = cgetg(3,t_COL);
    col2 = cgetg(3,t_COL);
    col1[1] = zero; col2[1] = lcopy(b);
    col1[2] = lcopy(c); col2[2] = zero;
    outmat[1] = (long) col1; outmat[2] = (long) col2;
  } else {
    outmat = gsqrt(x,0);
  }

  return gerepileupto(ltop,outmat);
}

GEN cprewrite(GEN quad) {
  long ltop = avma;
  GEN genout,msize,a,b,c,d;
  GEN negtrw,nmw,wpol,c0;
  long lgquad,j;

  if (typ(quad) == t_VEC) {
    lgquad = lg(quad);
    genout = cgetg(lgquad,t_VEC);
    for(j = 1; j < lgquad; j++) {
      genout[j] = (long) cprewrite((GEN) quad[j]);
    }
  } else if (typ(quad) == t_MAT) {
    msize = matsize(quad);
    if (cmpii((GEN) msize[1],gdeux) != 0 || cmpii((GEN) msize[2],gdeux) != 0) err(talker,"unexpected input %Z in cprewrite",quad);
    a = (GEN) ((GEN) quad[1])[1]; b = (GEN) ((GEN) quad[2])[1];
    c = (GEN) ((GEN) quad[1])[2]; d = (GEN) ((GEN) quad[2])[2];
    if (typ(c) == t_INT) {
      if (signe(c) == 0) c0 = gun;
      else c0 = c;
    } else if (typ(c) == t_PADIC) {
      if (signe((GEN) c[4]) == 0) c0 = gun;
      else c0 = c;
    } else c0 = c;
    negtrw = gdiv(gsub(d,a),c0); nmw = gneg(gdiv(b,c0));
    genout = cgetg(4,t_VEC);
    genout[1] = lcopy(d); genout[2] = lcopy(c);
    wpol = cgetg(5,t_POL);
    setsigne(wpol,1); setlgef(wpol,5); setvarn(wpol,0);
    wpol[2] = lcopy(nmw); wpol[3] = lcopy(negtrw); wpol[4] = un;
    genout[3] = (long) wpol;
  } else genout = gcopy(quad);

  return gerepileupto(ltop,genout);
}

GEN ellchangeofvars(GEN E1, GEN E2, long prec) {
  long ltop = avma,lbot;
  GEN u,usqr,r,s,t,outvec;

  usqr = gdiv(gmul((GEN) E1[11],(GEN) E2[10]),gmul((GEN) E1[10],(GEN) E2[11]));
  if (typ(usqr) == t_PADIC) u = cpsqrt(usqr); else u = gsqrt(usqr,prec);
  s = gdiv(gsub(gmul(u,(GEN) E2[1]),(GEN) E1[1]),gdeux);
  r = gsub(gdiv(gsub(gmul(usqr,(GEN) E2[2]),(GEN) E1[2]),stoi(3)),gdiv(gsub(gsqr((GEN) E1[1]),gmul(usqr,gsqr((GEN) E2[1]))),stoi(12)));
  t = gadd(gsub(gdiv(gsub(gmul(gmul(u,usqr),(GEN) E2[3]),(GEN) E1[3]),gdeux),gdiv(gsub(gmul(gmul(usqr,(GEN) E1[1]),(GEN) E2[2]),gmul((GEN) E1[1],(GEN) E1[2])),stoi(6))),gdiv(gsub(gmul((GEN) E1[1],gmul((GEN) E1[1],(GEN) E1[1])),gmul(usqr,gmul((GEN) E1[1],gsqr((GEN) E2[1])))),stoi(24)));

  outvec = cgetg(5,t_VEC);
  outvec[1] = lcopy(u); outvec[2] = lcopy(r); outvec[3] = lcopy(s); outvec[4] = lcopy(t);
  return gerepileupto(ltop,outvec);
}

GEN elltate(GEN E, GEN c) {
  long ltop = avma,ltop2,lbot2;
  long i,N;
  GEN q,Eq,X,Y,v;
  GEN Gi,t1,t3,t5,sa,s1,s3,s5,P1,P3,P5;
  GEN T1,T2,T3,T4,sX,sY;
  GEN u,r,s,t,usqrinv,usqrsqrinv,pt;

  if (typ(E) != t_VEC || lg(E) != 20) err(talker,"first argument should be elliptic curve");
  q = (GEN) E[17];
  if (typ(q) != t_PADIC) err(talker,"first argument should be p-adic elliptic curve");
  N = precp(q);

  /* compute s1,s3,s5 */
  /* this may be a waste if gsubst doesn't use Horner's rule */
  ltop2 = avma;
  t1 = cgetg((N+3),t_VEC); t3 = cgetg((N+3),t_VEC); t5 = cgetg((N+3),t_VEC);
  t1[1] = zero; t3[1] = zero; t5[1] = zero;
  for (i = 2; i <= N+2; i++) {
    Gi = stoi(i-1);
    t1[i] = (long) sumdiv(Gi);
    t3[i] = (long) gsumdivk(Gi,3);
    t5[i] = (long) gsumdivk(Gi,5);
  }
  P1 = gtopolyrev(t1,0); P3 = gtopolyrev(t3,0); P5 = gtopolyrev(t5,0);
  sa = cgetg(4,t_VEC);
  sa[1] = lsubst(P1,0,q); sa[2] = lsubst(P3,0,q); sa[3] = lsubst(P5,0,q);
  sa = gerepileupto(ltop2,sa);
  s1 = (GEN) sa[1]; s3 = (GEN) sa[2]; s5 = (GEN) sa[3];

  /* determine the curve E_q */
  Eq = cgetg(6,t_VEC);
  Eq[1] = un; Eq[2] = zero; Eq[3] = zero;
  Eq[4] = lmul(stoi(-5),s3);
  Eq[5] = ldiv(gadd(gmul(stoi(5),s3),gmul(stoi(7),s5)),stoi(-12));
  Eq = initell(Eq,0);

  /* compute the point on E_q */
  ltop2 = avma;
  sa = cgetg(3,t_VEC);
  T1 = gun;
  T2 = c;
  T3 = gsub(gun,T2);
  T4 = gsqr(T3);
  sX = gdiv(T2,T4);
  sY = gdiv(gsqr(T2),gmul(T3,T4));
  for (i = 1; i <= N+2; i++) {

    T1 = gpowgs(q,i);
    T2 = gmul(T1,c);
    T3 = gsub(gun,T2);
    T4 = gsqr(T3);
    sX = gadd(sX,gdiv(T2,T4));
    sY = gadd(sY,gdiv(gsqr(T2),gmul(T3,T4)));

    T1 = ginv(T1);
    T2 = gmul(T1,c);
    T3 = gsub(gun,T2);
    T4 = gsqr(T3);
    sX = gadd(sX,gdiv(T2,T4));
    sY = gadd(sY,gdiv(gsqr(T2),gmul(T3,T4)));
  }
  sa = cgetg(3,t_VEC);
  sa[1] = lcopy(sX); sa[2] = lcopy(sY);
  sa = gerepileupto(ltop2,sa);
  X = gsub((GEN) sa[1],gmul(gdeux,s1));
  Y = gadd((GEN) sa[2],s1);

  /* change variables and return */
  v = ellchangeofvars(Eq,E,0);
  u = (GEN) v[1]; r = (GEN) v[2]; s = (GEN) v[3]; t = (GEN) v[4];
  usqrinv = ginv(gsqr(u)); usqrsqrinv = gsqr(usqrinv);
  T1 = gsub(X,r);
  T2 = gmul(u,gadd(gsub(Y,gmul(s,X)),gsub(gmul(s,r),t)));
  pt = cgetg(3,t_VEC);
  pt[1] = lmul(T1,usqrinv); pt[2] = lmul(T2,usqrsqrinv);
  return gerepileupto(ltop,pt);
}